© ESO 2022

1. Introduction

In many astrophysical environments, radiation plays an important role in the system’s total energy or momentum budget. Some selected examples involve: solar convection simulations where radiation controls the heating and cooling of the photosphere (Stein & Nordlund 1998), wind outflows from massive stars where a strong radiation force lifts material off the stellar surface (Castor et al. 1975), evacuation of massive star-forming discs where radiative ablation may control the stellar upper mass limit (Kee & Kuiper 2019), accretion flows around (Jiang et al. 2019), and disc winds from (Proga & Kallman 2004) supermassive black holes, where radiation often is the dominant player controlling the energetics and dynamics.

The implementation of suitable radiation modules for dynamical simulations, however, it is a very difficult task in general (Mihalas & Mihalas 1984; Castor 2004) and, thus, a rather broad range of different approaches have been taken, depending on the specific application. For example, in the aforementioned convection simulations of solar-type stars, the heating and cooling terms are computed directly from solutions to the time-independent radiative transfer equation (RTE), using a sort of “frozen-in” approach (e.g., Stein & Nordlund 1998), and the effect on the momentum balance is neglected. Similar considerations apply for solar corona and wind simulations, however, here the radiative cooling term is often approximated by assuming an optically thin approximation (e.g., Schure et al. 2009). On the other hand, the winds of hot, massive stars are driven by momentum transfer from the radiation field to the gas. For such simulations, the force due to spectral lines is critical, and very elaborate radiation transport methods have been developed in the meantime for steady outflows (Sander et al. 2017; Sundqvist et al. 2019), whereas time-dependent simulations typically rely on various distribution-function approaches (Owocki et al. 1988; Sundqvist et al. 2018).

For general attempts to solve the full time-dependent radiation momentum and energy equations, two important complicating factors are: (i) how to obtain a suitable closure relation for the radiation equations and (ii) how to compute the radiation-matter interaction terms (the opacities) in a supersonic flow. Concerning the closure relation, it is possible to attempt to obtain closure by a variable Eddington tensor (VET) computed directly from the radiative transfer equation and using, for example, a short characteristics scheme (Jiang et al. 2012). Alternatively, various analytic closure relations are often applied, such as flux limited diffusion (FLD, see also Turner & Stone 2001; Krumholz et al. 2006), first moment (M1; Skinner & Ostriker 2013; Bloch et al. 2021), or a combination of both (Mignon-Risse et al. 2020). Concerning the radiation-matter interaction terms, the central issue is that typically the opacities are isotropic only in the frame that is co-moving with the fluid (the CMF). As such, it is only in this frame that suitable mean opacities can be defined without performing complicated double integrals of the specific intensity and opacity over frequency and angle. On the other hand, the radiative transfer equation itself is significantly more complicated in the CMF than in the laboratory frame (see discussions in Mihalas & Mihalas 1984; Castor 2004). As such, some trade-offs must typically be done, where one popular method is to compute the radiation variables in the laboratory frame and then via a first order expansion still treat the opacities in the CMF (Mihalas & Klein 1982; Lowrie et al. 1999; Jiang et al. 2012). A main drawback with this “mixed-frame” formulation, however, is that it is not suitable for computing line opacities in a supersonic flow, since the first-order expansion upon which it is based does not apply for the rapidly changing opacities of spectral lines.

In brief, MPI-AMRVAC¹ is a modern, MPI-parallelized computer code, aimed at solving partial differential equations (PDEs) on an adaptive, block structured quadtree/octree mesh (Xia et al. 2018; Keppens et al. 2021). Originally, the code focused on the equations of (magneto-) hydrodynamics, but recently the code has been expanded to solve general hyperbolic, parabolic, and elliptic PDEs. The code uses a finite volume method to solve hyperbolic advection equations from one to three dimensions and includes a large number of schemes and limiters for flux reconstruction, time discretization methods, and refinement strategies. In recent years, multiple new physics modules have been added to MPI-AMRVAC, such as descriptions of viscosity, conduction, and dust dynamics (Porth et al. 2014; Xia et al. 2018). One aspect that has been missing is a general description of the dynamical effects of a radiation field. Thus far, only the effects of optically thin radiative cooling have been implemented (Schure et al. 2009; van Marle & Keppens 2011). In this paper, we take a first step toward a more general description, implementing a non-equilibrium FLD method into MPI-AMRVAC. For many of our targeted applications line-opacity is crucial; as such we formulate the FLD equations in the CMF (e.g., Turner & Stone 2001) and also present a simplified way (based on Poniatowski et al. 2021) of accounting for supersonic line-opacities within the formalism. Recently, MPI-AMRVAC has been expanded to include the possibility to solve elliptic PDEs (Keppens et al. 2021) using a geometric multigrid library (Teunissen & Keppens 2019). This has been crucial for a proper treatment of the diffusive term in the FLD equations, ensuring a stable solution by applying an implicit method.

The paper is structured as follows: In Sect. 2, we introduce the radiative-hydrodynamics (RHD) equations and we describe the non-equilibrium FLD method for obtaining closure. The different aspects of the implementation in the MPI-AMRVAC code are explained in Sect. 3. Section 4 is devoted to benchmark tests for the newly developed code, and in Sect. 5, we present a first research application in the form of an optically thick radiation-driven wind outflow from a classical Wolf-Rayet star. Section 6 summarizes and discusses the paper as well as the general outlook.

2. RHD equations

Including source terms due to radiation, the hydrodynamical conservation equations of mass, momentum, and energy are:

$$\begin{array}{cc}& {\partial }_{\mathrm{t}}\rho +\mathrm{\nabla }·(\rho \mathit{v})=0,\hfill \end{array}$$(1)

$$\begin{array}{cc}& {\partial }_{\mathrm{t}}(\rho \mathit{v})+\mathrm{\nabla }·(\mathit{v}\rho \mathit{v}+p\mathit{I})={\mathit{f}}_{\mathrm{r}},\hfill \end{array}$$(2)

$$\begin{array}{cc}& {\partial }_{\mathrm{t}}e+\mathrm{\nabla }·(e\mathit{v}+p\mathit{v})=\mathit{v}·{\mathit{f}}_{\mathrm{r}}+\dot{q},\hfill \end{array}$$(3)

where ρ is the gas density, v is the gas velocity, e the total gas energy density (internal plus kinetic), p is the gas pressure, and I is the identity matrix. On the right-hand side, f_r is the radiation force density exerted on the matter and $\dot{q}$ the heating or cooling term due to the radiation. These equations are supplemented by the expression for gas energy, e, and the ideal gas law:

$$\begin{array}{cc}& e=\frac{p}{(\gamma -1)}+\frac{\rho v^2}{2},\hfill \end{array}$$(4)

$$\begin{array}{cc}& p=\frac{k_{\mathrm{B}}{T}_{\mathrm{g}}}{m_{\mathrm{p}}\mu }\rho .\hfill \end{array}$$(5)

Here, γ is the adiabatic index, k_B is the Boltzmann constant, T_g is the gas temperature, m_p is the proton mass, and μ is the mean molecular weight of the gas particles. As described below, the terms f_r and $\dot{q}$ can be computed by considering the energy and momentum equations for the radiation. We note further that the above equations neglect the effects of gravity and magnetic fields. However, the modular structure of MPI-AMRVAC enables to seamlessly integrate the new FLD module described in this paper in the main branch of the code in order to use it in physical situations that account for more physics.

2.1. Radiation energy and momentum equations

Conservative equations can be derived for the radiation field by taking angular moments of the radiative transfer equation for specific intensity, I_ν. However, when describing the transport of radiation using such moment equations, particular attention has to be paid to the reference frame of the radiation quantities. In an outside observer’s reference frame (equivalent to the Eulerian or laboratory frame of reference), the computation of $\dot{q}$ and f_r requires radiation-material interaction terms, where these involve an important angle-dependence induced by the Doppler shift. For moving fluids, this makes it necessary to always carry out double integrals over frequency and angle in order to obtain the correct coupling terms. These issues are avoided by transforming to a frame co-moving with the local velocity of the fluid. In this co-moving frame (CMF), the frequency-integrated radiation energy and momentum equations are (Mihalas & Mihalas 1984; Castor 2004):

$$\begin{array}{cc}& {\partial }_{\mathrm{t}}E+\mathrm{\nabla }·(E\mathit{v})+\mathrm{\nabla }·\mathit{F}+\mathit{P}:\mathrm{\nabla }\mathit{v}=-\dot{q},\hfill \end{array}$$(6)

$$\begin{array}{cc}& \frac{1}{c^2}({\partial }_{\mathrm{t}}\mathit{F}+\mathrm{\nabla }·(\mathit{F}\mathit{v}))+\mathrm{\nabla }·\mathit{P}=-{\mathit{f}}_{\mathrm{r}},\hfill \end{array}$$(7)

where E, F, and P are the frequency-integrated radiation energy density, flux vector, and pressure tensor, respectively, evaluated in the CMF. These radiation equations are now somewhat more complex than in the observer’s frame, as they contain extra terms stemming from the transformation (see Chapter 6 in Castor 2004 for the full transformation properties between the observer and co-moving frames). For example, the fourth term on the left-hand side of the energy Eq. (6) describes the dyadic product between the radiation pressure tensor and the gradient of the gas velocity vector. This term can be interpreted physically as the energy that leaves the radiation system when the radiation field accelerates and provides work, that is, it is the radiation work term. In optically thick supersonic media, this physical effect can become critical as it may diminish the radiation flux, and so also the radiation force, significantly. As such, it is sometimes called “photon tiring” (Owocki & Gayley 1997). On the other hand, all material interaction terms in the CMF are in most situations isotropic, which simplifies tremendously the evaluation of the right-hand-sides in Eqs. (6) and (7). These coupling terms in Eqs. (2), (3), (6) and (7) can now be written as:

$$\begin{array}{cc}& \dot{q}=c{\kappa }_{\mathrm{E}}\rho E-4\pi {\kappa }_{\mathrm{B}}\rho B,\hfill \end{array}$$(8)

$$\begin{array}{cc}& {\mathit{f}}_{\mathrm{r}}=\frac{\rho {\kappa }_{\mathrm{F}}\mathit{F}}{c},\hfill \end{array}$$(9)

