An asymptotically correct implicit–explicit time integration scheme for finite volume radiation-hydrodynamics

ABSTRACT Numerical radiation-hydrodynamics (RHD) for non-relativistic flows is a challenging problem because it encompasses processes acting over a very broad range of time-scales, and where the relative importance of these processes often varies by orders of magnitude across the computational domain. Here, we present a new implicit–explicit method for numerical RHD that has a number of desirable properties that have not previously been combined in a single method. Our scheme is based on moments and allows machine-precision conservation of energy and momentum, making it highly suitable for adaptive mesh refinement applications; it requires no more communication than hydrodynamics and includes no non-local iterative steps, making it highly suitable for massively parallel and Graphics Processing Unit (GPU)-based systems where communication is a bottleneck; and we show that it is asymptotically accurate in the streaming, static diffusion, and dynamic diffusion limits, including in the so-called asymptotic diffusion regime where the computational grid does not resolve the photon mean-free path. We implement our method in the GPU-accelerated RHD code quokka and show that it passes a wide range of numerical tests.


INTRODUCTION
Radiation-hydrodynamics (RHD) plays a significant role in astrophysics, influencing the evolution and energy distribution in various astrophysical systems or phenomena -stellar atmospheres (e.g.Mihalas 1978), planetary atmospheres (e.g.Zhang 2020), core-collapse supernovae (e.g.Skinner et al. 2016;Radice et al. 2018), star formation in a variety of environments (e.g.Thompson et al. 2005;Krumholz et al. 2009;Rosen et al. 2016;He et al. 2019;Menon et al. 2023), active galactic nuclei and jets (e.g.Davis & Tchekhovskoy 2020), and galactic outflows (e.g.Naab & Ostriker 2017;Zhang 2020).The dynamics of RHD systems vary substantially in a range of scales and physical conditions, parametrized by the typical optical depth and the sound speed of radiation-interacting gas, which determines how radiation is transported.While several numerical techniques exist to solve the RHD equations in various limits, developing methods that are accurate across all regimes and that run efficiently on modern, GPU-based architectures remains an ongoing challenge.
There are well-known difficulties associated with solving the RHD equations numerically for non-relativistic systems.One significant challenge arises from the huge difference in timescales associated with radiation and hydrodynamics -characterised by the speed of light  and gas flow speed , respectively.A second is that realistic RHD problems often contain a huge range of opacities, such that the photon mean free path may be comparable to the size of the entire simulation domain in some regions, while in others it may be so small ★ E-mail: Chongchong.He@anu.edu.au(CCH) as to be impossible to resolve at reasonable computational cost.The RHD equations are well-suited to hydrodynamics-like explicit solution methods in some regimes, but often the source terms coupling the gas and radiation are so stiff that an implicit method must be used to ensure stability.
To deal with these challenges, it is natural to solve both the transport and source terms in the RT equations implicitly so that the time step is not limited by the speed of light.This is the approach adopted by many RHD codes (e.g.Krumholz et al. 2007;Zhang et al. 2011;Jiang et al. 2012;Menon et al. 2022).However, this approach suffers from the need for an implicit treatment of the transport term, which is non-local.In multiple dimensions, implicit update of this term is usually accomplished by solving a sparse matrix system, and sparse matrix solvers show limited scalability and performance on modern massively parallel and GPU-accelerated architectures, where the high and unpredictable communication load they involve becomes a performance bottleneck.
Several authors have explored alternative approaches in which the source and transport terms are operator-split, with the former treated explicitly and the latter implicitly (e.g., Rosdahl & Teyssier 2015;Skinner et al. 2019;Wibking & Krumholz 2022).These schemes generally achieve better scaling and speed, particularly when running in parallel on large numbers of GPUs.One of the disadvantages of this approach is that the explicit treatment of the transport term requires a time step significantly smaller than the hydrodynamics one (since  ≫ ), but it is possible to alleviate this problem by subcycling the radiation transport step relative to hydrodynamics.Since the transport equations for radiation are much simpler than those for hydrodynamics, it is possible to carry out many radiation transport updates per hydrodynamic update at a comparable cost.If necessary the cost can be mitigated further using the reduced speed of light (RSLA) approximation (RSLA; Gnedin & Abel 2001;Skinner & Ostriker 2013).
A second challenge in the operator-split approach has received less attention in the literature, but is perhaps even more serious: the need to properly balance the transport and source terms across all asymptotic regimes.This balance is critical because the RT equation behaves as an advection equation in optically thin regimes, but transitions to a diffusion equation in optically thick regimes, and nearperfect cancellations between parts of the source and transport terms are responsible for this behaviour.Efforts to recover this property of the RHD system have thus far mostly involved ad hoc corrections to the Riemann solver or to the source terms to recover the correct asymptotic limit.For instance, Rosdahl & Teyssier (2015) add a "trapped photons" term to the source term to account for diffusion.Skinner et al. (2019) apply corrections to the wave speed in the Riemann solver, scaling down the characteristic speed of the radiation modes in optically thick regimes by a factor of 1/ √  cell , where  cell is the optical depth per cell, and , chosen empirically, is a small integer.A similar approach is used by Wibking & Krumholz (2022).Despite these efforts, such corrections are often only partly successful, and their accuracy across a wide range of parameter space has not been tested.Perhaps the most successful (in terms of accuracy) operator-split method published to date is the discontinuous Galerkin implicit-explicit (DG-IMEX) scheme proposed by McClarren et al. (2008), but as we discuss below even this scheme is not accurate in all RHD limits, and it has a number of other undesirable properties as well.
This situation motivates our goal of designing a method that achieves the best of both worlds: accuracy in all RHD regimes that is comparable to that achieved by fully-implicit methods, but without the need for poorly-scaling communication-intensive operations like sparse matrix solves.In this paper, we describe a method that achieves this goal using a novel time-integration scheme that recovers the proper asymptotic limit in the radiation diffusion regime without requiring non-local implicit updates.Our method is based on a convex-invariant, asymptotic-preserving IMEX approach that is second-order accurate in streaming limit, wherein the transport terms are handled explicitly, while the matter-radiation interacting part is treated implicitly and locally, eliminating the need for non-local implicit terms in iteration.
We begin in Section 2 by introducing the full set of two-moment RHD equations to be solved and deriving characteristic numbers and limiting behaviours for them, laying the foundation for our analysis.Then, in Section 3, we present the IMEX scheme.In Section 4, we derive some properties of our scheme and compare to alternative approaches.Finally, we present in Section 5 tests of our implementation, demonstrating its effectiveness and applicability.

