Issue 
A&A
Volume 657, January 2022



Article Number  A81  
Number of page(s)  14  
Section  Numerical methods and codes  
DOI  https://doi.org/10.1051/00046361/202141023  
Published online  14 January 2022 
Radiationhydrodynamics with MPIAMRVAC
Fluxlimited diffusion
^{1}
Instituut voor Sterrenkunde, KU Leuven, Celestijnenlaan 200D, 3001 Leuven, Belgium
email: nicolas.moens@kuleuven.be
^{2}
Centre for mathematical Plasma Astrophysics, Department of Mathematics, KU Leuven, Celestijnenlaan 200B, 3001 Leuven, Belgium
^{3}
Centrum Wiskunde and Informatica, PO Box 94079 1090 GB Amsterdam, The Netherlands
^{4}
Institut de Planétologie et d’Astrophysique de Grenoble (IPAG), Université Grenoble Alpes, 38058 Grenoble Cedex 9, France
Received:
8
April
2021
Accepted:
26
September
2021
Context. Radiation controls the dynamics and energetics of many astrophysical environments. To capture the coupling between the radiation and matter, however, is often a physically complex and computationally expensive endeavor.
Aims. We sought to develop a numerical tool to perform radiationhydrodynamics simulations in various configurations at an affordable cost.
Methods. We built upon the finite volume code MPIAMRVAC to solve the equations of hydrodynamics on multidimensional adaptive meshes and introduce a new module to handle the coupling with radiation. A nonequilibrium, fluxlimiting diffusion approximation was used to close the radiation momentum and energy equations. The timedependent radiation energy equation was then solved within a flexible framework, fully accounting for radiation forces and work terms and further allowing the user to adopt a variety of descriptions for the radiationmatter interaction terms (“opacities”).
Results. We validated the radiation module on a set of standard test cases for which different terms of the radiative energy equation predominate. As a preliminary application to a scientific case, we calculated spherically symmetric models of the radiationdriven and optically thick supersonic outflows from massive WolfRayet stars. This also demonstrates our code’s flexibility, as the illustrated simulation combines opacities typically used in static stellar structure models with a parametrized form for the enhanced lineopacity expected in supersonic flows.
Conclusions. This new module provides a convenient and versatile tool for performing multidimensional and highresolution radiativehydrodynamics simulations in optically thick environments with the MPIAMRVAC code. The code is ready to be used for a variety of astrophysical applications, where our first target is set to be multidimensional simulations of stellar outflows from WolfRayet stars.
Key words: radiation: dynamics / methods: numerical / stars: WolfRayet / hydrodynamics
© 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 starforming 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 solartype stars, the heating and cooling terms are computed directly from solutions to the timeindependent radiative transfer equation (RTE), using a sort of “frozenin” 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 timedependent simulations typically rely on various distributionfunction approaches (Owocki et al. 1988; Sundqvist et al. 2018).
For general attempts to solve the full timedependent 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 radiationmatter 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 (MignonRisse et al. 2020). Concerning the radiationmatter interaction terms, the central issue is that typically the opacities are isotropic only in the frame that is comoving 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 tradeoffs 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 “mixedframe” formulation, however, is that it is not suitable for computing line opacities in a supersonic flow, since the firstorder expansion upon which it is based does not apply for the rapidly changing opacities of spectral lines.
In brief, MPIAMRVAC^{1} is a modern, MPIparallelized 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 MPIAMRVAC, 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 nonequilibrium FLD method into MPIAMRVAC. For many of our targeted applications lineopacity 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 lineopacities within the formalism. Recently, MPIAMRVAC 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 radiativehydrodynamics (RHD) equations and we describe the nonequilibrium FLD method for obtaining closure. The different aspects of the implementation in the MPIAMRVAC 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 radiationdriven wind outflow from a classical WolfRayet 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:
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 righthand side, f_{r} is the radiation force density exerted on the matter and 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:
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 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 MPIAMRVAC 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 and f_{r} requires radiationmaterial interaction terms, where these involve an important angledependence 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 comoving with the local velocity of the fluid. In this comoving frame (CMF), the frequencyintegrated radiation energy and momentum equations are (Mihalas & Mihalas 1984; Castor 2004):
where E, F, and P are the frequencyintegrated 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 comoving frames). For example, the fourth term on the lefthand 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 righthandsides in Eqs. (6) and (7). These coupling terms in Eqs. (2), (3), (6) and (7) can now be written as:
where is the frequencyintegrated Planck function, with σ the StefanBoltzmann constant, and the opacities κ_{B} and κ_{E} are mass absorption coefficients measured in cm^{2} g^{−1}. Specifically, Eqs. (8) and (9) involve the Planck, energy density, and flux mean opacities:
which follows directly from considering the frequencydependent 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 . 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 (see Eq. (77) in Mihalas & Mihalas 1984, their page 336, and also the corresponding discussion on their page 472).
2.2. Nonequilibrium 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 multifrequency solutions to the radiative transfer equation in the CMF have thus far only been developed for 1D, steadystate 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 radiationhydrodynamics 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 socalled FLD approximation as our closure relation. Neglecting the first two terms in the radiation momentum equation, we have:
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):
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 meanfree 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 fluxlimiter λ suggested by Levermore & Pomraning (1981) (see also Turner & Stone 2001), expressed as:
and for the bridging law:
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 MPIAMRVAC code, the user can readily switch between different flux limiters like the Levermore (Levermore & Pomraning 1981) or the Minerbo (Minerbo 1978) prescriptions.
Using this fluxlimiter within Fick’s diffusion formulation, we can locally compute the comoving radiation flux according to:
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:
We also obtain the radiation pressure tensor and so we can omit Eq. (7) entirely. Here, 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 heatingcooling term , we further assume for the specific application in this paper that the Planck and energy density mean opacities are equal, κ_{P} = κ_{E} = κ, such that:
where we have reformulated E in terms of a radiation temperature, , where a_{r} = 4σ/c is the radiation constant. In the nonequilibrium 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 nonequilibrium FLD approximation consist of Eqs. (1)–(3) and (6). The conservative lefthand side parts of the hydrodynamic Eqs. (1)–(3) are entirely hyperbolic and can be solved using the existing variety of shock capturing, highresolution finite volume solvers in MPIAMRVAC (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 MPIAMRVAC for solving hyperbolic equations. Unless stated otherwise, in the test cases and applications presented in the following sections, we use the secondorder accurate shockcapturing total variation diminishing LaxFriedrichs (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 righthand 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 heatingcooling 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, 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 VACpreprocessor. 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 radiationMHD calculations.
3.1. Order of operations
Generally, the PDEs of the FLD system, Eqs. (1)–(3) and (6), can be written as:
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, are the entirely local source terms which are handled implicitly: cooling and heating. In our implementation in MPIAMRVAC, these four types of terms can be grouped in two operators: an implicit operator and an explicit operator:
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^{n} to is:
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 . This intermediate state is used in the local implicit heating and cooling terms as explained in Sect. 3.3 to retrieve another intermediate state . 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 secondorder 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^{n} to u^{n + 1} as follows:
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^{n}. 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^{n} to u^{n + 1} where G_{ex} computed from the midpoint state u^{n + 1/2} and the value of G_{im} is reused 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 CourantFriedrichsLewy (CFL) condition (Courant et al. 1928) is implemented. For radiationfree 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 radiativeacoustic 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 fivepoint stencil (Fornberg 1988). For instance, in a 2D setup, the xcomponent of the gradient is expressed as:
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:
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 fivestep 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:
Also, this source term is added explicitly in every cell. For the addition of Δt S_{ex} in the radiation energy density:
3.3. Heating and cooling
As mentioned above, our nonequilibrium 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 . 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 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}/2. Using (4), a pointimplicit, discretized formulation of the heating and cooling terms applied to the gas internal and radiation energy density are written as:
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 fourthdegree polynomial (Turner & Stone 2001):
where a_{1} = 4κσ(γ − 1)^{4}/ρ^{3}Δt and a_{2} = cκρΔt, with both κ and ρ evaluated at step n. To get the gas internal ϵ^{n + 1} energy at the next timestep, Halley’s rootfinding method (Press et al. 2007) is employed. This method is similar to NewtonRaphson, 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, ϵ^{n} + E^{n}). With the updated internal gas energy, the radiation energy density is given by:
In many astrophysical contexts, gases are radiationdominated (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 radiationdominated gas.
3.4. Radiation diffusion
For the diffusion term in Eq. (6), we make use of the newly developed geometric multigrid method library octreemg (Teunissen & Keppens 2019). This MPIparallelized library is capable of solving elliptic PDEs on 1D,2D or 3D Cartesian grids and is fully compatible with the blocktree AMR structure used in MPIAMRVAC. 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), ∂_{t}E + ∇ ⋅ F = 0, with F substituted using (14):
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 smallscale errors are smoothed on a finer grid. The multigrid library uses a standard GaussSeidel 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^{n} is the current radiation energy density and E^{n + 1} is the radiation energy density in the next time step:
As mentioned above, our methods are implemented in 1D, 2D, and 3D. In 2D, this equation is discretized using the following fivepoint stencil:
In addition to an initial state E^{n}, 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):
Locally, the diffusion constant in the cell centers is calculated from Eq. (18), where the fluxlimiter λ 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^{−5}. In an FMG cycle, the smoothing starts at the coarsest grid of the MGsolver, and then prolongs, smoothens and restricts to ever increasing finer grids. If the FMGcycle does not reach this tolerance, additional Vcycles are added until the convergence criterion is met. A Vcycle 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 nonradiative HD/MHD setup. These limitations in the multigrid solver arise from using pointwise 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 nontrivial analytic solutions to the FLD equations are scarce, we set up a number of testcases and then compare them to various semianalytic 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 nonequilibrium, there is a net heatingcooling term , 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^{−7} g cm^{−3}, opacity κ = 0.4 cm^{2} g^{−1}, adiabatic index γ = 5/3, and mean molecular weight μ = 0.6. The radiation energy density is initialized as E = 10^{12} erg cm^{−3}.
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^{−6} e_{equil} (‘x’ symbols in Fig. 1) and in the second test, it is set to 10^{2} 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 radiationdominated, 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 10^{2} 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 heatingcooling 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 loglog 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^{−12} s, which is regarded as a physically correct baseline; and once where Δt = 10^{−11} 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^{−11} s used here is greater than the initial cooling timescale, τ_{cool} ∼ e/(4πκρB)∼10^{−12} 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^{−5} s). In this case, the relative error in temperature converged to the requested tolerance after three time steps.
4.2. Radiationdominated shock
This test case here describes a steadystate radiationdominated 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 radiationmodified RankineHugoniot conditions, a left and righthand 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 lefthand side of the shock, gas density ρ_{l} = 10^{−2} g cm^{−3}, velocity v_{l} = 10^{9} cm s^{−1} and total (kinetic + internal) gas energy density e_{l} = 5.0 × 10^{15} erg cm^{−3}. The radiation energy density is set in equilibrium to E_{l} = 75.6 erg cm^{−3}. On the righthand side ρ_{r} = 6.85874 × 10^{−2} g cm^{−3}, v_{r} = 1.45 × 10^{8} cm s^{−1}, e_{r} = 1.93 × 10^{15} erg cm^{−3}, and E_{r} = 2.44 × 10^{16} erg cm^{−3}. The stability of the shock depends on having the initial values conform with the radiationmodified Rankine Hugoniot conditions. Empirically, at least five decimals for the righthand 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^{5} 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 radiationdominated shock after ten flow passing times. The profiles are normalized to the initial righthand side values, except for the velocity profile, which has been normalized to the initial lefthand 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 xaxis (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 highresolution 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:
where T_{0} = 10^{7} K, T_{1} = 2 × 10^{7} K, ρ_{0} = 1.2 g cm^{−3} 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^{2} g^{−1}. 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^{7} cm s^{−1}. 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 IMEXscheme. 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 timestepping 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 firstorder Euler IMEX scheme (marked ‘+’) and a secondorder 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 fluxlimited, 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} = 0.025 g cm^{−3} and constant opacity κ = 0.4 cm^{2} g^{−1}. 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:
This expression lights up the slab from the lefthand side, with E_{1} = 1.4 × 10^{11} erg cm^{−3}. There is no velocity and the gas energy everywhere is set to be in equilibrium with E_{0} = 1.4 × 10^{−11} erg cm^{−3}. 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^{−13} s, which is on the order of half a cellcrossing 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.
Spacetime 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 semianalytic 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_{a}e/(cE), with 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 energydominated and low when the gas is gas energydominated.
A slab of gas is considered with background values ρ_{0} = 3.216 × 10^{−9} g cm^{−3}, e_{0} = 26.02 × 10^{3} erg cm^{−3}, and E_{0} = 17.34 × 10^{3} erg cm^{−3}. The adiabatic index is γ = 5/3 leading to a Boltzman number of Bo = 10^{−3} and an energy density ratio of r = 0.1. The slab is perturbed on the lefthand side with the following boundary conditions to excite a traveling wave:
Here, the wavenumber k is chosen in such a way that for the constant opacity of κ = 0.4 cm^{2} g^{−1}, the optical depth across one wavelength is τ_{λ} = 10^{3}. The pulsation ω is related to k via the propagation speed of the RHD wave. The perturbation with density amplitude A_{ρ} = 0.01 ρ_{0} 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ρ_{0}) and A_{e} = A_{ρ}e_{0}/ρ_{0}. 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 rootfinding 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 fluxlimiter introduced in Sect. 2.2).
In Fig. 6, results are shown for a perturbation with the optical depth τ_{λ} = 10^{3} 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 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 nonoscillating 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 xaxis 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 bottomleft 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 WolfRayet wind outflow. Classical WolfRayet stars are massive stars that have evolved back to the blue side of the HertzsprungRussel 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 massloss rates. WolfRayet winds are believed to be accelerated by radiation (Sander & Vink 2020; Poniatowski et al. 2021), but due to the high massloss 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 setup 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 timedependent modeling of such WolfRayet outflows. Specifically, while the simulations by Poniatowski et al. (2021) assumed that the local radial flux is always set by L_{*}/(4πr^{2}), 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 timedependent 2D or 3D flows with local (and potentially nonradial) 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 WolfRayet simulation, we assume a fixed stellar mass M_{*} = 10 M_{⊙}, a gravity source term f_{g} = ρGM_{*}/r^{2}, a hydrostatic core lower boundary radius R_{*} = 1 R_{⊙}, and a stellar core luminosity of log_{10}(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 lineopacity expected in a supersonic flow:
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 Dopplershifted lines. At every time step, the CAKopacity is computed locally using a second order central difference derivative of the radial velocity with respect to radius (udDoula & Owocki 2002):
This combined opacity is then used in all source terms (radiative force, heating or cooling, photontiring) as well as in the computation of the diffusion coefficient.
In Eq. (50), represents the line opacity in the limit that all lines would be optically thin, and α represents the slope of the underlying assumed powerlaw 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 outerwind 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 and a socalled βvelocity law v(r) (see below). From L_{*} = 4πr^{2}F_{r} and the Eddington approximation 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^{4} 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:
where we set v_{∞} = 1000 km s^{−1} 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 radiationdriven 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:
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 3step scheme with a TVDLF solver and a minmod slope limiter. For this numerical setup, 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 abovedescribed 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 nonmonotonic 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 massloss rate for the FLD model is 1.87 × 10^{−5} M_{⊙} yr^{−1}, which again is very similar to the 1.64 × 10^{−5} M_{⊙} yr^{−1} found by Poniatowski et al. (2021).
Fig. 9.
Radial velocity profile of a spherically symmetric WolfRayet 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 WolfRayet 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 WolfRayet 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 lineabsorption 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 radiationpowered outflows by means of timedependent RHD modeling, this simulation illustrates the need for treating the enhanced lineopacity effect in supersonic flows. It also opens the door to the study of timevariable configurations. For instance, it is likely that, when run in a multidimensional setup, the lateral symmetry of these WolfRayet 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 radiationdominated supersonic environment.
6. Summary and perspectives
In this paper, we describe the implementation of a radiation module for the finite volume magnetohydrodynamics code MPIAMRVAC. 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 WolfRayet 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. Fluxlimiting 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 timedependent evolution equation for the radiative energy density is solved in the comoving frame to alleviate the angledependence 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 heatingcooling 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 MPIAMRVAC (Porth et al. 2014; Xia et al. 2018). Photontiring 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 GaussSeidel iterative relaxation method introduced in Teunissen & Keppens (2019). This module is fully compatible with the multidimensional blockbased adaptive mesh refinement at the basis of the domain decomposition strategy of MPIAMRVAC, which enables MPIparallelization 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 radiativelydominated 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 radiativehydrodynamics wave quantitatively matches the predicted behavior.
We next applied the FLD module to the launching of a radiativelydriven optically thick wind from the hydrostatic core of a WolfRayet star, using a superposition of the standard OPAL opacity tables used in hydrostatics and a simple parametrization of the significantly enhanced lineopacity 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 socalled “ironopacity bump” at T ≈ 2 × 10^{5} K), while the lineopacity 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 multidimensional and timevariable WolfRayet outflow; indeed, in order to make the purely 1D stellar outflow escape we had to make an ad hoc assumption that the lineforce in the outer wind is enhanced above the value expected for comparable Ostars in this region (by lowering the socalled CAKα parameter; see above and also discussion in Poniatowski et al. 2021). In a followup paper, we will extend this 1D WolfRayet 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 multidimensional 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 anisotropic diffusion coefficient, which might be of importance when treating lineofsight lineopacities in a multiD 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 MPIAMRVAC 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 WolfRayet outflows discussed above represent a research application that can be directly considered. Another target application for FLD regards “photontired” and very optically thick eruptive outflows from massive stars in their luminous blue variable phase (Owocki et al. 2019). Moreover, for the radiationdominated envelopes of massive stars in general, stellar models often lead to the notion that the radiative acceleration exceeds gravity at the socalled “ironopacity bump” mentioned above. It is thus possible that this, quite generally, might trigger turbulence in massivestar envelopes and atmospheres, which again might be characterized by coexisting 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 Ostars, for instance; and this strongly suggests the presence of supersonic velocities that are already in the photosphere (SimónDíaz et al. 2017).
In this respect, we plan to couple MPIAMRVAC to the 3D radiative transfer lineformation code by Hennicker et al. (2020) in order to compute postprocessed synthetic spectra directly from our dynamical simulations. In addition, this shortcharacteristics 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 MPIAMRVAC, 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 lineopacity 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łodowskaCurie 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 ERCADG 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
 Bloch, H., Tremblin, P., González, M., Padioleau, T., & Audit, E. 2021, A&A, 646, A123 [Google Scholar]
 Briggs, W. L., Henson, V. E., & Mccormick, S. F. 2000, A Multigrid Tutorial, 2nd edn. (Philadelphia: Society for Industrial and Applied Mathematics) [Google Scholar]
 Castor, J. I. 2004, Radiation Hydrodynamics (Cambridge: Cambridge University Press) [Google Scholar]
 Castor, J. I., Abbott, D. C., & Klein, R. I. 1975, ApJ, 195, 157 [Google Scholar]
 Courant, R., Friedrichs, K., & Lewy, H. 1928, Math. Ann., 100, 32 [Google Scholar]
 Crowther, P. A. 2007, ARA&A, 45, 177 [Google Scholar]
 Driessen, F. A., Sundqvist, J. O., & Wade, G. A. 2019, Proc. Int. Astron. Union, 346, 45 [Google Scholar]
 Fornberg, B. 1988, Math. Comput., 51, 699 [Google Scholar]
 Gábor, T., & Dušan, O. 1996, J. Comput. Phys., 128, 82 [Google Scholar]
 Hennicker, L., Puls, J., Kee, N. D., & Sundqvist, J. O. 2020, A&A, 633, A16 [Google Scholar]
 Hillier, D. J., & Miller, D. L. 1998, ApJ, 496, 407 [Google Scholar]
 Iglesias, C. A., & Rogers, F. J. 1996, ApJ, 464, 943 [Google Scholar]
 Jiang, Y.F., Stone, J. M., & Davis, S. W. 2012, ApJS, 199, 14 [Google Scholar]
 Jiang, Y.F., Cantiello, M., Bildsten, L., Quataert, E., & Blaes, O. 2015, ApJ, 813, 74 [Google Scholar]
 Jiang, Y.F., Stone, J. M., & Davis, S. W. 2019, ApJ, 880, 67 [Google Scholar]
 Kee, N. D., & Kuiper, R. 2019, MNRAS, 483, 4893 [Google Scholar]
 Kee, N. D., Owocki, S., & Sundqvist, J. O. 2016, MNRAS, 458, 2323 [Google Scholar]
 Keppens, R., Teunissen, J., Xia, C., & Porth, O. 2021, Comput. Math. Appl., 81, 316 [Google Scholar]
 Koren, B., & Vreugdenhil, C. 1993, Numerical Methods for AdvectionDiffusion Problems (Braunschweig: Vieweg) [Google Scholar]
 Krumholz, M. R., Klein, R. I., McKee, C. F., & Bolstad, J. 2006, ApJ, 667, 626 [Google Scholar]
 Lattimer, A. S., & Cranmer, S. R. 2021, ApJ, 910, 48 [Google Scholar]
 Levermore, G. D., & Pomraning, G. C. 1981, ApJ, 248, 321 [Google Scholar]
 Liu, X. D. 1994, J. Comput. Phys., 115, 200 [Google Scholar]
 Lowrie, R. B., Morel, J. E., & Hittinger, J. A. 1999, ApJ, 521, 432 [Google Scholar]
 MignonRisse, R., González, M., Commerçon, B., & Rosdahl, J. 2020, A&A, 635, A42 [Google Scholar]
 Mihalas, D., & Klein, R. I. 1982, J. Comput. Phys., 46, 97 [Google Scholar]
 Mihalas, D., & Mihalas, B. W. 1984, Foundations of Radiation Hydrodynamics (New York: Oxford University Press), 731 [Google Scholar]
 Minerbo, G. N. 1978, J. Quant. Spectr. Rad. Transf., 20, 541 [Google Scholar]
 Owocki, S. P., & Gayley, K. G. 1997, ASP Conf. Ser., 120, 121 [Google Scholar]
 Owocki, S. P., Castor, J. I., & Rybicki, G. B. 1988, ApJ, 335, 914 [Google Scholar]
 Owocki, S. P., Hirai, R., Podsiadlowski, P., & Schneider, F. R. N. 2019, MNRAS, 485, 988 [Google Scholar]
 Poniatowski, L. G., Sundqvist, J. O., Kee, N. D., et al. 2021, A&A, 647, A151 [Google Scholar]
 Porth, O., Xia, C., Hendrix, T., Moschou, S. P., & Keppens, R. 2014, ApJS, 214, 4 [Google Scholar]
 Press, W. H., Teukolsky, S. A., Vetterling, W. T., & Flannery, B. P. 2007, Numerical Recipes 3rd Edition: The Art of Scientific Computing, 3rd edn. (Cambridge: Cambridge University Press) [Google Scholar]
 Proga, D., & Kallman, T. R. 2004, ApJ, 616, 688 [Google Scholar]
 Puls, J., Springmann, U., & Lennon, M. 2000, A&AS, 141, 23 [Google Scholar]
 Puls, J., Najarro, F., Sundqvist, J. O., & Sen, K. 2020, A&A, 642, A172 [EDP Sciences] [Google Scholar]
 Roe, P. L. 1986, Annu. Rev. Fluid Mech., 18, 337 [Google Scholar]
 Sander, A. A. C., & Vink, J. S. 2020, MNRAS, 499, 873 [Google Scholar]
 Sander, A. A. C., Hamann, W. R., Todt, H., Hainich, R., & Shenar, T. 2017, A&A, 603, A86 [Google Scholar]
 Schure, K. M., Kosenko, D., Kaastra, J. S., Keppens, R., & Vink, J. 2009, A&A, 508, 751 [Google Scholar]
 SimónDíaz, S., Godart, M., Castro, N., et al. 2017, A&A, 597, A22 [Google Scholar]
 Skinner, M. A., & Ostriker, E. C. 2013, ApJS, 206, 21 [Google Scholar]
 Stein, R. F., & Nordlund, Å. 1998, ApJ, 499, 914 [Google Scholar]
 Sundqvist, J. O., & Owocki, S. P. 2013, MNRAS, 428, 1837 [Google Scholar]
 Sundqvist, J. O., Owocki, S. P., & Puls, J. 2018, A&A, 611, A17 [Google Scholar]
 Sundqvist, J. O., Björklund, R., Puls, J., & Najarro, F. 2019, A&A, 632, A126 [Google Scholar]
 Suresh, A., & Huynh, H. T. 1997, J. Comput. Phys., 136, 83 [Google Scholar]
 Teunissen, J., & Keppens, R. 2019, Comput. Phys. Commun., 245, 106866 [Google Scholar]
 Turner, N. J., & Stone, J. M. 2001, ApJS, 135, 95 [Google Scholar]
 udDoula, A., & Owocki, S. P. 2002, ApJ, 576, 413 [Google Scholar]
 van Marle, A. J., & Keppens, R. 2011, Comput. Fluids, 42, 44 [Google Scholar]
 Xia, C., Teunissen, J., Mellah, I. E., Chane, E., & Keppens, R. 2018, ApJS, 234, 30 [Google Scholar]
Appendix A: Pseudoplanar 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 Cartesianspherical pseudoplanar 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 rdirection of a spherical geometry. In 2D or 3D, the xdirection of a pseudo planar geometry plays the role of the rdirection 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 xdirection in Cartesian coordinates, which we note as ∇_{x} ⋅ (f), and in the rdirection in spherical coordinates, ∇_{r} ⋅ (f), respectively given by:
We can thus assume x ∼ r provided we account for a geometric source term . The conservation equation on a spherical grid can now be rewritten as the conservation equation on a Cartesian grid plus this geometric source term:
So, for the density, gas energy, and radiation energy:
For the evolution equation of a vectorlike conserved variable, such as momentum, the pseudoplanar 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:
as compared to a radial component:
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:
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 10^{2} 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 heatingcooling 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 radiationdominated shock after ten flow passing times. The profiles are normalized to the initial righthand side values, except for the velocity profile, which has been normalized to the initial lefthand 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 firstorder Euler IMEX scheme (marked ‘+’) and a secondorder 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.
Spacetime 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 xaxis 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 bottomleft corner. 

In the text 
Fig. 9.
Radial velocity profile of a spherically symmetric WolfRayet 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 WolfRayet 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 WolfRayet 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 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.