© ESO, 2015

1. Introduction

From Poynting-Robertson drag on dust grains to cosmological reionisation, and at every scale in between, there is no denying the importance of radiation in astrophysics. In particular, in accretion disks around young stars the effects of radiative cooling and stellar irradiation are crucial to determining the thermal and physical structure of the disk (e.g., flaring; Chiang & Goldreich 1997).

Until recently, models of these protoplanetary accretion disks that include radiation have been restricted to 1 + 1D dynamic or 2D and 3D static models (e.g., Gorti et al. 2009; Bruderer et al. 2010; Min et al. 2011). However, many phenomena, such as planetary migration (e.g., Kley et al. 2009), or the magneto-rotational instability (MRI; e.g., Flaig et al. 2010) are certainly dynamic and multidimensional, and so require time-dependent radiative transfer (RT) coupled to dynamics.

One approach is to self-consistently couple the equations of hydrodynamics (HD) and radiation transport. Although the combined subject of radiation hydrodynamics (RHD) is treated with great detail in many textbooks (e.g., Mihalas & Weibel Mihalas 1984), the development of practical methods for numerical modelling still remains an active area of astrophysical research.

A recent approach to RHD in accretion disks, made popular by Kuiper et al. (2010), is to combine flux-limited diffusion (FLD) with a simple and fast ray-tracing approach for the direct irradiation by stellar photons. In this algorithm, the direct stellar flux is followed until it is absorbed, and the resultant heat is added to the fluid. Any subsequent re-emission is handled as diffuse radiation using FLD. Alternatively, one could state that the stellar photons are followed until the surface of first absorption, whereupon they heat the fluid. Treating the stellar flux in this hybrid manner preserves the directional properties of the stellar radiation significantly better than with FLD alone, and the result can even be competitive with more complicated and expensive Monte Carlo methods (Kuiper et al. 2010; Kuiper & Klessen 2013). It is also well known that FLD is incapable of casting shadows (Hayes & Norman 2003), a potentially important effect for irradiated protoplanetary disks with puffed up inner rims (e.g., Dullemond et al. 2001), but which is possible with the hybrid irradiation + FLD algorithm.

Although this frequency-dependent algorithm is straightforward to implement and fast, the trade-off is that rays cannot travel in arbitrary directions: the irradiation flux is restricted to travelling parallel to coordinate axes. This naturally presents challenges for general circumstances, in particular scattering of radiation or multiple sources, but when applied to accretion disks using spherical coordinates, it is very efficient.

In this regard, there are certainly more accurate methods for radiation hydrodynamics available, including the variable Eddington tensor factor (e.g., Jiang et al. 2012), M1 closure (e.g., González et al. 2007), and Monte Carlo plus hydrodynamics methods (e.g., Haworth & Harries 2012). However, these methods are normally much more computationally expensive than the hybrid method discussed here.

It is also worth noting that the splitting of radiation into multiple components, whether diffusion, ray-tracing, stellar heating, or other, is not a new idea, and several examples in different contexts can be found in the literature (e.g., Wolfire & Cassinelli 1986; Murray et al. 1994; Edgar & Clarke 2003; Whalen & Norman 2006; Boley et al. 2007; Dobbs-Dixon et al. 2010).

We present a new radiation module for the AZEuS adaptive mesh refinement (AMR) magnetohydrodynamics (MHD) code (Ramsey et al. 2012), with a focus on protoplanetary disk (PPD) applications. We have combined frequency-dependent ray-tracing (Kuiper et al. 2010) with non-equilibrium (two-temperature) FLD (e.g., Bitsch et al. 2013; Kolb et al. 2013; Flock et al. 2013) and, in addition, coupled it to AMR (Berger & Colella 1989).

While there are other AMR FLD codes available in the astrophysical literature¹, including ORION (Krumholz et al. 2007), CRASH (van der Holst et al. 2011), CASTRO (Zhang et al. 2011), RAMSES (Commerçon et al. 2014), and FLASH (Klassen et al. 2014), AZEuS is the only AMR fluid code currently available which employs a fully staggered mesh, where the momentum and magnetic field are both face-centred quantities. Furthermore, as it is based on the ZEUS family of codes, it derives from one of the best documented, tested, and widely used codes in astrophysics.

The paper layout is as follows. In Sect. 2 we describe our algorithms for radiation hydrodynamics in the uniform grid case. In Sect. 3, we present out extensions to both FLD and irradiation for use with AMR. In Sect. 4, we perform tests to validate our algorithms. Finally, in Sect. 5, we summarise our results and present an outlook for future work.

2. Radiation hydrodynamics in AZEuS

AZEuS solves the time-dependent, frequency-integrated equations of RHD in the co-moving frame, to order , and assuming local thermodynamic equilibrium $\begin{array}{ccc}& & \frac{\mathit{\partial \rho }}{\mathit{\partial t}}\mathrm{+}\mathrm{\nabla }\mathrm{·}\left(\mathit{\rho }{v}\right)\mathrm{=}\mathrm{0}\mathrm{;}\\ & & \frac{\mathit{\partial }{s}}{\mathit{\partial t}}\mathrm{+}\mathrm{\nabla }\mathrm{·}\mathrm{(}{s}{{v}}^{\mathrm{)}}\mathrm{=}\mathrm{\nabla }\mathit{p}\mathrm{-}\mathrm{\nabla }\mathrm{·}\mathrm{Q}\mathrm{-}\mathit{\rho }\mathrm{\nabla \Phi }\mathrm{+}\mathrm{\nabla }\mathrm{·}\mathrm{T}\mathrm{+}\frac{\mathit{\rho }{\mathit{\kappa }}_{\mathrm{F}}}{\mathit{c}}{{F}}_{\mathrm{R}}\mathrm{;}\\ & & \frac{\mathit{\partial e}}{\mathit{\partial t}}\mathrm{+}\mathrm{\nabla }\mathrm{·}\mathrm{(}\mathit{e}{{v}}^{\mathrm{)}}\mathrm{=}\mathrm{-}\mathit{p}\mathrm{\nabla }\mathrm{·}{v}\mathrm{+}\mathit{c\rho }\mathrm{(}{\mathit{\kappa }}_{\mathrm{E}}{\mathit{E}}_{\mathrm{R}}\mathrm{-}{\mathit{\kappa }}_{\mathrm{P}}{\mathit{a}}_{\mathrm{R}}{\mathit{T}}^{\mathrm{4}}\mathrm{)}\\ & & \mathrm{+}\mathrm{T:\nabla }{v}\mathrm{-}\mathrm{Q:\nabla }{v}\mathrm{;}\\ & & \frac{\mathit{\partial }{\mathit{E}}_{\mathrm{R}}}{\mathit{\partial t}}\mathrm{+}\mathrm{\nabla }\mathrm{·}\mathrm{(}{\mathit{E}}_{\mathrm{R}}{{v}}^{\mathrm{)}}\mathrm{=}\mathrm{-}\mathit{c\rho }\mathrm{(}{\mathit{\kappa }}_{\mathrm{E}}{\mathit{E}}_{\mathrm{R}}\mathrm{-}{\mathit{\kappa }}_{\mathrm{P}}{\mathit{a}}_{\mathrm{R}}{\mathit{T}}^{\mathrm{4}}\mathrm{)}\mathrm{-}\mathrm{\nabla }\mathrm{·}{{F}}_{\mathrm{R}}\mathrm{-}\mathrm{P:\nabla }{v}\mathrm{;}\\ & & \frac{\mathit{\partial }{\mathit{E}}_{\mathrm{T}}}{\mathit{\partial t}}\mathrm{+}\mathrm{\nabla }\mathrm{·}(\mathrm{(}{\mathit{E}}_{\mathrm{T}}\mathrm{+}\mathit{p}\mathrm{)}{v}\mathrm{+}\mathrm{Q}\mathrm{·}{v})\mathrm{=}\mathrm{-}\mathrm{\nabla }\mathrm{·}{{F}}_{\mathrm{R}}\mathrm{-}\mathrm{P:\nabla }{v}\\ & & \hspace{0.17em} \mathrm{+}\mathrm{\nabla }\mathrm{·}\mathrm{(}\mathrm{T}\mathrm{·}{{v}}^{\mathrm{)}}\mathrm{+}\frac{\mathit{\rho }{\mathit{\kappa }}_{\mathrm{F}}}{\mathit{c}}{v}\mathrm{·}{{F}}_{\mathrm{R}}\mathit{,}\end{array}$where ρ is the mass density; v is the velocity; s = ρv is the momentum density; p is the thermal pressure; Q is the artificial viscous stress tensor (Von Neumann & Richtmyer 1950); Φ is the gravitational potential; $\mathrm{T}\mathrm{=}\mathit{\rho \nu }\mathrm{(}\mathrm{\nabla }{v}\mathrm{+}\mathrm{\nabla }{{v}}^{\mathit{T}}\mathrm{-}\frac{\mathrm{2}}{\mathrm{3}}\mathrm{I\nabla }\mathrm{·}{v}\mathrm{)}$ is the viscous stress tensor; ν is the kinematic viscosity; I is the identity matrix; e is the internal energy density; E_R is the radiation energy density; ${\mathit{E}}_{\mathrm{T}}\mathrm{=}\mathit{e}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}\mathit{\rho }{\mathit{v}}^{\mathrm{2}}\mathrm{+}{\mathit{E}}_{\mathrm{R}}$ is the total energy density; T is gas temperature; F_R is the radiation flux; P is the radiation pressure tensor; and κ_P, κ_E, and κ_F are the Planck, radiation energy, and radiation flux mean opacities. The radiation constant and the speed of light are a_R and c, respectively. Although different values for κ_P and κ_E are permitted in the algorithm, all the tests presented here assume that the radiation behaves like a black body, and therefore κ_P = κ_E.

The first closure relation applied to the above equation set is the ideal gas law: p = (γ − 1)e = ρC_VT, where γ is the ratio of specific heats, C_V = k/ (γ − 1)μm_H is the specific heat capacity at constant volume, k is Boltzmann’s constant, μ is the mean molecular weight, and m_H is the mass of hydrogen.

Magnetohydrodynamics is, of course, also available in AZEuS, but we restrict the discussion here to RHD for simplicity. AZEuS also includes the option to solve either the total energy equation or the internal energy equation. It should be noted, however, that the radiation solver uses the internal energy density. Furthermore, the total energy will generally not be conserved in the co-moving frame. For details on the MHD and AMR algorithms in AZEuS, including test problems, we refer the reader to Clarke (1996, 2010); Ramsey et al. (2012), and references therein.

2.1. Flux-limited diffusion

In order to close equation set (1)–(5), we adopt the FLD approximation, where the radiative flux is replaced by Fick’s law: ${{F}}_{\mathrm{R}}\mathrm{=}\mathrm{-}\mathit{\lambda }\mathrm{\left(}\mathit{R}\mathrm{\right)}\frac{\mathit{c}}{\mathit{\rho }{\mathit{\kappa }}_{\mathrm{R}}}\mathrm{\nabla }{\mathit{E}}_{\mathrm{R}}\mathrm{=}\mathrm{-}\mathit{D}\mathrm{\nabla }{\mathit{E}}_{\mathrm{R}}\mathrm{;}$(6)λ(R) is the so-called flux limiter, κ_R is the Rosseland mean opacity, and D is the diffusion coefficient. We assume that κ_R = κ_F, a common approximation in the diffusion regime, which preserves the correct energy and momentum transport.

By omitting a conservation equation for the radiative flux, and under the Eddington approximation (λ = 1/3), it is possible with Eq. (6) to generate an unphysically high flux (>cE_R), leading to an energy propagation speed that exceeds the speed of light. Flux limiters, by design, correct for this. One commonly used form for the flux limiter is given by (Levermore & Pomraning 1981) $\mathit{\lambda }\mathrm{\left(}\mathit{R}\mathrm{\right)}\mathrm{=}\frac{\mathrm{2}\mathrm{+}\mathit{R}}{\mathrm{6}\mathrm{+}\mathrm{3}\mathit{R}\mathrm{+}{\mathit{R}}^{\mathrm{2}}}\mathit{,}$(7)where $\mathit{R}\mathrm{=}\frac{\mathrm{\left|}\mathrm{\nabla }{\mathit{E}}_{\mathrm{R}}\mathrm{\right|}}{\mathit{\rho }{\mathit{\kappa }}_{\mathrm{R}}\mathrm{⟨}{\mathit{E}}_{\mathrm{R}}\mathrm{⟩}}\mathit{,}$(8)and ⟨⟩ denotes a spatially-averaged quantity. In Eq. (7), as R becomes large, λ → 1 /R, and the flux is limited to ≤ cE_R. Conversely, for small R, the flux limiter reduces to the Eddington approximation (λ = 1/3). Flux-limiters currently implemented in AZEuS include Levermore & Pomraning (1981), Minerbo (1978) and Kley (1989).