THE RHD EQUATIONS AND NUMERICAL METHODS FOR SOLVING THEM
We begin our analysis by describing the RHD system of equations in Section 2.1, non-dimensionalising it to obtain characteristic numbers and limiting regimes in Section 2.2, and then using those to gain insight into the challenges of designing numerical methods for RHD and to motivate our approach in Section 2.3.

The RHD system
Our method begins from the fundamental equations of RHD1 expressed in an inertial (lab) reference frame.These are where are vectors of the conserved quantities, the advection terms, and the source terms, respectively; in the equations above,  is the matter density,  is the matter velocity,  is the matter pressure,  gas is the gas total energy density,  is the radiation energy density,  is the radiation flux, P is the radiation pressure tensor, and ( 0 , ) is the radiation four-force.A significant advantage of working with labframe rather than comoving-frame radiation quantities is that these equations are manifestly conservative, a feature that the algorithm we describe below will preserve.As usual in the moment formulation, however, one must adopt a closure relation for the radiation pressure tensor P.There are a wide range of possible closures, and since our scheme is independent of this choice, we will not discuss closure relations further here.Our next step is to write out the radiation four-force in the mixedframe formulation, whereby we write the matter-radiation exchange coefficients in the comoving frame, where they are simplest, while all other quantities remain in the lab frame.We assume that the emitting matter is in local thermodynamic equilibrium, so its emissivity is proportional to the Planck function, and we neglect scattering.To order  2 / 2 the result, expressed in index notation and adopting the Einstein summation convention, is (Krumholz et al. 2007) Here a subscript 0 in  indicates the absorption coefficient is expressed in the comoving rather than the lab frame, and  0 ,  0 , and  0 are the comoving-frame Planck, energy, and flux mean absorption coefficients.The leading-order part of the time-like component of the mixed-frame radiation four-force  0 is the classical rate of radiation-matter energy exchange in the comoving frame, while the order / part combines the work done by the radiation force on matter (− 0     /) with a purely relativistic effect arising from the transformation of the opacity between the comoving and lab frames (2 0     /); the second-order terms are similarly relativistic effects arising from the frame transformation.In the space-like component   , the leading-order term is the radiation force in the comoving frame, while the remaining order / and  2 / 2 terms describe relativistic effects that can be interpreted as frame-dragging between matter and radiation.

Characteristic numbers and limiting regimes
We next derive characteristic dimensionless numbers for RHD systems and consider the limiting behaviour of various terms as we alter their relative sizes.We non-dimensionalise Equation 1 and Equation 2 following Lowrie et al. (1999).For the matter quantities, we let ℓ be the characteristic size of the system,  ∞ be the characteristic isothermal sound speed, and  ∞ be the characteristic density, and define Here the quantities with carets are dimensionless versions of the dimensional quantities.For the radiation quantities, we let  ∞ be the reference temperature and  ∞ be the reference length scale2 , and we define With these definitions, the radiation four-force becomes and the corresponding non-dimensionalised version of Equation 1 is where ∇ ≡ ℓ∇ and In these expressions, we have introduced three dimensionless quantities: These numbers represent, respectively, the dimensionless speed of light, the ratio of the system size to the photon mean free path (and thus is equal to the characteristic optical depth), and (up to factors of order unity) the ratio of radiation pressure to gas pressure at the characteristic temperature and sound speed of the system.It is therefore clear that an RHD system is determined by these three characteristic numbers; the first is always much greater than unity for a non-relativistic system, but the remaining two can be of any size.Now let us consider various limiting cases of the dimensionless numbers, focusing in particular on the evolution equations for the radiation quantities Ê and F (the last two entries in Equation 13); this will simplify our task since P does not appear in these equations.For numerical convenience, and since our goal here is insight rather than rigorous calculation, we will also at this point specialise to the case of grey material, which allows us to choose our scaling  ∞ such that χ0 = χ0 = χ0 = 1.In this case the non-dimensional radiation four-force reduces to For L ≪ 1, the system is optically thin and we are in the streaming regime.In this case it is clear that the largest term is ∇ • F Û , and so on a fluid flow timescale it is clear that the solution is simply ∇ • F = 0 and ∇ • P = 0. On the other hand, for L ≫ 1, we are in the diffusion limit, and it is clear that Ŝ Û , which is of order CL, is the largest term, and therefore on a fluid flow timescale to leading order we must have ( Ĝ0 , Ĝ) = (0, 0).This in turn requires that, to leading order, Ê = T4 and F = 0; though we have not shown it here, it is straightforward to show that in this limit we also have P = (1/3)I Ê, where I is the identity tensor (Mihalas & Mihalas 1984).Since the flux is zero to the leading order in this case, we must proceed to the next order to determine its value.To do so, we Taylor expand the radiation energy density, flux, and pressure tensor about the leading order solution: Ê = T4 + Ê(1) F = F(1) P = 1 3 where terms subscripted (1) are perturbations that we take to be small compared to the leading-order terms.We insert these expressions into Equation 13, Taylor expand, and linearise by dropping all terms that involve products of perturbed quantities.Which terms are at leading order after this procedure depends on the relative values of C and L.
If C ≫ L, known as the static diffusion limit, then the leading order surviving terms in the equation for the flux are Since both C and CL are large compared to unity, on a fluid flow timescale the terms proportional to these factors must cancel, and therefore we have This is the usual Fick's Law diffusion approximation.Armed with this leading order result, we can see that the relative sizes of the terms in the radiation evolution equations (/ t)( Ê, F) : C ∇ • ( F, P) : CL ( Ĝ0 , Ĝ) scale relative to one another as 1 : C/L : C/L.On the other hand, if we have C ≪ L, known as the dynamic diffusion case, then the leading non-zero terms are As in the previous case, since C and CL are large compared to unity, on a fluid flow timescale the two terms on the right hand side must cancel to leading order, and we instead have In this case the relative scalings of the terms in the radiation evolution equations (/ t) ( Ê, F) : C ∇ • ( F, P) : CL ( Ĝ0 , Ĝ) are 1 : 1 : 1, i.e., all terms are of equal order.