where $B={\int }_0^{\infty }{B}_{\nu }\text{d}\nu =\sigma T_{\text{g}}^4/\pi$ is the frequency-integrated Planck function, with σ the Stefan-Boltzmann constant, and the opacities κ_B and κ_E are mass absorption coefficients measured in cm² g⁻¹. Specifically, Eqs. (8) and (9) involve the Planck, energy density, and flux mean opacities:

$$\begin{array}{cc}& {\kappa }_{\mathrm{P}}\equiv \frac{\int B_{\nu }{\kappa }_{\nu }\mathrm{d}\nu }{B},\hfill \end{array}$$(10)

$$\begin{array}{cc}& {\kappa }_{\mathrm{E}}\equiv \frac{\int E_{\nu }{\kappa }_{\nu }\mathrm{d}\nu }{E},\hfill \end{array}$$(11)

$$\begin{array}{cc}& {\kappa }_{\mathrm{F},i}\equiv \frac{\int F_{\nu ,i}{\kappa }_{\nu }\mathrm{d}\nu }{F_i},\hfill \end{array}$$(12)

which follows directly from considering the frequency-dependent form of the coupling terms described by Eqs. (8) and (9). Here, F_i is the ith component of the flux vector. For the sake of simplicity, we assume an isotropic flux mean opacity throughout this paper, although the method can be expanded for anisotropic opacities as well. The quantities with subscript _ν are the frequency dependent radiation quantities (e.g., E_ν is the radiation energy density at frequency ν). We note that the above formulation does not necessarily mean the source functions must be Planckian. Indeed, the same expressions are also found for a model where the emission and extinction coefficients have both thermal absorption (“a”) and scattering (“s”) contributions, that is, when ${\kappa }_{\nu }={\kappa }_{\nu }^{\text{a}}+{\kappa }_{\nu }^{\text{s}}$. In such a situation, however, it is critical to keep in mind that although the flux mean then involves the total opacity κ_ν, the energy, and Planck means should be evaluated using only the thermal absorption part ${\kappa }_{\nu }^{\text{a}}$ (see Eq. (77) in Mihalas & Mihalas 1984, their page 336, and also the corresponding discussion on their page 472).

2.2. Non-equilibrium FLD closure relation

An additional relation5pt between P and E is needed to close the radiation moment Eqs. (6) and (7). In general, this relation must be obtained from full solutions of the frequency and angle dependent transfer equation in different directions. However, realistic multi-frequency solutions to the radiative transfer equation in the CMF have thus far only been developed for 1D, steady-state media with a monotonic velocity field (e.g., Hillier & Miller 1998; Puls et al. 2020). As such, analytic closure relations are often being used in practical radiation-hydrodynamics applications. Typically, these analytic relations recover the correct equilibrium limit in the optically thick limit and then apply some appropriate “bridging law” for extension into the opposite optically thin streaming limit.

In this paper, we apply the so-called FLD approximation as our closure relation. Neglecting the first two terms in the radiation momentum equation, we have:

$$\begin{array}{c}\hfill \mathrm{\nabla }·\mathit{P}=\mathrm{\nabla }·(\mathit{f}E)=-\frac{{\kappa }_{\mathrm{F}}\rho }{c}\mathit{F},\end{array}$$(13)

for the Eddington tensor f = P/E. This invites us to write the radiation flux from Fick’s diffusion law (see also Levermore & Pomraning 1981):

$$\begin{array}{c}\hfill \mathit{F}=-D\mathrm{\nabla }E,\end{array}$$(14)

with D a diffusion constant, which will depend on the local state of the gas and radiation field. It is important to again note here that this diffusion approximation is only applicable for the CMF quantities of the radiation flux and energy density (Mihalas & Mihalas 1984; Castor 2004).

In the limit of radiative diffusion, the energy density and pressure take their equilibrium values such that the scalar Eddington factor f = 1/3 and the Eddington tensor becomes f = fI with I the unit tensor. This is valid for very optically thick regions, where the photon mean-free path ℓ = 1/(ρκ_F) is small compared to the typical length scales over which the state variables vary. Letting this length scale be the radiation energy density scale height H_R = E/|∇E|, we require ℓ ≪ H_R in the radiative diffusion limit. This then results in a diffusion constant D = c/(3κ_Fρ). A basic issue with this diffusion approximation, however, is that the flux computed from it may exceed the physical limit |F|=cE in the opposite regime of freely streaming photons, where the diffusion constant will approach an arbitrarily large value as the local density and opacity approach zero. This suggests to introduce a bridging law that limits the flux in optically thin regions, while still recovering the optically thick limit. To this end, we apply the flux-limiter λ suggested by Levermore & Pomraning (1981) (see also Turner & Stone 2001), expressed as:

$$\begin{array}{cc}& f=\lambda +{\lambda }^2{R}^2,\hfill \end{array}$$(15)

$$\begin{array}{cc}& R=\frac{\ell }{H_{\mathrm{R}}}=\frac{|\mathrm{\nabla }E|}{\rho {\kappa }_{\mathrm{F}}E},\hfill \end{array}$$(16)

and for the bridging law:

$$\begin{array}{c}\hfill \lambda =\frac{2+R}{6+3R+R^2}·\end{array}$$(17)

Clearly this relation recovers the diffusion limit since ℓ ≪ H_R yields directly λ → 1/3 and so, f → 1/3. Similarly in the opposite free streaming limit λ → 1/R such that f → 1. We note that while we use this prescription throughout this paper, other variants have been suggested as well (e.g., Minerbo 1978); in the MPI-AMRVAC code, the user can readily switch between different flux limiters like the Levermore (Levermore & Pomraning 1981) or the Minerbo (Minerbo 1978) prescriptions.

Using this flux-limiter within Fick’s diffusion formulation, we can locally compute the co-moving radiation flux according to:

$$\begin{array}{c}\hfill \mathit{F}=\frac{-c\lambda }{{\kappa }_{\mathrm{F}}\rho }\mathrm{\nabla }E.\end{array}$$(18)

Here, it can be seen that in the thin limit, when ℓ ≫ H_R, λ → 1/R, consequently |F|→cE and causality is preserved. Thus, using the corresponding Eddington tensor suggested by Turner & Stone (2001), we have:

$$\begin{array}{c}\hfill \mathit{f}=\frac{1}{2}(1-f)\mathit{I}+\frac{1}{2}(3f-1)\widehat{\mathit{n}}\widehat{\mathit{n}}.\end{array}$$(19)

We also obtain the radiation pressure tensor and so we can omit Eq. (7) entirely. Here, $\widehat{\mathit{n}}=\mathrm{\nabla }E/|\mathrm{\nabla }E|$ is the unit vector in the direction of the gradient of the radiation energy density (i.e., that of the radiative flux). Thus, the only PDE that needs to be solved in order to advance the radiation subsystem is Eq. (6); this radiation energy equation describes the conservation of energy stored in the radiation field.

In summary, the radiation energy density, E, is integrated over time, the radiation flux, F, is then calculated locally using the analytic FLD approximation (18), and the pressure tensor, P, obtained locally from E using the analytic prescription for the Eddington tensor (19).

For the radiative heating-cooling term $\dot{q}$, we further assume for the specific application in this paper that the Planck and energy density mean opacities are equal, κ_P = κ_E = κ, such that:

$$\begin{array}{c}\hfill \dot{q}=\rho \kappa 4\sigma (T_{\mathrm{r}}^4-T_{\mathrm{g}}^4),\end{array}$$(20)

where we have reformulated E in terms of a radiation temperature, $E\equiv a_{\text{r}}{T}_{\text{r}}^4$, where a_r = 4σ/c is the radiation constant. In the non-equilibrium FLD method applied here, T_r can deviate from the gas temperature T_g, allowing for situations where the coupling between radiation and gas may be out of equilibrium.

3. Numerical implementation

The RHD system of PDEs in the non-equilibrium FLD approximation consist of Eqs. (1)–(3) and (6). The conservative left-hand side parts of the hydrodynamic Eqs. (1)–(3) are entirely hyperbolic and can be solved using the existing variety of shock capturing, high-resolution finite volume solvers in MPI-AMRVAC (Xia et al. 2018). However, when taking into account the coupling between the radiation field and its effects on the gas quantities, the system loses its purely hyperbolic property. This coupled system is now solved in an operator split manner, where different terms are added using different schemes, depending on the timescales of their effect. The advection term, ∇ ⋅ (Ev) in Eq. (6), is handled using the same solvers already available in MPI-AMRVAC for solving hyperbolic equations. Unless stated otherwise, in the test cases and applications presented in the following sections, we use the second-order accurate shock-capturing total variation diminishing Lax-Friedrichs (TVDLF; Gábor & Dušan 1996) scheme and a Koren slope limiter (Koren & Vreugdenhil 1993). Since the effects of the radiative force, f_r, and its work, v ⋅ f_r, contribute on the same timescale as the advection with the gas velocity (i.e., the dynamical timescale), the corresponding source terms in the right-hand side of the momentum and gas energy Eqs. (2) and (3), respectively, are added explicitly; this is the case with the photon tiring term P : ∇v in the radiation energy Eq. (6).

On the other hand, the heating-cooling and radiation diffusion take place on timescales τ_q and τ_diff, which can be, depending on the regime, several orders of magnitude shorter than the dynamical timescale. Consequently, to ensure numerical stability, their corresponding terms, $\dot{q}$ and ∇ ⋅ F, in the gas and radiation energy Eqs. (3) and (6) are computed following an implicit procedure. How these aforementioned source terms and methods are combined is explained in Sect. 3.1. The rest of Sect. 3 is then devoted to elaborating the different source terms individually.

The FLD module discussed below is implemented in Cartesian coordinates for 1D, 2D, and 3D setups, thanks to the VAC-preprocessor. For illustrative purposes, we only write out stencils corresponding to a 2D setup. An extension to 3D or reduction to a 1D setup is a matter of adding or subtracting an index in the equations below.

For now, in the application of the WR star in Sect. 5 a spherical correction is used to correct for the geometry. This is further described in Appendix A. Finally, the new module described in this paper can, in principle, also be directly applied to radiation-MHD calculations.

3.1. Order of operations