We note that, in Eqs. (6) and (8), the radiative flux F_R, the diffusion coefficients D, and the quantity R are all face-centred quantities while the product ρκ_R = σ_R is intrinsically zone-centred. Although, in principle, we could arithmetically average σ_R to the face, we instead adopt the surface formula of Howell & Greenough (2003) to improve results at sharp interfaces in optical depth. For example, in the 1-direction: ${\mathit{\sigma }}_{\mathrm{R}\mathit{,}\hspace{0.17em}\mathit{i}\mathrm{+}\mathrm{1}\mathit{/}\mathrm{2}}\mathrm{=}\begin{array}{c}\min \end{array}\left[\frac{{\mathit{\sigma }}_{\mathrm{R}\mathit{,}\hspace{0.17em}\mathit{i}}\mathrm{+}{\mathit{\sigma }}_{\mathrm{R}\mathit{,}\hspace{0.17em}\mathit{i}\mathrm{+}\mathrm{1}}}{\mathrm{2}}\mathit{,}\begin{array}{c}\max \end{array}\left(\frac{\mathrm{2}{\mathit{\sigma }}_{\mathrm{R}\mathit{,}\hspace{0.17em}\mathit{i}}\hspace{0.17em}{\mathit{\sigma }}_{\mathrm{R}\mathit{,}\hspace{0.17em}\mathit{i}\mathrm{+}\mathrm{1}}}{{\mathit{\sigma }}_{\mathrm{R}\mathit{,}\hspace{0.17em}\mathit{i}}\mathrm{+}{\mathit{\sigma }}_{\mathrm{R}\mathit{,}\hspace{0.17em}\mathit{i}\mathrm{+}\mathrm{1}}}\mathit{,}\frac{\mathrm{4}}{\mathrm{3}\mathrm{\Delta }{\mathit{x}}_{\mathrm{1}}}\right)\right]\mathrm{·}$(9)The radiation pressure tensor is assumed to have the form: $\mathrm{P}\mathrm{=}\mathrm{f}{\mathit{E}}_{\mathrm{R}}\mathit{,}$(10)where f is called the Eddington tensor, and is given by (Turner & Stone 2001; Hayes et al. 2006) $\begin{array}{ccc}& & \mathrm{f}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}\mathrm{(}\mathrm{1}\mathrm{-}\mathit{f}\mathrm{)}\mathrm{I}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}\mathrm{(}\mathrm{3}\mathit{f}\mathrm{-}\mathrm{1}\mathrm{)}\hat{{n}}\hat{{n}}\mathrm{;}\\ & & \mathit{f}\mathrm{=}\mathit{\lambda }\mathrm{+}{\mathit{\lambda }}^{\mathrm{2}}{\mathit{R}}^{\mathrm{2}}\mathit{,}\end{array}$and $\hat{{n}}\mathrm{=}\mathrm{\nabla }{\mathit{E}}_{\mathrm{R}}\mathit{/}\mathrm{\left|}\mathrm{\nabla }{\mathit{E}}_{\mathrm{R}}\mathrm{\right|}$.

In AZEuS, the coupled equations of radiation hydrodynamics are solved through the use of operator splitting, where the hydrodynamical terms are calculated explicitly, but the radiation terms are solved implicitly. The coupled gas-radiation system of equations that we wish to solve is given by $\begin{array}{ccc}& & \frac{\mathit{\partial e}}{\mathit{\partial t}}\mathrm{=}\mathit{c\rho }\mathrm{(}{\mathit{\kappa }}_{\mathrm{E}}{\mathit{E}}_{\mathrm{R}}\mathrm{-}{\mathit{\kappa }}_{\mathrm{P}}{\mathit{a}}_{\mathrm{R}}{\mathit{T}}^{\mathrm{4}}\mathrm{)}\mathrm{;}\\ & & \frac{\mathit{\partial }{\mathit{E}}_{\mathrm{R}}}{\mathit{\partial t}}\mathrm{=}\mathrm{-}\mathit{c\rho }\mathrm{(}{\mathit{\kappa }}_{\mathrm{E}}{\mathit{E}}_{\mathrm{R}}\mathrm{-}{\mathit{\kappa }}_{\mathrm{P}}{\mathit{a}}_{\mathrm{R}}{\mathit{T}}^{\mathrm{4}}\mathrm{)}\mathrm{+}\mathrm{\nabla }\mathrm{·}\left(\mathit{D}\mathrm{\nabla }{\mathit{E}}_{\mathrm{R}}\right)\mathrm{-}\mathrm{P}\mathrm{:}\mathrm{\nabla }{v}\mathit{,}\end{array}$and the advection, gravitational, compressive heating, and physical and artificial viscous terms which appear in Eqs. (1)–(5) are solved in the explicit part of the AZEuS algorithm.

Using the ideal gas law, we can substitute the temperature for the internal energy density to find $\frac{\mathit{\partial T}}{\mathit{\partial t}}\mathrm{=}\frac{\mathit{c}}{{\mathit{C}}_{\mathit{V}}}\mathrm{(}{\mathit{\kappa }}_{\mathrm{E}}{\mathit{E}}_{\mathrm{R}}\mathrm{-}{\mathit{\kappa }}_{\mathrm{P}}{\mathit{a}}_{\mathrm{R}}{\mathit{T}}^{\mathrm{4}}\mathrm{)}\mathit{.}$(15)We now write Eqs. (14) and (15) in difference form, $\begin{array}{ccc}& & \frac{{\mathit{T}}^{\mathit{n}\mathrm{+}\mathrm{1}}\mathrm{-}{\mathit{T}}^{\mathit{n}}}{\mathrm{\Delta }\mathit{t}}\mathrm{=}\frac{\mathit{c}}{{\mathit{C}}_{\mathit{V}}}\left[{\mathit{\kappa }}_{\mathrm{E}}^{\mathit{n}\mathrm{+}\mathrm{1}}{\mathit{E}}_{\mathrm{R}}^{\mathit{n}\mathrm{+}\mathrm{1}}\mathrm{-}{\mathit{\kappa }}_{\mathrm{P}}^{\mathit{n}\mathrm{+}\mathrm{1}}{\mathit{a}}_{\mathrm{R}}{\mathrm{(}{\mathit{T}}^{\mathit{n}\mathrm{+}\mathrm{1}}\mathrm{)}}^{\mathrm{4}}\right]\mathrm{;}\\ \hspace{0.17em}\begin{array}{c}\begin{array}{c}\end{array}\begin{array}{c}\end{array}\frac{{\mathit{E}}_{\mathrm{R}}^{\mathit{n}\mathrm{+}\mathrm{1}}\mathrm{-}{\mathit{E}}_{\mathrm{R}}^{\mathit{n}}}{\mathrm{\Delta }\mathit{t}}\mathrm{=}\mathrm{-}\mathit{c}{\mathit{\rho }}^{\mathit{n}}\left[{\mathit{\kappa }}_{\mathrm{E}}^{\mathit{n}\mathrm{+}\mathrm{1}}{\mathit{E}}_{\mathrm{R}}^{\mathit{n}\mathrm{+}\mathrm{1}}\mathrm{-}{\mathit{\kappa }}_{\mathrm{P}}^{\mathit{n}\mathrm{+}\mathrm{1}}{\mathit{a}}_{\mathrm{R}}{\mathrm{(}{\mathit{T}}^{\mathit{n}\mathrm{+}\mathrm{1}}\mathrm{)}}^{\mathrm{4}}\right]\\ \begin{array}{c}\end{array}\begin{array}{c}\end{array} \mathrm{+}\mathrm{\nabla }\mathrm{·}\mathrm{(}{\mathit{D}}^{\mathit{n}}\mathrm{\nabla }{{\mathit{E}}_{\mathrm{R}}^{\mathit{n}\mathrm{+}\mathrm{1}}}^{\mathrm{)}}\mathrm{-}\mathrm{\left(}{\mathrm{f}}^{\hspace{0.17em}\mathit{n}}{\mathit{E}}_{\mathrm{R}}^{\mathit{n}\mathrm{+}\mathrm{1}}\mathrm{\right)}\mathrm{:\nabla }{{v}}^{\mathit{n}}\mathrm{;}\end{array}& & \\ \hspace{0.17em}\begin{array}{c}\begin{array}{c}\end{array}\begin{array}{c}\end{array}\mathrm{\nabla }\mathrm{·}\mathrm{(}{\mathit{D}}^{\mathit{n}}\mathrm{\nabla }{{\mathit{E}}_{\mathrm{R}}^{\mathit{n}\mathrm{+}\mathrm{1}}}^{\mathrm{)}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{\Delta }{\mathit{x}}_{\mathrm{1}}}\left[{\mathit{D}}_{\mathit{i}\mathrm{+}\mathrm{1}}^{\mathit{n}}\frac{{\mathit{E}}_{\mathrm{R}\mathit{,}\hspace{0.17em}\mathit{i}\mathrm{+}\mathrm{1}}^{\mathit{n}\mathrm{+}\mathrm{1}}\mathrm{-}{\mathit{E}}_{\mathrm{R}\mathit{,}\hspace{0.17em}\mathit{i}}^{\mathit{n}\mathrm{+}\mathrm{1}}}{\mathrm{\Delta }{\mathit{x}}_{\mathrm{1}}}\mathrm{-}{\mathit{D}}_{\mathit{i}}^{\mathit{n}}\frac{{\mathit{E}}_{\mathrm{R}\mathit{,}\hspace{0.17em}\mathit{i}}^{\mathit{n}\mathrm{+}\mathrm{1}}\mathrm{-}{\mathit{E}}_{\mathrm{R}\mathit{,}\hspace{0.17em}\mathit{i}\mathrm{-}\mathrm{1}}^{\mathit{n}\mathrm{+}\mathrm{1}}}{\mathrm{\Delta }{\mathit{x}}_{\mathrm{1}}}\right]\\ \begin{array}{c}\end{array}\begin{array}{c}\end{array} \mathrm{+}\frac{\mathrm{1}}{\mathrm{\Delta }{\mathit{x}}_{\mathrm{2}}}\left[{\mathit{D}}_{\mathit{j}\mathrm{+}\mathrm{1}}^{\mathit{n}}\frac{{\mathit{E}}_{\mathrm{R}\mathit{,}\hspace{0.17em}\mathit{j}\mathrm{+}\mathrm{1}}^{\mathit{n}\mathrm{+}\mathrm{1}}\mathrm{-}{\mathit{E}}_{\mathrm{R}\mathit{,}\hspace{0.17em}\mathit{j}}^{\mathit{n}\mathrm{+}\mathrm{1}}}{\mathrm{\Delta }{\mathit{x}}_{\mathrm{2}}}\mathrm{-}{\mathit{D}}_{\mathit{j}}^{\mathit{n}}\frac{{\mathit{E}}_{\mathrm{R}\mathit{,}\hspace{0.17em}\mathit{j}}^{\mathit{n}\mathrm{+}\mathrm{1}}\mathrm{-}{\mathit{E}}_{\mathrm{R}\mathit{,}\hspace{0.17em}\mathit{j}\mathrm{-}\mathrm{1}}^{\mathit{n}\mathrm{+}\mathrm{1}}}{\mathrm{\Delta }{\mathit{x}}_{\mathrm{2}}}\right]\\ \begin{array}{c}\end{array}\begin{array}{c}\end{array} \mathrm{+}\frac{\mathrm{1}}{\mathrm{\Delta }{\mathit{x}}_{\mathrm{3}}}\left[{\mathit{D}}_{\mathit{k}\mathrm{+}\mathrm{1}}^{\mathit{n}}\frac{{\mathit{E}}_{\mathrm{R}\mathit{,}\hspace{0.17em}\mathit{k}\mathrm{+}\mathrm{1}}^{\mathit{n}\mathrm{+}\mathrm{1}}\mathrm{-}{\mathit{E}}_{\mathrm{R}\mathit{,}\hspace{0.17em}\mathit{k}}^{\mathit{n}\mathrm{+}\mathrm{1}}}{\mathrm{\Delta }{\mathit{x}}_{\mathrm{3}}}\mathrm{-}{\mathit{D}}_{\mathit{k}}^{\mathit{n}}\frac{{\mathit{E}}_{\mathrm{R}\mathit{,}\hspace{0.17em}\mathit{k}}^{\mathit{n}\mathrm{+}\mathrm{1}}\mathrm{-}{\mathit{E}}_{\mathrm{R}\mathit{,}\hspace{0.17em}\mathit{k}\mathrm{-}\mathrm{1}}^{\mathit{n}\mathrm{+}\mathrm{1}}}{\mathrm{\Delta }{\mathit{x}}_{\mathrm{3}}}\right]\end{array}\mathit{,}& & \end{array}$where n and n + 1 denote original and updated values within the implicit solve. We have assumed Cartesian coordinates here for brevity (x₁,x₂,x₃ = x,y,z); the curvilinear factors associated with cylindrical and spherical coordinates are presented in Appendix A of Ramsey et al. (2012). We also note that the diffusion coefficients (and therefore the flux limiter) and Eddington tensor f are time-lagged, which is a voluntary choice for stability and efficiency over accuracy.

As given, Eqs. (16)–(18) form a non-linear system of equations. If we now linearise the term in (T^{n + 1})⁴ by taking a Taylor expansion about Tⁿ (Commerçon et al. 2011), $\mathrm{(}{\mathit{T}}^{\mathit{n}\mathrm{+}\mathrm{1}}{\mathrm{)}}^{\mathrm{4}}\mathrm{\approx }\mathrm{4}\mathrm{(}{\mathit{T}}^{\mathit{n}}{\mathrm{)}}^{\mathrm{3}}{\mathit{T}}^{\mathit{n}\mathrm{+}\mathrm{1}}\mathrm{-}\mathrm{3}\mathrm{(}{\mathit{T}}^{\mathit{n}}{\mathrm{)}}^{\mathrm{4}}\mathit{,}$(19)then (T^{n + 1})⁴ can be calculated in terms of known quantities only, ${\mathit{T}}^{\mathit{n}\mathrm{+}\mathrm{1}}\mathrm{=}\frac{{\mathit{T}}^{\mathit{n}}\mathrm{+}\mathrm{3}{\mathit{a}}_{\mathrm{R}}\mathit{c}\frac{\mathrm{\Delta }\mathit{t}}{{\mathit{C}}_{\mathit{V}}}{\mathit{\kappa }}_{\mathrm{P}}^{\mathit{n}}\mathrm{(}{\mathit{T}}^{\mathit{n}}{\mathrm{)}}^{\mathrm{4}}\mathrm{+}\frac{\mathrm{\Delta }\mathit{t}}{{\mathit{C}}_{\mathit{V}}}\mathit{c}{\mathit{\kappa }}_{\mathrm{E}}^{\mathit{n}}{\mathit{E}}_{\mathrm{R}}^{\mathit{n}\mathrm{+}\mathrm{1}}}{\mathrm{1}\mathrm{+}\mathrm{4}{\mathit{a}}_{\mathrm{R}}\mathit{c}\frac{\mathrm{\Delta }\mathit{t}}{{\mathit{C}}_{\mathit{V}}}{\mathit{\kappa }}_{\mathrm{P}}^{\mathit{n}}\mathrm{(}{\mathit{T}}^{\mathit{n}}{\mathrm{)}}^{\mathrm{3}}}\mathit{,}$(20)thus eliminating the temperature equation from the radiation system of equations². We note that substituting Eqs. (19) and (20) into Eq. (17) would result in a linear system of equations (e.g., Commerçon et al. 2011) if it were not for the fact that we account for the generally non-linear temperature-dependence of the mean opacities, ${\mathit{\kappa }}_{\mathrm{P}}^{\mathit{n}\mathrm{+}\mathrm{1}}$ and ${\mathit{\kappa }}_{\mathrm{E}}^{\mathit{n}\mathrm{+}\mathrm{1}}$.

