I'll be implementing two changes to the Bondi Accretion module. The first will be a flag which will force sink particles to only accrete the *radial* momentum of gas (see previous blog post for motivation). The second is a check that the accretion volume does not contain any cells that violate the Jeans condition following Bondi accretion. If so, the excess gas will be removed.

Radial momentum accretion fix: description and equations

Here is a diagram defining the relevant variables. Note, in what follows primes (') will denote accreted quantities.

The momentum vector of each cell in the accretion volume is decomposed into its radial and transverse components (relative to the location of the sink particle):

Following accretion, only the radial component of the momentum will be accreted:

where

and

(i.e. fm = the fraction of mass to be left in the cell following accretion; the radial momentum will be accreted proportionally to mass).

Subtracting the first two equations gives us an equation for the update of the momentum vector:

which in cartesian coordinates is:

This last equation gives the accreted momentum in each cell within the accretion volume, for the case of radial momentum accretion.

• From discussion on Wednesday, see this page.

• For planet density in the high-flux case (assuming recombination-limited flows), see previous blog post.

Runs

Currently running small planet, low flux case (still recombination-limited, but looks like thanks to tidal forces the \tau=1 surface is staying well above the reset radius). Modified parameters (from first link above):

Also queued high flux, same size planet case. Python script​ was modified so that the planet profile is calculated assuming the totally-neutral \tau=1 surface is at the planet radius (though this is absolutely not the case for the new density). Modified parameters:

Code stuff

Got another soft crash on the small planet - need to investigate (race condition seems likely - Jonathan suggested layout_comms).

Running ionization-only test with radiation pressure implemented again to confirm and investigate results (second link above).

Further thoughts

Do we want to think about including a stellar wind at all? My thought: maybe for later simulations, but not for the parameter study, and probably not for the radiation pressure.

We had discussed possibly increasing the distance in y to capture more of the up- and downwind arms — may be nice for production runs at least (I wouldn't expect it to break anything, but could be worth testing first).

New Work

1. Plotted mass accretion rate using mass inside spheres of different radii around point particles.
2. Plotted mass budget with time.
3. Plotted center of mass of various components with time.

Results

​mass_results.pdf

Next steps

• Compute drag force.

I put a 1 solar mass object at the origin of the simulation box. There is no outflow from 0 to 6000 day. I list the physical quantities during the quiet phase below:

Radius of the object that has outflow: r=46.875R_{\odot}

Mass of the object: 1M_{\odot}

Density of the outflow: \rho=10^{-13}g\cdot cm^{-3}

Temperature of the object: T_{1}=30000K

After the quiet phase, there is an outflow. An important constraint is the momentum budget. It is related to the luminosity of the object by: \dot{p}\le\frac{L}{c}. Model ID 6 is the old simulation with only 500 days of quite phase.

I think the implied luminosity is still too high for the outflow. Also, it remains a question that what kind of driving mechanism is this outflow? Model 7 seems to be a nova outburst and model 6 seems to be a supernova….

Constants:

R_p = 1.5\times10^{10} \text{ cm}

\gamma = \frac53

M_p = 0.263 M_J

T_p = 3000 \text{ K}

Assume every photon entering a tube of dimensions (R₀-R_P, dx, dx) is absorbed by the time it reaches R_P, and all the hydrogen in the tube starts ionized (X=1). Using an adiabatic atmosphere, the density profile is given by

n(r) = n_p(R_{crit}(\frac1r - \frac1{R_0}))^{\frac1{\gamma-1}},

with R_{crit} = \frac{(\gamma-1)G M_p}{\gamma c_p^2 R_p}, R_0 = \frac{R_{crit}}{R_{crit}-1}, and n_p the number density at the surface of the planet (to be solved for).

Volumetric recombination rate is

\frac{dn_0}{dt} (X,r) = (X n_H(r))^2 \alpha,

\alpha = 2.59\times10^{-13} (at the wind temperature of 10⁴ K).

Total recombinations (per unit time) are then given by

dx^2 \int_{R_0}^{R_p} \frac{dn_0}{dt}(X,r) dr

Total photons entering the tube per unit time are given by

J_{phot} dx^2,

where J_{phot} is the photon number flux.

Solving the equality​ of these two quantities for n_p, find

n_p = 1.171\times10^{-17}\sqrt{J_{phot}}.

For the currently-running case, our flux of 2x10¹³ gives a mass density of \rho_p = 1.2\times10^{-11}, much larger than our current value of 1.32x10^-16 (but we're not recombination-limited in the current case). For the high-flux case, J_phot = 2x10¹⁷,

\rho_p = 5.23\times10^{-9}.

So I think I digested the Krumholz accretion paper sufficiently well and will attempt to distill its key points regarding angular momentum accretion onto a sub-grid sink particle, as well as how this may ultimately connect to outflows and jets launched from sink particles.

The Krumholz accretion algorithm is fairly sophisticated in that it calculates a modified Bondi accretion rate for rotating gas within the accretion volume, and then only removes gas from each grid point that is bound to the sink particle (drho = Mdot ~dt).

Currently, in AstroBEAR, the Bondi accretion routine also removes linear momentum and specific internal energy from each grid point based on the fraction of mass this drho constitutes. In doing so, each component of the linear momentum is reduced by the same fraction. This reduces the magnitude of the momentum vector in the ith cell, but not its orientation. In other words, this results in the loss of angular momentum within that zone, and thus over the kernel:

where L is the total angular momentum over the kernel, and the subscripts denote the given quantity in the ith cell. Note the radial vector points along a line connecting the cell center to the sink particle.

This loss in angular momentum is stored in the sink particle, which enables the accretion algorithm to effectively conserve mass, linear momentum, and angular momentum between the grid and the sink particle over an accretion step.

Now, onto the angular momentum problem of this routine…

Krumholz makes the argument (and after some thought I think I agree with him) that given the sink particle represents an object much smaller than the grid scale (i.e. think single star compared to protostellar core) the angular momentum of the sink particle is negligible to that of the gas on the scale of the accretion volume. For a rough idea, let's assume the accretion volume represents a protostellar core, than its angular momentum is approximately:

L_{core}=(100 Solar Masses)(1 pc)(km/s)

whereas the angular momentum of the sun is roughly:

L_{sun}=(1 Solar Masses)(10^{-8} pc)(km/s)

This argument (known as the 'angular momentum problem of star formation,' c.f. this ​review) implies single stars simply cannot acquire the enormous angular momentum of their host cores and clouds: else they would rotate faster than their break up speed.

Thus, Krumholz decomposes the linear momentum vectors of each cell into a parallel component and a transverse components (wrt to the line connecting the sink particle and cell center) and only accretes the radial component of the momentum. That is, since the angular momentum of the object the sink particle represents is negligible compared to that of the accretion volume, it should be ignored. By only accreting the radial component of the momentum vector the angular momentum of the gas is preserved (as well as the radial velocity).

Connection to outflows

If one does not accrete angular momentum from the grid, then one cannot inject angular momentum into the grid from sink particles upon the jet launching step. Thus, if we choose to adopt this way of looking at things, our injected protostellar jets will not be rotating on their own (but perhaps will spin up due to the bulk motion of gas within the accretion volume).

This contrasts with the other school of thought, namely Federrath's sink particle algorithms (which are currently implemented in the code). This method DOES accrete angular momentum from the grid (as described above for the Bondi accretion routine), and thus under those conditions, injecting a rotating jet back into the grid should (hypothetically) be self-consistent and conservative.

• Summarized important elements of radiation pressure papers (from Luca).

• Now that compute nodes are back up, can restart simulation (equilibrium state to restart from) from where it was killed.

• Registration for next semester — Geometric Methods in Fluids?

The Bondi accretion routine begins by calculating the Bondi radius of the sink particle's host cell using the following equation:

RB = GM/(cs²+vs²)

This equation gives the sonic surface of infalling gas onto a moving point particle.

Next the routine defines a region around the sink particle (on the finest AMR level) from which gas will be accreted. Within this region, the accuracy of the hydrodynamical solution is lost and so we want to minimize its size. However, we want to accurately-enough resolve inflowing material into this accretion region. To balance these two effects we choose a radius = 4dx. This is called r_acc in the code.

The amount of gas removed from the cells within r_acc is given by the accretion rate * the timestep. This accretion rate should be carefully chosen when the sonic radius lies within r_acc. This is because the accretion rate sets the resultant pressure in the kernel, and so if the accretion rate is too high (the resultant pressure too low) the hydrodynamical solution in cells surrounding the kernel could be degraded in the case of subsonic inflow.

In the other limit: R_BH > R_acc, then too high of an accretion rate shouldn't matter so much. The only issue I can foresee is if the accretion rate is too low, in which case gas would pile up within the kernel and eventually artificially fragment since new sink particles are not created within the accretion region of already existing particles. For this reason it might make most sense to allow sinks to form within the accretion regions of already existing sinks (if R_BH> R_acc) and then immediately merge (within a time-step).

Next the accretion rate is calculated using a semi-analytical expression found by Ruffert & Arnett (1994). See Krumholz 2004, eqn. 11. To calculate rho_inf in that equation, we calculate an average rho within the kernel using assigned weights to each cell within the kernel. These weights depend on the ratio R_BH : r_acc and effectively allow cells with higher inflow speeds to contribute more mass to the accreting sink particle. Correspondingly, these weights are then used to calculate the mass removed by any cell within r_acc during the accretion step.

Now this is all well and good if the gas is not rotating within the kernel. However, If it is rotating then the mass should not be accreted spherically symmetrically from all cells within the kernel. To calculate a correction due to rotation, each cell within the kernel is broken up into 8³ sub-cells each with 1/8³ the mass, momentum, and energy of the host cell. The minimum distance between the sink particle and these sub-cells is then calculated. If this distance is not less than dx/4, then that sub cell does not contribute to the accreted mass (it has too high an angular momentum to ever fall onto the surface of the sink). Additionally, if any of the cells within the accretion region is unbound, it does not contribute mass to the sink particle. Note, there is a maximum amount of mass any cell can contribute to the sink particle of 25%.

Now with regard to momentum changes to cells within the accretion radius, I was surprised to learn that the Krumholz's model for Bondi accretion reduces the radial momentum only (along a line between cell center and sink particle), in order to preserve the *radial velocity and angular momentum of the gas*. This is because the sink particle is envisioned as being so small in comparison to the grid size that its accreted angular momentum is negligible.