Radiation feedback from sink particles

The amount of thermal radiation produced in the grid is a function of temperature. Since sinks are a subgrid model they themselves do not have temperature (we are not sure how big the forming star is, how fast it is growing by contraction, etc., so there isn't an easy way of assigning the sub-grid object a 'temperature'). Thus, we are left to estimate the amount of radiation produced by the sink. We envision the sink particle as a protostar, and not yet producing radiation through fusion. Thus, the radiation our sinks produce comes from accretion.

Accretion Luminosity

For spherical symmetry, a gas parcel starting from rest and freely falling to the star from infinity will have its kinetic and gravitational energy balance at the stellar surface:

\frac{1}{2} m v_{ff}^2 = \frac{GmM_{*}}{R_{*}}

As material passes through the accretion shocks at the surface of the star, its kinetic energy is converted into heat that is then radiated away. For an accretion rate \dot{m}, the rate at which this heat is produced, or the luminosity L, is given by:

L = \frac{1}{2}\dot{m} v_{ff}^2 = \frac{G\dot{m}M_{*}}{R_{*}}

Since we do not track energy accretion onto the sink, we are left to assume gas that is accreted from the surrounding zones contribute to this accretion luminosity directly. Thus, the best we can do for tracking the energy released from infall is to calculate the RHS of this equation in the code and use it as an estimate of the true accretion luminosity. (By the way, this form of the accretion luminosity was shown to be a good approximation for our purposes ​here).

The accretion energy (L*dt) will then be distributed smoothly in a kernel surrounding the sink every time step. From there it will diffuse away from the sink via FLD radiative transfer. In this way, sinks will act as additional sources of radiation within the grid.

Tracking accretion luminosity in the code

At each time-step i, the luminosity is computed as:

L_i=\frac{G M_i dm_i}{R ~dt}

where M_i is the mass of the sink particle at i, dm_i is the total accreted mass for that time-step, G is the gravitational constant in computational units, R is the radius of the star, taken to be 1 solar radius by default, but modifiable by the user at run-time, and dt is the hydro time-step.

Kernel

To ensure that the sum of the differential energies over the kernel (i.e. volume of cells that we will be injecting the accretion energy into) equals the total accretion energy calculated, we have the following equation:

\sum \Delta E_i*dV_i=E ~~~~~~(1)

where E is the accretion energy from the time step, \Delta E_i is the differential amount of E to be distributed in the ith cell, and dV_i is the volume of the ith cell. As of now, the units don't balance in this equation. We then need to find a normalization constant that has units of 1/volume, which we will do next.

We want the amount of E in each cell to drop off smoothly with radius away from the sink. For this we choose a decaying exponential. Let,

\Delta E_i= k E e^{-r_i/\sigma} ~~~~~~~(2)

where \sigma is a scaling factor. In the code the exponential function used is such that it falls smoothly to zero at the boundary of the kernel, and covers 4 e-foldings. Now, to solve for the normalization constant, we insert (2) into (1):

\sum k E e^{-r_i/\sigma}dV_i=E

and solve for k:

k= \frac{1}{\sum e^{-r_i/\sigma} dV_i}

While the set of equations for the kernel is arbitrary, as we will see in the next section, this normalization constant allows us to easily feed into the source function a specific accretion energy (i.e. E/V), which is necessary for the code's solvers.

Work arrays and feeding L into the radiative source function

There is a subroutine in the code that is called by the radiative transfer module, 'apply kernel to work array'. This goes through and populates the work array index for each cell's accretion energy as:

info\%workarray(i,j,k,iaccretion\_energy= k*E*e^{-r/\sigma}/dt

This is an average luminosity over the hydro step. This is then fed into the source function as,

source(i,j,k,iaccretion\_energy)= info\%workarray(i,j,k,iaccretion\_energy)*dt_{rad}

where we now enter the total specific (recall k has units of 1/volume) accretion energy into the source function averaged over a radiative timestep. And now, we are ready to begin testing.

Tests

2D, fixed grid - radiating sink

Energy output

Different opacities

Marshak Waves

2D, fixed grid - collapsing clump