The right-hand sides of Eqs. (2) and (5) contain terms describing the effect of radiation pressure on the momentum and total energy densities, respectively, but are not included in the implicit radiation solve. They are instead treated as explicit source terms and, using Eq. (6), take the following forms: $\begin{array}{ccc}& & \frac{\mathit{\rho }{\mathit{\kappa }}_{\mathrm{F}}}{\mathit{c}}{{F}}_{\mathrm{R}}\mathrm{=}\mathrm{-}\mathit{\lambda }\mathrm{\nabla }{\mathit{E}}_{\mathrm{R}}\mathit{,}\\ & & \frac{\mathit{\rho }{\mathit{\kappa }}_{\mathrm{F}}}{\mathit{c}}{v}\mathrm{·}{{F}}_{\mathrm{R}}\mathrm{=}\mathrm{-}\mathit{\lambda }{v}\mathrm{·}\mathrm{\nabla }{\mathit{E}}_{\mathrm{R}}\mathit{.}\end{array}$In practice, since E_R is a zone-centred quantity, Eq. (21) is naturally face-centred and can therefore be applied to the momentum as is. Equation (22) is also naturally face-centred, although it is applied to the zone-centred total energy. As such, we adopt the following form (in Cartesian coordinates; cf. Zhang et al. 2011), $\begin{array}{ccc}& & \mathit{\lambda }{v}\mathrm{·}\mathrm{\nabla }{\mathit{E}}_{\mathrm{R}}\mathrm{=}{\mathit{\lambda }}_{\mathit{i,}\hspace{0.17em}\mathit{j,k}}(\hspace{0.17em}\frac{{\mathit{v}}_{\mathrm{1}\mathit{,}\hspace{0.17em}\mathit{i}\hspace{0.17em}\mathrm{+}\hspace{0.17em}\mathrm{1}}\mathrm{+}{\mathit{v}}_{\mathrm{1}\mathit{,}\hspace{0.17em}\mathit{i}}}{\mathrm{2}}\frac{\begin{array}{c}\mathrm{˜}\\ {\mathit{E}}_{\mathrm{R}\mathit{,}\hspace{0.17em}\mathit{i}\hspace{0.17em}\mathrm{+}\hspace{0.17em}\mathrm{1}}\end{array}\mathrm{-}\begin{array}{c}\mathrm{˜}\\ {\mathit{E}}_{\mathrm{R}\mathit{,}\hspace{0.17em}\mathit{i}}\end{array}}{\mathrm{\Delta }{\mathit{x}}_{\mathrm{1}}}\\ & & \mathrm{+}\frac{{\mathit{v}}_{\mathrm{2}\mathit{,}\hspace{0.17em}\mathit{j}\hspace{0.17em}\mathrm{+}\hspace{0.17em}\mathrm{1}}\mathrm{-}{\mathit{v}}_{\mathrm{2}\mathit{,}\hspace{0.17em}\mathit{j}}}{\mathrm{2}}{\frac{\begin{array}{c}\mathrm{˜}\\ {\mathit{E}}_{\mathrm{R}\mathit{,}\hspace{0.17em}\mathit{j}\hspace{0.17em}\mathrm{+}\hspace{0.17em}\mathrm{1}}\end{array}\mathrm{-}\begin{array}{c}\mathrm{˜}\\ {\mathit{E}}_{\mathrm{R}\mathit{,}\hspace{0.17em}\mathit{j}}\end{array}}{\mathrm{\Delta }{\mathit{x}}_{\mathrm{2}}}}^{}\\ & & \mathrm{+}\frac{{\mathit{v}}_{\mathrm{3}\mathit{,}\hspace{0.17em}\mathit{k}\hspace{0.17em}\mathrm{+}\hspace{0.17em}\mathrm{1}}\mathrm{+}{\mathit{v}}_{\mathrm{3}\mathit{,}\hspace{0.17em}\mathit{k}}}{\mathrm{2}}\frac{\begin{array}{c}\mathrm{˜}\\ {\mathit{E}}_{\mathrm{R}\mathit{,}\hspace{0.17em}\mathit{k}\hspace{0.17em}\mathrm{+}\hspace{0.17em}\mathrm{1}}\end{array}\mathrm{-}\begin{array}{c}\mathrm{˜}\\ {\mathit{E}}_{\mathrm{R}\mathit{,}\hspace{0.17em}\mathit{k}}\end{array}}{\mathrm{\Delta }{\mathit{x}}_{\mathrm{3}}})\mathit{,}\end{array}$(23)where $\begin{array}{c}{}_{\mathrm{˜}}\\ {\mathit{E}}_{\mathrm{R}}\end{array}{}_{\mathit{i}}$ denotes an upwinded interpolation of the radiation energy density in direction i (see, e.g., Clarke 1996).

The effects of radiation must also be included in the calculation of the time step. Following Krumholz et al. (2007), $\mathrm{\Delta }{\mathit{t}}_{\mathrm{RD}}^{-1}\mathrm{=}\frac{\mathrm{1}}{\mathrm{\Delta }{\mathit{x}}_{\min }}\sqrt{\frac{\mathrm{4}{\mathit{E}}_{\mathrm{R}}}{\mathrm{9}\mathit{\rho }}\left(\mathrm{1}\mathrm{-}{\mathrm{e}}^{\mathrm{-}{\mathit{\kappa }}_{\mathrm{R}}\mathit{\rho }\mathrm{\Delta }{\mathit{x}}_{\min }}\right)}$(24)is added to the existing time step calculation in quadrature to account for the (diffusion) radiation pressure. Finally, the time step can be postdictively limited based on the maximum desired fractional change of E_R within a single step. This is controlled via the user-defined tolerance dttoler (Table 1).

2.2. Irradiation

The implementation of stellar irradiation in AZEuS follows Kuiper et al. (2010); Bitsch et al. (2013), where the heating due to irradiation is added to the gas energy equation (Eq. (13)) as a source term: $\frac{\mathit{\partial e}}{\mathit{\partial t}}\mathrm{=}\mathit{c\rho }\mathrm{(}{\mathit{\kappa }}_{\mathrm{E}}{\mathit{E}}_{\mathrm{R}}\mathrm{-}{\mathit{\kappa }}_{\mathrm{P}}{\mathit{a}}_{\mathrm{R}}{\mathit{T}}^{\mathrm{4}}\mathrm{)}\mathrm{+}{\mathit{Q}}_{\mathrm{irr}}^{\mathrm{+}}\mathit{.}$(25)In regards to solving the radiation system of equations, this manifests as an additional term in Eq. (20): ${\mathit{T}}^{\mathit{n}\mathrm{+}\mathrm{1}}\mathrm{=}\frac{{\mathit{T}}^{\mathit{n}}\mathrm{+}\mathrm{3}{\mathit{a}}_{\mathrm{R}}\mathit{c}\frac{\mathrm{\Delta }\mathit{t}}{{\mathit{C}}_{\mathit{V}}}{\mathit{\kappa }}_{\mathrm{P}}^{\mathit{n}}\mathrm{(}{\mathit{T}}^{\mathit{n}}{\mathrm{)}}^{\mathrm{4}}\mathrm{+}\frac{\mathrm{\Delta }\mathit{t}}{{\mathit{C}}_{\mathit{V}}}\mathit{c}{\mathit{\kappa }}_{\mathrm{E}}^{\mathit{n}}{\mathit{E}}_{\mathrm{R}}^{\mathit{n}\mathrm{+}\mathrm{1}}\mathrm{+}\frac{\mathrm{\Delta }\mathit{t}}{\mathit{\rho }{\mathit{C}}_{\mathit{V}}}{\mathit{Q}}_{\mathrm{irr}}^{\mathrm{+}}}{\mathrm{1}\mathrm{+}\mathrm{4}{\mathit{a}}_{\mathrm{R}}\mathit{c}\frac{\mathrm{\Delta }\mathit{t}}{{\mathit{C}}_{\mathit{V}}}{\mathit{\kappa }}_{\mathrm{P}}^{\mathit{n}}\mathrm{(}{\mathit{T}}^{\mathit{n}}{\mathrm{)}}^{\mathrm{3}}}\mathrm{·}$(26)To include frequency-dependent transport of radiation in the irradiation algorithm, one simply sums over all frequency bins, $\begin{array}{ccc}{\mathit{Q}}_{\mathrm{irr}}^{\mathrm{+}}\mathrm{=}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{\infty }}{\mathit{Q}}_{\mathrm{irr}\mathit{,\nu }}^{\mathrm{+}}\hspace{0.17em}\mathrm{d}\mathit{\nu }\mathrm{\approx }\sum _{{\mathit{n}}^{\mathrm{\prime }}\mathrm{=}\mathrm{1}}^{\mathit{n}}{\mathit{Q}}_{\mathrm{irr}\mathit{,\nu }}^{\mathrm{+}}\mathrm{\Delta }\mathit{\nu }\mathrm{\left(}{\mathit{n}}^{\mathrm{\prime }}\mathrm{\right)}\mathit{,}& & \end{array}$(27)where n is the number of frequency bins, and Δν(n′) is the width of an individual frequency bin. In the context of accretion disks and spherical coordinates, if the source is assumed to be at r = 0, the heating rate is thus given by ${\mathit{Q}}_{\mathrm{irr}\mathit{,\nu }}^{\mathrm{+}}\mathrm{=}\{\begin{array}{c}\\ & \mathrm{if}\mathrm{\Delta }{\mathit{\tau }}_{\mathit{\nu }}\mathrm{\le }{\mathit{\tau }}_{\mathrm{thresh}}\mathrm{;}\\ & \mathrm{otherwise}\mathit{,}\end{array}$(28)where κ_{∗ ,ν} is the frequency-dependent opacity to stellar radiation, Δτ_ν is the change in optical depth across Δr, $\mathrm{\Delta }{\mathit{V}}_{\mathit{r}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{3}}\mathrm{(}{\mathit{r}}_{\mathit{i}\mathrm{+}\mathrm{1}\mathit{/}\mathrm{2}}^{\mathrm{3}}\mathrm{-}{\mathit{r}}_{\mathit{i}\mathrm{-}\mathrm{1}\mathit{/}\mathrm{2}}^{\mathrm{3}}\mathrm{)}$ is the radial volume element³, F_irr,ν is the attenuated stellar flux ${\mathit{F}}_{\mathrm{irr}\mathit{,\nu }}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}{\mathit{F}}_{\mathrm{irr}\mathit{,\nu }}\mathrm{\left(}{\mathit{R}}_{\mathrm{\ast }}\mathrm{\right)}\frac{{\mathit{R}}_{\mathrm{\ast }}^{\mathrm{2}}}{{\mathit{r}}^{\mathrm{2}}}{\mathit{e}}^{\mathrm{-}{\mathit{\tau }}_{\mathit{\nu }}\mathrm{\left(}\mathit{r}\mathrm{\right)}}\mathit{,}$(29)R_∗ is the stellar radius, ${\mathit{\tau }}_{\mathit{\nu }}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}{{}^{\mathrm{\int }}\mathit{R}_{\mathrm{\ast }}}_{\hspace{0.17em}\mathit{r}}^{}{\mathit{\kappa }}_{\mathrm{\ast }\mathit{,\nu }}\hspace{0.17em}\mathit{\rho }\hspace{0.17em}\mathrm{d}\mathit{r}$ is the radial optical depth, and F_irr,ν(R_∗) is the flux at the stellar surface. Although the two cases for the stellar heating rate in Eq. (28) are mathematically equivalent, when calculating differences on a discrete grid, the differential form is subject to numerical round-off errors at very low optical depths (Bruls et al. 1999), and thus we use the integral form when Δτ_ν ≤ τ_thresh, where typically τ_thresh = 10^-4.

Treating the stellar flux in this manner removes any explicit time-dependence, and thus it is valid only if the light travel time is shorter than the timescale for changes in the optical depth τ_ν(r). Conversely, this also means the calculation of the stellar flux at any location depends only on the instantaneous optical depth along a ray (and thus density and opacity), and can be decoupled from the hydrodynamics and diffuse radiation inside a single time step. As such, the irradiation algorithm is invoked before solving the radiation matrix, and its results are taken as constant during one combined hydrodynamic-radiation time step. To obtain a grey version of the irradiation algorithm, the frequency-dependence in all quantities is dropped, and κ_{∗ ,ν} is replaced by the Planck mean opacity at the stellar effective temperature, κ_P(T_∗).