Generally, the PDEs of the FLD system, Eqs. (1)–(3) and (6), can be written as:

$$\begin{array}{c}\hfill {\partial }_t\mathit{u}={\mathit{A}}_{\mathrm{adv}}(\mathit{u})+{\mathit{H}}_{\mathrm{diff}}(\mathit{u})+{\mathit{S}}_{\mathrm{ex}}(\mathit{u})+{\mathit{S}}_{\mathrm{im}}(\mathit{u}).\end{array}$$(21)

Here, u represents the vector of conservative variables [ρ, vρ, e, E]. A_adv(u) = − ∇ ⋅ [ρv, vρv + p, ev + pv, Ev] is the advection operator; H_diff(u) = [0, 0, 0, −∇⋅F] is the diffusion term; S_ex(u) = [0, f_r, v ⋅ f_r, −P : ∇v] is a collection of source terms which are to be handled explicitly: radiation force and work additions in momentum and total energy. Finally, ${\mathit{S}}_{\mathrm{im}}(\mathit{u})=[0,\mathbf{0},\dot{q},-\dot{q}]$ are the entirely local source terms which are handled implicitly: cooling and heating. In our implementation in MPI-AMRVAC, these four types of terms can be grouped in two operators: an implicit operator and an explicit operator:

$$\begin{array}{cc}& G_{\mathrm{ex}}(\mathit{u})={\mathit{A}}_{\mathrm{adv}}(\mathit{u})+{\mathit{S}}_{\mathrm{ex}}(\mathit{u})+{\mathit{S}}_{\mathrm{im}}(\mathit{u}),\hfill \end{array}$$(22)

$$\begin{array}{cc}& G_{\mathrm{im}}(\mathit{u})={\mathit{H}}_{\mathrm{diff}}(\mathit{u}).\hfill \end{array}$$(23)

We note that although G_ex is referred to as an explicit operator, it contains one local implicit step S_im(u) which is further discussed in Sect. 3.3. For the explicit operator G_ex(u) and for illustrative purposes, we only show the simplest formulation. The scheme employed to advance the vector of conservative variables from uⁿ to ${\mathit{u}}_{+}^{n+1}$ is:

$$\begin{array}{cc}& {\mathit{u}}_{\ast }^{n+1}={\mathit{u}}^n+\mathrm{\Delta }t{\mathit{S}}_{\mathrm{ex}}({\mathit{u}}^n),\hfill \end{array}$$(24)

$$\begin{array}{cc}& {\mathit{u}}_{\ast \ast }^{n+1}={\mathit{u}}_{\ast }^{n+1}+\mathrm{\Delta }t{\mathit{S}}_{\mathrm{im}}({\mathit{u}}_{\ast \ast }^{n+1},{\mathit{u}}_{\ast }^{n+1}),\hfill \end{array}$$(25)

$$\begin{array}{cc}& {\mathit{u}}_{+}^{n+1}={\mathit{u}}_{\ast \ast }^{n+1}+\mathrm{\Delta }t{\mathit{A}}_{\mathrm{adv}}({\mathit{u}}_{\ast \ast }^{n+1}),\hfill \end{array}$$(26)

where the states indexed with an asterisk or plus sign are intermediate states. As can be seen in Eqs. (24)–(26), the radiation force f_r, its work v ⋅ f_r and the photon tiring term P : ∇v are added as explained in Sect. 3.2 to retrieve a first intermediate state ${\mathit{u}}_{*}^{n+1}$. This intermediate state is used in the local implicit heating and cooling terms $\dot{q}$ as explained in Sect. 3.3 to retrieve another intermediate state ${\mathit{u}}_{\ast \ast }^{n+1}$. Finally, this updated state is used to calculate the advection fluxes of both the gas and radiation variables that are used in the Riemann solvers. This combination of terms described above can also be used when only advancing half a time step in, for example, a midpoint scheme, as discussed below.

The combination of G_ex and G_im takes place through an IMEX scheme. Multiple schemes of different orders are available; for illustrative purposes we show below only the second-order accurate Midpoint scheme. This is the scheme that was used for the test cases described in Sect. 4, unless stated otherwise. The midpoint scheme advances from uⁿ to u^n + 1 as follows:

$$\begin{array}{cc}& {\mathit{u}}_{+}^{n+1/2}={\mathit{u}}^n+\frac{1}{2}\mathrm{\Delta }t{G}_{\mathrm{ex}}({\mathit{u}}^n),\hfill \end{array}$$(27)

$$\begin{array}{cc}& {\mathit{u}}^{n+1/2}={\mathit{u}}_{+}^{n+1/2}+\frac{1}{2}\mathrm{\Delta }t{G}_{\mathrm{im}}({\mathit{u}}^{n+1/2}),\hfill \end{array}$$(28)

$$\begin{array}{cc}& {\mathit{u}}^{n+1}={\mathit{u}}^n+\mathrm{\Delta }t(G_{\mathrm{ex}}({\mathit{u}}^{n+1/2})+G_{\mathrm{im}}({\mathit{u}}^{n+1/2})),\hfill \end{array}$$(29)

where u^n + 1/2 is the state half a time step later. In the formulation above, (27) adds the G_ex term for half a time step, computed from uⁿ. Then, half of an implicit time step is added by advancing to u^n + 1/2 in (28). Finally, in (29), the full time step is computed by advancing uⁿ to u^n + 1 where G_ex computed from the midpoint state u^n + 1/2 and the value of G_im is re-used from (28). Here, in both the first and last operation (27) is performed following the scheme (24)–(26) for half a time step and a full time step, respectively.

In order to ensure numerical stability, an adapted Courant-Friedrichs-Lewy (CFL) condition (Courant et al. 1928) is implemented. For radiation-free hydrodynamics, a stable time step can be ensured by calculating the local maximum propagation speed of sound waves and limiting the time step in such a way that the distance traveled by a sound wave in one time step is smaller than the width of a cell. In an RHD setting, one has to take into account the propagation of radiative-acoustic waves when computing a constraint on the timestep. To this end, the radiation pressure is added to the gas pressure when calculating a time step constraint.

3.2. Radiation force and photon tiring

We next outline numerical aspects behind computing and adding the explicitly handled S_ex terms, namely, the radiation force (9), its work term, and the photon tiring term. Both the flux limiter λ and the radiation flux according to (18) depend on the gradient of E. This gradient is computed using a fourth order central difference with a five-point stencil (Fornberg 1988). For instance, in a 2D set-up, the x-component of the gradient is expressed as:

$$\begin{array}{c}\hfill {\left(\mathrm{\nabla }{E}_{i,j}\right)}_x=\frac{\frac{1}{12}{E}_{i-2,j}-\frac{2}{3}{E}_{i-1,j}+\frac{2}{3}{E}_{i+1,j}-\frac{1}{12}{E}_{i+2,j}}{\mathrm{\Delta }x}·\end{array}$$(30)

We use Eqs. (9) and (18) to add, in an explicit way, the source terms for the momentum and gas energy. In this way, the radiative force and its work are not purely local since they are derived from the radiative flux using Eq. (30) above.

The momentum and gas energy components of S_ex describe the radiation force and its work exerted on the gas. To add them explicitly for a time step of Δt, the following recipe is used:

$$\begin{array}{cc}& {\left(\rho \mathit{v}\right)}^{n+1}={\left(\rho \mathit{v}\right)}^n+\mathrm{\Delta }t{\rho }^n\frac{{\kappa }^n{\mathit{F}}^n}{c},\hfill \end{array}$$(31)

$$\begin{array}{cc}& e^{n+1}=e^n+\mathrm{\Delta }t\frac{{\kappa }^n{\rho }^n{\mathit{v}}^n·{\mathit{F}}^n}{c}·\hfill \end{array}$$(32)

Moreover, P  :  ∇v, which enters Eq. (6), is calculated by the dyadic product between the radiation energy pressure tensor and the gradient of the velocity vector. Using (15) and (19), the components of the radiation pressure tensor can be calculated. Here, again, the gradients of each of the components of the velocity field are calculated with a similar five-step stencil for improved accuracy. As an example, below we provide the representation of the photon tiring work term in the point (i, j) for a 2D setup:

$$\begin{array}{cc}\hfill {(\mathit{P}:\mathrm{\nabla }\mathit{v})}_{i,j}=& P_{i,j}^{x,x}{\left(\mathrm{\nabla }{(v_x)}_{i,j}\right)}_x+P_{i,j}^{x,y}{\left(\mathrm{\nabla }{(v_x)}_{i,j}\right)}_y\hfill \\ \hfill & +P_{i,j}^{y,y}{\left(\mathrm{\nabla }{(v_y)}_{i,j}\right)}_y+P_{i,j}^{y,x}{\left(\mathrm{\nabla }{(v_y)}_{i,j}\right)}_x.\hfill \end{array}$$(33)

Also, this source term is added explicitly in every cell. For the addition of Δt S_ex in the radiation energy density:

$$\begin{array}{c}\hfill E^{n+1}=E^n+\mathrm{\Delta }t({\mathit{P}}^n:\mathrm{\nabla }{\mathit{v}}^n).\end{array}$$(34)

3.3. Heating and cooling

As mentioned above, our non-equilibrium FLD description allows the gas temperature to be different from the radiation temperature. Locally, the gas temperature can be evaluated from the gas pressure and density using the ideal gas law (4). The radiation temperature is evaluated from $E=a_{\text{r}}{T}_{\text{r}}^4$. When both temperatures are equal, the system is in radiative equilibrium. However, when the radiation temperature is higher (lower) than the gas temperature, the gas will heat up (cool down) due to an energy exchange with the radiation field. This energy exchange is written as $\dot{q}$ in Eqs. (3) and (6). The timescale for heating and cooling is typically much shorter than a Courant time step. To ensure numerical stability over a time step, the gas heating and cooling terms are therefore added implicitly. This approach is based on the method described in Turner & Stone (2001). Since the gas cooling term depends only on the internal gas energy density ϵ (and is independent of kinetic energy), we first compute the internal gas energy density ϵ = e − ρv²/2. Using (4), a point-implicit, discretized formulation of the heating and cooling terms applied to the gas internal and radiation energy density are written as:

$$\begin{array}{cc}& {ϵ}^{n+1}={ϵ}^n+\mathrm{\Delta }t[c{\kappa }^n{\rho }^n{E}^{n+1}-4{\kappa }^n{\rho }^n\sigma {\left(\frac{m_{\mathrm{p}}\mu }{k_{\mathrm{B}}}(\gamma -1)\frac{{ϵ}^{n+1}}{{\rho }^n}\right)}^4],\hfill \end{array}$$(35)

$$\begin{array}{cc}& E^{n+1}=E^n+\mathrm{\Delta }t[-c{\kappa }^n{\rho }^n{E}^{n+1}+4{\kappa }^n{\rho }^n\sigma {\left(\frac{m_{\mathrm{p}}\mu }{k_{\mathrm{B}}}(\gamma -1)\frac{{ϵ}^{n+1}}{{\rho }^n}\right)}^4]·\hfill \end{array}$$(36)

The adiabatic index γ and mean molecular weight μ are constant for a given simulation as we do not take into account ionization effects. Solving these coupled equations comes down to first finding the root of the following fourth-degree polynomial (Turner & Stone 2001):

$$\begin{array}{c}\hfill {\left({ϵ}^{n+1}\right)}^4+\left(\frac{1+a_2}{a_1}\right){ϵ}^{n+1}-\left(\frac{1+a_2}{a_1}\right){ϵ}^n-\frac{a_2}{a_1}{E}^n=0,\end{array}$$(37)

where a₁ = 4κσ(γ − 1)⁴/ρ³Δt and a₂ = cκρΔt, with both κ and ρ evaluated at step n. To get the gas internal ϵ^n + 1 energy at the next timestep, Halley’s root-finding method (Press et al. 2007) is employed. This method is similar to Newton-Raphson, but it uses the second derivative for faster convergence. If this method does not reach a user set tolerance after a certain number of iterations (typically 100), we switch to a bisection scheme. Equation (36) only has one real positive root which has to be smaller than the sum of the internal and radiation energies and thus lies in the open interval (0, ϵⁿ + Eⁿ). With the updated internal gas energy, the radiation energy density is given by:

$$\begin{array}{c}\hfill E^{n+1}=\frac{a_1{\left({ϵ}^{n+1}\right)}^4+E^n}{1+a_2}·\end{array}$$(38)

In many astrophysical contexts, gases are radiation-dominated (r = E/4γϵ ≫ 1), meaning that the energy budget is controlled mainly by the radiation field. In such cases, it is crucial to have an accurate solution to Eq. (37), as small errors in a relatively small gas energy density as compared to the radiation energy could lead to a significant error in the gas temperature calculated by the ideal gas law. Hence, it is important to choose a sufficiently low tolerance. The RHD codes have previously been shown to sometimes be vulnerable to this problem in extreme conditions when the time step used to advance the system is several orders of magnitude larger than the time it takes the gas to reach thermal equilibrium (Jiang et al. 2012). In Sect. 4.1, the heating and cooling algorithm is tested for a highly radiation-dominated gas.

3.4. Radiation diffusion

For the diffusion term in Eq. (6), we make use of the newly developed geometric multigrid method library octree-mg (Teunissen & Keppens 2019). This MPI-parallelized library is capable of solving elliptic PDEs on 1D,2D or 3D Cartesian grids and is fully compatible with the block-tree AMR structure used in MPI-AMRVAC. The contribution by the diffusion term can be written in an operator split way by only considering the advection of the CMF radiation flux from Eq. (6), ∂_tE + ∇ ⋅ F = 0, with F substituted using (14):

$$\begin{array}{c}\hfill {\partial }_{\mathrm{t}}E+\mathrm{\nabla }·(-D\mathrm{\nabla }E)=0.\end{array}$$(39)

Due to the elliptic nature of the diffusion term and the difference between the gas and radiation dynamical timescales, it is computationally inefficient to use an explicit solver. Indeed, in many astrophysical contexts, such explicit solvers for radiative diffusion often require prohibitively short timesteps to be stable (typically several orders of magnitude smaller than that for the advection term). As such, an implicit method is used here instead. Geometric multigrid methods speed up the convergence rate of iterative relaxation by relaxing the PDE solution on a hierarchy of grids. In this way, large scale errors can be damped more efficiently on a coarser grid, while small-scale errors are smoothed on a finer grid. The multigrid library uses a standard Gauss-Seidel smoother over multiple grid levels to solve a Helmholtz equation.

The diffusion problem can be recast as a Helmholtz equation with a variable diffusion coefficient, where Eⁿ is the current radiation energy density and E^n + 1 is the radiation energy density in the next time step:

$$\begin{array}{c}\hfill E^{n+1}/\mathrm{\Delta }t-\mathrm{\nabla }\left(D^n\mathrm{\nabla }{E}^{n+1}\right)=E^n/\mathrm{\Delta }t.\end{array}$$(40)

As mentioned above, our methods are implemented in 1D, 2D, and 3D. In 2D, this equation is discretized using the following five-point stencil:

$$\begin{array}{cc}& [D_{i-1/2,j}^n(E_{i,j}^{n+1}-E_{i-1,j}^{n+1})-D_{i+1/2,j}^n(E_{i+1,j}^{n+1}-E_{i,j}^{n+1})]/\mathrm{\Delta }{x}^2\hfill \\ \hfill & +[D_{i,j-1/2}^n(E_{i,j}^{n+1}-E_{i,j-1}^{n+1})-D_{i,j+1/2}^n(E_{i,j+1}^{n+1}-E_{i,j}^{n+1})]/\mathrm{\Delta }{y}^2\hfill \\ \hfill & +E_{i,j}^{n+1}/\mathrm{\Delta }t=E_{i,j}^n/\mathrm{\Delta }t.\hfill \end{array}$$(41)

In addition to an initial state Eⁿ, in order to advance the radiation energy density, it is necessary to specify boundary conditions at every boundary. This can be done by either defining the value E^n + 1 in the ghost cells (Dirichlet conditions), by assuming a fixed gradient of E^n + 1 at the interface between the ghost cells and the numerical domain (Neumann conditions), or by a linear extrapolation into the ghost cells (continuous boundary conditions).

In (41), the diffusion coefficient is evaluated between two neighbouring cells using a harmonic mean as described by Teunissen & Keppens (2019):

$$\begin{array}{c}\hfill D_{i-1/2,j}^n=\frac{2D_{i,j}^n{D}_{i-1,j}^n}{D_{i,j}^n+D_{i-1,j}^n}·\end{array}$$(42)

Locally, the diffusion constant in the cell centers is calculated from Eq. (18), where the flux-limiter λ is computed according to (17). First, the diffusion solver will try a full multigrid (FMG) cycle to converge to a requested residual which for the tests described below was typically chosen to be 10⁻⁵. In an FMG cycle, the smoothing starts at the coarsest grid of the MG-solver, and then prolongs, smoothens and restricts to ever increasing finer grids. If the FMG-cycle does not reach this tolerance, additional V-cycles are added until the convergence criterion is met. A V-cycle starts at the finest grid, goes level by level to the coarsest grid and then back again to the finest grid. This process is elaborated in Teunissen & Keppens (2019).

Though fully compatible with the adaptive mesh refinement, this FLD module is currently limited to Cartesian settings due to the abilities of the multigrid solver used to solve the diffusion part of Eq. (6). In the future, we plan to expand this module to the cylindrical, polar and spherical meshes which are already available for the usual non-radiative HD/MHD setup. These limitations in the multigrid solver arise from using point-wise relaxation methods in smoothing the error in the solution. These relaxation methods are not efficient when used on spherical diffusion operators, as discussed in Teunissen & Keppens (2019) and Briggs et al. (2000).

4. Results: Test cases

In this section, we describe how we tested aspects of the newly developed FLD code in various regimes. Since non-trivial analytic solutions to the FLD equations are scarce, we set up a number of test-cases and then compare them to various semi-analytic solutions or predictions (see also, for example, Turner & Stone 2001; Krumholz et al. 2006). In Sect. 4.1, the testing of the implicit heating and cooling terms are described, for which the algorithm has been explained in Sect. 3.3. Both tests presented in Sects. 4.2 and 4.5 test the full system of equations, in a steady state and dynamic state, respectively. In Sect. 4.3, we focus on a test covering the Galilean invariance of the code. Finally, the test described in Sect. 4.4 examines the optically thin limit.

4.1. Heating and cooling

As a first test, we check the energy exchange by means of heating and cooling in a gas with zero velocity and a constant density. The 2D domain has periodic boundary conditions on all sides. The gas energy and radiation energy are initialized out of radiative equilibrium, so the gas temperature is not equal to the radiation temperature. Due to this non-equilibrium, there is a net heating-cooling term $\dot{q}$, which will relax the gas energy density to its equilibrium value and raise (or lower) the gas and radiation temperatures until they are equal. Since there is no velocity, there is no kinetic energy and the internal gas energy is equal to the total gas energy, ϵ = e.

We consider a gas with density ρ = 10⁻⁷ g cm⁻³, opacity κ = 0.4 cm² g⁻¹, adiabatic index γ = 5/3, and mean molecular weight μ = 0.6. The radiation energy density is initialized as E = 10¹² erg cm⁻³.

The equilibrium energy density e_equil for gas and radiation can be calculated from the conservation of the sum of the initial energies and the condition that there is no net cooling or heating when the system is in equilibrium. We consider two different uniform initial conditions for the gas energy density relative to this equilibrium energy density. In the first test, e is set to 10⁻⁶ e_equil (‘x’ symbols in Fig. 1) and in the second test, it is set to 10² e_equil (‘o’ symbols). Thereafter, we let the system relax to equilibrium and monitor the gas and radiation energy through time. When the gas reaches equilibrium, it is radiation-dominated, with r = 2180.

Fig. 1. Gas energy densities, as a function of time, of a gas that is relaxing to radiative equilibrium from two sets of initial conditions. Marked with a blue “o” is the relaxation where the initial gas energy is set as 102 times the expected equilibrium value. Marked with a red “x” is the initial gas energy is set as 10−6 times the equilibrium value. The dashed blue and red lines represent the same heating-cooling problem, but solved with a smaller numerical time step (see text). Since the time axis is logarithmic, the initial conditions at t = 0 are not visible.

