Research Notes

Brightness of spiral arms

Understanding generic hydrodynamical setups of disk-planet systems in Athena++

The biggest problem I have had is to develop a better understanding of frames of reference in which the equations are needed to be implemented in the code, say for example, in the MeshBlock::ProblemGenerator. The basic structure of a file looks like as follows

// Define initial conditions
void MeshBlock::ProblemGenerator(ParameterInput *pin) {
	// The idea now is to set the initial conditions at each grid-point
	// in the simulation domain. Initial conditions refer to the initial
	// density and momentum of the constituent fluids that make up the disk
	// In case of disk-planet interactions, the "planet" itself is not part of
	// initial conditions (even if it is assumed to start exerting its
	// gravitational influence at t=0), but rather in a separate source
	// function.
	
	// Assume the program has been configured for spherical coordinates.
	for (int k = ks; k <= ke; ++k) {
		// Notice the order of loops, the first index is `x3`, or in the case
		// of spherical coordinates, φ.
		Real x3 = pcoord->x3v(k); // φ
		// Also note the special variables offered by the `MeshBlock` class
		// in the form of `ks`, `ke`, `js`.. etc, providing the start and
		// end indices for us to loop over. Later when setting boundary conditions
		// we will see how handy they are when ghost cells also
		// get involved.
		for (int j = js; j <= je; ++j) {
			Real x2 = pcoord->x2v(j); // θ
			for (int i = is; i <= ie; ++i) {
				Real x3 = pcoord->x1v(i); // r
				
				// Now goes the logic to set the conditions. Again,
				// one must set the initial density and momentum at least.
				// Depending on the EOS, energy and pressure also need to be
				// set. One important thing to note is often in these problems
				// we assume locally isothermal or vertically isothermal EOS.
				// This needs to be set by the user explicitly as the 
				// `--eos=isothermal` configuration option while compiling the code
				// is only applicable to globally isothermal EOS. 
				
				// A simple way to understand this is through sound speed.
				// Globally isothermal implies `c_s` is constant everywhere.
				// Vertically isothermal implies `c_s` has some dependence in
				// radial direction.
			}
		}
	}
}

Understanding the planetary potential

Most model papers define gravity exerted by the planet by specifying its potential . Huang+25 uses the following potential

In the code, we have the 3 variables for each direction, depending on the choice of our coordinate system. For simplicity, let us work in spherical coordinates. is a scalar quantity, but the acceleration corresponding to this potential will be a vector. This acceleration , can be calculated as

Before we substitute in above equation it is nice to simplify it a bit. Let us define some more variables:

This reduces the potential to:

Now, let’s start by writing :

Similarly, the indirect term,

The third term has no dependence on so it vanishes. In the code, this is implemented as follows:

// In the function `PlanetaryGravity` in disk_planet_spherical.cpp and solar_system.cpp
// this d
distance_square(i) = SQR(x_dis(i)) + SQR(y_dis(i)) + SQR(z_dis(i)); 
// this is R - R_p cos (Δ) 
R_dis(i) = R_arr(i) - rad_planet[p] * std::cos(phi_arr(i) - phi_planet_move); // 
// this is GM_p / (\tilde{d^3}) (assuming `PlanetaryGravityOrder == 2`)
gravity(i) = (planet_gm / std::pow(distance_square(i) + SQR(rad_soft[p]), 1.5));
// this is the indirect term
indirect_acc_R(i) = planet_gm * std::cos(phi_arr(i) - phi_planet_move) / SQR(rad_planet[p]);

Similarly, one can derive expressions for and . The spherical polar components then can be obtained through simple co-ordinate transformation.

Implementing the Bae+25 potential

Bae+25 defines their potential as

when . For , there is no softening of the potential. So let’s proceed by calculating (the other two terms are same):

We can now now compare the terms on how to modify the code. Essentially we need to make the following replacement:

which translated to following in the code

Real d = std::sqrt(distance_square(i));
Real epsilon = rad_soft[p];
if (d <= epsilon) {
	Real x = d / epsilon;
	gravity(i) = planet_gm * (4.0 - 3.0 * x) / (epsilon * epsilon * epsilon);
} else {
	// No softening
	gravity(i) = planet_gm / (d * d * d);
}

The Hydro-RT iteration

For the solar-system models, we started with the Chiang & Goldreich (1997) mid-plane profile given by their Equation (14a)

At the location of Jupiter () this implies . We can derive the other parameters needed as follows:

To make scale-free, we can divide it by the Keplerian velocity

After putting in the values this gives . -value is simply , the exponent of the temperature profile. The surface density exponent, -value can subsequently be calculated. We start with following relation:

This is true for vertically isothermal disks. So we have

Ultimately, Chiang & Goldreich (1997) says to have , so we can put values of and find , which comes out to be .

Now, upon performing a radiative transfer run, we will have the dust temperature profile. We take dust temperature equal to gas temperature at midplane and get a new estimate to improve our initial MMSN profile (reason for this is LTE, but look at Appendix B of Speedie et al. (2022) for a more quantitative treatment). The above hydro iteration gave me a temperature of . First step is get the new , which is straightforward, we just need to put the values with temperature now as 43 K instead of 74 K. This gives new value of . To get the new value for , we need to fit the midplane temperature profile. I used np.polyfit to get , and thus .

Graph View