The radiative force density due to stellar irradiation is also considered, and is included as a source term in the momentum equation (Eq. (2)) akin to the radiation pressure term from FLD. As we consider only single-group (grey) hydrodynamics, the irradiative flux is integrated over frequency before the force density is calculated: ${\mathit{a}}_{\mathrm{irr}}\mathrm{=}\{\begin{array}{c}\\ & \mathrm{if}\mathrm{\Delta }\mathit{\tau }\mathrm{\le }{\mathit{\tau }}_{\mathrm{thresh}}\mathrm{;}\\ & \mathrm{otherwise}\mathit{.}\end{array}$(30)The above discussion on irradiation assumes spherical coordinates and a disk-like setting. However, the algorithm can be generalised to any curvilinear coordinate system and problem where the source irradiation propagates along coordinate axes. Astrophysical examples of this could include the collected radiation from massive stars at a large distance from a region of interest, e.g., a protoplanetary disk or a photodissociation region. Practically, the generalisation amounts to accounting for the appropriate geometric factors in Eqs. (28)–(30), and an example for Cartesian coordinates is presented in Sect. 4.7.

2.3. Solving the radiation matrix

We solve Eq. (17) for ${\mathit{E}}_{\mathrm{R}}^{\mathit{n}\mathrm{+}\mathrm{1}}$ via a globally convergent Newton-Raphson method (Press et al. 1992). To begin, we rewrite Eq. (17) in the form $\begin{array}{ccc}& & {\mathit{g}}_{\mathit{i,j,k}}\mathrm{=}{\mathit{E}}_{\mathrm{R}\mathit{,i,j,k}}^{\mathit{n}\mathrm{+}\mathrm{1}}\mathrm{-}{\mathit{E}}_{\mathrm{R}\mathit{,i,j,k}}^{\mathit{n}}\\ & & \mathrm{+}\mathrm{\Delta }\mathit{tc}{\mathit{\rho }}^{\mathit{n}}\mathrm{[}{\mathit{\kappa }}_{\mathrm{E}\mathit{,i,j,k}}^{\mathit{n}\mathrm{+}\mathrm{1}}{\mathit{E}}_{\mathrm{R}\mathit{,i,j,k}}^{\mathit{n}\mathrm{+}\mathrm{1}}\mathrm{-}{\mathit{\kappa }}_{\mathrm{P}\mathit{,i,j,k}}^{\mathit{n}\mathrm{+}\mathrm{1}}{\mathit{a}}_{\mathrm{R}}\mathrm{(}{\mathit{T}}_{\mathit{i,j,k}}^{\mathit{n}\mathrm{+}\mathrm{1}}{\mathrm{)}}^{\mathrm{4}}\mathrm{]}\\ & & \mathrm{-}\mathrm{\Delta }\mathit{t}\mathrm{\nabla }\mathrm{·}\mathrm{(}{\mathit{D}}^{\mathit{n}}\mathrm{\nabla }{{\mathit{E}}_{\mathrm{R}}^{\mathit{n}\mathrm{+}\mathrm{1}}}^{\mathrm{)}}\mathrm{-}\mathrm{\Delta }\mathit{t}\mathrm{\left(}{\mathrm{f}}^{\hspace{0.17em}\mathit{n}}{\mathit{E}}_{\mathrm{R}\mathit{,i,j,k}}^{\mathit{n}\mathrm{+}\mathrm{1}}\mathrm{\right)}\mathrm{:\nabla }{{v}}^{\mathit{n}}\mathit{,}\end{array}$(31)where ${\mathit{T}}_{\mathit{i,j,k}}^{\mathit{n}\mathrm{+}\mathrm{1}}$ is given by Eqs. (19) and (20) (or (26)if irradiation is included), fⁿ is given by Eq. (11), the radiative diffusion term is described by Eq. (18), and the mean opacities are assumed to be functions of T^{n + 1}. This equation forms a sparse linear system of N equations, where N is the total number of zones on the grid, and is of the form $\begin{array}{ccc}\mathrm{𝒥}\mathrm{\left(}{\mathit{E}}_{\mathrm{R}}\mathrm{\right)}\hspace{0.17em}\mathit{\delta }{\mathit{E}}_{\mathrm{R}}\mathrm{=}\mathrm{-}\mathit{g}\mathrm{\left(}{\mathit{E}}_{\mathrm{R}}\mathrm{\right)}\mathit{,}& & \end{array}$(32)where $\begin{array}{ccc}\mathrm{𝒥}\mathrm{\left(}{\mathit{E}}_{\mathrm{R}}\mathrm{\right)}\mathrm{=}\sum _{\mathit{n}\mathrm{=}\mathrm{1}}^{\mathit{N}}\frac{\mathit{\partial }{\mathit{g}}_{\mathit{i,j,k}}}{\mathit{\partial }{\mathit{E}}_{\mathrm{R}\mathit{,n}}}& & \end{array}$(33)is the Jacobian, and δE_R is the change in the radiation energy density. The desired solution of this system is g_i,j,k = 0.

Before attempting to solve Eq. (32), we first multiply each row of the matrix by the volume of its corresponding zone, ΔV = ΔV₁ΔV₂ΔV₃. This symmetrises the matrix (Hayes et al. 2006), and thus we only need to calculate and store the diagonal and super-diagonal (or sub-diagonal) elements of the Jacobian.

Although AZEuS includes tri-diagonal and successive over-relaxation (SOR) methods to solve Eq. (32), we instead primarily rely on the PETSc library (Balay et al. 2014) because of its flexibility and performance. Through experimentation, we find that the most generally reliable and efficient solver for our purposes is the iterative bi-conjugate gradient stabilised (BiCGSTAB) method plus a Jacobi (diagonal) preconditioner⁴, even though our radiation matrix is symmetric and permits the use of simpler solvers (such as the conjugate gradient method).

Convergence in the Newton-Raphson method is signalled when either $\mathrm{\left(}\max \mathrm{\right|}\mathit{\delta }{\mathit{E}}_{\mathrm{R}}^{\mathrm{\left(}\mathit{q}\mathrm{\right)}}\mathit{/}{\mathit{E}}_{\mathrm{R}}^{\mathit{n}\mathrm{+}\mathrm{1}\mathit{,}\mathrm{\left(}\mathit{q}\mathrm{\right)}}\mathrm{\left|}\mathrm{\right)}$ or $\mathrm{\left(}\max \mathrm{\right|}{\mathit{g}}_{\mathit{i,j,k}}^{\mathrm{\left(}\mathit{q}\mathrm{\right)}}\mathrm{\left|}\mathrm{\right)}$ falls below user-defined tolerances nrtoler and nrtolrhs, respectively, and where (q) denotes the number of Newton-Raphson steps taken. Furthermore, when using an iterative matrix solver, an additional user-defined tolerance is applied (slstol) to signal convergence with respect to the L₂-norm of the preconditioned residual (${\mathit{r}}^{\mathrm{\left(}\mathit{p}\mathrm{\right)}}\mathrm{=}\mathrm{-}\mathit{g}\mathrm{-}\mathrm{𝒥}\mathit{\delta }{\mathit{E}}_{\mathrm{R}}^{\mathrm{\left(}\mathit{p}\mathrm{\right)}}$, where (p) is the number of matrix solver iterations used).

At the end of the implicit solve, we calculate the change in the internal energy density (δe^{n + 1}) via the ideal gas law and, if solving the total energy equation, subsequently update the total energy density using $\mathit{\delta }{\mathit{E}}_{\mathrm{T}}^{\mathit{n}\mathrm{+}\mathrm{1}}\mathrm{=}\mathit{\delta }{\mathit{e}}^{\mathit{n}\mathrm{+}\mathrm{1}}\mathrm{+}\mathit{\delta }{\mathit{E}}_{\mathrm{R}}^{\mathit{n}\mathrm{+}\mathrm{1}}$.

3. Considerations for AMR

3.1. Flux-limited diffusion

We have chosen to implement AMR + FLD using the deferred synchronisation algorithm of Zhang et al. (2011). This approach derives from the algorithm of Howell & Greenough (2003), but does not require a simultaneous matrix solution over multiple levels of the AMR hierarchy for flux synchronisation, thus reducing the complexity of the algorithm and improving performance.

In the AMR algorithm, at the end of ξ time steps at level l + 1, where ξ is the refinement ratio and is typically equal to 2, zones at level l (which lie under zones at level l + 1) are conservatively overwritten using the higher resolution data. This introduces differences along the edges of level l + 1 grids because zones at level l were evolved using fluxes calculated on level l, but these are generally not the same as fluxes calculated at the same location on a grid at level l + 1. To restore consistency and, more importantly, conservation, fluxes calculated at level l that are co-spatial with the edges of grids at level l + 1 are replaced by the higher resolution fluxes. Operationally, these replacements are applied directly to the (M)HD variables at level l as so-called flux corrections before another time step is taken. This process of replacement and correction is often called restriction.

In the deferred synchronisation algorithm, radiative fluxes (F_R; Eq. (6)) calculated during an implicit radiation solve are saved along the edges of grids in the same way as fluxes for the (M)HD quantities are. While corrections to the (M)HD variables at level l are applied explicitly as described above, the radiative flux corrections are instead applied during the next implicit radiation solve as a source term on the right-hand side of the radiation matrix.

For example, in Cartesian coordinates and the 1-direction, the radiative flux correction for a coarse zone (I,J,K) at level l can be written in terms of a sum in time and space over the underlying fine radiative fluxes⁵$\mathit{\delta }{\mathit{F}}_{\mathrm{R}\mathit{,}\mathrm{1}}^{\mathit{N}\mathrm{+}\mathrm{1}}\mathrm{\left(}\mathit{I,J,K}\mathrm{\right)}\mathrm{=}{{\mathit{F}}^{\mathrm{\prime }}}_{\mathrm{R}\mathit{,}\mathrm{1}}^{\mathit{N}\mathrm{+}\mathrm{1}}\mathrm{\left(}\mathit{I,J,K}\mathrm{\right)}\mathrm{-}\sum _{\mathit{\beta ,\eta ,\tau }\mathrm{=}\mathrm{0}}^{\mathit{\xi }\mathrm{-}\mathrm{1}}{{\mathit{f}}^{\mathrm{\prime }}}_{\mathrm{R}\mathit{,}\mathrm{1}}^{\mathit{n}\mathrm{+}\mathit{\tau }\mathrm{+}\mathrm{1}}\mathrm{(}\mathit{i,j}\mathrm{+}\mathit{\beta ,k}\mathrm{+}\mathit{\eta }\mathrm{)}\mathit{,}$(34)where ${{\mathit{F}}^{\mathrm{\prime }}}_{\mathrm{R}\mathit{,}\mathrm{1}}^{\mathit{N}\mathrm{+}\mathrm{1}}$ is the coarse radiative 1-flux (with units F_RΔAΔt, where Δt and ΔA = Δx₂Δx₃ are the coarse time step and area element, respectively)⁶, and ${{\mathit{f}}^{\mathrm{\prime }}}_{\mathrm{R}\mathit{,}\mathrm{1}}^{\mathit{n}\mathrm{+}\mathit{\tau }\mathrm{+}\mathrm{1}}$ are the corresponding fine radiative 1-fluxes at level l + 1 (with units f_RδAδt, where δt and δA = δx₂δx₃ are the fine time step and area element, respectively).

Once the radiative flux corrections are known, Eq. (34) and the equivalent expressions in the 2- and 3-directions are added to the right-hand side of the matrix equation (Eq. (31)) to give $\begin{array}{ccc}& & {\mathit{g}}_{\mathit{i,j,k}}\mathrm{=}{\mathit{E}}_{\mathrm{R}}^{\mathit{n}\mathrm{+}\mathrm{1}}\mathrm{-}{\mathit{E}}_{\mathrm{R}}^{\mathit{n}}\mathrm{+}\mathrm{\Delta }\mathit{tc}{\mathit{\rho }}^{\mathit{n}}\mathrm{[}{\mathit{\kappa }}_{\mathrm{E}}^{\mathit{n}\mathrm{+}\mathrm{1}}{\mathit{E}}_{\mathrm{R}}^{\mathit{n}\mathrm{+}\mathrm{1}}\mathrm{-}{\mathit{\kappa }}_{\mathrm{P}}^{\mathit{n}\mathrm{+}\mathrm{1}}{\mathit{a}}_{\mathrm{R}}\mathrm{(}{\mathit{T}}^{\mathit{n}\mathrm{+}\mathrm{1}}{\mathrm{)}}^{\mathrm{4}}\mathrm{]}\\ & & \hspace{0.17em}\hspace{0.17em}\mathrm{-}\mathrm{\Delta }\mathit{t}\mathrm{\nabla }\mathrm{·}\mathrm{(}{\mathit{D}}^{\mathit{n}}\mathrm{\nabla }{{\mathit{E}}_{\mathrm{R}}^{\mathit{n}\mathrm{+}\mathrm{1}}}^{\mathrm{)}}\mathrm{-}\mathrm{\Delta }\mathit{t}\mathrm{\left(}{\mathrm{f}}^{\hspace{0.17em}\mathit{n}}{\mathit{E}}_{\mathrm{R}}^{\mathit{n}\mathrm{+}\mathrm{1}}\mathrm{\right)}\mathrm{:\nabla }{{v}}^{\mathit{n}}\\ & & \hspace{0.17em}\hspace{0.17em}\mathrm{+}\frac{\mathrm{1}}{\mathrm{\Delta }\mathit{V}\mathrm{\left(}\mathit{I,J,K}\mathrm{\right)}}\mathit{\delta }{\mathit{F}}_{\mathrm{R}}^{\mathit{N}\mathrm{+}\mathrm{1}}\mathrm{\left(}\mathit{I,J,K}\mathrm{\right)}\mathrm{;}\\ & & \mathit{\delta }{\mathit{F}}_{\mathrm{R}}^{\mathit{N}\mathrm{+}\mathrm{1}}\mathrm{=}\mathit{\delta }{\mathit{F}}_{\mathrm{R}\mathit{,}\mathrm{1}}^{\mathit{N}\mathrm{+}\mathrm{1}}\mathrm{+}\mathit{\delta }{\mathit{F}}_{\mathrm{R}\mathit{,}\mathrm{2}}^{\mathit{N}\mathrm{+}\mathrm{1}}\mathrm{+}\mathit{\delta }{\mathit{F}}_{\mathrm{R}\mathit{,}\mathrm{3}}^{\mathit{N}\mathrm{+}\mathrm{1}}\mathit{.}\end{array}$(35)The solution of the radiation matrix then proceeds as in Sect. 2.3.

In the absence of radiation, at the end of a time step at level l, and after the flux corrections have been applied to the hydrodynamic variables, all grids at levels l → l_max are co-temporal and consistent. As discussed in detail in Zhang et al. (2011), when applying the radiative flux corrections using deferred synchronisation, consistency no longer holds at the end of a time step. Although this does not affect results in practice, it does create three complications. First, whenever grids at level l′>l are adjusted, created, or destroyed (the regridding step) in order to maintain consistency, any existing radiative flux corrections must carry over to the modified grid structure for application during the next implicit radiation solve at level l. Second, when adjacent refined (l′>l) grids are present, radiative flux corrections should only be applied to zones on grids at level l where there are no overlying refined zones. To do so is to double correct the coarse zone (by overwriting and by flux correction). For the (M)HD variables, this never causes a problem because the overwriting takes place after the flux corrections, in contrast to the deferred synchronisation algorithm. Finally, checkpoints are created at the end of a level l = 1 time step when all levels are co-temporal, and the deferred radiative flux corrections must be included in the checkpoints to ensure consistency upon a restart.

3.2. Irradiation

Coupling the ray-tracing algorithm described in Sect. 2.2 to the AMR is relatively simple and straightforward. First, we interpolate the irradiative fluxes F_irr(I,J,K) from level l to any corresponding overlying boundary zones at level l + 1 with the requirement $\begin{array}{ccc}\hspace{0.17em}\begin{array}{c}\begin{array}{c}\end{array}{\mathit{F}}_{\mathrm{irr}}\mathrm{\left(}\mathit{I,}\hspace{0.17em}\mathit{J,K}\mathrm{\right)}\mathrm{\Delta }\mathit{A}\mathrm{\left(}\mathit{I,J,K}\mathrm{\right)}\mathrm{=}\sum _{\mathit{\beta ,\eta }\mathrm{=}\mathrm{0}}^{\mathit{\xi }\mathrm{-}\mathrm{1}}{\mathit{f}}_{\mathrm{irr}}\mathrm{(}\mathit{i,}\hspace{0.17em}\mathit{j}\mathrm{+}\mathit{\beta ,k}\mathrm{+}\mathit{\eta }\mathrm{)}\begin{array}{c}\end{array}\\ \begin{array}{c}\end{array}\mathrm{\times }\mathit{\delta A}\mathrm{(}\mathit{i,}\hspace{0.17em}\mathit{j}\mathrm{+}\mathit{\beta ,k}\mathrm{+}\mathit{\eta }\mathrm{)}\mathit{,}\end{array}& & \end{array}$(36)where F_irr and f_irr have the canonical definitions of flux (in contrast to F′_R,1 and f′_{R, 1} in Eq. (34)), and ΔA, δA are defined as before. If there is instead an adjacent or overlapping grid at the same level, data is then taken directly from this grid. This is the same procedure that is applied when interpolating other face-centred quantities over a surface area in AZEuS (e.g., the magnetic field; see Ramsey et al. 2012, for details), and is generally known as prolongation.

The interpolated values are then used as input (along with the local optical depth) to the ray-tracing algorithm to determine irradiation fluxes in the active zones of each grid at level l + 1⁷, followed by the irradiative heating. At the end of a time step at level l, the RHD variables are overwritten as described in Sect. 3.1, except for the irradiative fluxes. Since the irradiative flux at any point is globally determined, and is recalculated in the active zones at the beginning of every time step, it would be pointless to perform a restriction step. Rather, it is the consistency of the local optical depth Δτ, i.e., the density and opacity, that is important.

We emphasise that the constraint on the irradiation algorithm discussed earlier remains: the time scale for changes in the optical depth must be longer than the light travel time scale. Furthermore, when coupling the ray-tracing algorithm to AMR, this constraint extends to all levels. In other words, because a particular value of the irradiative flux can depend on Δτ from a different grid, and because it may have a different resolution, in order to maintain accuracy, changes in Δτ must occur on time scales comparable to or less than the smallest time step in the simulation.

4. Numerical tests

Below we present several test problems intended to validate different aspects of our radiation hydrodynamics algorithms. Unless otherwise noted, the BiCGSTAB linear solver and Jacobi preconditioner from PETSc is used with a relative tolerance of slstol = 10^-6, and the total energy equation is employed. A description of the most important user parameters, as well as the values used herein, can be found in Tables 1 and 2. For tests of the MHD and/or AMR algorithms in AZEuS, we direct the reader to Clarke (2010) and Ramsey et al. (2012).

4.1. Linear diffusion

The first test we present follows Commerçon et al. (2011, 2014); Kolb et al. (2013), and models a radiative energy pulse as it diffuses outward; all other physical modules are deactivated. Thus, the equation being solved is $\frac{\mathit{\partial }{\mathit{E}}_{\mathrm{R}}}{\mathit{\partial t}}\mathrm{-}\mathrm{\nabla }\mathrm{·}\left(\frac{\mathit{c}}{\mathrm{3}\mathit{\rho }{\mathit{\kappa }}_{\mathrm{R}}}\mathrm{\nabla }{\mathit{E}}_{\mathrm{R}}\right)\mathrm{=}\mathrm{0}\mathit{,}$(37)where λ = 1/3 (the Eddington approximation; no flux limiting), corresponding to an optically thick regime.

We set ρκ_R = 1, which then gives the following analytical solution to Eq. (37) in 1D, ${\mathit{E}}_{\mathrm{R}}\mathrm{\left(}\mathit{x,t}\mathrm{\right)}\mathrm{=}\begin{array}{c}\mathrm{˜}\\ {\mathit{E}}_{\mathrm{o}}\end{array}\frac{\mathrm{1}}{\mathrm{2}\sqrt{\mathit{D\pi t}}}\hspace{0.17em}{\mathrm{e}}^{\mathrm{-}{\mathit{x}}^{\mathrm{2}}\mathit{/}\mathrm{4}\mathit{Dt}}\mathit{,}$(38)where D = c/ (3ρκ_R) is the diffusion coefficient, and $\begin{array}{c}{}_{\mathrm{˜}}\\ {\mathit{E}}_{\mathrm{o}}\end{array}\mathrm{=}{\mathit{E}}_{\mathrm{R}}\mathrm{(}\mathit{x}\mathrm{=}\mathrm{0}\mathit{,t}\mathrm{=}\mathrm{0}\mathrm{)}\mathrm{\Delta }\mathit{x}\mathrm{=}{\mathrm{10}}^{\mathrm{5}}$ is the initial value of the energy pulse. Under these assumptions, this test is equivalent to linear thermal conduction (e.g., Mihalas & Weibel Mihalas 1984, Sect. 103). Similar to linear conduction, a diffusion coefficient that is independent of temperature and greater than zero everywhere results in an infinite signal propagation speed, and the energy pulse will diffuse at speeds greater than the speed of light.

A Cartesian domain of size x ∈ [ − 0.5,0.5] is employed, and the pulse is initialised on the coarsest level (l = 1) as a delta function at x = 0. In AZEuS, because the refinement ratio must be a power of two, and because interpolations do not introduce new extrema, the initial delta function thus becomes a finite width pulse at higher resolutions, but always satisfies $\begin{array}{c}{}_{\mathrm{˜}}\\ {\mathit{E}}_{\mathrm{o}}\end{array}\mathrm{=}{\mathit{E}}_{\mathrm{R}}\mathrm{\Delta }{\mathit{x}}_{\mathit{l}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\mathrm{1}}$. The coarsest level has 33 zones, a refinement ratio of 2 and 3 levels of refinement are used, giving an effective resolution of 264 zones. As in Commerçon et al. (2014), zones are flagged for refinement when the gradient in the radiation energy density exceeds 25%, and the time step on the coarsest level is maintained at 3.125 × 10^-15 s. For a summary of the parameters used in this test, and others, see Table 2.

Figure 1 shows the results at t = 2 × 10^-13 s, as well as the fractional relative difference from the analytical solution (Eq. (38)). The numerical result matches well with the analytical solution, and agrees with the results of Commerçon et al. (2014). The large differences toward the edges of the diffusion front are a result of the first-order backwards Euler implicit method; indeed, the relative differences in this case do not vary strongly with spatial resolution. It should be noted that, while the relative difference always remains quantitatively similar to a uniform grid simulation of the same effective resolution, the shape of the relative difference profile depends on the specific values of the refinement criteria (e.g., ibuff, geffcy, etc.).

4.2. Radiation-matter coupling test

This standard benchmark was originally presented in Turner & Stone (2001), and tests the coupling between radiation and fluid energies. A stationary, uniform combination of radiation and fluid is used, but it is initially out of equilibrium. The radiation energy density is initialised to 10¹² erg cm^-3, and is assumed dominant relative to the fluid energy, and therefore roughly constant in time. Under these assumptions, Eqs. (13)–(14) decouple, and the system simplifies to a single ordinary differential equation: $\frac{\mathrm{d}\mathit{e}}{\mathrm{d}\mathit{t}}\mathrm{=}\mathit{c}{\mathit{\kappa }}_{\mathrm{P}}\mathit{\rho }\mathrm{(}{\mathit{E}}_{\mathrm{R}}\mathrm{-}{\mathit{a}}_{\mathrm{R}}{\mathit{T}}^{\mathrm{4}}\mathrm{)}\mathit{.}$(39)Given a constant density and Planck opacity, plus the ideal gas law, this equation can be numerically integrated to determine a reference solution.

Two versions of the test are presented using very different initial internal energy densities: 10² erg cm^-3 and 10¹⁰ erg cm^-3, but the same radiation energy density (10¹² erg cm^-3). The density and Planck mean opacity are set to ρ = 10^-7 g cm^-3 and κ_P = 0.4 cm² g^-1, respectively, while the ratio of specific heats is set to γ = 5/3, and the mean molecular weight is set to μ = 0.6. To ensure good time resolution in the plotted results, we follow Kolb et al. (2013) and set our initial time step to 10^-20 s, and then allow it to increase by 5% each subsequent cycle until t = 10^-4 s is reached. Because the conditions should remain uniform, this test should be independent of numerical resolution, but for completeness, we employ a 1D grid with a resolution of 16 zones in Cartesian coordinates. Adaptive mesh refinement is not used for this test, neither is a flux limiter (λ = 1/3), and all other physical modules are deactivated.

Figure 2 plots the numerical results for the internal energy density as a function of time, the reference solution, and the relative difference between the two. For a time step growth rate of 5% per cycle, the relative difference between numerical and reference solutions has a maximum of ~10^-2, while the relative difference of the final solution fluctuates around 2 × 10^-7. We also find that the maximum relative difference decreases roughly linearly with decreasing growth factor, as would be expected for the first-order backwards Euler method used here. It is, however, worthwhile to note that, even with a growth rate of 50%, we still obtain a final solution to within a relative difference of 5 × 10^-7 .

4.3. Marshak waves

The non-equilibrium Marshak wave test is a time-dependent non-linear diffusion test for which there are analytic solutions (Su & Olson 1996). For this test, consider an initially zero temperature, homogeneous, static, and semi-infinite medium which is illuminated on one side by some external radiation flux F_inc. The material is assumed to have a constant opacity κ, but a specific heat capacity at constant volume given by C_V = αT³, where α is a parameter. With a specific heat capacity of this form, Eqs. (13), (14) become linear in E_R and T⁴.

Following Pomraning (1979), we adopt the following dimensionless variables: $\begin{array}{ccc}& & \mathit{x}\mathrm{\equiv }\sqrt{\mathrm{3}}\mathit{\kappa z}\mathrm{;}\\ & & \mathit{\tau }\mathrm{\equiv }\left(\frac{\mathrm{4}{\mathit{a}}_{\mathrm{R}}\mathit{c\kappa }}{\mathit{\alpha }}\right)\mathit{t}\mathrm{;}\\ & & \mathit{u}\mathrm{\left(}\mathit{x,\tau }\mathrm{\right)}\mathrm{\equiv }\frac{\mathit{c}}{\mathrm{4}}\left(\frac{{\mathit{E}}_{\mathrm{R}}\mathrm{\left(}\mathit{z,t}\mathrm{\right)}}{{\mathit{F}}_{\mathrm{inc}}}\right)\mathrm{;}\\ & & \mathit{v}\mathrm{\left(}\mathit{x,\tau }\mathrm{\right)}\mathrm{\equiv }\frac{\mathit{c}}{\mathrm{4}}\left(\frac{{\mathit{a}}_{\mathrm{R}}{\mathit{T}}^{\mathrm{4}}}{{\mathit{F}}_{\mathrm{inc}}}\right)\mathrm{·}\end{array}$The initial and boundary conditions, given in dimensionless variables, are $\begin{array}{ccc}& & \mathit{u}\mathrm{\left(}\mathrm{0}\mathit{,\tau }\mathrm{\right)}\mathrm{-}\frac{\mathrm{2}}{\sqrt{\mathrm{3}}}\frac{\mathit{\partial u}\mathrm{\left(}\mathrm{0}\mathit{,\tau }\mathrm{\right)}}{\mathit{\partial x}}\mathrm{=}\mathrm{1}\mathit{,}\\ & & \mathit{u}\mathrm{\left(}\mathrm{\infty }\mathit{,\tau }\mathrm{\right)}\mathrm{=}\mathit{u}\mathrm{\left(}\mathit{x,}\mathrm{0}\mathrm{\right)}\mathrm{=}\mathit{v}\mathrm{\left(}\mathit{x,}\mathrm{0}\mathrm{\right)}\mathrm{=}\mathrm{0.}\end{array}$The additional parameter, ϵ = 4a_R/α, is set to 0.1 to permit direct comparison with the solutions of Su & Olson (1996).

A specific heat capacity C_V ∝ T³ presents a complication for the FLD algorithm presented here. Equation (15), and consequently Eqs. (16) and (20), assume that C_V is independent of temperature, and thus solving the radiation matrix as described in Sect. 2.1 will give incorrect results for this test.

The gas energy equation is easily modified to correct for this. From Eq. (13) and the ideal gas law, we can write $\begin{array}{ccc}\frac{\mathit{\partial }\mathrm{\left(}{\mathit{C}}_{\mathit{V}}\mathit{T}\mathrm{\right)}}{\mathit{\partial t}}\mathrm{=}\mathit{\alpha }\frac{\mathit{\partial }\mathrm{\left(}{\mathit{T}}^{\mathrm{4}}\mathrm{\right)}}{\mathit{\partial t}}\mathrm{=}\mathit{c}\mathrm{(}{\mathit{\kappa }}_{\mathrm{E}}{\mathit{E}}_{\mathrm{R}}\mathrm{-}{\mathit{\kappa }}_{\mathrm{P}}{\mathit{a}}_{\mathrm{R}}{\mathit{T}}^{\mathrm{4}}\mathrm{)}\mathit{,}& & \end{array}$(46)from whence an expression for the updated temperature T^{n + 1} can be derived without the need to linearise any terms: $\mathrm{(}{\mathit{T}}^{\mathit{n}\mathrm{+}\mathrm{1}}{\mathrm{)}}^{\mathrm{4}}\mathrm{=}\frac{\mathit{\alpha }\mathrm{(}{\mathit{T}}^{\mathit{n}}{\mathrm{)}}^{\mathrm{4}}\mathrm{+}\mathit{c}\mathrm{\Delta }{\mathit{t}}^{\mathit{n}}{\mathit{\kappa }}_{\mathrm{E}}^{\mathit{n}}{\mathit{E}}_{\mathrm{R}}^{\mathit{n}\mathrm{+}\mathrm{1}}}{\mathit{\alpha }\mathrm{+}\mathit{c}\mathrm{\Delta }{\mathit{t}}^{\mathit{n}}{\mathit{\kappa }}_{\mathrm{P}}^{\mathit{n}}{\mathit{a}}_{\mathrm{R}}}\mathrm{·}$(47)Equation (47) replaces Eq. (20) in the implicit radiation solver.

Following Hayes et al. (2006), we model a 1D Cartesian domain with 200 zones and a length of 8 cm. The density is set to ρ = 1 g cm^-3, the opacity to $\mathit{\kappa }\mathrm{=}\mathrm{1}\mathit{/}\sqrt{\mathrm{3}}$, and the dimensionless time step to Δτ = 3 × 10^-4 (Zhang et al. 2011). As in the previous test, AMR is not used, λ = 1/3, and all other physical modules are deactivated.

Figure 3 plots the numerical results for the dimensionless radiation energy and internal energy densities as a function of x at two different times, τ = 0.01 and τ = 0.3. The corresponding benchmark data of Su & Olson (1996) is also plotted. The agreement between numerical and analytical data is very good.

4.4. Radiative shock waves

Radiative shock tubes, where the hydrodynamic shock structure is substantially modified by the presence of radiation (e.g., Mihalas & Weibel Mihalas 1984), are an important test for any RHD code. Determining analytical solutions for this problem is not possible in general, but Lowrie & Edwards (2008) have calculated piece-wise semi-analytic solutions for the non-equilibrium diffusion case, which are suitable for comparison with numerical results.

Radiative shock waves generally fall into two categories: sub- and super-critical. In the first case, radiation from the post-shock material diffuses upstream and heats the gas, although the temperature of this precursor remains significantly below the post-shock temperature. We note that gas and radiation are out of equilibrium in this region, and that its extent depends primarily on the opacity of the gas. Meanwhile, precursor and downstream states are connected by an embedded hydrodynamic shock and a subsequent relaxation region that eventually asymptotes to downstream values.

Following Zhang et al. (2011), we initialise a stationary sub-critical M = 2 shock tube in 1D, with domain −1000 ≤ x ≤ 500 cm, and two constant states separated by a discontinuity at x = 0. The left (upstream) state has the properties ρ₀ = 5.45887 × 10^-13 g cm^-3, T₀ = 100 K, and v₀ = 2.35435 × 10⁵ cm s^-1, while the right (downstream) state is given by ρ₁ = 1.24794 × 10^-12 g cm^-3, T₁ = 207.757 K, and v₁ = 1.02987 × 10⁵ cm s^-1. Gas and radiation are initially assumed to be in equilibrium (T = T_R; T_R = [E_R/a_R] ^1/4 is the radiation temperature). The Planck and Rosseland extinction coefficients are set to α_P = ρκ_P = 3.92664 × 10^-5 cm^-1, and α_R = ρκ_R = 0.848902 cm^-1, respectively, and we use γ = 5/3, μ = 1. Fixed boundary conditions are applied on either side of the domain with values given by the initial left and right states. Following Lowrie & Edwards (2008), we apply the Eddington approximation (λ = 1/3). For the AMR, the coarsest level has a resolution of 32 zones, and four levels of refinement are employed, resulting in an effective resolution of 512 zones. Zones are flagged based on gradients in the internal and total energy densities, at a threshold of 5%. As in Zhang et al. (2011), Commerçon et al. (2014), we find that the initial discontinuity must be shifted slightly to ensure the final steady-state position of the shock lies at x = 0. For the sub-critical case, an offset of −100 cm was required, but we note that the specific value of this shift is sensitive to the values of the artificial viscosity parameters (qcon, qlin), the Courant number, and the numerical resolution.

The left side of Fig. 4 plots the gas and radiation temperatures at t = 0.1 s for the sub-critical M = 2 case, along with the relative differences from the semi-analytic solution of Lowrie & Edwards (2008). With the exception of the region near the embedded shock, where the semi-analytical solution contains a true discontinuity, the agreement is excellent and differences are less than 1%.

In the case of a super-critical shock, the velocity of the upstream material is sufficiently high that the radiation from the shock has saturated and the maximum temperature in the precursor region now matches the temperature of the downstream (post-shock) region. The precursor remains and, notwithstanding the leading edge, the gas and radiation are nearly in equilibrium. Precursor and downstream regions are now connected via an embedded hydrodynamic shock and a Zel’Dovich spike (Mihalas & Weibel Mihalas 1984), before conditions very quickly relax to the downstream state.

We set up a M = 5 super-critical shock tube in a similar manner to the sub-critical case. A domain of −4000 ≤ x ≤ 4000 cm is used, with the left state given by ρ₀ = 5.45887 × 10^-13 g cm^-3, T₀ = 100 K, v₀ = 5.88588 × 10⁵ cm s^-1, and the right state by ρ₁ = 1.96405 × 10^-12 g cm^-3, T₁ = 855.720 K, and v₁ = 1.63592 × 10⁵ cm s^-1. As before, the initial discontinuity is shifted, in this case by −310 cm, to ensure that the final shock position is x = 0. The coarsest level has 256 zones, and four levels of refinement are used, giving an effective resolution of 4096 zones. All other conditions are identical to the sub-critical test.

The right side of Fig. 4 show our results at t = 0.1 s for a super-critical M = 5 shock tube, plus the relative differences with respect to the semi-analytic solution. Again, with the exception of the embedded shock region, the numerical and semi-analytic results are in good agreement.

4.5. Protoplanetary disk relaxation benchmarks

The primary science goal behind implementing FLD and irradiation in AZEuS is to model protoplanetary disks and their evolution. To this end, we now examine the structure of an idealised accretion disk which includes not only radiative cooling, viscous heating, and stellar irradiation, but also dynamics. Inspired by Kolb et al. (2013), we initialise a disk in axisymmetric (2.5D) spherical (r,θ) coordinates, with computational domain 0.4 ≤ r/a_Jup ≤ 2.5 and 90 ≤ θ ≤ 97°, where a_Jup is the semi-major axis of Jupiter. To facilitate comparison with Kolb et al. (2013), we use a uniform numerical resolution of r × θ = 256 × 32 zones (i.e., AMR is not employed). We have, however, confirmed with higher resolution models that the results presented here are converged.

The initial density profile is given by $\mathit{\rho }\mathrm{\left(}\mathit{r,\theta }\mathrm{\right)}\mathrm{=}{\mathit{\rho }}_{\mathrm{0}}\mathrm{\left(}\mathit{R}\mathrm{\right)}\exp \left(\frac{\sin \mathit{\theta }\mathrm{-}\mathrm{1}}{{\mathit{h}}^{\mathrm{2}}}\right)\mathit{,}$(48)where R = rsinθ is the cylindrical radius, h = H/R = 0.05 is the disk aspect ratio, H is the disk scale height, ${\mathit{\rho }}_{\mathrm{0}}\mathrm{\left(}\mathit{R}\mathrm{\right)}\mathrm{=}\mathrm{\Sigma }\mathrm{\left(}\mathit{R}\mathrm{\right)}\mathit{/}\sqrt{\mathrm{2}\mathit{\pi }}\mathit{H}\mathrm{\left(}\mathit{R}\mathrm{\right)}$ is the midplane density, Σ(R) = Σ₀(R/a_Jup)^{− 1/2} is the surface density, and Σ₀ is scaled such that the mass in the upper hemisphere of the disk is 0.005 M_⊙.

The initial pressure profile is given by $\mathit{p}\mathrm{=}\mathit{\rho }{\mathit{c}}_{\mathrm{iso}}^{\mathrm{2}}\mathrm{=}\mathit{H}{\mathrm{\Omega }}_{\mathrm{K}}\mathrm{=}\mathit{H}{\mathit{v}}_{\mathrm{K}}\mathit{/}\mathit{R}$, where c_iso is the isothermal sound speed, Ω_K is the Keplerian orbital frequency, and ${\mathit{v}}_{\mathrm{K}}\mathrm{=}\sqrt{\mathit{G}{\mathit{M}}_{\mathrm{\ast }}\mathit{/}\mathit{R}}$ is the Keplerian speed. The poloidal velocity is initially zero (v_r = v_θ = 0), but the disk is rotationally supported by a slightly sub-Keplerian toroidal velocity: ${\mathit{v}}_{\mathit{\phi }}\mathrm{=}\sqrt{\mathrm{1}\mathrm{-}\mathrm{2}{\mathit{h}}^{\mathrm{2}}}{\mathit{v}}_{\mathrm{K}}\mathit{.}$(49)Pressure and centrifugal forces are offset by a point source gravitational potential with M_∗ = 1 M_⊙, centred at r = 0. The temperature is given by the ideal gas law, with γ = 5/3, μ = 2.3.

Meanwhile, for the radiation we use a Rosseland mean opacity given by the piece-wise functions of Lin & Papaloizou (1985); κ_P = κ_R is assumed. We also assume radiation and gas are initially in thermal equilibrium. The flux limiter from Kley (1989), which is optimised for disks, is applied.

With respect to boundary conditions, the toroidal velocity is set to the local Keplerian value at the inner and outer radial boundaries, and reflecting elsewhere. All other hydrodynamic boundary conditions are set to reflecting. For the radiation energy density, we set the upper θ boundary to ${\mathit{E}}_{\mathrm{R}}\mathrm{=}{\mathit{a}}_{\mathrm{R}}{\mathit{T}}_{\mathrm{amb}}^{\mathrm{4}}$, with T_amb = 5 K, and reflecting elsewhere. This prescription allows the disk to radiatively cool unimpeded at the upper boundary.

Two different test cases of a radiative disk are considered here: first, a disk with internal viscous heating only, followed by a disk with both viscous and stellar irradiation heating. In both cases, a constant kinematic viscosity of ν = 10¹⁵ cm² s^-1 is applied (Kolb et al. 2013), and the system is evolved for 100 orbits at the radius of Jupiter (i.e., 1 orbit = t₀ = 3.732 × 10⁸ s).

In the first case, an equilibrium state should be reached when radiative cooling, with a rate given by ${\mathit{Q}}_{\mathrm{rad}}^{\mathrm{-}}\mathrm{=}\mathrm{2}\mathit{\sigma }{\mathit{T}}_{\mathrm{eff}}^{\mathrm{4}}\mathit{,}$(50)balances the viscous heating, given by ${\mathit{Q}}_{\mathrm{visc}}^{\mathrm{+}}\mathrm{=}\frac{\mathrm{9}}{\mathrm{4}}\mathrm{\Sigma }\mathit{\nu }{\mathrm{\Omega }}_{\mathrm{K}}^{\mathrm{2}}\mathit{,}$(51)where both rates are per unit area and take into account the two sides of the disk; σ is the Stefan-Boltzmann constant. Equating the two rates and solving for the effective temperature gives ${\mathit{T}}_{\mathrm{eff}}^{\mathrm{4}}\mathrm{\left(}\mathit{R}\mathrm{\right)}\mathrm{=}\frac{\mathrm{9}}{\mathrm{8}}\frac{\mathit{\nu }\mathrm{\Sigma }\mathrm{\left(}\mathit{R}\mathrm{\right)}{\mathrm{\Omega }}_{\mathrm{K}}^{\mathrm{2}}\mathrm{\left(}\mathit{R}\mathrm{\right)}}{\mathit{\sigma }}\mathit{,}$(52)which describes the temperature of the gas at the disk surface. To determine the temperature at the disk midplane, we first assume that all of the viscous heating occurs at the midplane, and then invoke radiative diffusion in the vertical direction to give (e.g., Armitage 2010): $\frac{\mathit{\partial }{\mathit{T}}^{\mathrm{4}}\mathrm{\left(}\mathit{z,R}\mathrm{\right)}}{\mathit{\partial z}}\mathrm{=}\mathrm{-}\frac{\mathrm{3}\mathit{\rho }\mathrm{\left(}\mathit{z,R}\mathrm{\right)}{\mathit{\kappa }}_{\mathrm{R}}\mathrm{\left(}\mathit{z,R}\mathrm{\right)}}{\mathrm{8}\mathit{\sigma }}{\mathit{Q}}_{\mathrm{visc}}^{\mathrm{+}}\mathit{.}$(53)Integrating from the disk surface to the disk midplane, taking into account the variation in density and opacity (i.e., using Lin & Papaloizou 1985 opacities), the midplane temperature can be expressed as $\begin{array}{ccc}{\mathit{T}}_{\mathrm{mid}\mathit{,}\mathrm{visc}}^{\mathrm{4}}\mathrm{\left(}\mathit{R}\mathrm{\right)}& \mathrm{=}& {\mathit{T}}_{\mathrm{eff}}^{\mathrm{4}}\mathrm{\left(}\mathit{R}\mathrm{\right)}\left(\mathrm{1}\mathrm{+}\frac{\mathrm{3}}{\mathrm{4}}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{\infty }}\mathit{\rho }\mathrm{\left(}\mathit{z,R}\mathrm{\right)}{\mathit{\kappa }}_{\mathrm{R}}\mathrm{\left(}\mathit{z,R}\mathrm{\right)}\mathrm{d}\mathit{z}\right)\\ & \mathrm{=}& {\mathit{T}}_{\mathrm{eff}}^{\mathrm{4}}\mathrm{\left(}\mathit{R}\mathrm{\right)}\left(\mathrm{1}\mathrm{+}\frac{\mathrm{3}}{\mathrm{4}}{\mathit{\tau }}_{\mathit{z}}\mathrm{\left(}\mathit{R}\mathrm{\right)}\right)\mathit{,}\end{array}$(54)where equilibrium between viscous heating and radiative cooling has been exploited, and τ_z(R) is the vertical optical depth.

