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 must come 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, kinetic energy will be converted into heat that is then radiated away. For an accretion rate \dot{m}, the rate at which this heat is produced, or the accretion 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 that the gas that is accreted from the surrounding zones contributes 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 energy in the code

At each time step i, the accretion energy is computed as:

E_i=\frac{G M dm}{R }

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

Kernel

After computing the accretion energy for the particle, this energy is smoothed over a kernel of cells, which later becomes a source term for the radiative transfer solver. To ensure that the sum of the differential energies over the kernel 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, \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, 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 E 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 accretion energy work array for each cell as:

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

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

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

Note, that the source is being populated by a specific accretion energy (recall k has units of 1/volume), averaged over a radiative time step.

Physical meaning of the accretion energy source term

-how it couples to the gas, how it couples to the radiation, some equations will go here.

Tests

2D, fixed grid - radiating sink

Energy output

Different opacities

Marshak Waves

2D, fixed grid - collapsing clump