In Fig. 1, the gas energy density is plotted through time on a log-log plot, together with the energy equilibrium value. For this test, it is important that the rate at which the temperature approaches its equilibrium is correct. To this end, we compare the previous simulation for two different fixed time steps: once where Δt = 10⁻¹² s, which is regarded as a physically correct baseline; and once where Δt = 10⁻¹¹ s. From the resulting curves, we can conclude that in the test where the gas energy density started out lower than the equilibrium energy and that it relaxes to the theoretically predicted value e_equil on a timescale which is independent of the used numerical time step. However, for the case where the gas energy is initiated greater than the equilibrium value, the heating lags behind for a shorter time step in the first couple of iterations. Later on, the correct equilibrium value is found on a correct timescale. We note that the time step of 10⁻¹¹ s used here is greater than the initial cooling timescale, τ_cool ∼ e/(4πκρB)∼10⁻¹² s.

Although the radiative energy density is left free to evolve, it remains essentially constant because its initialized value is much larger than the gas energy density. However, we verified that it matches with the equilibrium value for radiative energy. For both cases, at the end of the test, the relative difference between the gas and radiation temperatures is below 0.1%, and this value is steadily decreasing as the simulation time progresses. Although not plotted in the figure for reasons of clarity, similar agreements were obtained for relaxation from the same initial conditions, but when using a time step that is several orders of magnitude larger than the time it takes for the system to reach thermal equilibrium (Δt = 10⁻⁵ s). In this case, the relative error in temperature converged to the requested tolerance after three time steps.

4.2. Radiation-dominated shock

This test case here describes a steady-state radiation-dominated shock. Although there is no analytical solution for the exact shape of the shock front, there are approximations for its expected width (Mihalas & Mihalas 1984). Using radiation-modified Rankine-Hugoniot conditions, a left- and right-hand side gas state are calculated and used as the initial condition before they are relaxed to a steady state. With this test, we can assess how the code handles conserved quantities. Moreover, the discontinuity is also an ideal situation to test the adaptive mesh refinement (AMR).

On the left-hand side of the shock, gas density ρ_l = 10⁻² g cm⁻³, velocity v_l = 10⁹ cm s⁻¹ and total (kinetic + internal) gas energy density e_l = 5.0 × 10¹⁵ erg cm⁻³. The radiation energy density is set in equilibrium to E_l = 75.6 erg cm⁻³. On the right-hand side ρ_r = 6.85874 × 10⁻² g cm⁻³, v_r = 1.45 × 10⁸ cm s⁻¹, e_r = 1.93 × 10¹⁵ erg cm⁻³, and E_r = 2.44 × 10¹⁶ erg cm⁻³. The stability of the shock depends on having the initial values conform with the radiation-modified Rankine Hugoniot conditions. Empirically, at least five decimals for the right-hand side density were needed to converge to these profiles. We assume we are in the optically thick diffusion limit so that λ = 1/3. Both the radiation flux and radiation force will point upstream, widening the shock to a width which has been predicted to be w = D/v_l (Mihalas & Mihalas 1984; Turner & Stone 2001). We consider a fully ionised pure hydrogen gas with μ = 0.5 and γ = 7/5. The simulation box is 10⁵ cm wide and consists of 256 cells on the lowest AMR level (level 1). In Fig. 2, we show the relaxed shock after ten flow passing times.

Fig. 2. Density, velocity, gas pressure, and radiation energy profiles for a 1D simulation of a radiation-dominated shock after ten flow passing times. The profiles are normalized to the initial right-hand side values, except for the velocity profile, which has been normalized to the initial left-hand side value. The expected width of the shock is illustrated in grey.

As can be seen in Fig. 2, the shock width of the 1D simulation is conform with the predictions made by Mihalas & Mihalas (1984). We obtained the same results for both 2D and 3D simulations of the shock along the x-axis (performed, but not shown here). Moreover, we can check whether the solution is truly steady by comparing the momentum on both sides of the shock. From this, we retrieve a 0.2% relative error for a low resolution run on N_x = 64 without any AMR, a 0.1% relative error for a high-resolution run on N_x = 256 without any AMR and similarly a 0.1% relative error for a run with base resolution of N_x = 64 but an effective resolution of N_x = 256 by using three levels of refinement.

4.3. Galilean invariance

In this test, we check the Galilean invariance of the code during an advection diffusion problem. Again, there is no available analytic solution, but if the Galilean invariance is respected, the two profiles should keep the same shape. The same simulation was performed twice: once with and once without an initially constant background velocity field. If solved for correctly, the profiles of the conserved quantities in the two simulations will be the same but translated. Following Krumholz et al. (2006), we consider a slab of gas in total pressure equilibrium. In the center of the gas, there is a dip in density and gas pressure along with a corresponding bump in radiation energy density and pressure to preserve a constant total pressure. Once the simulation is begun, the radiation is diffused out of the dip, the total pressure equilibrium is lost, and the gas will be driven towards the center, where the gas pressure is lower. The initial conditions are given by:

$$\begin{array}{cc}& T=T_0+(T_1-T_0)exp\left(\frac{-x^2}{2w^2}\right),\hfill \end{array}$$(43)

$$\begin{array}{cc}& \rho ={\rho }_0\frac{T_0}{T}+\frac{a_{\mathrm{r}}\mu }{3k_{\mathrm{B}}}(\frac{T_0^4}{T}-T^3),\hfill \end{array}$$(44)

where T₀ = 10⁷ K, T₁ = 2 × 10⁷ K, ρ₀ = 1.2 g cm⁻³ and w = 24 cm. Both E and e are set in equilibrium to the above temperature profile. The mean molecular weight is taken to be μ = 2.33, which corresponds to a Helium abundance of 0.1. The adiabatic exponent γ = 5/3 and κ = 100 cm² g⁻¹. This simulation is run twice: a first time with zero background velocity and a second time with a constant background velocity of v = 5 × 10⁷ cm s⁻¹. In Fig. 3, we show both solutions after approximately 1 pulse width crossing times.

Fig. 3. Comparison of still and advected solutions in the Galilean invariance test. Top panel: density profiles with initial conditions (in dashed lines) are relaxed once with (marked with ‘o’) and once without (marked with ‘x’) an initial background velocity field. Bottom panel: same details but for the gas pressure (blue) and radiation energy density (red). The horizontal shift between the curves correspond to the expected rightward displacement.

As can be seen from Fig. 4, the relative difference between the advected and stagnant density profile peaks at 0.03% when using a second order Midpoint IMEX-scheme. However, when using a simpler first order Euler scheme, the relative difference in density goes up to 1.7%. Similar improvements can be seen for the relative differences in gas pressure and radiation energy density. This clearly shows the advantage of using higher order time-stepping schemes.

Fig. 4. Absolute value of the relative difference between the profiles with and the profiles without an initial background velocity field from the Galilean invariance test in Fig. 3, where the advected profile is translated back to the origin. We compare two different time integration schemes: a first-order Euler IMEX scheme (marked ‘+’) and a second-order Midpoint IMEX scheme (marked ‘⋅’). Different colors correspond to different quantities: gas density in black, gas pressure in blue, and radiation energy in red.

4.4. Optically thin limit

While all of the previous tests were situated in an optically thick regime, the following test assesses the workings of the FLD code in an optically thin situation. This is an important test, since if not flux-limited, the diffusion approximation breaks down in this regime.

Due to the flux limiter λ we should recover the free streaming flux |F|=cE when the optical depth approaches zero (ensuring the radiation energy density does not travel faster than the speed of light). In this test case, we start with a slab of gas with ρ₀ = 0.025 g cm⁻³ and constant opacity κ = 0.4 cm² g⁻¹. The numerical domain reaches from x = −0.5 cm to x = 1.5 cm and it is subdivided in N_x = 256 grid cells. This gives a total optical thickness of the slab of τ = 0.002. The initial radiation energy density is set according to an error function, centered around x = 0 with a width d = 0.05 cm:

$$\begin{array}{c}\hfill E(x)=E_0+\frac{1}{2}[1-\mathrm{erf}\left(\frac{x}{d}\right)]E_1.\end{array}$$(45)

This expression lights up the slab from the left-hand side, with E₁ = 1.4 × 10¹¹ erg cm⁻³. There is no velocity and the gas energy everywhere is set to be in equilibrium with E₀ = 1.4 × 10⁻¹¹ erg cm⁻³. For this setup, we only performed the radiation diffusion (thus the radiation force, heating, cooling, and photon tiring terms are not considered). Since the local sound speed is not relevant for the free streaming radiation field, the time step is set manually to Δt = 10⁻¹³ s, which is on the order of half a cell-crossing time for the speed of light. This time step is very stringent but is chosen only for the sake of illustrative purposes. In practice, larger time steps can be used.

The propagation speed of the radiation energy density is, as expected, limited by the speed of light as can be seen in Fig. 5. The initial profile has diffused a little, which is to be expected when executing what is practically an advection operator with a diffusion solver.

Fig. 5. Space-time diagram of the radiation energy front evolving in an optically thin slab of gas. In the optically thin limit, the FLD flux limiter λ prohibits radiation from moving faster than the speed of light. In this plot, the speed of light has been indicated by the black dotted lines. In addition to the colormap, the position of the radiation front has been evaluated at three arbitrarily values with the red contour lines. The total optical thickness of the slab is 0.002.

4.5. Linear RHD wave

When interacting with radiation, acoustic waves can be naturally damped by radiating the energy away. This process can be modeled with the FLD description of RHD. With the following RHD wave setup, the full system of equations can be tested against results of a semi-analytic perturbation relation. In addition to checking the different algorithms for adding the different source terms described in Sect. 3, this setup allows us to assess the numerical diffusion of the code. We also compare the 1D problem with a 2D setup which is symmetric along the first diagonal of the grid, which allows us to rule out any effects of grid anisotropy on the simulation outcome.

We test the code with dispersion relations that come out of an analytical linear perturbation analysis of the RHD system in the diffusion limit done by Mihalas & Mihalas (1984). This analysis shows that the damping length of a running wave is dependent on two dimensionless parameters. The Boltzman number Bo ∼ 4γc_ae/(cE), with $c_{\mathrm{a}}=\sqrt{\gamma p/\rho }$ the adiabatic sound speed, translates to how much radiation contributes to heat transfer. The ratio of energy densities r ∼ E/(4γϵ) is high when a gas is radiation energy-dominated and low when the gas is gas energy-dominated.