The top panel of Fig. 5 plots the midplane temperature as a function of radius in a viscously heated disk at different times. Also plotted (in black) is the midplane temperature predicted by Eq. (54) for the surface density and vertical optical depth from the t = 100t₀ model. As can be seen, by 50 orbits the temperature has effectively reached a steady state. The numerical results, however, do not agree with the analytical predictions of Eq. (54), with a factor of roughly 2^1/4 difference.

These differences can be attributed to two competing effects: First, the assumption that all of the viscous heating occurs at the midplane is false, and this will tend to decrease the disk midplane temperature. Second, in Eq. (53), the Eddington approximation is implicitly assumed (λ = 1/3), even though we employ the flux limiter of Kley (1989) in these models. Near the midplane, even though the Eddington approximation does roughly hold, it certainly does not in the upper layers of the disk (where λ ≪ 1/3). This limits the rate at which energy diffuses out of the disk, and will cause an increase in the midplane temperature. Unfortunately, analytically accounting for both of these effects is not straightforward. Although the numerical results do not agree with the simple analytical model, they are nevertheless in good agreement with previous numerical studies (Kolb et al. 2013).

In the second test case, heating by stellar irradiation is added to the model. We consider only frequency-independent (grey) irradiation by a star with effective temperature T_∗ = 5000 K and radius R_∗ = 3 R_⊙. The flux at the stellar surface is then given by ${\mathit{F}}_{\mathrm{\ast }}\mathrm{\left(}{\mathit{R}}_{\mathrm{\ast }}\mathrm{\right)}\mathrm{=}\mathit{\sigma }{\mathit{T}}_{\mathrm{\ast }}^{\mathrm{4}}$. Following Bitsch et al. (2013), the stellar opacity is set to κ_∗ = 0.1κ_P. We also assume that the disk continues inside the inner radial boundary by attenuating the stellar flux in the boundary zones using values of κ_P and ρ from the first active zone, effectively shielding the disk midplane (but not the upper layers) from direct stellar heating.

