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’s 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);
}