A slab of gas is considered with background values ρ₀ = 3.216 × 10⁻⁹ g cm⁻³, e₀ = 26.02 × 10³ erg cm⁻³, and E₀ = 17.34 × 10³ erg cm⁻³. The adiabatic index is γ = 5/3 leading to a Boltzman number of Bo = 10⁻³ and an energy density ratio of r = 0.1. The slab is perturbed on the left-hand side with the following boundary conditions to excite a traveling wave:

$$\begin{array}{cc}& \rho ={\rho }_0+A_{\rho }sin(kx-\omega t),\hfill \end{array}$$(46)

$$\begin{array}{cc}& \mathit{v}=A_{v}sin(kx-\omega t),\hfill \end{array}$$(47)

$$\begin{array}{cc}& e=e_0+A_{\mathrm{e}}sin(kx-\omega t).\hfill \end{array}$$(48)

Here, the wavenumber k is chosen in such a way that for the constant opacity of κ = 0.4 cm² g⁻¹, the optical depth across one wavelength is τ_λ = 10³. The pulsation ω is related to k via the propagation speed of the RHD wave. The perturbation with density amplitude A_ρ = 0.01 ρ₀ will travel along through the gas, but it will be damped by the effects of the radiation field. Amplitudes for the velocity and gas energy perturbation are set by A_v = A_ρω/(kρ₀) and A_e = A_ρe₀/ρ₀. Gas energy in the oscillation will be transferred to the radiation field by means of the cooling mechanism. After this, the radiation energy will leak outward and leave the system due to diffusion. The dispersion relation given by Mihalas & Mihalas (1984) can be solved with a standard root-finding algorithm to provide a theoretical dampening length for the oscillation. For the numbers provided here, the predicted dampening length is 8.16λ, with λ the wavelength of the induced acoustic oscillation (not to be confused with the flux-limiter introduced in Sect. 2.2).

In Fig. 6, results are shown for a perturbation with the optical depth τ_λ = 10³ in a 1D setup. For the 2D version, the same problem is set up along the line of the first diagonal. Now, the wave is driven at any point where x + y < λ. Instead of (kx − ωt), the sine in the driving conditions now take $(k(x+\mathit{y})/\sqrt{2}-\omega t)$ as an argument. Additionally, boundary conditions are copied from their diagonal neighbours for all conserved quantities, such that they correspond to the correct phase in the diagonally traversing wave. For the left and top boundary, the conserved quantities in cell (i, j) are copied from their bottom right neighbour: u_i, j = u_{i + 1, j − 1}. For the right and bottom boundary, they are copied from the top left neighbour: u_i, j = u_{i − 1, j + 1}.

Fig. 6. Comparison between numerical and analytical solutions of the RHD wave. From top to bottom: the density perturbation, velocity, gas energy perturbation, and radiation energy perturbation after 40 oscillation periods, for an RHD wave where the optical thickness of a wavelength is equal to 1000. For the top three panels, this is compared to the analytic profile (black dotted lines) that is calculated as described by Mihalas & Mihalas (1984). The wave is driven from the grey shaded area on the left side of the plot.

The solution depends heavily on the flux limiter that is used in the approximate Riemann solver. More diffusive schemes such as minmod (Roe 1986) or Koren (Koren & Vreugdenhil 1993) give a shorter dampening length, while more advanced, higher order schemes, such as a fifth order weighted non-oscillating scheme (Liu 1994) or those that are monotonicity preserving (Suresh & Huynh 1997), are better at approaching the correct signal. Overall, Figs. 6 and 7 show that the numerical solution is a very good match to the first order perturbation relation, both in wavelength and dampening length. As seen from Fig. 8, the solution stays symmetric along the first diagonal, as expected.

Fig. 7. Comparison between 1D and 2D numerical solutions of the RHD wave. Top panel: density perturbation profile of an RHD wave after 40 oscillation periods. In black, the analytically predicted profile, in blue results from a 1D simulation and in red results from a 2D simulation with the wave vector parallel to the first diagonal. The x-axis is here chosen in the direction of the wave vector. Bottom panel: respective residuals calculated as |ρ − ρanalytic|/ρanalytic, where the analytic density is calculated from perturbation theory.

Fig. 8. Map of density perturbation in an RHD wave along the first diagonal (black dotted line) of a 2D plane, after 40 oscillation periods. The wave is driven from the grey area in the bottom-left corner.

5. First research application: Wolf Rayet wind

As a final example of our code’s applicability, we performed a 1D simulation of a supersonic, optically thick Wolf-Rayet wind outflow. Classical Wolf-Rayet stars are massive stars that have evolved back to the blue side of the Hertzsprung-Russel diagram, after shedding their outer hydrogen layers thus exposing a helium core (see review by Crowther 2007). Observationally they are known to have supersonic wind outflows characterized by high mass-loss rates. Wolf-Rayet winds are believed to be accelerated by radiation (Sander & Vink 2020; Poniatowski et al. 2021), but due to the high mass-loss rate their hydrostatic surface lies deep within optically thick layers. This provides an interesting first research application for the FLD code, testing the effects of a strong radiation force in an optically thick and highly supersonic environment.

The simulation set-up is based on the recent work by Poniatowski et al. (2021), but now using the FLD method described in previous sections to compute the radiation force and energy balance. This then allows for a more flexible (and complete) approach for time-dependent modeling of such Wolf-Rayet outflows. Specifically, while the simulations by Poniatowski et al. (2021) assumed that the local radial flux is always set by L_*/(4πr²), where L_* is the stellar core luminosity, in the FLD method presented here fluxes (and radiation work terms, neglected in Poniatowski et al. 2021) are computed directly from the evolving energy density. As such, in contrast to the Poniatowski et al. (2021) model, the FLD method presented here could be readily extended to time-dependent 2D or 3D flows with local (and potentially non-radial) fluxes and forces, as further discussed below. Also, a spherically symmetric wind is assumed, hence, the Cartesian formulation of the RHD equations has to be corrected for spherical geometry. For this, we follow the recipe by Sundqvist et al. (2018), as further explained in the appendix.

As in Poniatowski et al. (2021), for this first 1D Wolf-Rayet simulation, we assume a fixed stellar mass M_* = 10 M_⊙, a gravity source term f_g = ρGM_*/r², a hydrostatic core lower boundary radius R_* = 1 R_⊙, and a stellar core luminosity of log₁₀(L_*/L_⊙) = 5.416.

5.1. Opacities

Since the WR outflow is initiated by the radiation force, a key feature in this model regards the applied opacities. To this end, we follow Poniatowski et al. (2021) and use a superposition of equilibrium opacities computed in the static limit and a simple parametrized form for the large enhancement of line-opacity expected in a supersonic flow:

$$\begin{array}{c}\hfill \kappa ={\kappa }^{\mathrm{OPAL}}(\rho ,T_{\mathrm{g}})+{\kappa }^{\mathrm{CAK}}(\rho ,\mathrm{d}v/\mathrm{d}r).\end{array}$$(49)

Here κ^OPAL, which represents the opacities computed for static media, is taken from the tabulations by Iglesias & Rogers (1996) and κ^CAK uses a variant of the parametrization first introduced by Castor et al. (1975) (‘CAK’) to represent the accumulative effect of Doppler-shifted lines. At every time step, the CAK-opacity is computed locally using a second order central difference derivative of the radial velocity with respect to radius (ud-Doula & Owocki 2002):

$$\begin{array}{c}\hfill {\kappa }_i^{\mathrm{CAK}}=\frac{{\kappa }_{\mathrm{e}}\overline{Q}}{1-\alpha (r)}{\left(\frac{1}{c{\kappa }_{\mathrm{e}}\overline{Q}{\rho }_i}\left|\frac{v_{i-1}-v_{i+1}}{r_{i-1}-r_{i+1}}\right|\right)}^{\alpha (r)}·\end{array}$$(50)

This combined opacity is then used in all source terms (radiative force, heating or cooling, photon-tiring) as well as in the computation of the diffusion coefficient.

In Eq. (50), ${\kappa }_{\mathrm{e}}\overline{Q}$ represents the line opacity in the limit that all lines would be optically thin, and α represents the slope of the underlying assumed power-law distribution of lines. In general, these parameters should be derived from excitation and ionization calculations using full line lists (Puls et al. 2000; Lattimer & Cranmer 2021). In this first application, however, we assume the same set of parameters as in Poniatowski et al. (2021), which means we also here have introduced a radial variation of α. In the inner wind, for 0 < 1 − R_*/r < 0.4, the exponent is set constant at α = 0.66. In the outer wind, for 0.65 < 1 − R_*/r < 1, α = 0.5, to ensure a steady outflow. In the transition region, where 0.4 < 1 − R_*/r < 0.65, α decreases linearly as a function of 1 − R_*/r from 0.66 to 0.5. As discussed in that paper, this assumed variation makes the radial outflow stable against fallback by ensuring an outer-wind radiation force strong enough to accelerate the gas towards infinity.

5.2. Initial and boundary conditions

The initial conditions for the density and radial velocity in this simulation are derived from a constant mass loss rate $\dot{M}=4\pi r^2\rho {\mathit{v}}_{\mathrm{r}}$ and a so-called β-velocity law v(r) (see below). From L_* = 4πr²F_r and the Eddington approximation $F_{\mathrm{r}}=\frac{c}{3{\kappa }_0\rho }\mathrm{\nabla }E$ with a constant opacity (set to the electron scattering value for a fully ionised helium plasma), we obtain initial conditions for E by integrating the radiation energy density inward from the outermost point of the simulation where we assume a known floor radiation temperature T(r_max) = 5 × 10⁴ K. Finally, the ideal gas law is used to compute the gas pressure everywhere by assuming thermal equilibrium with the radiation. In putting this together, we obtain the following for our initial conditions:

$$\begin{array}{cc}& \rho (r)=\frac{\dot{M}}{4\pi r^{2}v},\hfill \end{array}$$(51)