With the addition of irradiation, the disk midplane temperature should reach a new equilibrium state, given by the combined balance of viscous heating, irradiation, and radiative cooling. Under the assumption that the stellar irradiation is only absorbed at a single height in the disk, the temperature of this layer should be given by the equilibrium between stellar heating and radiative cooling. The stellar heating rate per unit area is given by Eqs. (28)and (29): $\begin{array}{ccc}{\mathit{Q}}_{\mathrm{irr}}^{\mathrm{+}}& \mathrm{=}& \mathrm{2}\mathrm{\nabla }\mathrm{·}\left({\mathit{F}}_{\mathrm{\ast }}\mathrm{\left(}{\mathit{R}}_{\mathrm{\ast }}\mathrm{\right)}\frac{{\mathit{R}}_{\mathrm{\ast }}^{\mathrm{2}}}{{\mathit{r}}^{\mathrm{2}}}{\mathrm{e}}^{\mathrm{-}{\mathit{\tau }}_{\mathrm{\ast }}\mathrm{\left(}\mathit{r}\mathrm{\right)}}\right)\\ & \mathrm{=}& \mathrm{2}{\mathit{F}}_{\mathrm{\ast }}\mathrm{\left(}{\mathit{R}}_{\mathrm{\ast }}\mathrm{\right)}\frac{{\mathit{R}}_{\mathrm{\ast }}^{\mathrm{2}}}{{\mathit{r}}^{\mathrm{2}}}{\mathrm{e}}^{\mathrm{-}{\mathit{\tau }}_{\mathrm{\ast }}\mathrm{\left(}\mathit{r}\mathrm{\right)}}\mathrm{(}\mathrm{1}\mathrm{-}{\mathrm{e}}^{\mathrm{-}\mathrm{\Delta }{\mathit{\tau }}_{\mathrm{\ast }}\mathrm{\left(}\mathit{r}\mathrm{\right)}}\mathrm{)}\mathit{,}\end{array}$where Δτ_∗(r) = ρ(r)κ_∗(r)Δr, and accounting for both sides of the disk. Equating this to the radiative cooling rate (Eq. (50)) then gives a temperature of ${\mathit{T}}_{\mathrm{surf}\mathit{,}\mathrm{irr}}^{\mathrm{4}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}{\mathit{T}}_{\mathrm{\ast }}^{\mathrm{4}}\frac{{\mathit{R}}_{\mathrm{\ast }}^{\mathrm{2}}}{{\mathit{r}}^{\mathrm{2}}}{\mathrm{e}}^{\mathrm{-}{\mathit{\tau }}_{\mathrm{\ast }}\mathrm{\left(}\mathit{r}\mathrm{\right)}}\mathrm{(}\mathrm{1}\mathrm{-}{\mathrm{e}}^{\mathrm{-}\mathrm{\Delta }{\mathit{\tau }}_{\mathrm{\ast }}\mathrm{\left(}\mathit{r}\mathrm{\right)}}\mathrm{)}\mathit{.}$(57)With the further assumption that, in the presence of irradiation heating only, the disk will be locally isothermal with a temperature profile given by Eq. (57), then the midplane temperature of a disk including both viscous heating and irradiation can be estimated from Eqs. (54) and (57) as ${\mathit{T}}_{\mathrm{mid}\mathit{,}\mathrm{visc}\mathrm{+}\mathrm{irr}}^{\mathrm{4}}\mathrm{\left(}\mathit{R}\mathrm{\right)}\mathrm{=}{\mathit{T}}_{\mathrm{mid}\mathit{,}\mathrm{visc}}^{\mathrm{4}}\mathrm{\left(}\mathit{R}\mathrm{\right)}\mathrm{+}{\mathit{T}}_{\mathrm{surf}\mathit{,}\mathrm{irr}}^{\mathrm{4}}\mathrm{\left(}\mathit{R}\mathrm{\right)}\mathit{.}$(58)The bottom panel of Fig. 5 shows the midplane temperature from the model including both viscous heating and stellar irradiation. Also plotted (in black) is the midplane temperature predicted by Eq. (58).