Design considerations for numerical methods
While the exploration of the limiting cases here is not new (e.g., Mihalas & Mihalas 1984;Lowrie et al. 1999;Krumholz et al. 2007), revisiting it allows us to make some important observations about design considerations for numerical methods.First, which terms are large on a fluid flow timescale, and relative to each other, changes from one RHD regime to another -in the streaming limit the transport terms proportional to C dominate, in the static diffusion regime these terms come into balance with the source terms and both are at order C/L ≫ 1, while in the dynamic diffusion regime all terms, including the time derivative, become of order unity.Therefore if an RHD scheme is to correctly capture the limiting behaviour in each of these regimes, it cannot rely on any particular assumptions about the relative orderings of these terms, and must be able to cope with situations where each of them is both dominant and sub-dominant.
A second consideration, which we have already introduced in Section 1, is that in the diffusion regime obtaining the correct solution depends on the dominant terms cancelling to high-order.That is, we have seen that the source terms are naturally of size CL, but in static diffusion the leading order terms cancel and so the dominant nonvanishing term is smaller by a factor L 2 ≫ 1, while for dynamic diffusion it is smaller by a factor CL ≫ 1.Similarly, the natural sizes of the transport terms are C, but this is reduced by a factor L in static diffusion, and by a factor C in dynamic diffusion.This creates problems for schemes where the transport and source terms, or parts of the source terms, are operator-split, because in an operator split scheme it is difficult to recover the proper near-exact cancellations between the various terms.Overcoming this problem will be our primary objective.
Before describing our proposed solution, however, we pause to discuss another possible approach that is somewhat close to ours in spirit: the DG-IMEX method of McClarren et al. (2008).This method uses a temporal discretisation that accurately captures cancellations in the static diffusion limit while avoiding non-local implicit solves, and McClarren et al. show that it passes a number of tests that other operator-split methods fail.However, the price for this is high: the method requires a careful and complex reconstruction scheme that is heavy in terms of both computation and memory usage -indeed, the method requires eight degrees of freedom per cell, and therefore requires eight times as much memory as the finite volume scheme we present below, a major problem for GPU-based computations where memory is at a premium.Additionally, in the full radiationhydrodynamics context (not considered by McClarren et al.), coupling a DG radiation solver to a finite volume hydrodynamics code requires careful consideration of how the fluid internal energy is mapped from the finite volume grid to the DG nodes in order to maintain the asymptotic-preserving property.While existing methods for this mapping numerically manifest the correct asymptotic diffusion solution, they have not been subjected to a rigorous asymptotic analysis such as the one we present below for our scheme, and thus their ability to produce correct asymptotic behaviour over all parameter regimes remains unproven (Bolding et al. 2017).Moreover, extending a DG scheme to adaptive mesh refinement (AMR) would be a substantial challenge, whereas the finite volume scheme we propose is fully compatible with existing AMR frameworks.Finally, McClarren et al. focus exclusively on the static diffusion regime, and their scheme is not easily extensible to either the streaming or dynamic diffusion limits -the former because the method relies on a specific closure relation that is a poor approximation for streaming radiation, and the latter because the scheme does not include the velocity-dependent terms that become order unity in the dynamic diffusion limit.Our scheme, by contrast, applies to all RHD regimes.

A NEW ASYMPTOTIC-PRESERVING SCHEME FOR RADIATION HYDRODYNAMICS
We now proceed to describe our new numerical method.We first describe our overall strategy for the full RHD system in Section 3.1, then the IMEX scheme we use for the radiation subsystem in Section 3.2, and finally our method for carrying out each IMEX stage in Section 3.3.

Overall time stepping strategy
We solve the system formed by Equation 1 using an operator split approach consisting of two major steps.In the first step, we advance the hydrodynamic transport subsystem using an explicit method.In the second step, we update the radiation transport subsystem and matterradiation coupling terms using a mixed implicit-explicit (IMEX) method.
The hydrodynamic transport subsystem consists of the partial differential equations (PDEs) Our scheme for radiation does not depend on the numerical method used to solve this subsystem; for the purposes of our implementation in Quokka (Wibking & Krumholz 2022), which we use for all the tests below, we adopt a semidiscrete approach, discretising the spatial variables onto a grid while initially keeping the time variable continuous, thereby transforming the partial differential equations into a large set of ODEs.These ODEs are then integrated in time utilizing the second-order accurate, strong stability preserving Runge-Kutta method (RK2-SSP; Shu & Osher 1988).For an in-depth explanation of the method we direct readers to Wibking & Krumholz (2022).
In the present work, our emphasis is on the second major step, the radiation update and matter-radiation coupling.These steps are subcycled with respect to the hydrodynamic step, since they require smaller time steps.

