Issue 
A&A
Volume 574, February 2015



Article Number  A81  
Number of page(s)  15  
Section  Numerical methods and codes  
DOI  https://doi.org/10.1051/00046361/201424954  
Published online  29 January 2015 
Radiation hydrodynamics including irradiation and adaptive mesh refinement with AZEuS
I. Methods
Universität Heidelberg, Zentrum für Astronomie, Institut für Theoretische Astrophysik, AlbertÜberleStr. 2, 69120 Heidelberg, Germany
email: ramsey@uniheidelberg.de
Received: 10 September 2014
Accepted: 12 November 2014
Aims. The importance of radiation to the physical structure of protoplanetary disks cannot be understated. However, protoplanetary disks evolve with time, and so to understand disk evolution and by association, disk structure, one should solve the combined and timedependent equations of radiation hydrodynamics.
Methods. We implement a new implicit radiation solver in the AZEuS adaptive mesh refinement magnetohydrodynamics fluid code. Based on a hybrid approach that combines frequencydependent raytracing for stellar irradiation with nonequilibrium flux limited diffusion, we solve the equations of radiation hydrodynamics while preserving the directionality of the stellar irradiation. The implementation permits simulations in Cartesian, cylindrical, and spherical coordinates, on both uniform and adaptive grids.
Results. We present several hydrostatic and hydrodynamic radiation tests which validate our implementation on uniform and adaptive grids as appropriate, including benchmarks specifically designed for protoplanetary disks. Our results demonstrate that the combination of a hybrid radiation algorithm with AZEuS is an effective tool for radiation hydrodynamics studies, and produces results which are competitive with other astrophysical radiation hydrodynamics codes.
Key words: accretion, accretion disks / protoplanetary disks / hydrodynamics / methods: numerical / radiative transfer
© ESO, 2015
1. Introduction
From PoyntingRobertson 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 magnetorotational instability (MRI; e.g., Flaig et al. 2010) are certainly dynamic and multidimensional, and so require timedependent radiative transfer (RT) coupled to dynamics.
One approach is to selfconsistently 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 fluxlimited diffusion (FLD) with a simple and fast raytracing 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 reemission 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 frequencydependent algorithm is straightforward to implement and fast, the tradeoff 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, raytracing, 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; DobbsDixon 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 frequencydependent raytracing (Kuiper et al. 2010) with nonequilibrium (twotemperature) 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^{1}, 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 facecentred 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 timedependent, frequencyintegrated equations of RHD in the comoving frame, to order , and assuming local thermodynamic equilibrium 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; 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; 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_{V}T, 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 comoving 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. Fluxlimited diffusion
In order to close equation set (1)–(5), we adopt the FLD approximation, where the radiative flux is replaced by Fick’s law: (6)λ(R) is the socalled 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) (7)where (8)and ⟨⟩ denotes a spatiallyaveraged 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). Fluxlimiters 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 facecentred quantities while the product ρκ_{R} = σ_{R} is intrinsically zonecentred. 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 1direction: (9)The radiation pressure tensor is assumed to have the form: (10)where f is called the Eddington tensor, and is given by (Turner & Stone 2001; Hayes et al. 2006) and .
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 gasradiation system of equations that we wish to solve is given by 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 (15)We now write Eqs. (14) and (15) in difference form, where n and n + 1 denote original and updated values within the implicit solve. We have assumed Cartesian coordinates here for brevity (x_{1},x_{2},x_{3} = 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 timelagged, which is a voluntary choice for stability and efficiency over accuracy.
As given, Eqs. (16)–(18) form a nonlinear system of equations. If we now linearise the term in (T^{n + 1})^{4} by taking a Taylor expansion about T^{n} (Commerçon et al. 2011), (19)then (T^{n + 1})^{4} can be calculated in terms of known quantities only, (20)thus eliminating the temperature equation from the radiation system of equations^{2}. 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 nonlinear temperaturedependence of the mean opacities, and .
The righthand 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: In practice, since E_{R} is a zonecentred quantity, Eq. (21) is naturally facecentred and can therefore be applied to the momentum as is. Equation (22) is also naturally facecentred, although it is applied to the zonecentred total energy. As such, we adopt the following form (in Cartesian coordinates; cf. Zhang et al. 2011), (23)where 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), (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 userdefined 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: (25)In regards to solving the radiation system of equations, this manifests as an additional term in Eq. (20): (26)To include frequencydependent transport of radiation in the irradiation algorithm, one simply sums over all frequency bins, (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 (28)where κ_{∗ ,ν} is the frequencydependent opacity to stellar radiation, Δτ_{ν} is the change in optical depth across Δr, is the radial volume element^{3}, F_{irr,ν} is the attenuated stellar flux (29)R_{∗} is the stellar radius, 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 roundoff 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 timedependence, 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 hydrodynamicradiation time step. To obtain a grey version of the irradiation algorithm, the frequencydependence 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 singlegroup (grey) hydrodynamics, the irradiative flux is integrated over frequency before the force density is calculated: (30)The above discussion on irradiation assumes spherical coordinates and a disklike 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 via a globally convergent NewtonRaphson method (Press et al. 1992). To begin, we rewrite Eq. (17) in the form (31)where is given by Eqs. (19) and (20) (or (26)if irradiation is included), f^{n} 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 (32)where (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_{1}ΔV_{2}ΔV_{3}. This symmetrises the matrix (Hayes et al. 2006), and thus we only need to calculate and store the diagonal and superdiagonal (or subdiagonal) elements of the Jacobian.
Although AZEuS includes tridiagonal and successive overrelaxation (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 biconjugate gradient stabilised (BiCGSTAB) method plus a Jacobi (diagonal) preconditioner^{4}, even though our radiation matrix is symmetric and permits the use of simpler solvers (such as the conjugate gradient method).
Convergence in the NewtonRaphson method is signalled when either or falls below userdefined tolerances nrtoler and nrtolrhs, respectively, and where (q) denotes the number of NewtonRaphson steps taken. Furthermore, when using an iterative matrix solver, an additional userdefined tolerance is applied (slstol) to signal convergence with respect to the L_{2}norm of the preconditioned residual (, 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 .
3. Considerations for AMR
3.1. Fluxlimited 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 cospatial 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 socalled 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 righthand side of the radiation matrix.
For example, in Cartesian coordinates and the 1direction, 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^{5}(34)where is the coarse radiative 1flux (with units F_{R}ΔAΔt, where Δt and ΔA = Δx_{2}Δx_{3} are the coarse time step and area element, respectively)^{6}, and are the corresponding fine radiative 1fluxes at level l + 1 (with units f_{R}δAδt, where δt and δA = δx_{2}δx_{3} 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 3directions are added to the righthand side of the matrix equation (Eq. (31)) to give (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 cotemporal 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 cotemporal, and the deferred radiative flux corrections must be included in the checkpoints to ensure consistency upon a restart.
3.2. Irradiation
Coupling the raytracing 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 (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 facecentred 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.
Description of important userset parameters in AZEuS.
The interpolated values are then used as input (along with the local optical depth) to the raytracing algorithm to determine irradiation fluxes in the active zones of each grid at level l + 1^{7}, 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 raytracing 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 (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, (38)where D = c/ (3ρκ_{R}) is the diffusion coefficient, and 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 . 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.
Fig. 1
Linear diffusion test. Top: radiation energy density at t = 2 × 10^{13} s (black), and corresponding analytical solution (red). Bottom: relative difference between analytical and numerical solutions. Successively darker shading indicates higher levels of refinement (from l = 1 → 4). 
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 firstorder 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. Radiationmatter 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^{12} 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: (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^{2} erg cm^{3} and 10^{10} erg cm^{3}, but the same radiation energy density (10^{12} erg cm^{3}). The density and Planck mean opacity are set to ρ = 10^{7} g cm^{3} and κ_{P} = 0.4 cm^{2} 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.
Fig. 2
Radiationmatter coupling test. Top: internal energy density as a function of time given two different initial states: 10^{2} erg cm^{3} (triangles) and 10^{10} erg cm^{3} (circles). The reference solution is also plotted (red). Bottom: relative difference between numerical and reference solutions. 
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 firstorder 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 nonequilibrium Marshak wave test is a timedependent nonlinear diffusion test for which there are analytic solutions (Su & Olson 1996). For this test, consider an initially zero temperature, homogeneous, static, and semiinfinite 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^{3}, where α is a parameter. With a specific heat capacity of this form, Eqs. (13), (14) become linear in E_{R} and T^{4}.
Following Pomraning (1979), we adopt the following dimensionless variables: The initial and boundary conditions, given in dimensionless variables, are 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^{3} 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 (46)from whence an expression for the updated temperature T^{n + 1} can be derived without the need to linearise any terms: (47)Equation (47) replaces Eq. (20) in the implicit radiation solver.
Fig. 3
Marshak wave test. The dimensionless radiation energy density u and gas energy density v are plotted at times τ = 0.01 and τ = 0.3. Numerical results are shown as curves, while reference data are shown as symbols (squares: u; triangles: v). 
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 , 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 piecewise semianalytic solutions for the nonequilibrium diffusion case, which are suitable for comparison with numerical results.
Fig. 4
Radiative shock wave tests. Top: gas (red) and radiation (green) temperatures in M = 2 subcritical (left) and M = 5 supercritical (right) radiative shocks at t = 0.1 s, plus corresponding semianalytic solutions (black). Bottom: relative difference between numerical and semianalytical solutions. As before, successively darker shading indicates higher levels of refinement (from l = 1 → 5). 
Radiative shock waves generally fall into two categories: sub and supercritical. In the first case, radiation from the postshock material diffuses upstream and heats the gas, although the temperature of this precursor remains significantly below the postshock 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 subcritical 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 ρ_{0} = 5.45887 × 10^{13} g cm^{3}, T_{0} = 100 K, and v_{0} = 2.35435 × 10^{5} cm s^{1}, while the right (downstream) state is given by ρ_{1} = 1.24794 × 10^{12} g cm^{3}, T_{1} = 207.757 K, and v_{1} = 1.02987 × 10^{5} 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 steadystate position of the shock lies at x = 0. For the subcritical 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 subcritical M = 2 case, along with the relative differences from the semianalytic solution of Lowrie & Edwards (2008). With the exception of the region near the embedded shock, where the semianalytical solution contains a true discontinuity, the agreement is excellent and differences are less than 1%.
In the case of a supercritical 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 (postshock) 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 supercritical shock tube in a similar manner to the subcritical case. A domain of −4000 ≤ x ≤ 4000 cm is used, with the left state given by ρ_{0} = 5.45887 × 10^{13} g cm^{3}, T_{0} = 100 K, v_{0} = 5.88588 × 10^{5} cm s^{1}, and the right state by ρ_{1} = 1.96405 × 10^{12} g cm^{3}, T_{1} = 855.720 K, and v_{1} = 1.63592 × 10^{5} 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 subcritical test.
The right side of Fig. 4 show our results at t = 0.1 s for a supercritical M = 5 shock tube, plus the relative differences with respect to the semianalytic solution. Again, with the exception of the embedded shock region, the numerical and semianalytic 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 semimajor 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 (48)where R = rsinθ is the cylindrical radius, h = H/R = 0.05 is the disk aspect ratio, H is the disk scale height, is the midplane density, Σ(R) = Σ_{0}(R/a_{Jup})^{− 1/2} is the surface density, and Σ_{0} is scaled such that the mass in the upper hemisphere of the disk is 0.005 M_{⊙}.
The initial pressure profile is given by , where c_{iso} is the isothermal sound speed, Ω_{K} is the Keplerian orbital frequency, and is the Keplerian speed. The poloidal velocity is initially zero (v_{r} = v_{θ} = 0), but the disk is rotationally supported by a slightly subKeplerian toroidal velocity: (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 piecewise 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 , 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^{15} cm^{2} 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_{0} = 3.732 × 10^{8} s).
In the first case, an equilibrium state should be reached when radiative cooling, with a rate given by (50)balances the viscous heating, given by (51)where both rates are per unit area and take into account the two sides of the disk; σ is the StefanBoltzmann constant. Equating the two rates and solving for the effective temperature gives (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): (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 (54)where equilibrium between viscous heating and radiative cooling has been exploited, and τ_{z}(R) is the vertical optical depth.
Fig. 5
Midplane gas temperatures in the PPD relaxation tests. Top: temperature as a function of radius at different times (t_{0} = 3.732 × 10^{8} s) for the model including radiative cooling and viscous heating. Bottom: midplane temperature for the model, which also includes stellar irradiation. Temperatures predicted by Eqs. ((54); top) and (58; bottom) are shown in black. 
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_{0} 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.
Fig. 6
2D gas temperature structure at t = 100t_{0} (filled colour contours) in the PPD relaxation test including radiative cooling, viscous heating, and stellar irradiation. 
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 frequencyindependent (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 . 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): 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 (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 (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_{0}. 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_{0}. 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 raytracing 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 (59)where h(R) = z_{0} (R/R_{0})^{1.125} is the scale height, R_{0} = 500 AU, z_{0} = 125 AU, and ρ_{0} is used to scale the disk mass accordingly. The temperature is assumed to be initially locally isothermal with profile T(R) = T_{0}(R/ 1000 AU)^{− 1/3}, where T_{0} = 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 dusttogas 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 rdirection, 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, zerogradient 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.
Fig. 7
Midplane gas temperatures in the low and moderate optical depth static disk RT benchmarks (following Pascucci et al. 2004). Top: temperature as a function of radius as determined by RADMC as well as both grey and frequencydependent results as determined by AZEuS. Bottom: relative difference between AZEuS and RADMC results. We note that the RADMC and AZEuS grid points are not cospatial, and so the RADMC data is linearly interpolated to the location of AZEuS zones as necessary. 
We calculate models using both FLD + frequencydependent 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^{8}. 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.
Parameters and results for the static disk RT benchmarks.
Figure 7 plots the results for disks of mass M_{dust} = 1.1 × 10^{7} and 1.1 × 10^{4}M_{⊙}, 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 frequencydependent (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 frequencydependent 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)(60)where r_{0} = 100 AU, z_{0} = 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.
Fig. 8
Midplane gas temperatures in the high and extreme optical depth static disk RT benchmarks (following Pinte et al. 2009). Top: temperature as a function of radius as determined by RADMC as well as both grey and frequencydependent results as determined by AZEuS. Bottom: relative difference between AZEuS and RADMC results. 
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^{3}, 1.22 × 10^{4}, and 1.22 × 10^{6} (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 setup used previously does not satisfy this criterion. Therefore, we instead adopt the following setup:
τ_{810 nm} = 1.22 × 10^{3}:: there are 20 ratioed zones in the rdirection 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^{4}:: 29 ratioed zones are placed in the rdirection 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^{6}:: 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^{3} 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^{5}M_{⊙}. 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 frequencydependent optical depths. In order to improve upon these results, we would need to adopt at least a frequencydependent (i.e., multigroup) FLD approach.
The most extreme case, τ_{810 nm} = 1.22 × 10^{6}, 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 NewtonRaphson 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 + raytracing 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 (61)where r = (x/a)^{2} + (y/b)^{2} ≤ 1, a = 0.01 cm, b = 0.06 cm, and ρ_{1} = 10ρ_{0}. The ambient medium is initialised to a gas temperature of T_{0} = 290 K and a density ρ_{0} = 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 (62)where α_{0} = 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_{0}. 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%.
Fig. 9
Steadystate (T_{0} = 290 K) shadow test results at t = 3.3 × 10^{7} s, corresponding to 10^{4} lightcrossing times in the ambient medium. The source of irradiation is located at the left boundary. Top: gas temperature. Middle: radiation temperature. Bottom: total radiation temperature (diffuse + irradiation). Dashed lines denote AMR grid boundaries. 
In this setup, 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 + raytracing algorithm at t = 3.3 × 10^{7} s, or 10^{4} lightcrossing 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 (, where ), 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_{0} = 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 lightcrossing 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.
Fig. 10
Results for the evaporating shadow test (T_{0} = 103530 K) at t = 1.67 × 10^{6} s. The irradiation source is located at the left boundary. Top: uniform grid results. Bottom: AMR results at the same time and effective resolution as the uniform grid. Left: gas density. Middle: gas temperature. Right: total radiation temperature (diffuse + irradiation). Dashed lines denote AMR grid boundaries. 
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 steadystate 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 AMRFLD 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 multilevel solver (e.g., Howell & Greenough 2003).
5. Summary and discussion
We have presented a new implementation for RHD in the AZEuS AMRMHD fluid code, combining twotemperature (nonequilibrium) FLD and frequencydependent stellar irradiation using the simple and fast raytracing algorithm of Kuiper et al. (2010). The radiation subsystem is operator split from the hydrodynamics, and is solved implicitly using a globally convergent NewtonRaphson method. Meanwhile, we employ the flexible and freelyavailable PETSc library (Balay et al. 2014) to solve the radiation matrix inside the NewtonRaphson 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 raytracing 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 AMRRHD grid codes, even for very challenging disk radiative transfer problems^{9}.
The driving motivation behind implementing RHD in AZEuS is to study photoevaporation in PPDs. The combination of FLD and frequencydependent 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 twofluid 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 farUV powered (i.e., at energies 6 <hν< 13.6 eV) photoevaporation in a timedependent and selfconsistent manner.
Nestedgrid, or static mesh refinement FLD codes are also available, e.g., Tomida et al. (2013).
An alternative approach would be to eliminate the internal energy equation via a partial LU decomposition (e.g., Hayes et al. 2006; Tomida et al. 2013).
For a list of the matrix solvers and preconditioners available in PETSc, see http://www.mcs.anl.gov/petsc/documentation/linearsolvertable.html
For additional information on the code, including results from additional numerical tests, we direct the reader to the code website: http://www.ica.smu.ca/azeus/
Acknowledgments
This work is supported by DFG grant DU 414/91. 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 PhotoEvaporation in Astrophysical Systems (June 2013) where part of this work was carried out.
References
 Armitage, P. J. 2010, Astrophysics of Planet Formation (Cambridge, UK: Cambridge University Press) [Google Scholar]
 Balay, S., Abhyankar, S., Adams, M. F., et al. 2014, PETSc Web page, http://www.mcs.anl.gov/petsc [Google Scholar]
 Berger, M. J., & Colella, P. 1989, J. Comp. Phys., 82, 64 [Google Scholar]
 Bitsch, B., Crida, A., Morbidelli, A., Kley, W., & DobbsDixon, I. 2013, A&A, 549, A124 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Boley, A. C., Durisen, R. H., Nordlund, Å., & Lord, J. 2007, ApJ, 665, 1254 [NASA ADS] [CrossRef] [Google Scholar]
 Bruderer, S., Benz, A. O., Stäuber, P., & Doty, S. D. 2010, ApJ, 720, 1432 [NASA ADS] [CrossRef] [Google Scholar]
 Bruls, J. H. M. J., Vollmöller, P., & Schüssler, M. 1999, A&A, 348, 233 [NASA ADS] [Google Scholar]
 Chiang, E. I., & Goldreich, P. 1997, ApJ, 490, 368 [NASA ADS] [CrossRef] [Google Scholar]
 Clarke, D. A. 1996, ApJ, 457, 291 [NASA ADS] [CrossRef] [Google Scholar]
 Clarke, D. A. 2010, ApJS, 187, 119 [NASA ADS] [CrossRef] [Google Scholar]
 Commerçon, B., Teyssier, R., Audit, E., Hennebelle, P., & Chabrier, G. 2011, A&A, 529, A35 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Commerçon, B., Debout, V., & Teyssier, R. 2014, A&A, 563, A11 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 DobbsDixon, I., Cumming, A., & Lin, D. N. C. 2010, ApJ, 710, 1395 [NASA ADS] [CrossRef] [Google Scholar]
 Draine, B. T., & Lee, H. M. 1984, ApJ, 285, 89 [Google Scholar]
 Dullemond, C. P., & Dominik, C. 2004, A&A, 417, 159 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Dullemond, C. P., Dominik, C., & Natta, A. 2001, ApJ, 560, 957 [Google Scholar]
 Edgar, R., & Clarke, C. 2003, MNRAS, 338, 962 [NASA ADS] [CrossRef] [Google Scholar]
 Flaig, M., Kley, W., & Kissmann, R. 2010, MNRAS, 409, 1297 [NASA ADS] [CrossRef] [Google Scholar]
 Flock, M., Fromang, S., González, M., & Commerçon, B. 2013, A&A, 560, A43 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 González, M., Audit, E., & Huynh, P. 2007, A&A, 464, 429 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gorti, U., Dullemond, C. P., & Hollenbach, D. 2009, ApJ, 705, 1237 [NASA ADS] [CrossRef] [Google Scholar]
 Haworth, T. J., & Harries, T. J. 2012, MNRAS, 420, 562 [Google Scholar]
 Hayes, J. C., & Norman, M. L. 2003, ApJS, 147, 197 [Google Scholar]
 Hayes, J. C., Norman, M. L., Fiedler, R. A., et al. 2006, ApJS, 165, 188 [NASA ADS] [CrossRef] [Google Scholar]
 Howell, L. H., & Greenough, J. A. 2003, J. Comp. Phys., 184, 53 [NASA ADS] [CrossRef] [Google Scholar]
 Jiang, Y.F., Stone, J. M., & Davis, S. W. 2012, ApJS, 199, 14 [NASA ADS] [CrossRef] [Google Scholar]
 Kamp, I., & Dullemond, C. P. 2004, ApJ, 615, 991 [NASA ADS] [CrossRef] [Google Scholar]
 Klassen, M., Kuiper, R., Pudritz, R. E., et al. 2014, ApJ, 797, 4 [NASA ADS] [CrossRef] [Google Scholar]
 Kley, W. 1989, A&A, 208, 98 [NASA ADS] [Google Scholar]
 Kley, W., Bitsch, B., & Klahr, H. 2009, A&A, 506, 971 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kolb, S. M., Stute, M., Kley, W., & Mignone, A. 2013, A&A, 559, A80 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Krumholz, M. R., Klein, R. I., McKee, C. F., & Bolstad, J. 2007, ApJ, 667, 626 [NASA ADS] [CrossRef] [Google Scholar]
 Kuiper, R., & Klessen, R. S. 2013, A&A, 555, A7 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kuiper, R., Klahr, H., Dullemond, C., Kley, W., & Henning, T. 2010, A&A, 511, A81 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Levermore, C. D., & Pomraning, G. C. 1981, ApJ, 248, 321 [NASA ADS] [CrossRef] [Google Scholar]
 Lin, D. N. C., & Papaloizou, J. 1985, in Protostars and Planets II, eds. D. C. Black, & M. S. Matthews, 981 [Google Scholar]
 Lowrie, R. B., & Edwards, J. D. 2008, Shock Waves, 18, 129 [NASA ADS] [CrossRef] [Google Scholar]
 Mihalas, D., & Weibel Mihalas, B. 1984, Foundations of radiation hydrodynamics (New York: Oxford University Press) [Google Scholar]
 Min, M., Dullemond, C. P., Kama, M., & Dominik, C. 2011, Icarus, 212, 416 [NASA ADS] [CrossRef] [Google Scholar]
 Minerbo, G. N. 1978, J. Quant. Spec. Radiat. Transf., 20, 541 [NASA ADS] [CrossRef] [Google Scholar]
 Murray, S. D., Castor, J. I., Klein, R. I., & McKee, C. F. 1994, ApJ, 435, 631 [NASA ADS] [CrossRef] [Google Scholar]
 Pascucci, I., Wolf, S., Steinacker, J., et al. 2004, A&A, 417, 793 [Google Scholar]
 Pinte, C., Harries, T. J., Min, M., et al. 2009, A&A, 498, 967 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Pomraning, G. C. 1979, J. Quant. Spec. Radiat. Transf., 21, 249 [NASA ADS] [CrossRef] [Google Scholar]
 Press, W. H., Teukolsky, S. A., Vetterling, W. T., & Flannery, B. P. 1992, Numerical recipes in FORTRAN. The art of scientific computing, 2nd edn. (Cambridge: University Press) [Google Scholar]
 Ramsey, J. P., Clarke, D. A., & Men’shchikov, A. B. 2012, ApJS, 199, 13 [NASA ADS] [CrossRef] [Google Scholar]
 Su, B., & Olson, G. L. 1996, J. Quant. Spec. Radiat. Transf., 56, 337 [Google Scholar]
 Tomida, K., Tomisaka, K., Matsumoto, T., et al. 2013, ApJ, 763, 6 [Google Scholar]
 Turner, N. J., & Stone, J. M. 2001, ApJS, 135, 95 [Google Scholar]
 van der Holst, B., Tóth, G., Sokolov, I. V., et al. 2011, ApJS, 194, 23 [NASA ADS] [CrossRef] [Google Scholar]
 Von Neumann, J., & Richtmyer, R. D. 1950, J. Appl. Phys., 21, 232 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Weingartner, J. C., & Draine, B. T. 2001, ApJ, 548, 296 [NASA ADS] [CrossRef] [Google Scholar]
 Whalen, D., & Norman, M. L. 2006, ApJS, 162, 281 [NASA ADS] [CrossRef] [Google Scholar]
 Wolfire, M. G., & Cassinelli, J. P. 1986, ApJ, 310, 207 [Google Scholar]
 Zhang, W., Howell, L., Almgren, A., Burrows, A., & Bell, J. 2011, ApJS, 196, 20 [NASA ADS] [CrossRef] [Google Scholar]
All Tables
All Figures
Fig. 1
Linear diffusion test. Top: radiation energy density at t = 2 × 10^{13} s (black), and corresponding analytical solution (red). Bottom: relative difference between analytical and numerical solutions. Successively darker shading indicates higher levels of refinement (from l = 1 → 4). 

In the text 
Fig. 2
Radiationmatter coupling test. Top: internal energy density as a function of time given two different initial states: 10^{2} erg cm^{3} (triangles) and 10^{10} erg cm^{3} (circles). The reference solution is also plotted (red). Bottom: relative difference between numerical and reference solutions. 

In the text 
Fig. 3
Marshak wave test. The dimensionless radiation energy density u and gas energy density v are plotted at times τ = 0.01 and τ = 0.3. Numerical results are shown as curves, while reference data are shown as symbols (squares: u; triangles: v). 

In the text 
Fig. 4
Radiative shock wave tests. Top: gas (red) and radiation (green) temperatures in M = 2 subcritical (left) and M = 5 supercritical (right) radiative shocks at t = 0.1 s, plus corresponding semianalytic solutions (black). Bottom: relative difference between numerical and semianalytical solutions. As before, successively darker shading indicates higher levels of refinement (from l = 1 → 5). 

In the text 
Fig. 5
Midplane gas temperatures in the PPD relaxation tests. Top: temperature as a function of radius at different times (t_{0} = 3.732 × 10^{8} s) for the model including radiative cooling and viscous heating. Bottom: midplane temperature for the model, which also includes stellar irradiation. Temperatures predicted by Eqs. ((54); top) and (58; bottom) are shown in black. 

In the text 
Fig. 6
2D gas temperature structure at t = 100t_{0} (filled colour contours) in the PPD relaxation test including radiative cooling, viscous heating, and stellar irradiation. 

In the text 
Fig. 7
Midplane gas temperatures in the low and moderate optical depth static disk RT benchmarks (following Pascucci et al. 2004). Top: temperature as a function of radius as determined by RADMC as well as both grey and frequencydependent results as determined by AZEuS. Bottom: relative difference between AZEuS and RADMC results. We note that the RADMC and AZEuS grid points are not cospatial, and so the RADMC data is linearly interpolated to the location of AZEuS zones as necessary. 

In the text 
Fig. 8
Midplane gas temperatures in the high and extreme optical depth static disk RT benchmarks (following Pinte et al. 2009). Top: temperature as a function of radius as determined by RADMC as well as both grey and frequencydependent results as determined by AZEuS. Bottom: relative difference between AZEuS and RADMC results. 

In the text 
Fig. 9
Steadystate (T_{0} = 290 K) shadow test results at t = 3.3 × 10^{7} s, corresponding to 10^{4} lightcrossing times in the ambient medium. The source of irradiation is located at the left boundary. Top: gas temperature. Middle: radiation temperature. Bottom: total radiation temperature (diffuse + irradiation). Dashed lines denote AMR grid boundaries. 

In the text 
Fig. 10
Results for the evaporating shadow test (T_{0} = 103530 K) at t = 1.67 × 10^{6} s. The irradiation source is located at the left boundary. Top: uniform grid results. Bottom: AMR results at the same time and effective resolution as the uniform grid. Left: gas density. Middle: gas temperature. Right: total radiation temperature (diffuse + irradiation). Dashed lines denote AMR grid boundaries. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.