As in the purely viscous heating model, the temperature reaches a steady state by t = 50t₀. The disk is flat (H_p/R ≃ constant; H_p is the height where most of the stellar irradiation is absorbed), and virtually all of the stellar irradiation is absorbed at low radii, leading to higher temperatures at low R and a temperature profile in the outer disk that is dominated by viscous heating. The similar midplane temperatures in the two panels of Fig. 5 at large radii is evidence of this.

While the numerical results generally agree with the temperatures predicted by Eq. (58) at very low radii, they otherwise disagree. Some of the disagreement is certainly due to the effects discussed above, but these cannot entirely explain the observed differences. Figure 6 shows the 2D temperature structure of the disk model including stellar heating at t = 100t₀. Clearly, the disk is not locally isothermal as assumed above. Furthermore, Eq. (58) does not consider any (radial or vertical) diffusion of the stellar heating. Taking these in concert with the effects discussed for the first case, the differences between numerical and analytical results in Fig. 5 are not surprising.

4.6. Static disk radiative transfer benchmarks

One of the motivations behind using a hybrid ray-tracing plus FLD algorithm is to significantly improve the results for stellar irradiated disks relative to pure FLD. We present tests following Kuiper et al. (2010), Kuiper & Klessen (2013) that serve as a comparison of the hybrid algorithm to more exact radiative transfer algorithms, as applied to protoplanetary disks. More specifically, we simulate the static disks described in Pascucci et al. (2004), Pinte et al. (2009), and compare the results to those produced with the RADMC Monte Carlo radiative transfer code (Dullemond & Dominik 2004).

The first set of tests follows Pascucci et al. (2004), where the density structure of a disk in cylindrical coordinates is given by $\mathit{\rho }\mathrm{\left(}\mathit{R,z}\mathrm{\right)}\mathrm{=}{\mathit{\rho }}_{\mathrm{0}}{\left(\frac{\mathit{R}}{{\mathit{R}}_{\mathrm{0}}}\right)}^{-1}\exp \left(\mathrm{-}\frac{\mathit{\pi }}{\mathrm{4}}{\left(\frac{\mathit{z}}{\mathit{h}\mathrm{\left(}\mathit{R}\mathrm{\right)}}\right)}^{\mathrm{2}}\right)\mathit{,}$(59)where h(R) = z₀ (R/R₀)^1.125 is the scale height, R₀ = 500 AU, z₀ = 125 AU, and ρ₀ is used to scale the disk mass accordingly. The temperature is assumed to be initially locally isothermal with profile T(R) = T₀(R/ 1000 AU)^{− 1/3}, where T₀ = 14.7 K, and is related to the gas energy by the ideal gas law (with γ = 5/3 and μ = 1). The diffuse radiation is also assumed to initially be in equilibrium with the gas, although this will not generally be the case for the converged results.

For the stellar irradiation, we use a central source with effective temperature T_∗ = 5800 K and stellar radius R_∗ = R_⊙. Opacities are given by Draine & Lee (1984) for astronomical silicates, with a grain radius of 0.12μm and dust bulk density of 3.6 g cm^-3. The opacity table used has 61 frequency bins between 0.12μm and 2000μm (Pascucci et al. 2004). A constant dust-to-gas ratio of 0.01 is also assumed. For the diffuse radiation component, the Planck and Rosseland mean opacities are calculated by integrating the dust opacities over frequency using the canonical definitions for each mean opacity.

We apply spherical coordinates for these tests, with a domain extending from r_in = 1 ≤ r ≤ 1000 AU and from π/ 4 ≤ θ ≤ 3π/ 4, with 128 logarithmically spaced zones in the r-direction, and 180 uniformly spaced zones in the θ-direction. As the models are static, all physical modules other than irradiation and FLD are deactivated. Therefore, the only boundary conditions we need are for the gas and radiation energies. For the gas temperature, zero-gradient conditions are applied at all boundaries except for the outer radial boundary, where a value of 14.7 K is maintained at all times. The same is true for the radiation energy density, except that vacuum boundary conditions are applied at the inner radial boundary. It is also important to note that, in contrast to Sect. 4.5, the inner edge of the disk is sharp and there is no material present inside r_in. Furthermore, scattering is neglected in both the AZEuS and RADMC models.

We calculate models using both FLD + frequency-dependent stellar irradiation (with frequency bins given by the opacity data) and FLD + grey stellar irradiation, where a stellar Planck mean opacity is calculated using the effective stellar temperature and opacity table. In order to consistently obtain solutions for these tests, we use a direct LU method (rather than BiCGSTAB) to solve the radiation matrix⁸. Although slow, it has proven to be the most reliable method for these benchmarks, particularly at the extremely high optical depths of the Pinte et al. (2009) disks presented below. As it is a direct solver, the parameter slstol does not apply.

Figure 7 plots the results for disks of mass M_dust = 1.1 × 10^-7 and 1.1 × 10^-4M_⊙, which correspond to optical depths of τ_{550 nm} = 0.1 and 100, respectively (see Table 3, Pascucci et al. 2004). To obtain these results, we iterate on the temperature structure until the relative change between successive calls to the radiation solver decreases below 10^-4 (cf. Kuiper et al. 2010).

At low optical depths (Fig. 7, left), diffusion of radiation is unimportant, and the temperature structure of the disk is determined by stellar irradiation. In this case, the AZEuS results match very well with the corresponding RADMC models. Furthermore, as the disk is optically thin in all frequency bins, grey and frequency-dependent (freq.-dep.) models give very similar answers. At moderate optical depths (Fig. 7, right), the differences between grey and freq.-dep. results are suddenly quite evident. In this particular situation, the use of the Planck mean opacity in the grey model either over- or underestimates the stellar opacity at ultraviolet or infrared wavelengths, respectively, and as a result, the photon penetration depth varies significantly with frequency. Clearly, the τ_{550 nm} = 100 disk is an excellent example of when frequency-dependent treatment of stellar irradiation matters. Meanwhile, Kuiper & Klessen (2013) have already demonstrated that using purely FLD for these tests gives results that are significantly less accurate than even the combination of grey irradiation and FLD.

The relative differences between AZEuS and RADMC results are summarised in Table 3. In addition, the AZEuS results are generally in good agreement with Kuiper & Klessen (2013). Differences are seen at large radii, but this may be due to the assumed initial temperature profile, which undergoes only minimal evolution in the outer disk.