Implicit-explicit method for the radiation subsystem
Using the same method of lines approach for the radiation subsystem and radiation-matter coupling terms as for hydrodynamics, we define     as the volume-averaged radiation energy density in cell   , and similarly for all other variables, and express the radiation subsystem for each cell as where, dropping the    subscript from this point forward for convenience, we define The characteristic time scales associated with the transport term, , which consists of the divergence of radiation flux, ∇ • , and the divergence of the radiation pressure tensor ∇ • P, are C −1 , which is potentially fast compared to hydrodynamics, but in the diffusion regime is much longer than the timescale (CL) −1 associated with the source term  .3Consequently, we opt to evolve the transport terms using an explicit method.This choice is strategic; the explicit update obviates the need for global communication across the computational domain.This advantage becomes particularly pronounced on GPUs, where inter-device communication often represents the primary bottleneck in performance.In problems where C is so large that this requires infeasibly-many explicit steps, we may elect to solve the approximate RSLA equations instead (Gnedin & Abel 2001;Skinner & Ostriker 2013;Wibking & Krumholz 2022), for which we instead have where ĉ is the reduced speed of light, chosen so to be ≪  but still much greater than any hydrodynamic speed.This approximation reduces the size of the transport term from C to ( ĉ/)C, and thus allows larger time steps.For generality in what follows we will write our update scheme using the RSLA equations, but these can be reduced to the exact equations of RHD simply by setting ĉ = .
The short time scales associated with the matter-radiation coupling terms,  , require an implicit treatment.Notably, as these terms do not have spatial derivatives, they allow for the independent update of each cell within the computational domain.Therefore, in the whole radiation update, there is no non-local implicit update, effectively eliminating the requirement for additional inter-domain communication beyond what is standard in a pure hydrodynamic update.
To ensure that our choice to operator-split between the source and transport terms in this manner does not compromise accuracy, we utilize the asymptotic-preserving IMEX PD-ARS integrator (Chu et al. 2019), a choice motivated by its proven effectiveness in preserving the diffusion limit while maintaining second-order accuracy and stability in the streaming limit. 4The IMEX PD-ARS integrator can be characterized by its double Butcher tableau where  is a free parameter in the range [0, 0.5); in our implementation we adopt  = 0.When expressed in equations, the scheme to advance the system by a time Δ is where the superscript () indicates the state at the start of the radiation update (but after the operator-split hydrodynamic update), ( + 1) indicates the state at the end of the radiation update, and ( + 1/2) indicates an intermediate stage.
One key feature of this update scheme is that the transport and source terms appear symmetrically at each of the two stages, so that cancellations can be captured properly.A second key feature is that, on modern architectures where communication is expensive, this scheme is only marginally more costly than an update like RK2-SSP, because the only excess work it requires is an extra iterative solve at the first stage (Equation 28).Crucially, this extra solve is purely local, because only the local source terms, rather than the non-local transport terms, must be iterated.Thus there is no extra communication in this scheme relative to RK2-SSP.

Update strategy for IMEX stages
At each stage of the integrator above, the right-hand side contains both explicit terms -those that are evaluated at state () during the first IMEX stage (Equation 28) and at states () or ( + 1/2) during the second stage (Equation 29) and implicit terms that are evaluated at state ( + 1/2) during the first stage and state ( + 1) during the second stage.We therefore begin each stage by evaluating the explicit terms.This in turn requires that we evaluate the transport terms  , which are always explicit in keeping with our overall strategy.We do this using a Godunov method.Specifically, we compute the fluxes for  and  between cells by using PPM reconstruction to obtain the states at the cell edges and then using a Harten, Lax, van Leer (HLL) Riemann solver (Harten et al. 1983) to compute the fluxes from these states.The details of this step are elaborated in Wibking & Krumholz (2022), and hence we will not repeat them here.However, we modify the scheme in an important way compared to the method presented in Wibking & Krumholz: that scheme, and all previous comparable explicit treatments of radiation, have been forced to invoke an ad hoc correction to the wavespeeds computed in the Riemann solver in an attempt to recover the diffusion limit.With our new time-stepping scheme, no such modification is required, and we can simply use the uncorrected HLL fluxes in our update.We discuss this issue further in Section 4.3.
After carrying out the explicit part of the update, we are left with the implicit part.For both stages of the IMEX integrator we can express this stage in the generic form for the energies, and for the momenta.Here terms with superscript () denote the state after applying the explicit terms -for example during the first stage (Equation 28) we have  ( ) =  () − Δ ( ĉ/)∇ •  () -and terms with superscript ( + 1) denote the final state for which we are attempting to solve; the factor  is unity during the first stage and 1/2 during the second stage.
In each cell this is a system of eight equations in eight unknowns -radiation and gas energy, three components of gas velocity, and three components of radiation flux -that must be solved simultaneously.While we could do so straightforwardly using a single Newton-Raphson iteration scheme or similar, it is more efficient to separate the iteration procedure into an inner stage where we solve the energy equations while freezing  and , and a second, outer stage where update the flux and gas velocity and then if necessary go back to the inner stage and recompute  and  gas using the updated values of  and .The reason this is more efficient is that it allows the inner iteration stage to consider only two variables rather than eight, and in most cases the work terms proportional to  and  are small in the energy equations, so  and  gas change little to none as a result of the update to  and  and the whole procedure converges in a single or at most a few outer iterations.
The inner stage consists of updating the energy quantities using a modified version of the iteration scheme initially proposed by Howell & Greenough (2003) and modified by Wibking & Krumholz (2022).We solve Equation 30 and Equation 31for  (+1) gas and  (+1) by performing Newton-Raphson iteration on the non-linear implicit system where and  is an optional term to include, for example, the addition of radiation by stellar sources.Note that , P, and all the variables that can depend on  gas - 0 ,  0 ,  0 , and  -carry the superscripts ( + 1) which are omitted here for the sake of brevity.By contrast, we freeze  and  at their values at the start of the inner stage as noted above.
A single Newton-Raphson iteration consists of solving the linearised equations where  is the set of variables to be updated, Δ is the change in these variables during this iteration, () is the vector whose zero we wish to find, and J is the Jacobian matrix of ().Instead of taking  = ( gas , ), as in the original Howell & Greenough (2003) scheme, we use  = ( gas , ) as the base variables over which to iterate; we find that the system generally converges in fewer iterations using this basis, likely because at high optical depth the system almost immediately converges to the solution  = 0, and thus the remaining iterations are effectively on  gas alone The Jacobian in this basis is We have omitted the  2 / 2 terms and assumed  (  /  )/ = 0 in the calculation of the Jacobian for simplicity, but this simplification only changes the rate of convergence; it does not affect the converged solution.
After solving Equation 37 for Δ, we update  ←  + Δ; we do this repeatedly until the system converges, as determined by the condition where is the total radiation and material energy at the beginning of the timestep accounted for reduced speed of light.We set the relative tolerance  = 10 −13 by default.Once this Newton-Raphson system converges, we have the updated gas total energy  (+1) gas and we can compute the updated radiation energy as (41) Note that, for ĉ =  and  = 0, this procedure ensures that our scheme conserves total energy to machine precision regardless of the level of accuracy with which we have iterated the equations to convergence.
We then proceed to the outer stage of the iteration where we solve the flux and momentum update equations, Equation 32 and Equation 33, with the updated gas temperature, opacity, and radiation energy.To order /, the solution is straightforward: To order  2 / 2 , in cases where  0 −  0 = 0, the solution is When  0 ≠  0 and to order  2 / 2 , the solution is slightly more complex because all three components of  appear in   .In this case,  (+1)  is the solution of a set of three linear equations; these are straightforward to solve analytically, but the resulting expressions are somewhat lengthy and so we omit them here for brevity.Lastly, following Equation 32, we update the gas momentum via This update also ensures momentum conservation to machine precision.After we update the gas momentum, we also recalculate the gas's internal energy (which in Quokka we track separately because we implement a dual energy formalism) by subtracting the updated kinetic energy from the updated total energy.As previously indicated, the gas velocity  and radiation flux  we use in the inner stage of the iteration are lagged.This can cause significant inaccuracies at high optical depths when the velocitydependent terms are non-negligible.To eliminate errors like this and render this scheme fully implicit, we now repeat the inner iteration using the updated values of  and , and compute new estimates for  and , repeating this procedure until either the relative change in the value of the terms proportional to  and  in  (Equation 36) from one outer iteration to the next is below 10 −13 or the absolute change is below 10 −13 .Except in the dynamic diffusion limit, where the velocity-dependent terms are at the same order as all other terms, this iterative process typically terminates after just one iteration, and in all the tests we present below, and for all the test problems presented in Wibking & Krumholz (2022), we never require more than a handful of outer iterations.Thus the cost is modest.This completes accounting for all terms in the radiation four-force, thus completing radiation-matter coupling.