$$\begin{array}{cc}& v(r)=v_{\infty }{(1-\frac{R_{\ast }}{r})}^{\beta },\hfill \end{array}$$(52)

$$\begin{array}{cc}& E(r)=a_{\mathrm{r}}{(T(r_{\max }))}^4+{\displaystyle {\int }_{r_{\max }}^r\frac{3{\kappa }_0\rho L_{\ast }}{4\pi r^{2}c}\mathrm{d}r,}\hfill \end{array}$$(53)

$$\begin{array}{cc}& p(r)=\frac{k_{\mathrm{B}}\rho (r)}{m_{\mathrm{p}}\mu }{\left(\frac{E(r)}{a_{\mathrm{r}}}\right)}^{\frac{1}{4}},\hfill \end{array}$$(54)

where we set v_∞ = 1000 km s⁻¹ and β = 1/2 for the initial velocity field. We note, however, that it is clear that the relaxed steady state will not resemble this simple β-law model; the conditions above just provide a good setup for initiating the actual simulation.

At the lower boundary, which is subsonic and so bound to the star, we follow the basic setup of previous radiation-driven wind simulations and fix the gas density while letting the velocity float (Sundqvist & Owocki 2013; Driessen et al. 2019; Poniatowski et al. 2021). However, here we additionally need to set the boundary condition for the radiation energy density, E. Instead, we use the fixed stellar luminosity to obtain the gradient of the radiation energy, using the diffusion coefficient calculated at the previous time step. A simple finite difference then gives for the lower boundary energy density:

$$\begin{array}{c}\hfill E_{i-1}^n=E_i^n+\frac{L_{\ast }}{4\pi r_i^2{D}_i^{n-1}}\mathrm{\Delta }r.\end{array}$$(55)

Finally, the gas energy in the ghost cells is set to be in equilibrium with the radiation field.

At the supersonic outer boundary, the density, momentum, and gas energy are extrapolated. The radiation energy density is set here by first computing the local optical depth and from this the radiation temperature, with the optical depth obtained by analytic inward integration from r = ∞, assuming the wind has reached its asymptotic velocity at the outer boundary and again that the opacity outside this is a constant that is set by electron scattering.

The simulation is run on a 1D Cartesian grid which stretches from 1R_* to 11R_*. The FLD solver module is only constructed for such Cartesian geometry. A 1D stellar outflow, however, is a spherically symmetric problem. For this reason, a correction term for the spherical divergence is added to all conservation equations following Sundqvist et al. (2018), as explained further in the appendix. To increase the resolution near the base (important to resolve the subsonic region), a constant refinement is used in the first 2 stellar radii. On the finest level, the 10R_* are resolved by 2048 grid cells. For this simulation we use a 3-step scheme with a TVDLF solver and a minmod slope limiter. For this numerical set-up, the full simulation discussed below takes about 20 min on a single laptop core to evolve over 20 flow passing times.

5.3. Relaxed profiles

The above-described simulations are run until they reach a relaxed steady state. Figure 9 shows a resulting relaxed velocity profile calculated with our new FLD module, comparing this with the velocity profile calculated by Poniatowski et al. (2021). The main feature, which is the stagnated and non-monotonic velocity profile, is similar between the two methods. Both profiles reach the same terminal velocity, however, the bump in the FLD profile has a slightly lower velocity than that from Poniatowski et al. (2021). The resulting stable mass-loss rate for the FLD model is 1.87 × 10⁻⁵ M_⊙ yr⁻¹, which again is very similar to the 1.64 × 10⁻⁵ M_⊙ yr⁻¹ found by Poniatowski et al. (2021).

Fig. 9. Radial velocity profile of a spherically symmetric Wolf-Rayet outflow in a relaxed state. In red, the profile is computed with the 1D FLD method while in black, it is computed following the alternative method by Poniatowski et al. (2021).

Figure 10 shows a similar comparison, but for the radiation temperature structure, where qualitatively the models match at the boundary density, in the region experiencing “wind blanketing” from the additional CAK force, and in the outer wind. Finally, in Fig. 11, the optical depth through the stellar wind is computed. Here, again, the spherically corrected total opacity was used (see Eq. (8) in Poniatowski et al. 2021). As seen from this figure, the simulation spans a wide domain, from the thick core at an optical depth τ_sc ∼ 20 to the optically thin outer wind where τ_sc ∼ 0.005. The energy budget in the wind everywhere is mainly dominated by the radiation, the ratio of energies r rises from ∼10 near the hydrostatic core to ∼300 near the outer edge of the simulation. However, due to the high velocities reached in the outer wind, the gas kinetic energy is even greater than both the radiation and gas internal energies by two orders of magnitude.

Fig. 10. Temperature profile of a spherically symmetric Wolf-Rayet outflow in a relaxed state. In red we show the profile is computed with the 1D FLD method. In black we show it computed following the alternative method by Poniatowski et al. (2021).

Fig. 11. Spherically corrected optical depth of a spherically symmetric Wolf-Rayet outflow in a relaxed state, calculated using the total opacity. In red, the profile is plotted as performed with the 1D FLD method and in black as computed by Poniatowski et al. (2021). In the latter, this is used to calculate the temperature profile.

In the deepest layers of the simulation, near the hydrostatic core, acceleration is essentially ensured by the OPAL opacity while in the region transiting towards the optical photosphere and beyond, the CAK opacity dominates. It highlights the complementary role played by the components: while the OPAL opacity lifts up the material from the dense and hot inner regions, resonant line-absorption not only prevents the flow from falling back but also provides it with additional momentum. When the gas reaches the photosphere, the outflow is already highly supersonic. Overall, the wind launching and final escape is thus made possible thanks to the joint action of both opacities.

In the interest of understanding radiation-powered outflows by means of time-dependent RHD modeling, this simulation illustrates the need for treating the enhanced line-opacity effect in supersonic flows. It also opens the door to the study of time-variable configurations. For instance, it is likely that, when run in a multi-dimensional setup, the lateral symmetry of these Wolf-Rayet models will be broken (see also discussion in Poniatowski et al. 2021). As discussed in the next section, this allows us to study structure formation in a radiation-dominated supersonic environment.

6. Summary and perspectives

In this paper, we describe the implementation of a radiation module for the finite volume magneto-hydrodynamics code MPI-AMRVAC. We validated it with a set of classic benchmark tests and applied it to a more realistic setup, with wind launching and mass loss in the supersonic, expanding atmospheres of Wolf-Rayet stars. The coupling between matter and radiation is performed in the diffusion approximation which provides a closure relation that enables us to deduce the radiative flux and radiative pressure from the energy density of the radiation field. Flux-limiting is applied in order to retrieve the free streaming limit in the optically thin regime, while smoothly transiting to a fully diffusive behavior in highly optically thick environments. The time-dependent evolution equation for the radiative energy density is solved in the co-moving frame to alleviate the angle-dependence of emission, absorption, and scattering induced by the Doppler effect. By default, local thermodynamical equilibrium is not assumed and heat exchanges between matter and the radiative field are accounted for. The opacities which enter the formalism (i.e., the energy, Planck, and flux means) can be prescribed a priori or dynamically computed based on hydrodynamical quantities such as the gas density, temperature, and velocity gradient. Radiative feedback on the ambient gas is ensured by the radiative force in the conservation of gas momentum and by the heating-cooling terms in the conservation of gas energy. In the radiative energy equation, the advection term is handled thanks to the high order approximate Riemann solvers already available in MPI-AMRVAC (Porth et al. 2014; Xia et al. 2018). Photon-tiring is added as an explicit source term while heating and cooling are computed in an implicit way. The diffusive term is treated with the multigrid solver based on a Gauss-Seidel iterative relaxation method introduced in Teunissen & Keppens (2019). This module is fully compatible with the multi-dimensional block-based adaptive mesh refinement at the basis of the domain decomposition strategy of MPI-AMRVAC, which enables MPI-parallelization up to an arbitrarily high number of cores. It performs well on a variety of test cases. Precursors and realistic shock thickness are retrieved in 1D setups of radiatively-dominated shocks. In optically thin environments, front shocks propagate at a speed very close to the speed of light. Galilean invariance is respected and linear damping of a radiative-hydrodynamics wave quantitatively matches the predicted behavior.

We next applied the FLD module to the launching of a radiatively-driven optically thick wind from the hydrostatic core of a Wolf-Rayet star, using a superposition of the standard OPAL opacity tables used in hydrostatics and a simple parametrization of the significantly enhanced line-opacity expected in a supersonic outflow. In agreement with the results obtained by Poniatowski et al. (2021), we find the OPAL opacity to be decisive in the deep and optically thick layers of the star (initiating the supersonic outflow from the so-called “iron-opacity bump” at T ≈ 2 × 10⁵ K), while the line-opacity mechanism takes over in the outer wind, preventing the flow from falling back by bringing the outflow above the local escape speed.

We note, however, that it is far from obvious that this is what would actually take place in a multi-dimensional and time-variable Wolf-Rayet outflow; indeed, in order to make the purely 1D stellar outflow escape we had to make an ad hoc assumption that the line-force in the outer wind is enhanced above the value expected for comparable O-stars in this region (by lowering the so-called CAK-α parameter; see above and also discussion in Poniatowski et al. 2021). In a follow-up paper, we will extend this 1D Wolf-Rayet model to 2D and 3D, in order to investigate the properties of the significant wind structure formation and time variability that presumably will occur if we instead assume more realistic conditions and, thus, also allow for gas that starts to fall back upon the stellar core. The FLD code developed here is ideally suited for this project, as it is fast enough for such a multi-dimensional application while simultaneously accounting for the potential feedback from the structures on the radiative fluxes and forces.