The second set of (higher optical depth) disk benchmarks performed have a density structure as described in Pinte et al. (2009)$\mathit{\rho }\mathrm{\left(}\mathit{R,z}\mathrm{\right)}\mathrm{=}{\mathit{\rho }}_{\mathrm{0}}{\left(\frac{\mathit{R}}{{\mathit{R}}_{\mathrm{0}}}\right)}^{-2.625}\exp \left(\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}{\left(\frac{\mathit{z}}{\mathit{h}\mathrm{\left(}\mathit{R}\mathrm{\right)}}\right)}^{\mathrm{2}}\right)\mathit{,}$(60)where r₀ = 100 AU, z₀ = 10 AU, and h(R) is as before. Meanwhile, for the stellar component, an effective temperature of T_∗ = 4000 K and a stellar radius of R_∗ = 2 R_⊙ are used. The opacity table used is taken from Weingartner & Draine (2001) for astronomical silicates of radius 1 μm and bulk grain densities of 3.5 g cm^-3. This opacity table contains 100 wavelength bins between 0.1 and 3000 μm.

The disks extend from r = r_in = 0.1 AU to 400 AU, and from θ = π/ 4 to 3π/ 4, and are scaled to midplane optical depths of τ_{810 nm} = 1.22 × 10³, 1.22 × 10⁴, and 1.22 × 10⁶ (Table 3; Pinte et al. 2009). To obtain the correct temperature structure, it is critical to numerically resolve strong gradients in optical depth, and in particular transitions from optically thin to thick. Unfortunately, because of the extreme optical depths of the Pinte et al. (2009) disks, the simple grid set-up used previously does not satisfy this criterion. Therefore, we instead adopt the following set-up:

τ_{810 nm} = 1.22 × 10³:: there are 20 ratioed zones in the r-direction from 0.1 to 0.15 AU, with a constant ratio of Δr(i + 1)/Δr(i) = 1.28, resulting in a minimum zone size of Δr = 1.01 × 10^-4 AU at r_in. From 0.15 to 400 AU, 100 logarithmically spaced zones are employed. In the θ-direction, we continue to use 180 uniform zones between π/ 4 and 3π/ 4.

τ_{810 nm} = 1.22 × 10⁴:: 29 ratioed zones are placed in the r-direction between 0.1 to 0.15 AU with Δr(i + 1)/Δr(i) = 1.22, followed by 100 logarithmically spaced zones. The minimum zone size is Δr = 3.5 × 10^-5 AU. In the θ-direction, there are 20 ratioed zones between π/ 4 and 3π/ 8, followed by 130 uniform zones from 3π/ 8 to 5π/ 8, and subsequently another 20 ratioed zones. The minimum Δθ in the ratioed zones is set to match smoothly with the uniform zones.

τ_{810 nm} = 1.22 × 10⁶:: 60 ratioed zones from 0.1 to 0.15 AU, with Δr(i + 1)/Δr(i) = 1.31 (Δr_min = 1.4 × 10^-9 AU), followed by 100 logarithmically spaced zones. In the θ-direction, we employ 50 ratioed zones between π/ 4 and 7π/ 10, 500 uniform zones in the middle, followed by another 50 zones between 13π/ 20 and 3π/ 4.

We have chosen not to use AMR for these tests because the extreme resolution required at the inner edge of the disk clashes with how grids are refined in the AMR framework. For example, in the τ_{810 nm} = 1.22 × 10³ case, nine levels of refinement with a refinement ratio of 4 would be required to obtain the same resolution at the inner edge, resulting in many more zones than needed using the above prescription. The situation worsens considerably for even higher optical depths.

Figure 8 plots the results for the Pinte et al. (2009) benchmarks with masses of M_dust = 3 × 10^-8, 3 × 10^-7, and 3 × 10^-5M_⊙. In these tests, the optical depths are sufficiently high that the disk is opaque to stellar photons at all frequencies, and the grey and freq.-dep. results become virtually identical.

For all the tests in this set, it can be seen that the AZEuS results overestimate the midplane temperature between (r − r_in) ≃ 0.1−10 AU by up to ~ 50%. This region corresponds to where the grey vertical optical depth transitions through one, and thus we are seeing a similar effect to the τ_{550 nm} = 100 disk with regard to the importance of frequency-dependent optical depths. In order to improve upon these results, we would need to adopt at least a frequency-dependent (i.e., multi-group) FLD approach.

The most extreme case, τ_{810 nm} = 1.22 × 10⁶, is a particularly challenging problem for any radiative transfer code. For example, Pinte et al. (2009) report that some of the codes included in the comparison did not converge for this particular test. Indeed, from the right panel of Fig. 8, it can be seen that RADMC has difficulty converging to a solution, even after we increase the number of photon packets by a factor of 10 relative to the other tests. For the AZEuS results, we find that setting nrtolrhs= 10^-14 is necessary in order to consistently obtain convergence during the Newton-Raphson iterations. Furthermore, it takes >2 orders of magnitude more wall clock time to reach convergence than for the lower opacity tests. While the right panel of Fig. 8 presents the results when the relative change in temperature drops below 10^-4 (as before), we have continued running the models and find that the dip in the midplane temperature between 0.1–10 AU disappears with ~30% more iterations. That said, the AZEuS results remain in reasonable agreement with the RADMC results, except (as before) for the vertical optical depth transition region, and are also in reasonably good agreement with Kuiper & Klessen (2013).

4.7. Shadow test

One of the continual criticisms against FLD is its inability to cast shadows (e.g., Hayes & Norman 2003). By not solving a conservation equation for the radiative flux, and by using Eq. (6) as a closure relation, the radiative flux follows gradients in E_R, and results in radiation propagating around corners. Here, we present a test derived from Hayes & Norman (2003) and Jiang et al. (2012) to demonstrate how the hybrid FLD + ray-tracing algorithm can alleviate this weakness in certain contexts. Additionally, in the second part of this test, we demonstrate that the hybrid algorithm is suitable for dynamical studies of photoevaporation, i.e., thermal winds driven by radiation.

Inspired by Jiang et al. (2012), we consider a 2D Cartesian domain of size −0.5 ≤ x ≤ 0.5 cm and −0.25 ≤ y ≤ 0.25 cm. An overdense ellipse is centred at the origin with a density profile prescribed by $\begin{array}{ccc}\mathit{\rho }\mathrm{\left(}\mathit{x,y}\mathrm{\right)}\mathrm{=}{\mathit{\rho }}_{\mathrm{0}}\mathrm{+}\frac{{\mathit{\rho }}_{\mathrm{1}}\mathrm{-}{\mathit{\rho }}_{\mathrm{0}}}{\mathrm{1}\mathrm{+}\exp \mathrm{\left(}\mathrm{10}\mathrm{\right(}\mathit{r}\mathrm{-}\mathrm{1}\mathrm{\left)}\mathrm{\right)}}\mathit{,}& & \end{array}$(61)where r = (x/a)² + (y/b)² ≤ 1, a = 0.01 cm, b = 0.06 cm, and ρ₁ = 10ρ₀. The ambient medium is initialised to a gas temperature of T₀ = 290 K and a density ρ₀ = 1 g cm^-3. The entire domain is initially in pressure equilibrium, and the diffuse radiation is assumed to be in equilibrium with the gas. The gas opacity is given by $\mathit{\alpha }\mathrm{=}\mathit{\kappa \rho }\mathrm{=}{\mathit{\alpha }}_{\mathrm{0}}{\left(\frac{\mathit{T}}{{\mathit{T}}_{\mathrm{0}}}\right)}^{-3.5}{\left(\frac{\mathit{\rho }}{{\mathit{\rho }}_{\mathrm{0}}}\right)}^{\mathrm{2}}\mathit{,}$(62)where α₀ = 1, and we have assumed that κ_P = κ_R. We apply outflow boundary conditions for all variables and sides of the box, with the exception that we set a uniform and constant irradiation source across the left boundary (at x = −0.5 cm) characterised by a black body with effective temperature T_irr = 6T₀. Equation (62) is used for both the diffuse radiation and the irradiation. For the grid, we employ two levels of refinement above the base grid, giving an effective resolution of 512 × 256 zones. Refinement is based on gradients in the gas temperature and radiative energy density, both with thresholds of 5%.

In this set-up, the ambient medium is optically thin, and the photon mean free path is equal to the length of the domain. Meanwhile, the clump is optically thick, with an interior mean free path of ≃ 3.2 × 10^-6 cm. As a result, the clump will absorb the irradiation incident upon it, and should cast a shadow behind it. Figure 9 shows the results using the hybrid FLD + ray-tracing algorithm at t = 3.3 × 10^-7 s, or 10⁴ light-crossing times in the ambient medium. Given these initial conditions, the dynamics are entirely negligible, and the solution quickly reaches a steady state. The bottom panel of Fig. 9 shows the total radiation temperature ($\mathrm{=}\mathrm{[}{\mathit{T}}_{\mathrm{rad}}^{\mathrm{4}}\mathrm{+}{\mathit{T}}_{\mathrm{irr}}^{\mathrm{4}}{\mathrm{]}}^{\mathrm{1}\mathit{/}\mathrm{4}}$, where ${\mathit{F}}_{\mathrm{irr}}\mathrm{=}\mathit{\sigma }{\mathit{T}}_{\mathrm{irr}}^{\mathrm{4}}$), and demonstrates not only the ability of the hybrid algorithm to cast shadows, but also the ability of the AMR to handle irradiation across changes in resolution.

To make things more interesting, and following Jiang et al. (2012), we increase the initial ambient temperature by a factor of 357 (⇒ T₀ = 103 530 K) while keeping the remainder of the initial conditions the same. In this case, heating due to irradiation (Eq. (28)integrated over frequency) is no longer negligible, and heats the left side of the clump. The irradiative heating establishes a pressure gradient and subsequent thermal wind which removes material from the clump over time (i.e., it photoevaporates). Figure 10 shows the results of this process at t = 1.67 × 10^-6 s (5000 light-crossing times in the initial ambient medium) for both uniform and AMR grids, demonstrating that the AMR algorithms for irradiation and FLD are capable of reproducing the uniform solution.

Although the uniform and AMR results are qualitatively very similar, there are slight differences in the shape and position of the evaporation surface. The shape and location of this surface is determined by the local irradiation flux integrated over many time steps, and the value of this flux is naturally sensitive to the column of material along the ray. Any differences in the density along the ray between uniform and AMR models could then, over time, lead to a disagreement in the shape or position of the evaporation surface. Densities are continually being interpolated in the AMR simulation to fill boundary zones for existing refined grids, or when a new refined grid is created, and these will not generally be equivalent to the densities at the same time and location in the uniform grid model. It is therefore not surprising to find differences between uniform and AMR models, but it is reassuring that they remain small.

Finally, a cautionary statement: in order to produce the results seen in the bottom panel of Fig. 10, we increased the number of AMR buffer zones (ibuff) while decreasing the grid efficiency parameter (geffcy) relative to the steady-state case (see Table 2). For a test problem such as this, where the radiation mean free path is large in much of the domain, a nearly global solution at any given location may be necessary in order to obtain the correct solution. As realised by other authors (e.g., Commerçon et al. 2014), this can be problematic for any AMR-FLD algorithm which calculates a local solution and later applies corrections, such as the deferred synchronisation algorithm used here. To minimise this problem, we have intentionally relaxed the values for geffcy and ibuff to favour a few large grids over many small grids. The alternative would be instead to use a more complex, expensive, and global multi-level solver (e.g., Howell & Greenough 2003).

5. Summary and discussion

We have presented a new implementation for RHD in the AZEuS AMR-MHD fluid code, combining two-temperature (non-equilibrium) FLD and frequency-dependent stellar irradiation using the simple and fast ray-tracing algorithm of Kuiper et al. (2010). The radiation subsystem is operator split from the hydrodynamics, and is solved implicitly using a globally convergent Newton-Raphson method. Meanwhile, we employ the flexible and freely-available PETSc library (Balay et al. 2014) to solve the radiation matrix inside the Newton-Raphson iterations. Both FLD (via the deferred synchronisation algorithm of Zhang et al. 2011) and stellar irradiation are available for use with AMR and in curvilinear coordinates.

The RHD implementation in AZEuS inherits the general limitations of the FLD approximation, including difficulties with transitions from optically thick to thin, and the loss of directional information for the radiative flux. However, in certain contexts, the combination with ray-tracing mitigates these shortcomings; the resultant hybrid algorithm is even capable of casting shadows. Furthermore, FLD is typically much less computationally expensive than other, more accurate, radiative transfer methods (e.g., VET, Jiang et al. 2012; M1 moment method, González et al. 2007; Monte Carlo methods, Haworth & Harries 2012).

We have presented several benchmarks to validate our methods and demonstrated the usefulness of the code not only for RHD simulations of PPDs, but also more general astrophysical contexts. These tests also demonstrate that AZEuS is competitive with other available AMR-RHD grid codes, even for very challenging disk radiative transfer problems⁹.

The driving motivation behind implementing RHD in AZEuS is to study photoevaporation in PPDs. The combination of FLD and frequency-dependent stellar irradiation presented here provides us with the necessary foundations to do so. In an upcoming paper (Ramsey et al., in prep.), we expand the code further to include a two-fluid approximation to enable the decoupling of dust and gas temperatures as is expected in the upper layers of PPDs (e.g., Kamp & Dullemond 2004). We also implement a simplified chemistry module, providing us with all the tools we need to examine far-UV powered (i.e., at energies 6 <hν< 13.6 eV) photoevaporation in a time-dependent and self-consistent manner.

---

Acknowledgments

This work is supported by DFG grant DU 414/9-1. We would like to thank the referee for a timely and constructive report. J.P.R. would like to thank Rolf Kuiper for helpful technical discussions on irradiation and the tests of Sect. 4.6 during the early stages of this research, as well as Bertram Bitsch for useful discussions on what developed into the tests of Sect. 4.5. J.P.R. would also like to acknowledge the NORDITA program on Photo-Evaporation in Astrophysical Systems (June 2013) where part of this work was carried out.

References

All Tables

All Figures