PROPERTIES OF THE SCHEME IN THE DIFFUSION LIMIT
Before proceeding to numerical tests of the scheme we have described, we first present an analysis of its behaviour in the so-called asymptotic diffusion regime, where the photon mean free path is not resolved by the computational grid, in order to demonstrate directly why it succeeds in capturing this limit while other schemes fail.

The asymptotic diffusion limit of the discrete IMEX equations
We begin our analysis by recalling the results from Section 2.2, which are that in the static diffusion limit for grey material with χ0 = χ0 = χ0 = 1, to leading order the radiation energy in non-dimensionalised variables is and the radiation flux is Inserting these limits into the evolution equations for matter and radiation energy that we solve during the radiation update (Equation 13; i.e., omitting changes in the matter energy due to fluid processes), we have and thus the evolution equation for the total matter plus radiation energy reduces to the usual radiation-diffusion form,5 For a numerical method for thermal radiative transfer to be "asymptotic preserving", it must give a valid discretization of Equation 49and enforce the conditions in Equation 45and Equation 46 when C ≫ L ≫ 1.With an asymptotic preserving method, it is possible to use cells that are optically thick and still obtain accurate solutions of radiative transfer.
To verify that our IMEX scheme satisfies this condition, we begin by writing down the two steps of the IMEX update for the radiation energy (Equation 28 and Equation 29), again using nondimensionalised variables and assuming grey material.For simplicity we will adopt v = 0 and ĉ =  as well.This gives followed by From Equation 50we can solve for  (+1/2) and expand to first order in L −1 to obtain Similarly, performing the same operation on Equation 51, one can show to leading order in L −1 that We have therefore established the our scheme enforces Equation 45at both stages of the IMEX update.
We next examine the two IMEX stages for the radiation flux,6 Equation 54 implies and if we again take the limit C ≫ L ≫ 1, to leading order we have Applying the same procedure to Equation 55 yields leading-order terms Thus the leading term of the radiation flux reduces the form given by Equation 46, simply lagged by a half-step.This nonetheless means that our scheme obeys this constraint.Finally, we write down the two IMEX stages for the matter energy update Adding Equation 50to Equation 59 and using Equation 57 we get where the superscript ( − 1/2) indicates the intermediate state of the previous time step.Similarly, adding Equation 51to Equation 60and using Equation 58gives The combination of Equation 61and Equation 62 represents a valid discretisation of the diffusion equation Equation 49, albeit one using a temperature that is lagged by half a step, thus proving our IMEX PD-ARS scheme preserves the asymptotic diffusion limit.
It is worth noting that schemes such as the one described by Equation 61 and Equation 62, while they represent valid discretisations of the diffusion equation, can be unstable depending on how the spatial derivatives are evaluated.In particular, Radice et al. (2018) point out that, in a scheme where computation of F in the diffusion limit effectively reduces to computing a centred difference on Ê, and in turn evaluating ∇ • F in the equations above reduces to evaluating a centred difference on F, the resulting scheme has the property that the solutions in even-and odd-numbered cells are decoupled, i.e., the solution in cell  depends only on the states in cells  − 2 and  + 2, not  − 1 and  + 1; this in turn can give rise to an even-odd instability where numerical oscillations with a period of two cells are not damped and can grow large.Radice et al. propose a method to suppress this instability.In our tests, while we observe faint hints of the instability in our scheme, these appear only in the dynamic diffusion regime and only in tests at very low resolution (e.g., < 64 cells per linear dimension).In all other problems the instability, if it exists at all, is imperceptibly small.We therefore do not use Radice et al.'s correction for any of the tests we present below.However, we have implemented it in Quokka, and allow users to enable it via a compile-time option should it prove useful at some point in the future.

