Non relativistic Radiation Transport

Governing equations

Two separate radiation transport modules (enabled with -nr_radiation or -implicit_radiation) are supported, which solve the time dependent evolution equation for specific intensities alone discrete angles. The equations and numerical algorithms are described in Jiang (2021, ApJS, 253, 49), Jiang (ApJS, 2022, 263, 4) as well as an earlier paper Jiang et al. (2014, ApJS, 213, 7).

We describe the radiation field using specific intensity $I_f(t,k,j,i,n)$, which is a function of time $t$ , spatial locations $k,j,i$, angular direction $\hat{n}$ and frequency group $f$. It is defined as frequency integrated specific intensity over each frequency group $I_f\equiv \int_{\nu_f}^{\nu_{f+1}} I_{\nu} d\nu.$ If there one frequency group, the integral will be from $0$ to $\infty$ and it will be the gray approximation.

The modules solve the following radiation transport equation $$\frac{\partial I_f}{\partial t}+c\symbf{n}\cdot\nabla I_f=c(\eta_f-\chi_c I_f).$$ Here $I_f$, emissivity $\eta_f$ and opacity $\chi_f$ are all defined in the lab frame. We transform them to the co-moving frame and source terms take the format $$\rho(\kappa_s+\kappa_{af})(J_{0,f} - I_{0,f}) + \rho\kappa_{pf} \epsilon_{0,f} - \rho\kappa_{pf,e} J_{0,f}.$$ Here the co-moving frame opacities $\rho\kappa_s$, $\rho\kappa_{af}$, $\rho\kappa_{pf}$ and $\rho\kappa_{pf,e}$ correspond to $\sigma_s$, $\sigma_a$, $\sigma_p$ and $\sigma_{pe}$ arrays in the code. In the case of one frequency group, the emissivity $\epsilon_{0,f}$ will be just blackbody emission based on local gas temperature $a_rT^4/4\pi$.

The source terms are always solved implicitly. When -radiation is enabled, we also solve the transport term $c\symbf{n}\cdot\nabla I_f$ explicitly so that the time step is limited by the speed of light. When -implicit_radiation is enabled, we solve the whole equation implicitly so that the time step is not limited by the speed of light. This typically requires iterations to converge.

Units

The two modules solve the radiation transport equation using a dimensionless system. The user needs to choose a temperature unit $T_0$, density unit $\rho_0$ as well as length unit $L_0$. For a constant mean molecular weight $\mu$, the gas pressure and energy density units are $P_0\equiv k_B \rho_0 T_0/(\mu m_p)$. Gas velocity is in unit of $v_0=\sqrt{k_B T_0/(\mu m_p)}$. Therefore dimensionless speed of light is $\mathcal{C}=c/v_0$. Specific intensities are in the unit of $a_rT_0^4$ and the commonly used radiation energy density is $E_r=\int Id\Omega$ as we take $4\pi$ to be 1. Radiation flux $F_r$ in the output file is defined as $F_r=\int \symbf{n} I d\Omega$ and it is in the unit of $ca_rT_0^4$. Therefore, to convert $E_r$ and $F_r$ generated by the code to cgs, we just need to multiply them by $a_rT_0^4$ and $ca_rT_0^4$ respectively. Because radiation and gas are in different unit system, when we sum radiation and gas pressure together, we need to use $P_g+P_r * \mathcal{P}$, where $\mathcal{P}\equiv a_rT_0^4/P_0$, and the total pressure will be in the unit of $P_0$. Similarly, when we sum gas and radiation momentum together, we need to use $\rho v +F_r \mathcal{P}/\mathcal{C}$ and the total momentum will be in the unit of $\rho_0 v_0$.

Features and Limitations

• The two modules only works for Newtonian flows.
• Can work with shearing box, but no orbital advection
• Can work with SMR/AMR and all coordinate systems in 1D/2D/3D
• The source terms assume gamma law equation of state for gas

Configuration

• use -nr_radiation to turn on the explicit radiation transport module
• use -implicit_radiation to turn on the implicit radiation transport module
• They cannot be turned on simultaneously

Input File

All the parameters for the radiation module are set in the <radiation> block. Here is an example for the explicit radiation module with two frequency groups and pre-calculated $\mathcal{P}$ and $\mathcal{Crat}$

    <radiation>
    nmu            = 4      # default angular system with 80 angles per cell in 3D
    prat           = 10     # $a_rT_0^4/P_0$
    crat           = 1000   # $c/v_0$
    n_frequency    = 2      # two frequency groups
    frequency_min  = -10    # boundary of the two groups at $\nu=10 k_BT_0/h$

Another example with units specified in the input file

    <radiation>
    angle_flag       = 1        # using the angular system that rotates with local coordinate
    nzeta            = 5        # number of angles from 0 to pi/2 along polar direction
    npsi             = 6        # number of angles from 0 to pi along toroidal direction
    unit             = 1        # adopt the units in the input file
    T_unit           = 5.e4     # T_0=5.e4 K
    density_unit     = 1.e-10   # \rho_0=1.e-10 g/cm^3
    length_unit      = 7.078e16 # L_0=7.078e16 cm
    molecular_weight = 0.62     # mean molecular weight

Another example to use the implicit radiation transport module

    <radiation>
    angle_flag       = 1        # using the angular system that rotates with local coordinate
    nzeta            = 5        # number of angles from 0 to pi/2 along polar direction
    npsi             = 6        # number of angles from 0 to pi along toroidal direction
    unit             = 1        # adopt the units in the input file
    T_unit           = 5.e4     # T_0=5.e4 K
    density_unit     = 1.e-10   # \rho_0=1.e-10 g/cm^3
    length_unit      = 7.078e16 # L_0=7.078e16 cm
    molecular_weight = 0.62     # mean molecular weight
    nlimit           = 50       # maximum number of iterations to try
    error_limit      = 1.e-6    # tolerance level of iteration

Check out the regression tests in nr_radiation and implicit_radiation for more examples.

Boundary condition

Boundary condition for specific intensities share the same boundary condition flag as hydro variables. For example, periodic, outflow, reflecting can be used. In order to use different boundary conditions for radiation and hydro variables, user defined boundary condition will be used. For radiation, it will be the function EnrollUserRadBoundaryFunction.

User defined Opacity

All the absorption and scattering opacities can be calculated based on local gas quantities. This is done with a user defined opacity function, which can be enrolled using the function EnrollOpacityFunction provided by the radiation class. Opacities in the ghost cells also need to be calculated in this function.

Here is one example. First, define a function

void StarOpacity(MeshBlock *pmb, AthenaArray<Real> &prim);

Then, Inside the function void MeshBlock::InitUserMeshBlockData(ParameterInput *pin), enroll the opacity function as

pnrrad->EnrollOpacityFunction(StarOpacity);

Inside the Opacity function StarOpacity(MeshBlock *pmb, AthenaArray<Real> &prim), the user needs to set the four opacity arrays pmb->pnrrad->sigma_s(k,j,i,fre), pmb->pnrrad->sigma_a(k,j,i,fre),pmb->pnrrad->sigma_p(k,j,i,fre)and pmb->pnrrad->sigma_pe(k,j,i,fre),which are functions of spatial location k,j,i and frequency fre. These opacity arrays should have unit of 1/length. The array sigma_s+sigma_a should be the Rosseland mean opacity while sigma_s is the scattering component. sigma_p and sigma_pe are typically set to the Planck mean opacity.