An Implicit Finite Volume Scheme to Solve the Time Dependent Radiation Transport Equation Based on Discrete Ordinates
DDraft version February 5, 2021
Typeset using L A TEX twocolumn style in AASTeX63
An Implicit Finite Volume Scheme to Solve the Time Dependent Radiation Transport EquationBased on Discrete Ordinates
Yan-Fei Jiang( 姜 燕 飞 ) Center for Computational Astrophysics,Flatiron Institute,New York, NY 10010, USA
ABSTRACTWe describe a new algorithm to solve the time dependent, frequency integrated radiation transport(RT) equation implicitly, which is coupled to an explicit solver for equations of magnetohydrodynam-ics (MHD) using
Athena++ . The radiation filed is represented by specific intensities along discreterays, which are evolved using a conservative finite volume approach for both cartesian and curvilinearcoordinate systems. All the terms for spatial transport of photons and interactions between gas andradiation are calculated implicitly together. An efficient Jacobi-like iteration scheme is used to solvethe implicit equations. This removes any time step constrain due to the speed of light in RT. Weevolve the specific intensities in the lab frame to simplify the transport step. The lab-frame specificintensities are transformed to the co-moving frame via Lorentz transformation when the source term iscalculated. Therefore, the scheme does not need any expansion in terms of v/c . The radiation energyand momentum source terms for the gas are calculated via direct quadrature in the angular space. Thetime step for the whole scheme is determined by the normal Courant—Friedrichs—Lewy condition inthe MHD module. We provide a variety of test problems for this algorithm including both opticallythick and thin regimes, and for both gas and radiation pressure dominated flows to demonstrate itsaccuracy and efficiency.
Keywords:
Computational astronomy — radiative transfer INTRODUCTIONThermal properties of the gas are determined by emis-sion, absorption and scattering of photons as well astransport of the radiation field through the gas. In cer-tain conditions, for example the flow around compactobjects and in massive stars, the intense radiation fieldcan be the dominant source of pressure force via themomentum exchange between photons and gas, whichcontrols the dynamics of the system. Solving the radi-ation transport (RT) equation together with magneto-hydrodynamic (MHD) simulations becomes a very im-portant task in numerical simulations of various astro-physical systems (Teyssier 2015, and references therein).RT equation is typically much more expensive to solvecompared with the cost of MHD simulations. The spe-cific intensity, which is the fundamental quantity to de-
Corresponding author: Yan-Fei Jiangyjiang@flatironinstitute.org scribe the radiation field, is a function of time, spatiallocations, propagating directions and frequency. Evenadopting the gray approximation (frequency indepen-dent), resolving the angular distribution of photons isstill harder to do compared with MHD. In addition, pho-tons travel with the speed of light in the optically thinregime, which is much faster than the flow speed of manynon-relativistic systems. While in the very opticallythick limit, photon mean free path can be much smallerthan any other length scale of interest. Capturing differ-ent limits of the RT equation self-consistently while stillevolving the MHD equations efficiently can be challeng-ing. For these reasons, different assumptions have beenmade to simplify the RT equation. One type of the sim-plified approach is to take the moments of specific inten-sity and make certain assumptions to close the momentequations. Flux-limited diffusion (FLD) is a widely usedmethod (e.g., Levermore & Pomraning 1981; Turner &Stone 2001; Krumholz et al. 2007; van der Holst et al.2011), which assumes that the mean radiation energydensity follows the diffusion equation with a modified a r X i v : . [ a s t r o - ph . I M ] F e b Jiang diffusion coefficient to limit the effective photon speedto be smaller than the speed of light. This assumptionis well known to have issues in the optically thin regime.But more important, in the regime with intermediateoptical depth when it is unclear whether diffusion ap-proximation can be made or not, simulations adoptingthe FLD approximation can result in misleading conclu-sions (Krumholz & Thompson 2013), which will not beknown until the results are compared with simulationswith more accurate methods (Davis et al. 2014; Tsang &Milosavljevi´c 2018; Smith et al. 2020). Another increas-ingly popular alternative is the M1 method (Dubroca& Feugeas 1999; Pons et al. 2000; Gonz´alez et al. 2007;Skinner & Ostriker 2013; S¸adowski et al. 2013; McKin-ney et al. 2013; Skinner et al. 2019), which evolves thezeroth and first moments of specific intensity with an as-sumed closure relation based on local radiation energydensity and flux. It is often claimed that the M1 ap-proach is accurate in both optically thin and opticallythick regimes because it can handle the shadow cast byone beam of radiation field. This is definitely incorrectas the assumed closure relation in M1 can easily fail inother cases. For example, it will merge multiple beamsof photons into one since radiation field is treated asanother fluid in this method.Other methods have been developed to solve the equa-tion for specific intensities directly without assumingany closure relation. Monte Carlo based RT approachhas been explored widely (Whitney 2011, and referencestherein), which is able to handle various radiative pro-cess accurately. Some novel schemes have also beendeveloped to reduce the noise and computational costof this type of method in different regimes (Densmoreet al. 2007; Steinacker et al. 2013; Roth & Kasen 2015;Tsang & Milosavljevi´c 2018; Foucart 2018; Smith et al.2020; Ryan & Dolence 2020). Another type of widelyused method is based on discrete ordinates, or S n , whichsolves the transport equation along discrete angles. Par-ticularly, formal solution to the RT equation can be usedas a closure to the radiation moment equation so that itcan remove the limitations of FLD and M1. This vari-able Eddington tensor (VET) scheme has been used fora wide range of applications (Stone et al. 1992; Hayes& Norman 2003; Gehmeyr & Mihalas 1994; Jiang et al.2012; Davis et al. 2012; Asahina et al. 2020). The key ofthis approach is that the full time dependent evolutionof radiation MHD is handled by the radiation momentequation, while the closure is given by the formal solu-tion to the time independent transport equation. Eventhough the RT equation is solved twice in different ways,it is often assumed that it will still be more efficient com-pared with solving the full time dependent RT equation directly. In this paper, we show that it is actually possi-ble to solve the full time dependent RT equation in a fi-nite volume approach efficiently with the computationalcost comparable or even better than VET in some cases.The finite volume approach also makes it much easier touse in curvlinear coordinates, which can be complicatedfor the VET scheme.The transport term in the RT equation is easier tosolve in the laboratory frame while the interaction termsbetween radiation and gas are greatly simplified in theco-moving frame of the fluid. For non-relativistic sys-tems where the flow velocity v is much smaller thanthe speed of light c , RT equation is typically simplifiedby transforming the source terms from the co-movingframe to the laboratory frame only to certain orders of v/c (Mihalas & Mihalas 1984; Mihalas & Klein 1982).However, it has also been pointed out by Lowrie et al.(1999) that expansion only to the first order of v/c is notenough. Some second order v/c terms are needed to getthe correct asymptotic limit in different regimes, partic-ularly the dynamic diffusion regime (Krumholz et al.2007). We want to point out that these expansionscan be avoided completely if the specific intensities aresolved directly. We can apply the Lorentz transforma-tion to connect the radiation source terms between dif-ferent frames. Even for non-relativistic systems, thisapproach is simple, accurate for all regimes and doesnot cause too much additional computational cost orcomplexity.When the typical flow speed in the system is close tospeed of light, for example around compact objects, itis more efficient to solve the transport term of the RTequation explicitly by limiting the time step based on thespeed of light (Jiang et al. 2014; S¸adowski et al. 2013;McKinney et al. 2013; Anninos & Fragile 2020; Asahinaet al. 2020). However, for a wide range of applicationswith non-relativistic systems, it is necessary to solve theMHD equation explicitly while the RT equation is solvedimplicitly so that the time step is not limited by thespeed of light. This is the approach adopted by mostalgorithms when the radiation moment equations aresolved. When specific intensities are solved directly, thetime dependent term is typically neglected. There havebeen several attempts to solve the full time dependentRT equations implicitly with mixed results (Stone & Mi-halas 1992; Sumiyoshi & Yamada 2012). When the timestep is much larger than the light cross time per cell, thechallenge is to keep a good balance between the trans-port term and the source term in all the asymptoticregimes. This is because the RT equation is a simpleadvection equation in the optically thin regime and itbecomes a diffusion equation in the optically thick case. mplicit RT Algorithm EQUATIONSThe radiation MHD equations we solve are similar towhat are used by Jiang et al. (2014). However, we donot expand the source terms in orders of v/c . We willdescribe the RT equations and how they are coupled tothe MHD equations in the following subsections.2.1.
Equations for Radiation Transport
We describe the radiation field using frequency inte-grated specific intensity I in the lab frame, which is afunction of time t , spatial location ( x, y, z ) and angulardirection as defined by the unit vector n . The trans-port term is greatly simplified in the lab frame whilethe source terms describing the interactions between ra-diation and gas are convenient in the co-moving frameof the fluid. Therefore, we also adopt this mixed frameapproach (Lowrie et al. 1999; Jiang et al. 2012, 2014).The RT equation we solve is (Mihalas & Klein 1982;Mihalas & Mihalas 1984) ∂I∂t + c n · ∇ I = c ( η − χI ) , (1)where c is the speed of light. The emissivity η and opac-ity χ are formally defined in the lab frame in the aboveequation. But the whole source term is calculated in theco-moving frame via the approach below. Traditionally,in order to get the expression for the source term, expan-sion to the first order of O ( v/c ) is usually used (Mihalas& Klein 1982), with some additional second order terms(Lowrie et al. 1999; Krumholz et al. 2007; Jiang et al.2014). Here we apply Lorentz transformation to the spe-cific intensities to avoid the need of expansion in termsof v/c completely. Although this may not sound neces-sary for non-relativistic systems, this approach does notcause additional computational cost and it results in amuch cleaner calculation of the source terms.The frequency integrated specific intensities in thelab frame I ( n ) is related to the co-moving frame val-ues I ( n (cid:48) ) via the Lorentz transformation as (Mihalas & Mihalas 1984) I ( n (cid:48) ) = γ (1 − n · v /c ) I ( n ) ≡ Γ ( n , v ) I ( n ) , (2)where v is the flow velocity, γ ≡ / (cid:112) − v /c is thecorresponding Lorentz factor and n (cid:48) is the angle in theco-moving frame given by n (cid:48) = 1 γ (1 − n · v /c ) (cid:20) n − γ v c (cid:18) − γγ + 1 n · v c (cid:19)(cid:21) . (3)We have defined Γ( n , v ) ≡ γ (1 − n · v /c ) for conve-nience. We multiply the left hand sides of the sourceterm with Γ to get ∂I ∂t = c Γ [ η − χ I ] . (4)Here the frequency integrated emissivity ( η ) and opac-ity ( χ ) in the co-moving frame are related to η, χ as η = Γ η, χ = Γ − χ. (5)Under the assumption of local thermal equilibrium, wecan write the frequency integrated source term in theco-moving frame as ∂I ∂t = c Γ [ ρκ s ( J − I )+ ρκ a (cid:18) a r T π − I (cid:19) + ρκ δP (cid:18) a r T π − J (cid:19)(cid:21) = c Γ [ ρ ( κ s + κ a ) ( J − I )+ ρ ( κ a + κ δP ) (cid:18) a r T π − J (cid:19)(cid:21) , (6)and the corresponding equation for gas internal energy E g : ∂E g ∂t = − cρ ( κ a + κ δP ) (cid:0) a r T − πJ (cid:1) . (7)Here a r is the radiation constant and ρ, T are gas den-sity and temperature, which should be co-moving framequantities in principle. When the RT algorithm is cou-pled to a Newtonian MHD solver, we do not distinguishgas quantities as well as time in the lab and co-movingframes. The angular averaged mean intensity in the co-moving frame J is defined as J ≡ π (cid:90) I d Ω , (8)where Ω is the angular element in the co-movingframe. The specific scattering opacity κ s (with unitof cm /g ) is assumed to be isotropic and does notchange the mean photon energy such as electron scat-tering. The Rosseland mean absorption opacity is κ a Jiang and we take κ δP to be the difference between thePlanck mean and Rosseland mean opacity. This is de-signed to ensure that when we integrate the zeroth andfirst moment of the co-moving frame specific intensityover the angular space, we get the energy and mo-mentum source terms ρ ( κ a + κ δP ) (cid:0) a r T − πJ (cid:1) and − ρ ( κ s + κ a ) (4 π H ) with H ≡ (cid:82) I n (cid:48) d Ω / (4 π ). Thisis consistent with the Rosseland and Planck mean opac-ities normally used for radiation moment equations. Wemultiply both sides of equation 6 with Γ − to get thesource term for lab-frame specific intensities.2.2. Equations of Radiation MHD
Photons are coupled to the gas via momentum andenergy exchange. Since the total energy and momentumare conserved in the lab frame, we integrate the sourceterms for lab frame specific intensities over the angles toget the full radiation MHD equations as ∂ρ∂t + ∇ · ( ρ v ) = 0 ,∂ ( ρ v ) ∂t + ∇ · ( ρ vv − BB + P ∗ ) = − S r ( P ) ,∂E∂t + ∇ · [( E + P ∗ ) v − B ( B · v )] = − S r ( E ) ,∂ B ∂t − ∇ × ( v × B ) = 0 ,∂I∂t + c n · ∇ I = cS I ,S I ≡ Γ − [ ρ ( κ s + κ a ) ( J − I )+ ρ ( κ a + κ δP ) (cid:18) a r T π − J (cid:19) ] ,S r ( E ) ≡ πc (cid:90) S I d Ω , S r ( P ) ≡ π (cid:90) n S I d Ω . (9)Here B is the magnetic field, P ∗ ≡ ( P + B / I (with I the unit tensor), and the magnetic permeability µ istaken to be 1 in this unit system. The total gas energydensity is E = E g + 12 ρv + B . (10)We assume ideal gas so that the gas internal energy isrelated to gas pressure via the adiabatic index γ g as E g = P/ ( γ g −
1) for γ g (cid:54) = 1. The gas temperature iscalculated with T = P/ ( R ideal ρ ), where R ideal is theideal gas constant. THE IMPLICIT SOLVER FOR THE RTEQUATION We will use the normal explicit Godunov scheme withconstrained transport in
Athena++ to evolve the MHDpart of equations 9. The time step is limited by the CFLcondition based on the maximum flow velocity, soundspeed and Alfv´en speed. For the RT equation, we willsolve all the terms implicitly to remove the time stepconstrain due to speed of light. We will describe how wediscretize and solve the RT equation in this section.3.1.
Spatial and Angular Discretization
The radiation field shares the same spatial grid asused for the gas quantities in
Athena++ . We label thethree coordinates as x, y, z , which can be any orthogo-nal curvilinear systems such as cartesian, cylindrical orspherical polar coordinates. The total number of gridpoints along the three directions are N x , N y , N z respec-tively. Notice that our finite volume scheme does notneed the grid spacing to be uniform. What we needis the volume of each cell V i,j,k for i ∈ [0 , N x − , j ∈ [0 , N y − , k ∈ [0 , N z −
1] as well as surface areas of eachcell facing the three directions A x , A y , A z . All the spe-cific intensities are defined at the volume center of eachcell.We adopt the same angular discretization scheme asused by Davis et al. (2012) and Jiang et al. (2014) to beour default choice . The angular unit vectors n are de-fined with respect to a global cartesian coordinate andthey do not change with spatial locations. The threecomponents of each angle n = ( µ x , µ y , µ z ), as well asthe quadrature weight w n are determined based on thesame algorithm developed by Bruls et al. (1999) follow-ing the original quadrature principle described in Carl-son (1963). This set of angles ensures that the followingrelations are satisfied numerically N − (cid:88) n =0 w n = 1 , N − (cid:88) n =0 µ x w n = N − (cid:88) n =0 µ y w n = N − (cid:88) n =0 µ z w n = 0 , N − (cid:88) n =0 µ x w n = N − (cid:88) n =0 µ y w n = N − (cid:88) n =0 µ z w n = 13 , (11)where N is the total number of angles.Notice that these angles are always defined in the labframe. The angles in the co-moving frame are deter-mined based on equation 3 while the quadrature weightin the co-moving frame is w (cid:48) n = Γ − w n . For infinitenumber of angles, it can be shown that (cid:82) Γ − d Ω = 1 Other angular discretization approaches can also be used for spe-cific applications as shown in Section 3.2.4 mplicit RT Algorithm w (cid:48) n as w (cid:48) n ≡ Γ − w n (cid:80) N − n =0 Γ − w n . (12)Therefore (cid:80) N − w (cid:48) n = 1 is always satisfied. We find itis necessary to correct this angular truncation error toavoid accumulative energy error when flow velocity isnot zero. The commonly used radiation energy density,radiation flux and radiation pressure in the lab and co-moving frames are just derived quantities from specificintensities as E r = 4 π N − (cid:88) n =0 I n w n , E r, = 4 π N − (cid:88) n =0 I ,n w (cid:48) n , F r = 4 πc N − (cid:88) n =0 I n n w n , F r, = 4 πc N − (cid:88) n =0 I ,n n (cid:48) w (cid:48) n , P r = 4 π N − (cid:88) n =0 I n nn w n , P r, = 4 π N − (cid:88) n =0 I ,n n (cid:48) n (cid:48) w (cid:48) n . (13)3.2. Discretized Equation for RT
The radiation energy source term can be very stiff inthe optically thick regions while the radiation momen-tum source term is typically not (Sekora & Stone 2010;Jiang et al. 2012). Therefore, when we solve the RTequation, we will assume the flow velocity and gas den-sity do not change. But we evolve the gas temperaturewith the radiation variables together to ensure that cor-rect thermal equilibrium state is achieved. As higherorder spatial reconstruction without flux limiter cannotguarantee total variation diminishing, which is necessaryto avoid oscillatory behavior around local maximums orminimums, we will only use first order spatial recon-struction for our algorithm. The commonly used fluxlimiters will make the transport term non-linear, whichis very hard to solve implicitly. The algorithm will alsobe first order accurate in time for RT to avoid oscillatorybehavior during the temporal evolution of gas tempera-ture and radiation energy density when the source termis very stiff.For given initial conditions I mn , ρ m , T m , v m at timestep m , the specific intensities I m +1 n with n ∈ [0 , N − m + 1 are calculated based on the following N + 1 equations I m +1 n − I mn + ∆ tc n · ∇ I m +1 n =∆ tc Γ − n [ ρ m κ s (cid:0) J m +10 − I m +10 ,n (cid:1) + ρ m κ a (cid:32) a r (cid:0) T m +1 (cid:1) π − I m +10 ,n (cid:33) + ρ m κ δP (cid:18) a r ( T m +1 ) π − J m +10 (cid:19) ] ,ρ m R ideal γ g − (cid:0) T m +1 − T m (cid:1) = − ∆ tcρ m ( κ a + κ δP ) (cid:2) a r ( T m +1 ) − πJ m +10 (cid:3) , (14)where I m +10 ,n = Γ n I m +1 n , J m +10 = N − (cid:88) n =0 w (cid:48) n I m +10 ,n , (15)and Γ n is calculated using v m for each angle. The opac-ities ( κ s , κ a and κ δP ) can be functions of local gas prop-erties in principle. But we will assume they keep thevalues at time step m in the above equations. This is aset of fully coupled equations for gas temperature andspecific intensities at all angles (via the terms with J m +10 and spatial locations (via the transport term n · ∇ I m +1 n ).These equations are only linear with respect to I m +1 n butnot for T m +1 due to the ( T m +1 ) term. All the trans-port and source terms are calculated using the variablesat time step m + 1 to make the scheme unconditionalstable for any time step ∆ t . We will use an iterativeapproach to solve these equations as described below.3.2.1. The Transport Term
Following Jiang et al. (2014), we can rewrite the trans-port term in the above equations as c ∇ · (cid:0) n I m +1 n (cid:1) withour default angular discretization because n does notchange with spatial locations. Then this term becomesthe conservative format and can be calculated withthe standard finite volume approach. For interface at i − / , j, k , we simply take the left and right states for I m +1 n to be I m +1 n ( i − , j, k ) and I m +1 n ( i, j, k ) respec-tively. The same approach is used for other directions.Special care is needed to determine the surface fluxesfor each specific intensity. If we simply take the up-wind fluxes as in Stone & Mihalas (1992), although totalenergy is conserved, numerical diffusion can overcomethe physical diffusion in the very optically thick regime.This can be easily demonstrated with only two anglesin one dimension. If we assume the two angles have µ x, > µ x, <
0, the upwind fluxes at i − / cµ x, I ( i −
1) and cµ x, I ( i ). Similarly, the fluxes at i + 1 / cµ x, I ( i ) and cµ x, I ( i + Jiang w I + w I at cell i is then determined by cw µ x, I ( i +1)+ cw µ x, I ( i ) − cw µ x, I ( i ) − cw µ x, I ( i − F r ( i + 1) − F r ( i −
1) = cw µ x, I ( i +1)+ cw µ x, I ( i +1) − cw µ x, I ( i − − cw µ x, I ( i − | I ( i + 1) − I ( i + 1) | can bemuch smaller than | I ( i + 1) − I ( i ) | . This will cause thenumerical flux due to the truncation errors to be muchlarger than the physical flux due to photon diffusion,which results in completely wrong numerical solutions.We have designed a modified HLLE flux to overcome theissue.The dynamical diffusion regime also needs to betreated carefully as discussed in Jiang et al. (2014). Thenumerical truncation error due to finite flow velocity canalso overcome the physical diffusion when the opticaldepth per cell is very large particularly with first orderspatial reconstruction. In principle, these numericallyerrors can be minimized with second or third order spa-tial reconstructions. However, it is not an option for ourimplicit scheme. Instead, we find a specially designedHLLE type flux can work very well in all the regimes.We rewrite the transport term as c ∇ · (cid:0) n I m +1 n (cid:1) = ∇ · (cid:2) ( c n − f v m ) I m +1 n (cid:3) + ∇ · (cid:0) f v m I m +1 n (cid:1) . (16)Here f is a function of optical depth per cell τ c ≡ ρ m ( κ a + κ s ) ∆ L with ∆ L to be a constant factor timesthe cell size. The function is chosen as f = 1 − exp( − τ c ) . (17)This is designed to separate the advection term whenthe photon mean free path is not resolved ( f ≈ f ≈ ∇ · (cid:0) f v m I m +1 (cid:1) explicitly, which means c ∇ · (cid:0) n I m +1 n (cid:1) ≈ ∇ · (cid:2) ( c n − f v m ) I m +1 n (cid:3) + ∇ · ( f v m I mn ) . (18)We use second order spatial reconstruction for I mn to getthe values at the cell faces and then upwind flux basedon f v m to calculate the term ∇ · ( f v m I mn ).For the interface at ( i − / , j, k ), we take the leftand right states for each specific intensity to be I n,L = I m +1 n ( i − , j, k ) and I n,R = I m +1 n ( i, j, k ). Then we cal-culate the flux F n ( i − /
2) at the interface for each angle based on a modified HLLE formula (Appendix A). Thekey modification is to make sure the flux is the upwindvalue in the optically thin case and it gets close to thecentered value in the optically thick case based on theparameter τ c . Similarly, we calculate the flux F n ( i +1 / i +1 / , j, k ) as well as fluxes along otherdirections F n ( j − / , F n ( j + 1 / , F n ( k − / , F n ( k +1 / ∇ · (cid:2) ( c n − f v m ) I m +1 n (cid:3) = 1 V i,j,k [ A x ( i + 1 / F n ( i + 1 / − A x ( i − / F n ( i − / A y ( j + 1 / F n ( j + 1 / − A y ( j − / F n ( j − / A z ( k + 1 / F n ( k + 1 / − A z ( k − / F n ( k − /
2) ] . (19)Reorganizing all the terms results in the following equa-tion that we will apply our iterative solver g I m +1 n + g I m +1 n ( i −
1) + g I m +1 n ( i + 1)+ g I m +1 n ( j −
1) + g I m +1 n ( j + 1)+ g I m +1 n ( k −
1) + g I m +1 n ( k + 1)= I mn − ∆ t ∇ · ( f v m I mn )+∆ tc Γ − n [ ρ m κ s (cid:0) J m +10 − I m +10 ,n (cid:1) + ρ m κ a (cid:32) a r (cid:0) T m +1 (cid:1) π − I m +10 ,n (cid:33) + ρ m κ δP (cid:18) a r ( T m +1 ) π − J m +10 (cid:19) ] ,ρ m R ideal γ g − (cid:0) T m +1 − T m (cid:1) = − ∆ tcρ m ( κ a + κ δP ) (cid:2) a r ( T m +1 ) − πJ m +10 (cid:3) . (20)Expressions for the coefficients g to g are given in theAppendix A. Notice that their values depend on the lo-cal optical depth τ c and they only need to be calculatedonce for each time step.3.2.2. The Iterative Solver
We need to solve the N + 1 equations (20) simulta-neously to get the solutions for I m +1 n ( n ∈ [0 , N − T m +1 . These equations are only linear for I m +1 n .These terms can be rewritten as a sparse matrix withsize of N × N x × N y × N z . If the specific intensities areordered in a 1D array such that I m +1 n ( i, j, k ) is locatedat kN y N x N + jN x N + iN + n , the matrix has the specialstructure that the diagonal direction are blocks with size N × N , which couple all the angles at the same spatiallocations. All the other elements are only non-zero at lo-cations corresponding to ( n, i ± , j, k ), ( n, i, j ± , k ) and mplicit RT Algorithm n, i, j, k ± T m +1 are non-linearbut T m +1 at different spatial locations are completelydecoupled.This motivates us to solve these equations iterativelysimilar to the Jacobi approach used for linear systems.For each cell, we take all the specific intensities at differ-ent spatial locations to be the values from the last step.Then we go through each cell and solve the equationsfor all the specific intensities at once. We update theoff-diagonal terms with the updated specific intensitiesand continue this process until changes of the specific in-tensities between neighboring steps are below a certainthreshold.Specifically, we solve the following equations itera-tively g I m +1 n,l + g I m +1 n,l − ( i −
1) + g I m +1 n,l − ( i + 1)+ g I m +1 n,l − ( j −
1) + g I m +1 n,l − ( j + 1)+ g I m +1 n,l − ( k −
1) + g I m +1 n,l − ( k + 1)= I mn − ∆ t ∇ · ( f v m I mn )+∆ tc Γ − n [ ρ m κ s (cid:16) J m +10 ,l − I m +10 ,n,l (cid:17) + ρ m κ a (cid:32) a r (cid:0) T m +1 l (cid:1) π − I m +10 ,n,l (cid:33) + ρ m κ δP (cid:18) a r ( T m +1 l ) π − J m +10 ,l (cid:19) ] ,ρ m R ideal γ g − (cid:0) T m +1 l − T m (cid:1) = − ∆ tcρ m ( κ a + κ δP ) (cid:104) a r ( T m +1 l ) − πJ m +10 ,l (cid:105) , (21)where l represents each step during the iterative process.At the beginning, we need an initial guess for specificintensities at the whole grid for all the angles. This isusually taken to be the solution from the last time stepor the initial condition. Then we go through all the gridpoints to solve the above equation for all the angles ateach cell together in the following way.For each step during the iteration, since all the vari-ables I m +1 n,l − ( i − , I m +1 n,l − ( i + 1) , I m +1 n,l − ( j − , I m +1 n,l − ( j +1)) , I m +1 n,l − ( k − , I m +1 n,l − ( k +1) are already known, we canreorganize all the equations we need to solve at each cellas g I m +1 n,l = I c + ∆ tc Γ − n [ ρ m κ s (cid:16) J m +10 ,l − I m +10 ,n,l (cid:17) + ρ m κ a (cid:32) a r (cid:0) T m +1 l (cid:1) π − I m +10 ,n,l (cid:33) + ρ m κ δP (cid:18) a r ( T m +1 l ) π − J m +10 ,l (cid:19) ] , (22) ρ m R ideal γ g − (cid:0) T m +1 l − T m (cid:1) = − ∆ tcρ m ( κ a + κ δP ) (cid:104) a r ( T m +1 l ) − πJ m +10 ,l (cid:105) , (23)where I c includes all the other terms that are alreadyknown before step l . We multiply the left and rightsides of the equation for each I m +1 n,l with Γ n to get theequation for co-moving frame specific intensities as g I m +10 ,n,l = Γ n I c + ∆ tc Γ n [ ρ m κ s (cid:16) J m +10 ,l − I m +10 ,n,l (cid:17) + ρ m κ a (cid:32) a r (cid:0) T m +1 l (cid:1) π − I m +10 ,n,l (cid:33) + ρ m κ δP (cid:18) a r ( T m +1 l ) π − J m +10 ,l (cid:19) ] . (24)We sum up the equations with the appropriate weight w (cid:48) n for each angle and get g J m +10 ,l = (cid:88) n w (cid:48) n Γ n I c + ∆ tc Γ n [+ ρ m κ a (cid:32) a r (cid:0) T m +1 l (cid:1) π − J m +10 ,l (cid:33) + ρ m κ δP (cid:18) a r ( T m +1 l ) π − J m +10 ,l (cid:19) ] . (25)The scattering term drops out because it does notchange the mean photon energy. Together with equa-tion 23, we can get the solution for J m +10 ,l and T m +1 l byperforming a root finding for a fourth order polynomial.Then specific intensity for each angle I m +10 ,n,l can then becalculated from equation 24 independently, which canbe converted to I m +1 n,l directly. In this way, the non-linearity of the source term, as well as the coupling be-tween different angles, all come down to a four orderpolynomial solver, which is trivial and independent ofthe number of angles. Therefore, the whole cost for eachstep during the iteration is linearly proportional to thenumber of angles. We continue the process for each gridand stop the iteration when∆ I ≡ (cid:88) n,i,j,k | I m +1 n,l − I m +1 n,l − | / (cid:88) n,i,j,k | I m +1 n,l | (26)is smaller than a threshold, typically 10 − .3.2.3. The Radiation Source Terms
Since total energy and momentum are conserved in thelab frame, we use the conservation laws to determine thesource terms we need to add to the gas. The radiationenergy and momentum at the beginning of the time stepis E mr = 4 π N − (cid:88) n =0 I mn w n , F mr = 4 πc N − (cid:88) n =0 n I mn w n . (27) Jiang
At the end of the time step, they can be calculated as E m +1 r = 4 π N − (cid:88) n =0 I m +1 n w n , F m +1 r = 4 πc N − (cid:88) n =0 n I m +1 n w n . (28)Part of the change is due to the transport term, whichcan also be calculated as∆ E r = 4 π N − (cid:88) n =0 { ∆ t ∇ · (cid:2) ( c n − f v m ) I m +1 n (cid:3) + ∆ t ∇ · ( f v m I mn ) } w n , ∆ F r = 4 πc N − (cid:88) n =0 { ∆ t ∇ · (cid:2) ( c n − f v m ) I m +1 n (cid:3) + ∆ t ∇ · ( f v m I mn ) } n w n . (29)The source terms are then given by S r ( E ) = E m +1 r − E mr − ∆ E r S r ( P ) = 1 c (cid:0) F m +1 r − F mr − ∆ F r (cid:1) . (30)3.2.4. Additional Terms with Spherical Polar AngularSystem
The choice of the angular system is not unique. Ifthe angles vary with spatial locations, additional termswill show up when we convert n · ∇ I to ∇ · ( n I ). TheRT equation with general angular systems is describedin Davis & Gammie (2020). One special case that isuseful for non-relativistic system is the angular systemin spherical polar coordinate, which allows us to usethis scheme in 1D or 2D spherical polar coordinates.Notice that our default angular system as described inSection 3.1 can also be used in principle. The alternativeangular system we describe here can be more convenientand efficient for certain applications.In spherical polar coordinate, we define a set of angleswith respect to the local coordinate ( r, θ, φ ) in each cell.We divide the longitudinal angles ψ uniformly from 0 to2 π . For the polar angle ζ , it is uniformly divided in cos ζ for ζ ∈ [0 , π ]. These angles point to different directionsat different spatial locations, which introduce additionalterms in the transport term as (Davis & Gammie 2020) ∂I∂t + c n · ∇ I = ∂I∂t + c ∇ · ( n I ) − r sin ζ ∂ (cid:0) sin ζI (cid:1) ∂ζ − sin ζ cos θr sin θ ∂ (sin ψI ) ∂ψ . (31)All the other terms are unchanged. The numericalscheme described in the previous sections remains the same except that we need to calculate the change of I due to the new terms.The new terms represent transport of specific intensi-ties in the angular space, which we can calculate implic-itly as I m +1 n − I mn ∆ t − r sin ζ n +1 / I m +1 n +1 / − sin ζ n − / I m +1 n − / cos ζ n − / − cos ζ n +1 / − sin ζ n cos θr sin θ sin ψ n +1 / I m +1 n +1 / − sin ψ n − / I m +1 n − / ψ n +1 / − ψ n − / = 0 . (32)Here the variables with subscripts n − / n + 1 / I m +1 n − / and I m +1 n +1 / .During each iteration, we modify the coefficient g inequation 21 to include additional contributions for I m +1 n,l from this term. For all the other specific intensities withdifferent angles in the above equation, we use valuesfrom the last iteration, which are added to the righthand side of equation 21. We find it is necessary toinclude all the terms during the iterative process in orderto maintain a balance between all the transport termsand source terms.3.3. Algorithm for the Full Radiation MHD Equation
The implicit solver is coupled to the two step van-Leer MHD solver in
Athena++ (Stone et al. 2020). Wewill not repeat details of the Godunov scheme with con-strained transport for MHD here. Instead, we just out-line the steps of the full algorithm here.(a) Determine the time step ∆ t based on the flowspeed, sound speed, Alfv´en speed and cell sizeswith CFL number 0 . t/ Athena++ .(c) Evolve the radiation field for half time step ∆ t/ S r ( E ) , S r ( P )(equation 30). Add these terms to gas momentumand total energy.(d) Apply the MHD solver for a full time step ∆ t usingthe partially updated variables at ∆ t/ t and add the energy as well asmomentum source terms to the gas. mplicit RT Algorithm t/t E r , a r T E r a r T . . . . . . t/t E r , a r T E r a r T Figure 1.
Histories of E r and a r T when gas and radiationrelax to the thermal equilibrium state from static and uni-form initial conditions. The left and right panels show caseswith different relative values between a r T and E r . The timeunit is t = 1 / ( cρκ a ) in the two plots. The dotted blue linesindicate the analytic solution during the thermal equilibriumstate. 4. TESTS OF THE ALGORITHMThe algorithm we developed here can in principle workin a wide range of parameter space in terms of opticaldepth and ratio between radiation pressure and gas pres-sure, since no approximation is made when we solve thetransport equation. When the time step is made suffi-ciently small, the algorithm can also be reduced to theexplicit scheme as described by Jiang et al. (2014). How-ever, in practice, convergent properties for our iterativescheme can be very problem dependent. We will repeatmost of the tests done by Jiang et al. (2012) and Jianget al. (2014) to explore the accuracy and efficiency ofour algorithm. Unless otherwise specified, we will as-sume the temperature, density and length units to be T , ρ , L respectively. The velocity unit v is chosento be the isothermal sound speed corresponding to thetemperature T and the time unit is given by L /v .Radiation energy density is scaled to a r T and radia-tion flux is scaled to ca r T . Solutions given in the di-mensionless format can be converted to any unit systemgiven the two dimensionless parameters C ≡ c/v and P = a r T / ( ρ R ideal T ) (Jiang et al. 2012).4.1. Thermal Equilibrium in a uniform medium
To test our treatment of the source term, we setupa uniform box in 2D covering the region [0 , × [0 , ρ = 1 andgas temperature T = 1. The radiation field is initializedisotropically with the mean energy density E r = 100. .
00 0 .
05 0 .
10 0 .
15 0 .
20 0 . πµ x I − . − . − . . . . . π µ y I Figure 2.
Angular distribution of the lab frame specificintensities in the steady state for a moving uniform gas asdescribed in Section 4.1. The black and blue lines connect theorigin and the points (4 πµ x I, πµ y I ) for the twelve angles wehave in the numerical solution. All the angles representedby the blue lines have µ z = 0 .
333 while the other anglesrepresented by the black lines have µ z = 0 . µ x , µ y ) with the same corresponding µ z respectively. We fix the two dimensionless parameters P and C to be1 and 100 respectively. We use 32 grid points for bothdirection. A constant absorption opacity ρκ a = 100 isadopted for the test. The typical thermalization timescale is then 1 / ( C ρκ a ) = 10 − while the time step weuse is roughly 0 . E r − a r T so that thereis no oscillatory behavior as a function of time, whichcan happen for high order time integrators when thetime step is much larger than the thermalization timescale. As another test, we change the initial conditionto be T = 100 , E r = 1 and ρκ a = 1 so that the gas ismuch hotter than the equilibrium state in this case. Wecan still get the correct equilibrium state even thoughthere is a huge diffrence between a r T and E r initiallyas shown in the right panel of Figure 1.Our treatment of the source term is very similar tothe scheme as described in Jiang et al. (2014), where theflux divergence term is calculated explicitly. However,there are non-trivial differences for the term n · ∇ I even0 Jiang
Figure 3.
Spatial distribution of mean radiation energydensity E r after two beams of photons propagate through the2D box with zero opacity. The two beams are injected from x = ± . , y = − µ x = ± . , µ y = 0 . for the case with a uniform spatial distribution. Whenthis term is calculated using the variables from the lasttime step, it is guaranteed to be zero to the level ofround-off error as there is no spatial gradient at eachstep. This is not the case for the fully implicit scheme.During each cycle of the iteration when the gradientis calculated for each cell, variables in the neighboringcells are taken from last iterative step while we use thevalues from this step for the current cell. Before thesystem reaches the equilibrium state, the flux divergenceterm is actually not zero. After the solution converges,we recover the same solution as given by the explicitscheme with the accuracy set by the tolerance level. TheJacobi like iterative scheme we adopt also guaranteesthat all the cells are treated in the same way during eachiteration. Therefore, we do not generate any artificialspatial gradient during the iteration.We have also tested the velocity dependent terms bygiving the gas a uniform horizontal velocity v x = 3.Specific intensities are initialized isotropically in the labframe with the mean energy density E r = 1. The initialgas temperature and density are T = 1 , ρ = 1 respec- tively. We choose ρκ a = 1 , P = 1 and C = 100 for thistest and use 12 angles per cell. The system settles downto a steady state with gas velocity reduced to v x = 2 . F r,x is very close to v x ( E r + P r,xx ), which isthe value given by the first order expansion in terms of v/c in radiation moment equations (Lowrie et al. 1999)since v x / C is only 3% in this test. During this pro-cess, the change of kinetic energy is converted to ra-diation energy and lab frame E r = 1 .
13 in the steadystate. The specific intensities can also be compared withthe analytical values directly, which are shown in Fig-ure 2. In steady state, the radiation field is isotropicin the co-moving frame. Angular distribution of spe-cific intensities in the lab frame is then determined bythe Lorentz transformation (equation 2). The numericalsolution agrees with the analytical values very well.Notice that the radiation field is not isotropic in thelab frame and the Eddington tensor is P r,xx /E r = 0 . / I in the lab frame.4.2. Crossing Beams in Vacuum
Our iterative scheme is typically easier to convergewhen the source term dominates over the transportterm, because this means the diagonal coefficients in theresulting matrix of equation 20 are larger than the off-diagonal one. The most difficult scenario is the trans-port of photons through a medium with zero opacity,where all the source terms are zero and we only iter-ative over the transport term. To show how the al-gorithm works in this case, we inject two beams ofphotons from the bottom of a box covering the region[ − . , . × [ − , x = ± . , y = −
2, the specificintensities with µ x = ± . , µ y = 0 .
577 are set to be0 . x direction and outflowat the top. Opacity coefficients of the gas are set to 0and the gas does not evolve in this test. The spatialresolution is 64 ×
256 for the horizontal and vertical di-rections and we use 4 angles per cell. The time step is4 × − , which is also the light cross time across the ver-tical direction of the box. During the first step, it takes100 iterations to reach the relative error 10 − , which de-clines very slowly with more iterations. After one step,the two beams have crossed the whole simulation box.In about 10 time steps, the system reaches steady state mplicit RT Algorithm E r is shown in the leftpanel of Figure 3. Compared with the same test done byJiang et al. (2014), the solution is very similar to the casewith first order spatial reconstruction but more diffusivecompared with the solution using second order spatialreconstruction. This is expected as our implicit algo-rithm only uses first order spatial reconstruction. Theoverall cost of the implicit algorithm is also comparableto the cost if the transport equation is solved explicitlyfor this case, because the number of iterations we needfor convergence is now comparable to the ratio betweenspeed of light and sound speed. Another reason that thistest requires a lot of iterations for convergence is thatthe initial condition is far away from the steady statesolution. The injected photons at the bottom boundaryneed to be communicated with the whole simulation boxduring the iterations. However, multigrid technique cansignificantly speed up the convergence in this case be-cause the information at the boundary can propagationto the domain much faster at lower resolution levels.Our algorithm is a finite volume scheme, which evolvesthe specific intensities using fluxes at cell faces as othercell centered variables such as density in Athena++ .Therefore, the existing mesh refinement module in thecode (Stone et al. 2020) can be directly used for our RTalgorithm. One example is shown in the right panel ofFigure 3. The setup is the same as the left panel exceptthat we refine the region [ − . , . × [ − ,
1] with twolevels of refinements. Photons can cross the refinementboundaries smoothly without any numerical artifact.4.3.
Shadow Test
One widely used test to demonstrate the differencebetween solving the angular resolved RT and radiationmoment equations based on simplified assumptions ofthe closure relation is to capture the shadow cast byan optically thick cloud with beams of photons(Hayes& Norman 2003; Gonz´alez et al. 2007). Diffusion likeapproximation will always send photons to the regionwhere radiation energy density gradient exists and thuscannot keep the shadow. Other approximations thattreat the radiation field as a fluid (such as M1, Gonz´alezet al. 2007; McKinney et al. 2013) have their own caveatssuch as artifically merging photons coming from differentdirections. Our implicit scheme should get the correctshadows as we still solve the specific intensities alongdifferent directions independently.We used the same setup as in Jiang et al.(2012) and Jiang et al. (2014) for direct compar-ison. The simulation domain covers the cartesianbox [ − . , . × [ − . , .
3] with resolution 512 × . . . . − . − . . . . . − . − . − . − . . . . . . . . − . − . − . E r . . . . . . . . f xx . . . . . . . . f yy Figure 4.
The top panel shows the distribution of radiationenergy density E r after two beams of photons pass an opti-cally thick cloud, which is located at the center of the boxas described in section 4.3. The middle and bottom panelsshow the x − x and y − y components of Eddington tensor f xx and f yy respectively. tion. The gas density has the distribution ρ ( x, y ) =1 + 9 / (cid:2) (cid:2) (cid:0) ( x/ . + ( y/ . − (cid:1)(cid:3)(cid:3) whilethe temperature T is set to achieve a constant gas pres-sure P = 1 for the whole box. We use pure absorp-tion opacity κ a = ρT − . so that the total optical depthacross the box in the hot gas around the cloud is only1, while optical depth across the cloud is more than10 . For the left boundary, specific intensities alongthe directions ± ◦ with respect to the horizontal axishave the fixed value 103 .
13 but all the others are set0. For the right boundary, outgoing specific intensitiesare just copied from the last activate zones to the ghostzones while incoming specific intensities are 0. Peri-odic boundary condition is used for the up and bottomboundaries. We choose the dimensionless speed of lightto be C = 1 . × , which does not affect the steadystate solution. All specific intensities are initialized tobe 0 and we keep the gas quantities fixed for the shadowtest. It takes about 1000 iterations to reach a relative2 Jiang . . . . − . − . . . . . − . − . − . − . . . . . . . . − . − . − . ρ Figure 5.
Density evolution of an optically thick cloud ir-radiated by two beams of photons from the left boundary.The three snapshots are shown at times 0 . t , . t , . t from top to bottom respectively, where t is the sound cross-ing time across the whole box horizontally in the low densitygas. error of 10 − with three time steps and the radiationquantities are very close to the steady state values acrossthe whole simulation box. The spatial distributions ofmean radiation energy density E r and components ofthe Eddington tensor f xx and f yy are shown in Figure4. Our implicit algorithm gets the umbra and penumbracorrectly and the results are very similar to the solutionsreturned by the explicit scheme as shown in Jiang et al.(2014) as well as the VET method as described by Jianget al. (2012).One great advantage of our implicit algorithm is thatit can also capture the time dependent evolution of theradiation field efficiently even though the speed of lightis 1 . × times the typical sound speed. To demon-strate that our algorithm can also capture the full dy-namical evolution of the cloud, we use the same setup asin the shadow test but evolve the hydrodynamic quan-tities with the radiation field together. We choose thetemperature scale such that P = 10 − so that the ra- diation pressure of the injected photons from the leftboundary is about 10 − of the initial gas pressure. No-tice that the static shadow test is independent of thisratio. The parameter is also chosen to match the sameexperiment done by Jiang et al. (2012). Outflow bound-ary condition is used for the hydrodynamic variables atthe left and right boundaries while periodic boundary isused for the top and bottom. The typical sound crossingtime scale across the box is t = 1 while the light cross-ing time in the low density gas is ≈ . × − t . Thephoton diffusion time scale across the cloud is compara-ble to t . Since the effective temperature of the incom-ing radiation field is larger than the temperature of thecold gas by a factor of 30, the gas will be heated up inthe time scale of 2 × − t . Density distributions of thecloud at three different times t = 0 . t , . t , . t areshown in Figure 5. The cloud is heated up at the surfaceof the left side initially. The heated gas drives shocksinto the cloud and some gas is also evaporated simulta-neously. The shock propagates inside the cloud and itgets destructed eventually. The overall evolution is verysimilar to the results shown by Jiang et al. (2012) andProga et al. (2014), despite some small differences of thedetailed structures. Compared with the VET scheme,there is a subtle difference on effectively what equationsare solved. The VET scheme typically calculates the Ed-dington tensor at the beginning of each time step, whichis then held fixed when the radiation moment equationsand hydrodynamic equations are evolved for one timestep. The algorithm developed here solves the specificintensities implicitly with the same time evolution asthe hydrodynamic variables. If we integrate the specificintensities over the angular space, effectively we evolvethe radiation moments using the Eddington tensor de-termined at the end of each time step consistent withthe implicit update of the radiation field. This will notmake any difference for steady state solutions or if thetime step is sufficiently small. However, when the timestep is much larger than the light crossing time for eachcell, change of the Eddington tensor due to time evolu-tion of the radiation field within one time step may notbe negligible. Future studies are needed to quantify thedifference.4.4. Static and Dynamic Diffusion
Photon diffusion in the optically thick regime, which iseasy to solve with the classical diffusion approximationby design, turns out to require careful treatment whenthe full time dependent transport equation is solved.This is because in the very optically thick regime, thenumerical diffusion associated with the standard HLLEsolver can be much larger than the physical diffusive flux mplicit RT Algorithm − . − . . . . x − − − − − E r √ Dt = 0 . √ Dt = 0 . √ Dt = 0 . Figure 6.
Profiles of radiation energy density E r at threedifferent times corresponding to √ Dt = 0 . , . , .
22 froman initial Gaussian profile as indicated by the dotted blackline. The diffusion coefficient is D = C / (3 ρκ s ) = 8 . × − for this test. The solid lines are numerical solutions while thedashed lines are analytical solutions to the diffusion equation. (see the Appendix of Jiang et al. 2013). To demonstratethat our treatment of the transport term as described inSection 3.2.1 is sufficient to capture the photon diffusioncorrectly, we set up a radiation field with mean energydensity E r ( x ) = exp (cid:0) − x (cid:1) in the box x ∈ ( − . , . | x | > .
5, we set E r = exp ( − − ,
1) and 2 an-gles per cell so that the Eddington tensor is always 1 / ρκ s = 4 × andwe do not evolve the gas quantities for this test. Wechoose the dimensionless speed of light to be C = 10.Evolution of E r can be described by the solution to thediffusion equation as (Sekora & Stone 2010; Jiang et al.2014) E r ( x, t ) = 1(160 Dt + 1) / exp (cid:18) − x Dt + 1 (cid:19) , (33)where D ≡ C / (3 ρκ s ) is the diffusion coefficient. Thespecific intensities are initialized isotropically with themean energy density matching the initial profile of E r .Profiles of E r at three different times from the numerical − − x − − − − E r t = 4 t = 8 t = 16 Figure 7.
Diffusion of photons in a moving background asdescribed in section 4.4. The gas has a constant scatter-ing opacity ρκ s = 4 × , velocity v = 1 while the speedof light is 10 in this test. The colored solid lines are nu-merical solutions at times t = 4 , ,
16 respectively, while thecorresponding analytical solutions are indicated by the blackdashed lines. The dotted line at x = 0 is the initial profile. solution are shown in Figure 6, which match the analyti-cal solution very well. Notice that because we do not setthe boundary condition to match the diffusion solutionexactly, we can only compare the region roughly √ Dt away from the boundary. We have also tried cases withlarger opacity and our algorithm can still get the solu-tion correctly. However, we need to use a larger value ofthe parameter τ c (equation 17) to get the same accuracywith a larger κ s .Another important regime is dynamic diffusion, wherethe photons are advected with the flow while diffusionhappens. We use a similar setup as for the static diffu-sion case but extend the simulation domain to ( − , v = 1. The dimen-sionless speed of light is also increased to C = 1000 forthis test so that the typical diffusion speed is not neg-ligible compared with the advection speed. Profiles ofthe mean radiation energy densities at time t = 4 , , Jiang solutions can also be estimated by replacing the position x in equation 33 with x − vt to the first order of v/c ,which is 10 − in this test. The analytical solutions areshown as the black dashed lines in Figure 7 and theymatch the numerical solution very well. − − − τ − − − − E r (cid:15) = 10 − (cid:15) = 10 − (cid:15) = 0 . (cid:15) = 10 − (cid:15) = 10 − Figure 8.
Profiles of radiation energy density as a func-tion of optical depth for the atmosphere test as describedin Section 4.5. The solid lines are numerical solutions whilethe dashed black lines are corresponding analytical solutions.The ratio between absorption and total opacity (cid:15) varies from0 . − as indicated for each group of lines in the plot. Non-LTE Atmosphere
The test of a uniform temperature, scattering domi-nated atmosphere done by Davis et al. (2012) is useful toshow how our iterative scheme performs when there is alarge dynamical range of radiation energy density withboth optically thin and optically thick regions. Partic-ularly, we can compare the performance of our iterativescheme with the classical short characteristic methodfor the problem that the classical method is efficient tosolve. We use 1280 grid points covering a simulationdomain x ∈ [ − ,
10] in 1D. Gas density follows the ex-ponential profile ρ = 10 − exp(10 − x ) and gas tempera-ture is fixed to be T i = 1. The ratio between absorptionand total opacity is (cid:15) , which is taken to be a constant inthe whole simulation domain. Therefore, the scattering Iterations − − ∆ I / I c/v = 10 c/v = 10 c/v = 100 c/v = 10 c/v = 10 Figure 9.
Change of the relative error for specific intensitiesat x = 10 of the atmosphere test described in Section 4.5 as afunction of the total number of iterations. This is measuredfor the case with destruction fraction (cid:15) = 10 − . Differentlines represent cases with three different values of speed oflight we used compared with the fiducial sound speed as in-dicated in the legend. and absorption opacity are κ a = (cid:15) and κ s = 1 − (cid:15) . Thegas quantities are fixed to be the initial conditions withzero velocity and does not evolve. Since our algorithmalways solves the gas internal energy and radiation fieldtogether, we simply reset the gas quantities to be theinitial profiles at the end of each step. We use one angleper octant for this test so that the Eddington tensor isalways 1 /
3. At the bottom boundary, we set the specificintensities to be the isotropic value a r T i / (4 π ). At thetop boundary, the incoming specific intensity is set tobe 0 while the outgoing specific intensity is copied fromthe last active zone to the ghost zone. Specific inten-sities inside the simulation domain is initialized to be a r T i / (4 π ).The steady state profile of radiation energy density forthis atmosphere setup is E r = a r T i (cid:34) − exp (cid:0) −√ (cid:15)τ (cid:1) √ (cid:15) (cid:35) , (34) mplicit RT Algorithm τ ≡ (cid:82) x ρdx is the total optical depth integratedfrom x = 10. Upward and downward specific intensi-ties are given by I + = (1 / π ) [ E r + (1 / µ x ) ∂E r /∂τ ] and I − = (1 / π ) [ E r − (1 / µ x ) ∂E r /∂τ ] respectively. Thesteady state numerical solutions for (cid:15) varying from 0 . − are shown as the solid lines in Figure 8, whichagree with the analytical solutions very well. Our algo-rithm requires more iterations to reach the same level ofaccuracy when (cid:15) decreases. The relative error is alsolarger in the optically thin surface and the error in-creases with reducing (cid:15) for a fixed spatial resolution. Allthese properties are very similar to results given by theshort characteristic method for this test problem (Daviset al. 2012). We have also tried the test in a 3D domainwith the atmosphere profile aligned with the x axis. Weuse periodic boundary condition for y and z directionsin this case. The steady state solutions are identical tothe 1D case and the rate of convergence is also the same.The overall cost is just proportional to the total numberof grid points.To quantify the rate of convergence for our algorithmand compare with the short characteristic method, wepick the (cid:15) = 10 − case and measure the relative error∆ I/I of the specific intensities at x = 10 after each it-eration. Here we calculate ∆ I as the difference betweenthe solution at each iterative step and the analytical so-lution. We only focus on the location at x = 10 becausethe difference between the initial condition and the an-alytical solution is the largest there. It is also the mostoptically thin region and takes more iterations to con-verge. The criterion is different from what Davis et al.(2012) used because we do not iterate over the sourceterm. But both criteria serve the same purpose. Noticethat the steady state solution does not depend on thespeed of light. However, our algorithm always solvesthe time dependent transport equation and relaxes tothe steady state solution with the time step ∆ t fixed bythe gas sound speed v . We have tried three differentvalues of c/v and the changes of ∆ I/I with number ofiterations are shown in Figure 9. When c/v is small,although the implicit scheme is easier to converge pertime step, it actually takes more total number of iter-ations to get the steady state solution. When c/v islarger than 10 so that c ∆ t is larger than the lengthof the simulation box, the convergence rate is indepen-dent of c/v anymore. Our implicit algorithm also takesmore iterations to converge with this initial conditioncompared with the short characteristic method. Thistest demonstrates that the implicit scheme is not opti-mal just to find steady state solution as we are solvingthe full time dependent RT, although we can still getthe correct answer. 4.6. Homogeneous Sphere
Solving the RT equations in 1D spherical polar coor-dinate can be useful for many systems that are close tospherical symmetric (such as stars), particularly with fu-ture development of multiple frequency groups. Spher-ical polar angular systems as described in section 3.2.4need to be used in this case because of the assumedsymmetry. One useful test for this case is the homoge-neous sphere problem, which is widely used to test thenumerical scheme for radiation and neutrino transport(Abdikamalov et al. 2012; Radice et al. 2013). For thistest, we setup the gas with a constant density ρ = 1and temperature T = 1 from radius 0.05 to R = 1. Be-tween radius 1 and the outer boundary, which is chosento be 7, the density is 10 − and temperature is 10 − . Weassume absorption opacity ρ κ a = 10 for radius smallerthan 1 and the opacity is 0 in other region. The gas isheld with zero velocity and does not evolve. The unitscan be chosen to match any system and do not affectthe solution to the transport equation. We set the spe-cific intensity to be isotropic with the value a r T / (4 π )at the inner boundary. For the outer boundary, we re-quire the outgoing specific intensities to be continuouswhile the incoming specific intensities to be 0. We use1000 uniformly distributed grid points for the spatialresolution and 40 angles per cell for the specific inten-sities. This setup covers both the optically thick and asharp transition to the optically thin region, which canbe challenging if the transport equation is not solved ac-curately. The specific intensities are initialized to be theisotropic value a r T / (4 π ) in the whole grid. We choosethe speed of light to be 100 times the gas sound speedand the code finds the steady state solution in roughly20 iterations. Notice that the steady state solution doesnot depend on the speed of light.This setup also has an analytical solution, which wecan compare with. The specific intensities as a functionof radius r and angles n are (Smit et al. 1997) I ( r, n ) = a r T π [1 − exp ( ρ κ a s ( r, n r ))] , (35)where n r is the radial component of the unit vector n and the function s ( r, n r ) is defined as s ( r, n r ) = rn r + R g ( r, n r ) if r < R R g ( r, n r ) if r ≥ R & (cid:113) − ( R /r ) ≤ n r ≤
10 otherwise , (36)6 Jiang with the function g ( r, n r ) to be g ( r, n r ) = (cid:115) − (cid:18) rR (cid:19) (1 − n r ) . (37)Moments of the radiation field ( E r , F r , P r ) can be cal-culated directly via angular quadrature of I ( r, n ).Radial profiles of E r , F r and P r /E r from the numer-ical solution are shown as the solid lines in Figure 10,which agree with the analytical solutions very well. Par-ticularly, our numerical algorithm can capture the sharptransition at r = 1 without showing any artificial nu-merical effect. The Eddington tensor P r /E r is 1 / / r = 1. At large distance, the sphere becomes a pointsource effectively and only outgoing specific intensitiespointing radially are non-zero. Therefore, the Edding-ton tensor approaches 1. The slight difference betweenthe numerical and analytical solutions for P r /E r is dueto the finite angular resolution to resolve the angulardomain (cid:113) − ( R /r ) ≤ n r ≤ r . − − E r / a r T , P r / E r E r /a r T P r /E r r − − − − F r / c a r T Figure 10.
Profiles of radiation energy density E r , Ed-dington tensor P r /E r (the top panel) and radiation flux F r (the bottom panel) from the test of homogeneous sphere asdescribed in section 4.6. The solid lines are from the numer-ical solution while the analytical solutions are given by thedashed lines. Figure 11.
Convergence of the L1 error (cid:15) for the linear wavetests as a function of spatial resolution N . The dots withdifferent colors are for different optical depth per wavelengthas labeled in the plot. The left and right panels are forbackground states with P = 0 .
01 and 10 respectively. . . . . R e ( ω ) / k a r T /P = 0 . a r T /P = 1 a r T /P = 1010 − − τ a − − − − I m ( ω ) Figure 12.
Comparison of the propagation speed (the toppanel) and damping rate (the bottom panel) for the radia-tion modified linear waves between the numerical calculatedvalues (the dots) and analytical solutions (the solid lines).The lines with black, blue and red colors are for backgroundstates with P = 0 . , mplicit RT Algorithm Convergence of Linear Waves
Properties of convergence for the full radiation MHDalgorithm can be checked with linear wave tests. Thelinearized equations we solve here are identical to whatare evolved by Jiang et al. (2014). In order to comparedirectly, we adopt the same setup here. The backgroundgas has uniform density ρ = 1, and temperature T = 1with adiabatic index γ = 5 /
3. We use one angle peroctant with mean radiation energy density E r = 1 sothat the Eddington tensor is always 1 / I as used in theanalytical solutions. We fix the dimensionless speed oflight to be C = 10 and only use a constant pure ab-sorption opacity. We use a 2D simulation box with size L x = 1 , L y = 1 and change the spatial resolution from32 ×
32 to 512 × (cid:15) be-tween the numerical solutions and the analytical solu-tions over the whole simulation box for each resolutionand check how the errors change with resolution. Wehave also tried the same test in 3D and get the sameresults, although these 3D tests are more expensive torun. We initialize the calculation with eigenmodes ofthe radiation modified sound waves as given by the Ap-pendix B of Jiang et al. (2012) and stop the calculationafter the waves propagate one period along x direction.The two key dimensionless parameters that control theproperties of radiation modifieid sound waves are opticaldepth per wavelength length τ a and the ratio betweenradiation pressure and gas pressure P in the backgroundstate after C is fixed. Figure 11 shows the change of (cid:15) with resolution N for three different values τ a and twodifferent values of P . When the error is dominated bythe hydrodynamic part as in the optically thin regime,small radiation pressure or low resolution, the error de-creases with increasing resolution as N . When the er-ror is dominated by the RT module, it only shows firstorder convergence as expected. Part of the truncationerror in RT is due to different upwind directions for spe-cific intensities propagating along opposite directions asexplained in section 3.2.1. For the case with the highestresolution N = 512, we find that increasing τ c (equation17) in the HLLE flux by a factor of 2 from our defaultvalue can help decrease the L1 error.We can also measure the numerically calculated prop-agation speed and damping rate of the linear waves withthe resolution N = 256, which are shown in Figure 12 for τ a varying from 10 − to 10 with three different valuesof P in each case. When the optical depth is small, ouralgorithm is able to capture both the adiabatic soundwave with P = 0 .
01 and isothermal sound wave with P = 10. When τ a is larger than 10 so that radiation andgas are tightly coupled, we are also able to capture the radiation driven acoustic modes correctly. The resultsare the same as shown in Figure 10 of Jiang et al. (2014)as the same parameters are used. But the time step inthis algorithm is a factor of 10 larger compared with thecase when the time step is limited by the speed of light.4.8. Radiation Shocks
Another challenging test for the full algorithm is theradiation modified shocks, which check the numericalscheme in the non-linear regime. Particularly, we cansee how good the source terms balance the transportterms for both the RT and hydrodynamic part, which iscrucial to get the structures of the shock correctly, eventhough the time step is much larger than the light cross-ing time per cell. Since our algorithm has the full angu-lar distribution of the radiation field, we will comparewith the semi-analytical solutions provided by Fergusonet al. (2017), which gets the radiation shock solutionsby solving the time independent RT and hydrodynamicequations numerically without making any assumptionon the closure relation. These are the same radiationshock solutions that are used to test the explicit RTalgorithm in Jiang et al. (2014). . . . . T r . . . . T − − − × F r . . . ρ − . − .
025 0 .
000 0 .
025 0 . x . . . . f xx − . − .
025 0 .
000 0 .
025 0 . x . . . v Figure 13.
Structures of the sub-critical radiation modi-fied shock with upstream Mach number M = 2. The panelsshow steady state profiles of radiation temperature T r , gastemperature T , radiation flux F r , density ρ , horizontal com-ponent of Eddington tensor f xx and flow velocity v . Theblack lines are from our numerical solution in steady statewhile the dashed red lines are semi-analytical solutions asdescribed in section 4.8. We choose a unit system such that the upstream adi-abatic sound speed is 1 and the dimensionless speed oflight C = 1 . × . The upstream gas has densityand temperature ρ = 1 and T = 0 . Jiang T r T − − F r ρ − . − .
025 0 .
000 0 .
025 0 . x . . . f xx − . − .
025 0 .
000 0 .
025 0 . x v Figure 14.
The same as Figure 13 but for the case of super-critical shock with Mach number M = 3. dex γ = 5 /
3. We change upstream flow velocity for caseswith different Mach numbers. The dimensionless param-eter P is 7 . × − , which is roughly the ratio betweenradiation pressure and gas pressure in the upstream. Wechoose a constant absorption opacity ρκ a = 577 . a r T . For the right boundary,we simply copy the gas quantities from the last activezone to the ghost zones but we find it is necessary tokeep the gradients of specific intensities to be continu-ous across the boundary.Profiles of various quantities after the numerical so-lution reaches steady state are shown in Figure 13 forthe case with upstream flow velocity v = 2. This isthe case with a sub-critical shock (Mihalas & Mihalas1984; Sincell et al. 1999; Lowrie & Edwards 2008), wherethe gas temperature has a Zel’dovich spike (Zel’Dovich& Raizer 1967) while the radiation temperature is con-tinuous. For the sub-critical case, the gas temperaturebefore the spike is smaller than the downstream temper-ature. Particularly, the Eddington tensor is smaller than1 / / v = 3 is shown inFigure 14. This is the case with a super-critical shock,which means the gas temperature at the upstream sideof the spike is the same as the downstream tempera-ture. The downstream gas is hotter with an increasedratio between radiation pressure and gas pressure. Ouralgorithm is still able to keep the shock structures verywell as demonstrated in Figure 14 by comparing our nu-merical solution with the semi-analytical solution.4.9. Radiation Driven Wind
The simulation of radiation driven dusty wind men-tioned in the introduction is a useful test to demonstratethe difference between different approaches to solve theRT equation. This is a test that requires solving thefull radiation hydrodynamic equations and the simula-tion results can be different significantly if the RT equa-tion is not solved accurately. In order to compare withthe results based on VET method directly, we use thesame setup as in the run T3 F0 . − , L × [0 , L , wherethe length scale L = 2 . × cm is the fiducial scaleheight. The spatial resolution is 512 × T = 82K with the corresponding fiducial ve-locity v = 5 . × cm / s and the fiducial time scale t = 3 . × s. We adopt the same Rosseland andPlanck mean dust opacities as in Davis et al. (2014): κ a = 0 . (cid:18) T (cid:19) cm g − ,κ aP = 0 . (cid:18) T (cid:19) cm g − . (38)Notice that κ δP = κ aP − κ a is what we need in equa-tion 6. We initialize the gas with a constant density ρ = 7 . × − g cm − for 9 < y/L <
10 and den-sity floor 10 − ρ for all the other regions. The gas hasa uniform temperature T = T and a constant gravita-tional acceleration g = 1 . × − cm s − along the − y direction. The radiation field is initialized to give E r = a r T and F r = 0 . c g/κ a ( T ) through the wholesimulation box. We assume upward specific intensitieshave the same value at each cell and the same assump-tion is made for downward specific intensities. The sameapproach is used to set the bottom boundary conditionwith the same F r but E r = 9 . a r T . This is designedto maintain a constant vertical radiation flux F r with mplicit RT Algorithm τ = 3. At the top bound-ary, we set incoming specific intensities to be 0 and justcopy the outgoing specific intensities from the last acti-vate zones to the ghost zones. For gas quantities, we usereflecting boundary condition at the bottom and out-flow boundary condition at the top. We use randomperturbation on density with amplitude 12 .
5% to seedthe instability. The ratio between speed of light and thefiducial sound speed is C = 5 . × for this test. Figure 15.
Density distribution for two snapshots from thetest of radiation driven dusty wind as described in section4.9. The left panel shows the development of radiation drivenRayleigh-Taylor instability at time t = 35 t , while the rightpanel shows the fully non-linear stage at time t = 80 t . Here t is the typical sound cross time across one scale height. Density distribution at two snapshots t = 35 t and80 t are shown in Figure 15. The gas shell is pushedaway by the radiation force initially and Rayleigh-Taylorinstability has developed at t = 35 t . The high densityfilaments fall back temporarily but they get acceleratedagain later on. At time t = 80 t , the majority of the gasis still flowing outward. At each snapshot, we calculatethe box averaged Eddington ratio f E , volume ( τ V ) andflux weighted ( τ F ) optical depths as f E = (cid:104) ρκ a F r (cid:105) cg (cid:104) ρ (cid:105) ,τ V = 1000 L (cid:104) κ a ρ (cid:105) ,τ F = 1000 L (cid:104) κ a ρF r,y (cid:105)(cid:104) F r,y (cid:105) , (39)where (cid:104)·(cid:105) represents volume average over the whole sim-ulation box. Histories of f E , τ V and τ F /τ V are shownin Figure 16. The Eddington ratio dips below 1 around . . . f E τ V t/t . . . . τ F / τ V Figure 16.
Histories of various volume averaged quantitiesfrom the test of radiation driven dusty wind (section 4.9).The top panel is the volume averaged Eddington ratio f E while the Rosseland mean optical depth across the horizontaldirection of the simulation box is shown in the middle panel.The bottom panel shows the ratio between the flux weighted( τ F ) and volume averaged ( τ V ) optical depth. t due to the Rayleigh-Taylor instability but returnsto 1 quickly. The timescale when the dip shows up andthe evolution history are basically the same as Figure 6in Davis et al. (2014) despite some small difference at theearly phase, which is likely due to slightly different initialcondition. This demonstrates that the two completelydifferent algorithms capture the same growth rate andnonlinear structure of the radiation driven Rayleigh-Taylor instability. The Eddington ratio in our simu-lation eventually becomes larger than 1 and the gas isblown away from the top boundary. The 2D densitystructures shown in Figure 15 can also be directly com-pared with Figure 4 in Davis et al. (2014) and they sharevery similar structures.It is also interesting to compare the performance ofour algorithm for this test with the performance of VETscheme. The time step for our algorithm is 6 × − t and we use 50 iterations per time step to achieve a rel-ative accuracy of 10 − except the first few time steps,where 200 iterations are used per time step to pass theinitial transient. With 256 MPI ranks using skylake0 Jiang nodes, we are able to run roughly 2 × steps to reach120 t within about four hours of wall clock time. Thecost per time step is comparable to the cost of the VETmethod used by Davis et al. (2014). However, Daviset al. (2014) used a time step that was a factor of 10smaller to enable convergence, which caused the overallcost significantly larger than the cost of our new schemedeveloped here. Our implicit scheme also reduces thecost by a factor of ≈ if we would solve the problemusing explicit scheme by limiting the time step with thespeed of light. DISCUSSIONS AND CONCLUSIONSWe have developed a novel finite volume algorithmthat solves the full time dependent radiation magneto-hydrodynamic equations based on the discrete ordinatesapproach. The hydrodynamic equations are solved usingthe explicit MHD integrator in
Athena++ while the RTequation is solved fully implicitly. The gas and radiationare coupled via both energy and momentum exchange.The velocity dependent terms in the RT equation arehandled via Lorentz transformation so that no expan-sion in terms of v/c is needed. The time step of thewhole algorithm is only determined by the normal CFLcondition for the MHD module. We have carried out anextensive set of tests to demonstrate the accuracy andefficiency of the algorithm for both gas pressure and ra-diation pressure dominated flows with a wide range ofoptical depth. The finite volume approach also allowsthe algorithm to be easily used in both Cartesian andSpherical polar coordinate systems as demonstrated inour tests.The difference between the algorithm developed hereand radiation hydrodynamic schemes that adopt sim-plified assumption of closure relations for the radiationmoment equations (such as FLD and M1) are clear. Wedo not make any assumptions on the closure relationas we do not need it. Therefore, we do not subject tothe well known artifacts of FLD and M1, although ouralgorithm is generally more expensive. Our algorithmalso subjects to the well known ray effects as many othermethods based on discrete ordinates (Castor 2004). Thealgorithm is also different from the classical short char-acteristic method as well as the VET approach. We aresolving the full time dependent transport equation forspecific intensities directly, instead of the time indepen-dent transport equation, which is normally adopted forshort characteristic method. It is often argued that fornon-relativistic systems, the light crossing time is muchshorter than the local dynamical time scale and there-fore it is reasonable to drop the time dependent termfor the radiation field. This is widely used to estimate the heating and cooling rate of the gas for the modelingof stellar atmosphere. This is true for the low opticaldepth regime with negligible radiation pressure. How-ever, when the thermalization time scale is reduced withincreasing optical depth, this approach can be numeri-cally unstable even for gas pressure dominated regime(Davis et al. 2012). There are also systems where thetypical sound speed is much smaller than the speed oflight while radiation pressure is still larger than the gaspressure, such as the envelope of massive stars. Our al-gorithm can overcome the limitations of the short char-acteristic method and it can be suitable for all theseregimes. The VET method can also be used in principle,where a short characteristic module is only used to pro-vide a closure relation. As discussed in Section 4.3, thereis still a subtle difference compared with the algorithmdeveloped here in that the VET approach calculates theEddington tensor at the lab frame but neglecting theflow velocity. The Eddingon tensor is also calculated atthe beginning of each time step assuming that the radia-tion field has reached steady state with the gas. Futurestudies will be needed to investigate the consequencesof these differences for applications. Another advantageof our algorithm compared with VET is that the radi-ation field is entirely described one set of variables (thespecific intensities), while VET and Monte Carlo basedclosure schemes typically need to evolve two completelydifferent sets of radiation variables to get the Eddingtontensor and determine the time evolution of the radiationmoments. It can be an issue to keep the two sets of vari-ables to be consistent with each other (Park et al. 2012),which we do not need to worry about. The algorithm wehave developed here shares many common features withthe neutrino transport scheme developed by Sumiyoshi& Yamada (2012) and Nagakura et al. (2014). The maindifferences are the modified HLLE solver we have devel-oped here and the approach to solve the radiation sourceterms, where we do not linearize the equations.The computational cost of this algorithm can be veryproblem dependent, because the convergent propertiescan change dramatically depending on the opacity andhow quickly the radiation field is changing within onetime step. The cost to perform one iteration is smallerthan the cost of one step for the explicit scheme. Theoverall cost is roughly proportional to the number of it-erations we need for convergence. The algorithm alsoshows very good parallel scaling with domain decom-position because unlike the short characteristic module,there is no requirement to perform the iteration with anyparticular order. This means that each domain can per-form the iterations independently and they only need toexchange boundary condition after each iteration. How- mplicit RT Algorithm
Athena++ with ourRT module is our immediate next step.The current algorithm solves the frequency integratedtransport equation as the first step. It can be easily ex- tended to the frequency dependent case using the multi-group approach (Vaytet et al. 2011), which has beenwidely used to study neutrino transport (Just et al.2015; Skinner et al. 2019). The key for this extensionis to choose appropriate frequency bins, which can bein the lab frame or the fluid rest frame, or both (Na-gakura et al. 2014). All these options will be exploredwith future development.ACKNOWLEDGEMENTSThe author thanks James Stone, Shane Davis, Zhao-huan Zhu and the anonymous referee for valuable com-ments that have improved the draft. The Center forComputational Astrophysics at the Flatiron Institute issupported by the Simons Foundation.REFERENCES
Abdikamalov, E., Burrows, A., Ott, C. D., et al. 2012, ApJ,755, 111, doi: 10.1088/0004-637X/755/2/111Anninos, P., & Fragile, P. C. 2020, ApJ, 900, 71,doi: 10.3847/1538-4357/abab9cAsahina, Y., Takahashi, H. R., & Ohsuga, K. 2020, ApJ,901, 96, doi: 10.3847/1538-4357/abaf51Bjørgen, J. P., & Leenaarts, J. 2017, A&A, 599, A118,doi: 10.1051/0004-6361/201630237Bruls, J. H. M. J., Vollm¨oller, P., & Sch¨ussler, M. 1999,A&A, 348, 233Carlson, B. G. 1963, Methods in Computational Physics,Vol. 1, Statistical Physics, ed. B. Alder, S. Fernbach &M. RotenbergCastor, J. I. 2004, Radiation Hydrodynamics, ed. Castor,J. I.Davis, S. W., & Gammie, C. F. 2020, ApJ, 888, 94,doi: 10.3847/1538-4357/ab5950Davis, S. W., Jiang, Y.-F., Stone, J. M., & Murray, N.2014, ApJ, 796, 107, doi: 10.1088/0004-637X/796/2/107Davis, S. W., Stone, J. M., & Jiang, Y.-F. 2012, ApJS, 199,9, doi: 10.1088/0067-0049/199/1/9Densmore, J. D., Urbatsch, T. J., Evans, T. M., & Buksas,M. W. 2007, Journal of Computational Physics, 222, 485,doi: 10.1016/j.jcp.2006.07.031Dubroca, B., & Feugeas, J. 1999, Academie des SciencesParis Comptes Rendus Serie Sciences Mathematiques,329, 915, doi: 10.1016/S0764-4442(00)87499-6Fabiani Bendicho, P., Trujillo Bueno, J., & Auer, L. 1997,A&A, 324, 161Ferguson, J. M., Morel, J. E., & Lowrie, R. B. 2017, HighEnergy Density Physics, 23, 95,doi: 10.1016/j.hedp.2017.02.010 Foucart, F. 2018, MNRAS, 475, 4186,doi: 10.1093/mnras/sty108Gehmeyr, M., & Mihalas, D. 1994, Physica D NonlinearPhenomena, 77, 320, doi: 10.1016/0167-2789(94)90143-0Gonz´alez, M., Audit, E., & Huynh, P. 2007, A&A, 464, 429,doi: 10.1051/0004-6361:20065486Hayes, J. C., & Norman, M. L. 2003, ApJS, 147, 197,doi: 10.1086/374658Jiang, Y.-F., Stone, J. M., & Davis, S. W. 2012, ApJS, 199,14, doi: 10.1088/0067-0049/199/1/14—. 2013, ApJ, 767, 148, doi: 10.1088/0004-637X/767/2/148—. 2014, ApJS, 213, 7, doi: 10.1088/0067-0049/213/1/7Just, O., Obergaulinger, M., & Janka, H. T. 2015, MNRAS,453, 3386, doi: 10.1093/mnras/stv1892Krumholz, M. R., Klein, R. I., McKee, C. F., & Bolstad, J.2007, ApJ, 667, 626, doi: 10.1086/520791Krumholz, M. R., & Thompson, T. A. 2013, MNRAS, 434,2329, doi: 10.1093/mnras/stt1174Levermore, C. D., & Pomraning, G. C. 1981, ApJ, 248, 321,doi: 10.1086/159157Lowrie, R. B., & Edwards, J. D. 2008, Shock Waves, 18,129, doi: 10.1007/s00193-008-0143-0Lowrie, R. B., Morel, J. E., & Hittinger, J. A. 1999, ApJ,521, 432, doi: 10.1086/307515McKinney, J. C., Tchekhovskoy, A., Sadowski, A., &Narayan, R. 2013, ArXiv e-prints.https://arxiv.org/abs/1312.6127Mihalas, D., & Klein, R. I. 1982, Journal of ComputationalPhysics, 46, 97, doi: 10.1016/0021-9991(82)90007-9Mihalas, D., & Mihalas, B. W. 1984, Foundations ofradiation hydrodynamics, ed. Mihalas, D. & Mihalas,B. W. Jiang
Nagakura, H., Sumiyoshi, K., & Yamada, S. 2014, ApJS,214, 16, doi: 10.1088/0067-0049/214/2/16Park, H., Knoll, D. A., Rauenzahn, R. M., Wollaber, A. B.,& Densmore, J. D. 2012, Transport Theory andStatistical Physics, 41, 284,doi: 10.1080/00411450.2012.671224Pons, J. A., Ib´a˜nez, J. M., & Miralles, J. A. 2000, MNRAS,317, 550, doi: 10.1046/j.1365-8711.2000.03679.xProga, D., Jiang, Y.-F., Davis, S. W., Stone, J. M., &Smith, D. 2014, ApJ, 780, 51,doi: 10.1088/0004-637X/780/1/51Radice, D., Abdikamalov, E., Rezzolla, L., & Ott, C. D.2013, Journal of Computational Physics, 242, 648,doi: 10.1016/j.jcp.2013.01.048Roth, N., & Kasen, D. 2015, ApJS, 217, 9,doi: 10.1088/0067-0049/217/1/9Ryan, B. R., & Dolence, J. C. 2020, ApJ, 891, 118,doi: 10.3847/1538-4357/ab75e1S¸adowski, A., Narayan, R., Tchekhovskoy, A., & Zhu, Y.2013, MNRAS, 429, 3533, doi: 10.1093/mnras/sts632Sekora, M. D., & Stone, J. M. 2010, Journal ofComputational Physics, 229, 6819,doi: 10.1016/j.jcp.2010.05.024Sincell, M. W., Gehmeyr, M., & Mihalas, D. 1999, ShockWaves, 9, 391, doi: 10.1007/s001930050169Skinner, M. A., Dolence, J. C., Burrows, A., Radice, D., &Vartanyan, D. 2019, ApJS, 241, 7,doi: 10.3847/1538-4365/ab007fSkinner, M. A., & Ostriker, E. C. 2013, ApJS, 206, 21,doi: 10.1088/0067-0049/206/2/21 Smit, J. M., Cernohorsky, J., & Dullemond, C. P. 1997,A&A, 325, 203Smith, A., Kannan, R., Tsang, B. T. H., Vogelsberger, M.,& Pakmor, R. 2020, ApJ, 905, 27,doi: 10.3847/1538-4357/abc47eSteinacker, J., Baes, M., & Gordon, K. D. 2013, ARA&A,51, 63, doi: 10.1146/annurev-astro-082812-141042Stone, J. M., & Mihalas, D. 1992, Journal of ComputationalPhysics, 100, 402, doi: 10.1016/0021-9991(92)90246-UStone, J. M., Mihalas, D., & Norman, M. L. 1992, ApJS,80, 819, doi: 10.1086/191682Stone, J. M., Tomida, K., White, C. J., & Felker, K. G.2020, ApJS, 249, 4, doi: 10.3847/1538-4365/ab929bSumiyoshi, K., & Yamada, S. 2012, ApJS, 199, 17,doi: 10.1088/0067-0049/199/1/17Teyssier, R. 2015, ARA&A, 53, 325,doi: 10.1146/annurev-astro-082214-122309Tsang, B. T. H., & Milosavljevi´c, M. 2018, MNRAS, 478,4142, doi: 10.1093/mnras/sty1217Turner, N. J., & Stone, J. M. 2001, ApJS, 135, 95,doi: 10.1086/321779van der Holst, B., T´oth, G., Sokolov, I. V., et al. 2011,ApJS, 194, 23, doi: 10.1088/0067-0049/194/2/23Vaytet, N. M. H., Audit, E., Dubroca, B., & Delahaye, F.2011, JQSRT, 112, 1323, doi: 10.1016/j.jqsrt.2011.01.027Whitney, B. A. 2011, Bulletin of the Astronomical Societyof India, 39, 101. https://arxiv.org/abs/1104.4990Zel’Dovich, Y. B., & Raizer, Y. P. 1967, Physics of shockwaves and high-temperature hydrodynamic phenomena,ed. Zel’Dovich, Y. B. & Raizer, Y. P. mplicit RT Algorithm A. COEFFICIENTS OF THE DISCRETIZED RT EQUATIONIn our algorithm, the RT equation is solved fully implicitly as (equation 14 and section 3.2.1) I m +1 n + ∆ t ∇ · (cid:2) ( c n − f v m ) I m +1 n (cid:3) = I mn − ∆ t ∇ · ( f v m I mn )+ ∆ tc Γ − n (cid:34) ρ m κ s (cid:0) J m +10 − I m +10 ,n (cid:1) + ρ m κ a (cid:32) a r (cid:0) T m +1 (cid:1) π − I m +10 ,n (cid:33) + ρ m κ δP (cid:32) a r (cid:0) T m +1 (cid:1) π − J m +10 (cid:33)(cid:35) . (A1)At each cell ( i, j, k ), the left hand side is calculated as I m +1 n + ∆ tV i,j,k [ A x ( i + 1 / F n ( i + 1 / − A x ( i − / F n ( i − / A y ( j + 1 / F n ( j + 1 / − A y ( j − / F n ( j − / A z ( k + 1 / F n ( k + 1 / − A z ( k − / F n ( k − / , (A2)where A x , A y , A z are the cell areas at the corresponding surfaces while V i,j,k is the cell volume. The surface flux iscalculated via the modified HLLE solver (taking F n ( i − /
2) as an example) F n ( i − /
2) = S + i − / S + i − / − S − i − / [ cµ x − f v mx ( i − / I m +1 n ( i − − S − i − / S + i − / − S − i − / [ cµ x − f v mx ( i − / I m +1 n ( i )+ S + i − / S − i − / S + i − / − S − i − / (cid:2) I m +1 n ( i ) − I m +1 n ( i − (cid:3) . (A3)Here the maximum ( S + i − / ) and minimum ( S − i − / ) signal speeds are defined as S + i − / = cµ x (cid:112) [1 − exp ( − τ c )] /τ c if µ x > − cµ x (cid:112) [1 − exp ( − τ c )] /τ c if µ x < , (A4) S − i − / = − cµ x (cid:112) [1 − exp ( − τ c )] /τ c if µ x > cµ x (cid:112) [1 − exp ( − τ c )] /τ c if µ x < . (A5)The optical depth per cell is defined as τ c ≡ α [ ρ m ( i −
1) + ρ m ( i )] [ κ a ( i −
1) + κ a ( i ) + κ s ( i −
1) + κ s ( i )] ∆ x , where ∆ x is the local cell size and α is a free parameter with the default value chosen to be 5. The value of α is chosen tominimize numerical diffusion while still keep the numerical scheme stable. For all the test problems we have shown,the solution is independent of α as long as it is large enough. The signal speeds as well as the HLLE flux are chosento preserve the asymptotic limit such that when τ c → ∞ , S + i − / → c | µ x | /τ c , S − i − / → − c | µ x | /τ c and F n ( i − / →
12 [ cµ x − f v mx ( i − / (cid:2) I m +1 n ( i −
1) + I m +1 n ( i ) (cid:3) − c | µ x | τ c (cid:2) I m +1 n ( i ) − I m +1 n ( i − (cid:3) . (A6)In the optically thin limit such that τ c → F n ( i − /
2) is reduced to the simple upwind flux. Similarly, surfacefluxes from other directions as well the signal speeds S + j − / , S − j − / , S + j +1 / , S − j +1 / , S + k − / , S − k − / , S + k +1 / , S − k +1 / canbe calculated.Replacing the expressions for the fluxes at cell faces, we get the discretized equation g I m +1 n + g I m +1 n ( i −
1) + g I m +1 n ( i + 1) + g I m +1 n ( j − g I m +1 n ( j + 1) + g I m +1 n ( k −
1) + g I m +1 n ( k + 1)= I mn − ∆ t ∇ · ( f v m I mn )+ ∆ tc Γ − n (cid:34) ρ m κ s (cid:0) J m +10 − I m +10 ,n (cid:1) + ρ m κ a (cid:32) a r (cid:0) T m +1 (cid:1) π − I m +10 ,n (cid:33) + ρ m κ δP (cid:32) a r (cid:0) T m +1 (cid:1) π − J m +10 (cid:33)(cid:35) , (A7)4 Jiang where the coefficients are given by g = 1 − A x ( i − / tV i,j,k S − i − / S + i − / − S − i − / (cid:104) − cµ x + f v mx ( i − /
2) + S + i − / (cid:105) + A x ( i + 1 / tV i,j,k S + i +1 / S + i +1 / − S − i +1 / (cid:104) cµ x − f v mx ( i + 1 / − S − i +1 / (cid:105) − A y ( j − / tV i,j,k S − j − / S + j − / − S − j − / (cid:104) − cµ y + f v my ( j − /
2) + S + j − / (cid:105) + A y ( j + 1 / tV i,j,k S + j +1 / S + j +1 / − S − j +1 / (cid:104) cµ y − f v my ( j + 1 / − S − j +1 / (cid:105) − A z ( k − / tV i,j,k S − k − / S + k − / − S − k − / (cid:104) − cµ z + f v mz ( k − /
2) + S + k − / (cid:105) + A z ( k + 1 / tV i,j,k S + k +1 / S + k +1 / − S − k +1 / (cid:104) cµ z − f v mz ( k + 1 / − S − k +1 / (cid:105) ,g = − A x ( i − / tV i,j,k S + i − / S + i − / − S − i − / (cid:104) cµ x − f v mx ( i − / − S − i − / (cid:105) ,g = A x ( i + 1 / tV i,j,k S − i +1 / S + i +1 / − S − i +1 / (cid:104) − cµ x + f v mx ( i + 1 /
2) + S + i +1 / (cid:105) ,g = − A y ( j − / tV i,j,k S + j − / S + j − / − S − j − / (cid:104) cµ y − f v my ( j − / − S − j − / (cid:105) ,g = A y ( j + 1 / tV i,j,k S − j +1 / S + j +1 / − S − j +1 / (cid:104) − cµ y + f v my ( j + 1 /
2) + S + j +1 / (cid:105) ,g = − A z ( k − / tV i,j,k S + k − / S + k − / − S − k − / (cid:104) cµ z − f v mz ( k − / − S − k − / (cid:105) ,g = A z ( k + 1 / tV i,j,k S − k +1 / S + k +1 / − S − k +1 / (cid:104) − cµ z + f v mz ( k + 1 /
2) + S + k +1 / (cid:105) . (A8)The convergent properties for our iterative scheme depends strongly on the values of g . If g is large and positivefor all I m +1 n , the diagonal term will dominate over the off-diagonal terms and the iteration can converge quickly.Otherwise, the convergence may be slow and it can even diverge sometimes. To improve the robustness of our iterativescheme, we have an alternative procedure to perform the iteration compared with equation 21. The expression for g includes terms that can be positive or negative depending on the sign of µ x , µ y , µ z . For each angle, we separate g into two parts g = g (cid:48) + g (cid:48)(cid:48) , where g (cid:48) includes all the terms that are positive while all the negative terms are includedin g (cid:48)(cid:48) . At each iterative step l , we replace g I m +1 n,l in equation 21 with g (cid:48) I m +1 n,l + g (cid:48)(cid:48) I m +1 n,l −1