Comparison to schemes with purely explicit intermediate stages
It is instructive at this point to repeat the analysis we have just performed for the IMEX discretisation for the RK2-SSP discretisation used in Wibking & Krumholz (2022), since this will let us see why this scheme does not successfully capture the asymptotic diffusion regime.While our analysis will be specific to this particular time stepping approach, we will see that the results generalise straightforwardly to any scheme where the intermediate stage is fully explicit and includes only the transport terms, and thus to other schemes such as those proposed by Rosdahl & Teyssier (2015) and Skinner et al. (2019).The RK2-SSP scheme uses the time update  (+1/2) =  () + Δ (  () , ) (63) (  () , ) + 1 2 (  (+1/2) ,  + Δ) Compared to Equation 28and Equation 29, the primary difference is in the treatment of the fast terms  ; in the IMEX scheme these appear at both stages of the update, while in the RK2-SSP scheme they appear only at the final stage.This will prove to be crucial in what follows.
Again adopting the limit C ≫ L ≫ 1 and setting P = (1/3)I, as appropriate for static diffusion, we can write down the leading-order terms in the two stages of this update for the radiation energy as For the radiation flux, we obtain Finally, if we write down the update for the gas energy (which is non-trivial only for the second stage, since  = 0 for the gas energy) and add it to that for the radiation energy, Substituting Equation 67and Equation 68 into Equation 69, we obtain Thus we see that there are two modes of radiation diffusion in this numerical scheme: one is physical diffusion with a diffusion coefficient of C/3L, and the other is numerical diffusion with a coefficient of C 2 Δt/3.This numerical mode will dominate for any time step Δt ≳ 1/CL, or, in dimensional terms, Δ >  ∞ /, i.e., whenever the time step is large enough that we do not resolve the light-crossing time of a photon mean free path.In the asymptotic diffusion regime, where the photon mean free path is smaller than the size of a cell, this means that numerical diffusion dominates any time that the time step is larger than the light crossing time of a cell.In practice, the time step is always much larger than this -since otherwise one might as well use a fully explicit scheme -which explains why discretisations such as RK2-SSP fail in the asymptotic diffusion regime.
Comparing Equation 70to Equation 62, we see that the numerical diffusion mode is removed in the IMEX discretisation, and by comparing the calculations leading up to these equations we can also understand why.The numerical diffusion term in the RK2-SSP update originates in a term that appears at the intermediate stage of the flux update (Equation 67).This term does not appear to leading order in the intermediate stage of the IMEX update (Equation 57) because it is overwhelmed by the source term, which is a factor of L larger; it is a cancellation within the source term that forces the radiation flux to the correct value for diffusion.Thus the IMEX update winds up with an estimate for the intermediate-state flux that is of order 1/L independent of the time step, while the RK2-SSP update, because it ignores the order CL source term but retains the order C transport term during the intermediate stage, obtains a flux estimate that is of order CΔt instead.This overestimate is what gives rise to the artificial numerical diffusion of energy.An important conclusion to draw from this argument is that the failing in the RK2-SSP scheme for this problem is not specific to that scheme, but is instead generic to any update scheme that contains a stage that includes only the transport terms and not the source terms.Such a scheme will always overestimate the transport in diffusion-regime problems where the source term is responsible for suppressing them.

On modifications to the radiation wave speed
The problem we have identified in the RK2-SSP scheme is not new; indeed, multiple authors have pointed out that the HLL solver applied to the radiation subsystem, in a scheme where the source and transport terms are operator-split, yields fluxes that fail to preserve the asymptotic diffusion limit (Lowrie & Morel 2001;Audit et al. 2002;Jiang et al. 2013;Skinner et al. 2019;Wibking & Krumholz 2022).In an attempt to circumvent this problem, these authors have proposed a range of corrections to the energy fluxes; for instance, Skinner et al. (2019) suggest where  , are the characteristic left and right wavespeeds, respectively, whose values absent the RSLA are given by  = ± √︁  , where  is the component of the Eddington tensor along that direction.In this expression,  is an empirical correction factor that smoothly transitions from 1 in the streaming limit to 1/ in the optically thick limit, where  is an empirical estimate of the optical depth, usually computed as the optical depth across a handful of computational cells.This correction implies an effective wave speed of In light of the preceding discussion, we can see that this reduction in the flux is effectively a correction that attempts to reduce the numerical diffusion mode in Equation 70, by in turn forcing the overestimated intermediate time-flux (Equation 67) back toward the solution that would have been obtained by retaining rather than dropping the source term.However, the accuracy of this fix, particularly in the regime of intermediate optical depth, is not known.
In our IMEX PD-ARS method, we use the wavespeed  without any correction in the Riemann solver.No correction is necessary because, by retaining the source term, we automatically recover the correct flux in the diffusion limit (Equation 57), and our scheme correctly and smoothly goes between the optically thin limit, where the transport term dominates, and the optically thick limit, where the source term does.Indeed, a corollary of this analysis is that, in the highly optically thick regime where the source term is dominant, we need not even obey the CFL condition for the radiation in order to obtain the correct answer.To understand why this is the case, we notice that Equation 62 is an effective temporal discretisation of the diffusion equation which has a diffusion coefficient of Ddiff = C/(3L).The effective speed of the diffusion in time interval t is approximately To better understand the magnitude of Ŝdiff , we reintroduce dimensions to the variables and get where CFL ≡ /Δ is the radiation CFL number.This evaluation of diffusion speed aligns with the wave speed correction 'hack' (Equation 72) when CFL = 1.Our IMEX scheme's asymptotic preserving nature, as demonstrated, allows the CFL number to exceed 1 while still accurately capturing the diffusion limit.This is exemplified in the subsequent discussion on the asymptotic Marshak wave test (Section 5.2), where our numeric scheme accurately catches the position of the diffusion front for radiation CFL numbers up to ∼ 10.This observation is in concordance with Jiang et al. (2013), who, through a linear analysis of the diffusion equations, showed that if one calculates the effective propagation speed with a wavelength of 10 cells, the numerical diffusion is small enough not to affect the solution.This compatibility with the findings from linear wave analysis of diffusion equations not only validates our scheme's robustness in limiting numerical diffusion effectively but also underscores the versatility of the IMEX PD-ARS scheme in accommodating higher CFL numbers without compromising the asymptotic preserving property.