To this end, the FLD module has also further been tested with an an-isotropic diffusion coefficient, which might be of importance when treating line-of-sight line-opacities in a multi-D medium (e.g., Kee et al. 2016). In this formalism, the diffusion constant becomes a diagonal tensor, with the diagonal elements representing the diffusion constant in the direction of each grid line. For simplicity, however, we did not include this aspect explicitly in the paper; the reader and user can readily transform the corresponding notation in Sect. 3. More generally, the extension of MPI-AMRVAC toward general radiation(-magneto)-hydrodynamics provides us with a powerful tool that is suitable for a range of astrophysical applications. The FLD method is a first important step for this, and the Wolf-Rayet outflows discussed above represent a research application that can be directly considered. Another target application for FLD regards “photon-tired” and very optically thick eruptive outflows from massive stars in their luminous blue variable phase (Owocki et al. 2019). Moreover, for the radiation-dominated envelopes of massive stars in general, stellar models often lead to the notion that the radiative acceleration exceeds gravity at the so-called “iron-opacity bump” mentioned above. It is thus possible that this, quite generally, might trigger turbulence in massive-star envelopes and atmospheres, which again might be characterized by co-existing regions of upflows and downflows (see Jiang et al. 2015, for some promising first simulation results). In turn, this might then provide a natural explanation for the very broad photospheric absorption lines typically observed for O-stars, for instance; and this strongly suggests the presence of supersonic velocities that are already in the photosphere (Simón-Díaz et al. 2017).

In this respect, we plan to couple MPI-AMRVAC to the 3D radiative transfer line-formation code by Hennicker et al. (2020) in order to compute post-processed synthetic spectra directly from our dynamical simulations. In addition, this short-characteristics code will provide the basis for another key component of our planned future work, namely, an extension of the FLD method presented here toward full radiative transfer within MPI-AMRVAC, where the Eddington tensor can be computed from the actual RTE instead of an analytic closure relation. Here, an important aspect is related to the careful evaluation of the analytic closure relation applied in the FLD method for various regimes, as well as critical testing (and extension) of the simple line-opacity formalism for supersonic flows described in the previous section.

Acknowledgments

NM and JS acknowledge support by the Belgian Research Foundation Flanders (FWO) Odysseus program under grant number G0H9218N. LP and JS acknowledge support from the KU Leuven C1 grant MAESTRO C16/17/007. IEM has received funding from the Research Foundation Flanders (FWO), from the European Union’s Horizon 2020 research and innovation program under the Marie Skłodowska-Curie grant agreement No 665501 and from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No 863412). JT was supported by postdoctoral fellowship 12Q6117N from Research Foundation – Flanders (FWO). RK received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No. 833251 PROMINENT ERC-ADG 2018). The resources and services used in this work were provided by the VSC (Flemish Supercomputer Center), funded by the Research Foundation – Flanders (FWO) and the Flemish Government. We thank the referee for their constructive feedback on the manuscript.

References

Appendix A: Pseudo-planar correction

Since the multigrid method implemented for the FLD module is not capable of solving the Helmholtz equation on spherical meshes, the full system of PDEs (1), (2), (3), and (6) is solved on a Cartesian grid. This means that for spherical problems, the advection terms have to be modified for spherical fluxes. In the type of simulation presented in Section 5, the calculations are therefore done on a hybrid Cartesian-spherical pseudo-planar grid as presented by Sundqvist et al. (2018). This allows for 1D, 2D, and even 3D simulations of radially extended systems on a Cartesian grid. In a 1D setting, the pseudo planar geometry is equivalent to the r-direction of a spherical geometry. In 2D or 3D, the x-direction of a pseudo planar geometry plays the role of the r-direction of a spherical geometry. If we neglect curvature effects, that is, if the lateral extension of the slab is small compared to its extension along x, then the fluxes along the lateral y and z directions do not require any correction. To illustrate the method, we consider scalar conservation equations such as for (1), (3) and (6). The divergence of the vector f_u in the conservation equation for the conserved quantity u contains a term with partial derivatives in the x-direction in Cartesian coordinates, which we note as ∇_x ⋅ (f), and in the r-direction in spherical coordinates, ∇_r ⋅ (f), respectively given by:

$$\begin{array}{cc}\hfill {\mathrm{\nabla }}_x·{\mathit{f}}_u& ={\partial }_x{f}_{u,x}\hfill \end{array}$$(A.1)

$$\begin{array}{cc}\hfill {\mathrm{\nabla }}_r·{\mathit{f}}_u& =\frac{1}{r^2}{\partial }_r(r^2{f}_{u,r})={\partial }_r{f}_{u,r}+\frac{2f_{u,r}}{r}.\hfill \end{array}$$(A.2)

We can thus assume x ∼ r provided we account for a geometric source term $S_u^g=-2f_{u,r}/r$. The conservation equation on a spherical grid can now be re-written as the conservation equation on a Cartesian grid plus this geometric source term:

$$\begin{array}{c}\hfill {\partial }_{t}u+{\mathrm{\nabla }}_x·{\mathit{f}}_u=S_u^g.\end{array}$$(A.3)

So, for the density, gas energy, and radiation energy:

$$\begin{array}{cc}\hfill S_{\rho }^g& =-\frac{2\rho v_x}{x}\hfill \end{array}$$(A.4)

$$\begin{array}{cc}\hfill S_e^g& =-\frac{2(e+p)v_x}{x}\hfill \end{array}$$(A.5)

$$\begin{array}{cc}\hfill S_E^g& =-\frac{2(E{v}_x+F_x)}{x}·\hfill \end{array}$$(A.6)

For the evolution equation of a vector-like conserved variable, such as momentum, the pseudo-planar correction is different as we work with the divergence of a tensor instead of a vector. In the pseudo planar approach, θ ∼ π/2 since we are working near the equatorial plane of the spherical coordinate system. The coordinates y and z are locally equivalent to rθ and rϕ. The correction term can then be calculated for each of the spatial components of the momentum equation. Due to how the divergence of a tensor is defined, the correction term is different for lateral components:

$$\begin{array}{c}\hfill S_{\rho v_y}^g=-\frac{3(\rho v_x{v}_y)}{x};S_{\rho v_z}^g=-\frac{3(\rho v_x{v}_z)}{x},\end{array}$$(A.7)

as compared to a radial component:

$$\begin{array}{c}\hfill S_{\rho v_x}^g=-\frac{2\rho v_x^2}{x}+\frac{2\rho v_y^2}{x}+\frac{2\rho v_z^2}{x}·\end{array}$$(A.8)

Finally, the radiation work term in equation (6) also features a divergence operator; thus, this one also needs correcting. This additional geometric correction source term in the radiation energy equation is:

$$\begin{array}{c}\hfill S_E^{g,\mathrm{rad}.\mathrm{work}}=\frac{2v_x{\mathit{P}}_{\mathit{xx}}}{x}·\end{array}$$(A.9)

All Figures

Fig. 1. Gas energy densities, as a function of time, of a gas that is relaxing to radiative equilibrium from two sets of initial conditions. Marked with a blue “o” is the relaxation where the initial gas energy is set as 102 times the expected equilibrium value. Marked with a red “x” is the initial gas energy is set as 10−6 times the equilibrium value. The dashed blue and red lines represent the same heating-cooling problem, but solved with a smaller numerical time step (see text). Since the time axis is logarithmic, the initial conditions at t = 0 are not visible.

In the text

	Fig. 2. Density, velocity, gas pressure, and radiation energy profiles for a 1D simulation of a radiation-dominated shock after ten flow passing times. The profiles are normalized to the initial right-hand side values, except for the velocity profile, which has been normalized to the initial left-hand side value. The expected width of the shock is illustrated in grey.
In the text

Fig. 3. Comparison of still and advected solutions in the Galilean invariance test. Top panel: density profiles with initial conditions (in dashed lines) are relaxed once with (marked with ‘o’) and once without (marked with ‘x’) an initial background velocity field. Bottom panel: same details but for the gas pressure (blue) and radiation energy density (red). The horizontal shift between the curves correspond to the expected rightward displacement.

In the text

Fig. 4. Absolute value of the relative difference between the profiles with and the profiles without an initial background velocity field from the Galilean invariance test in Fig. 3, where the advected profile is translated back to the origin. We compare two different time integration schemes: a first-order Euler IMEX scheme (marked ‘+’) and a second-order Midpoint IMEX scheme (marked ‘⋅’). Different colors correspond to different quantities: gas density in black, gas pressure in blue, and radiation energy in red.

In the text

Fig. 5. Space-time diagram of the radiation energy front evolving in an optically thin slab of gas. In the optically thin limit, the FLD flux limiter λ prohibits radiation from moving faster than the speed of light. In this plot, the speed of light has been indicated by the black dotted lines. In addition to the colormap, the position of the radiation front has been evaluated at three arbitrarily values with the red contour lines. The total optical thickness of the slab is 0.002.

In the text

Fig. 6. Comparison between numerical and analytical solutions of the RHD wave. From top to bottom: the density perturbation, velocity, gas energy perturbation, and radiation energy perturbation after 40 oscillation periods, for an RHD wave where the optical thickness of a wavelength is equal to 1000. For the top three panels, this is compared to the analytic profile (black dotted lines) that is calculated as described by Mihalas & Mihalas (1984). The wave is driven from the grey shaded area on the left side of the plot.

In the text

Fig. 7. Comparison between 1D and 2D numerical solutions of the RHD wave. Top panel: density perturbation profile of an RHD wave after 40 oscillation periods. In black, the analytically predicted profile, in blue results from a 1D simulation and in red results from a 2D simulation with the wave vector parallel to the first diagonal. The x-axis is here chosen in the direction of the wave vector. Bottom panel: respective residuals calculated as |ρ − ρanalytic|/ρanalytic, where the analytic density is calculated from perturbation theory.

In the text

	Fig. 8. Map of density perturbation in an RHD wave along the first diagonal (black dotted line) of a 2D plane, after 40 oscillation periods. The wave is driven from the grey area in the bottom-left corner.
In the text

	Fig. 9. Radial velocity profile of a spherically symmetric Wolf-Rayet outflow in a relaxed state. In red, the profile is computed with the 1D FLD method while in black, it is computed following the alternative method by Poniatowski et al. (2021).
In the text

	Fig. 10. Temperature profile of a spherically symmetric Wolf-Rayet outflow in a relaxed state. In red we show the profile is computed with the 1D FLD method. In black we show it computed following the alternative method by Poniatowski et al. (2021).
In the text

	Fig. 11. Spherically corrected optical depth of a spherically symmetric Wolf-Rayet outflow in a relaxed state, calculated using the total opacity. In red, the profile is plotted as performed with the 1D FLD method and in black as computed by Poniatowski et al. (2021). In the latter, this is used to calculate the temperature profile.
In the text