NUMERICAL TESTS
The original Quokka paper (Wibking & Krumholz 2022) introduced a series of tests to validate the accuracy and convergence capabilities of Quokka, including tests of the hydrodynamic subsystem, the radiation transport, and coupled RHD.This new scheme passes all these tests, so we will not repeat them here.We instead introduce additional tests for our implementation of the IMEX PD-ARS scheme in Quokka that are specifically designed to test the code's ability to preserve the asymptotic diffusion limit.The full source code and outputs for all tests are available in the Quokka github repository (see Data Availability statement for details).

Non-equilibrium Radiation Shock
We begin our tests with the classical grey non-equilibrium radiative shock test described by Zel'dovich & Raizer (2012).Radiation can modify the structure of a shock because it diffuses and interacts with matter.Lowrie & Edwards (2008) have found a semi-analytic exact solution of radiative shocks with grey -non-equilibrium diffusion.
Using their parameters with an upstream Mach number of M = 3, we obtain a subcritical radiation shock whose temperature jumps discontinuously at the shock interface.We present our numerical calculation of this problem and compare them to the solutions of Lowrie & Edwards (2008).opacities to  = 577 cm 2 g −1 and mean molecular weight to  =  H .We use an adiabatic EOS with an adiabatic index to  = 5/3.The shock is simulated in a 1D region with  ∈ [0, 0.01575] cm resolved with 512 grids.The problem is set in the rest frame of the shock initialised at  = 0, and the conditions on the left and right sides of the shock are uniform, with densities, temperatures, and velocities given by   = 5.69g cm −3 ,   = 2.18 × 10 6 ,   = 5.19 × 10 7 cm s −1 , and   = 17.1 g cm −3 ,   = 7.98 × 10 6 ,   = 1.73 × 10 7 cm s −1 , respectively.In order to exactly match the assumptions used in the semi-analytic solution, we use the Eddington approximation, P = (1/3)I, to calculate the radiation pressure tensor.Following Skinner et al. ( 2019), we use a reduced speed of light ĉ = 10(  +  , ), where  , is the adiabatic sound speed of the left-side state.We use a CFL number of 0.4 and evolve until  = 10 −9 s.We show the resulting temperature profile in Figure 1.The agreement between our numerical calculation and the semi-analytic solution is excellent.The  1 norm of the relative error of the gas temperature is 0.38 per cent, which is as good as the solution of Skinner et al. (2019).

The Asymptotic Marshak wave problem
We test the code's ability to accurately capture radiation in the optically thick limit via the diffusive Marshak-wave problem.This problem consists of a semi-infinite medium of material with a variable absorption coefficient  = 300 cm −1 /(  /keV) 3 .We input incoming isotropic radiation on the boundary at  = 0 from a 1 keV temperature source.We use an initial condition of  () =  0 = 10 −3 keV and  () =    4 0 for our test problem.Figure 2 shows our numerical results compared to the semianalytic solution for the diffusive Marshak-wave problem, alongside a comparison to the solution using the SSP-RK2 scheme.The numerical calculation uses a spatial grid of 60 cells across the domain (0, 0.66 cm), with each cell spanning 3 × 10 9 mean-free paths at the minimum temperature and 3 mean-free paths at the high-est temperature.The heat capacity in this problem is constant at 3 × 10 15 erg cm −3 keV −1 .We compare this solution with the semianalytic equilibrium-diffusion solution from Zel 'dovich & Raizer (2002).We find a relative L1-norm error at about 4.5%.The simulation results presented here are obtained using a CFL number of 0.9, but the simulation runs with a CFL number up to ∼ 10, implying a time step roughly 10 times the largest permissible time step in the streaming limit.

Advecting radiation pulse in the static diffusion limit
In this test, introduced by Krumholz et al. (2007), we simulate the advection of a radiation pulse by the gas motion in the diffusion regime.We set the initial gas and radiation temperatures equal in this test, and the initial temperature and density profiles are with  0 = 10 7 K,  1 = 2 × 10 7 K,  0 = 1.2 g cm −3 ,  = 24 cm, and  = 2.33  = 3.9 × 10 −24 g.The radiation pressure is estimated via the Eddington approximation.The opacity of the gas is set at  0 =  0 =  0 = 100 cm 2 g −1 .The system is initially in pressure balance.If there were no radiation diffusion, the system would be in an equilibrium between the gas pressure and radiation pressure.Because of radiation diffusion, the balance is lost and the gas moves.
We solve the problem numerically in two different frames, one in the lab frame and the other in the comoving frame, and compare the results in Figure 3.In the comoving frame run, the velocity is 0 in the beginning everywhere.In the lab frame run, the initial velocity is  0 = 10 6 cm s −1 .In both runs, the optical depth across the pulse is  =  = 2.9×10 3 , and in the lab-frame run  ≡  0 / = 3.3×10 −5 , giving  ≈ 0.1 and placing this problem in the static diffusion limit.The computational domain in both runs is a 1D region of (−512 cm, 512 cm) with periodic boundaries, and the grid consists of 512 uniform cells.We show the density, temperature, and velocity profiles from both runs at  = 4.8 × 10 −5 s in Figure 3.The results of the lab-frame run have been shifted in space by  0  for comparison.The density, temperature, and velocity profiles from both runs are almost identical, demonstrating the accuracy of our scheme in capturing radiation advection.We also show the relative difference from both runs in Figure 4.The relative difference between the non-advecting and advecting cases is below 0.03% anywhere.The relative difference between the temperature of the gas and radiation is below 3 × 10 −9 .We compare these results with that of Krumholz et al. (2007) and Zhang et al. ( 2011) and find good agreement.

Advecting radiation pulse in the dynamic diffusion limit
In this test, we rerun the advecting radiation pulse problem in the dynamic diffusion regime.Compared to the static diffusion case, the advection speed is increased to 3 × 10 7 cm s −1 and the opacity is increased to 500 cm 2 g −1 .Numerically,  = 10 −3 ,  = 1.4 × 10 4 , and  = 14, placing this problem in the dynamic diffusion limit.We increase the number of cells to 1024 to reduce the magnitude of odd-even decoupling instability, and evolve the simulation to  = 4.8 × 10 −5 s.The relative difference between the non-advecting and advecting cases is below 0.06 per cent anywhere (Figure 5), demonstrating the excellent accuracy of our scheme in capturing CFL number = 0.9 exact numerical Figure 2. Comparing numerical solutions to the Marshak-wave problem in the asymptotic diffusion limit from the SSP-RK2 scheme (left) and IMEX PD-ARS scheme we introduce here (right).The solid lines are the semi-analytic solutions and the dots are the nodal values of the numerical solution.Both simulations run with Δ = 0.011 cm and CFL = 0.9.The IMEX scheme accurately captures the diffusion of radiation while in the SSP-RK2 scheme solution, the wavefront propagates about two times too fast.radiation advection in the dynamic diffusion regime.Although there is no analytic solution to this problem, we can estimate the expected pulse width via the diffusion equation The expected distance that the pulse has diffused at time  is approximately The temperature profile of our numerical calculation Figure 5 agrees well with this analytic estimation.

Advecting uniform medium in the dynamic diffusion regime
We test the code's ability to incorporate relativistic corrections up to order / or order  2 / 2 with a 1-D test problem of advecting a uniform medium in the dynamic diffusion regime.In this test, the gas is initially moving at a uniform velocity of  0 = 0.01 and has a uniform density.We configure the opacity and grid size such that each grid's optical depth is 10 5 .The radiation is initially in thermal equilibrium with the medium at the comoving frame.
In the lab frame, to order / the initial radiation quantities are and we can show that the space-like component of the radiation four-force, Equation 4, to order / vanishes: To order  2 / 2 the initial radiation quantities are and we can show that both the time-like and space-like component of the radiation four-force vanishes: The key feature of this test is that, in the lab frame, the radiation energy is out of equilibrium with the thermal radiation.However, this excess is offset by the work done to the radiation by the matter.One can show that our backwards Euler scheme, Equation 42 and Equation 43, preserves this equilibrium, leading to an equilibrium state ( (+1) =  () and  (+1) =  () ) that aligns precisely with the expected outcomes of this test (Figure 6), achieving precision up to machine accuracy even in the dynamic diffusion limit (Table 1).Conversely, the source term from Howell & Greenough (2003) and Wibking & Krumholz (2022), defined as ∝ (4/ −   ), fails to attain this balance.

CONCLUSION
We have presented a novel mixed implicit-explicit (IMEX) time integration scheme for finite-volume RHD that is second-order accurate in the streaming limit and accurately preserves the asymptotic diffusion limit.Our method uniquely combines the robustness of local implicit methods in handling stiff source terms with the scalability and simplicity of explicit methods for advective transport, making it particularly suitable for adaptive mesh refinement and massively parallel, GPU-accelerated architectures.This scheme addresses a critical gap in current RHD solvers by ensuring asymptotic accuracy in the streaming, static diffusion, and dynamic diffusion limits, including in the asymptotic diffusion regime where the photon mean free path is much smaller than the computational grid, without necessitating non-local implicit steps or ad hoc adjustments to the radiation flux based on the optical depth.We have implemented our algorithm in the GPU-accelerated AMR RHD code Quokka (Wibking & Krumholz 2022).Quokka is open-source and a link to download the source code is in the DATA AVAILABILITY section.
We have verified our algorithm using a variety of established quantitative tests, including the non-equilibrium radiation shock test, the asymptotic Marshak wave test, and the advecting radiation pulse test in the static and dynamic diffusion limit.These tests demonstrate the scheme's capability to recover the correct asymptotic limits.
Looking forward, we envision several areas for further development.Extending our scheme to include more complex radiation transport models, such as multi-frequency radiation transfer and nonthermal emission, could enhance its applicability to a broader range of astrophysical problems.Additionally, integrating our RHD method time ct = 1000 radiation (numerical) radiation (exact) gas (numerical) gas (exact) Figure 6.Advecting uniform medium test in the dynamic diffusion limit ( = 0.01,   = 10 3 ).The numerical results match the analytic solution to the limits of machine precision.

Figure 1 .
Figure 1.Comparing Radiation and gas temperature from the numerical calculation (solid lines) with the exact steady-state solution (dots) in a subscritical radiation shock with  = 3.

Figure 3 .
Figure 3. Advecting radiation pulse tests in the static diffusion regime at  = 4.8 × 10 −5 s using the RK2-SSP scheme (left) versus our newly introduced IMEX scheme (right).The IMEX scheme accurately captures the radiation diffusion, agreeing well with the results presented in Krumholz et al. (2007) and Zhang et al. (2011).

Figure 4 .Figure 5 .
Figure 4. Top: Relative errors in density (dashed lines) and temperature (solid lines) between the advected and unadvected runs.Bottom: Relative difference between gas and radiation temperature in the advected run.

Table 1 .
Relative errors of the matter temperature in numerical calculation with respect to theoretical solution in the advecting uniform medium tests.Col. 3: the / order used in the source terms.