A GENERAL HYBRID RADIATION TRANSPORT SCHEME FOR STAR FORMATION SIMULATIONS ON AN ADAPTIVE GRID

We use the publicly available flash high-performance general application physics code, currently in its fourth major version (Fryxell et al. 2000; Dubey et al. 2009). The code is modular, and has physics capabilities that now include 2T radiation hydrodynamics (our contribution), magnetohydrodynamics, multi-group flux-limited diffusion, self-gravity, and a variety of options for the equation of state.

The code has been expanded to include sink particles, originally implemented in version 2.5 (Federrath et al. 2010), and then later ported to version 4.0 (Safranek-Shrader et al. 2012). Multigroup flux-limited diffusion for radiation hydrodynamics was added in version 4.0 for the study of high-energy density physics (HEDP), such as radiative shocks and laser energy deposition experiments. Despite these original motivations, the code is general and can be applied to astrophysical scenarios.

A time-independent raytracer with hybrid characteristics was added in version 2.5 by Rijkhorst et al. (2006) and then modified for the study of collapse calculations by Peters et al. (2010a). The raytracer was used to calculate photoelectric heating and photoionization rates, as well as heating by optically thin nonionizing radiation, to study H ii regions and ionization feedback in molecular clouds (Peters et al. 2010a, 2010b, 2010c, 2011, 2012; Klassen et al. 2012). However, this approach is not suitable to propagate non-ionizing radiation in optically thick clouds.

We ported this raytracer into the present version of the flash code in order to calculate the irradiation of gas and dust by point sources within the domain and solve the radiation hydrodynamic equations.

flash solves the fluid equations on an Eulerian mesh with adaptive mesh refinement (AMR) using the piecewise parabolic method (Colella & Woodward 1984). The method is well-suited for dealing with shocks. The diffusion equation is solved implicitly via the generalized minimal residual (GMRES) method. The refinement criterion used by flash is an error estimator based on Löhner (1987), which is a modified second derivative normalized by the average of the gradient over one computational cell. It has the advantage of being dimensionless and local, so one has the flexibility to choose the grid variable upon which to refine.

In this paper we describe the implementation of a hybrid radiation transfer scheme, involving raytracing coupled to a flux-limited diffusion (FLD) solver. The raytracer finds the flux each cell receives from all point sources (e.g., stars) in the simulation while the FLD solver evolves the diffuse radiation field. The combined radiation field is used to update the matter temperature.

In the sections below we recapitulate the basic theory and how these equations are implemented in flash.

If the radiation fields are assumed to have a blackbody spectrum, then in the absence of magnetic and gravitational fields, and assuming local thermodynamic equilibrium (LTE), the frequency- and angle-integrated radiation hydrodynamic equations in the comoving frame are (Turner & Stone 2001; Mihalas & Mihalas 1984)

$$\begin{eqnarray} &\frac{D \rho }{Dt} + \rho \boldsymbol{\nabla }\boldsymbol{\cdot } \boldsymbol{v} = 0 \end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray} &\rho \frac{Dv}{Dt} = {-} \boldsymbol{\nabla }p + \frac{1}{c}\kappa _P \rho \boldsymbol{F}_r \end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray} &\rho \frac{D}{Dt}\left(\frac{e}{\rho }\right) + p \boldsymbol{\nabla }\boldsymbol{\cdot } \boldsymbol{v} = - \kappa _P \rho (4 \pi B - c E_r) \end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray} &\hspace{-8pt}\rho \frac{D}{Dt}\left(\frac{E_r}{\rho }\right) + \boldsymbol{\nabla }\boldsymbol{\cdot } \boldsymbol{F}_r + \boldsymbol{\nabla }\boldsymbol{v} \boldsymbol{:} \boldsymbol{P} = \kappa _P \rho (4 \pi B - c E_r), \end{eqnarray} \tag{ 7 }$$

where $D/Dt \equiv \partial /\partial t + \boldsymbol{v \cdot \nabla }$ is the convective derivative. The fluid variables $\rho, e, \boldsymbol{v}$, and p are the matter density, internal energy density, fluid velocity, and scalar pressure, respectively, while $E_r, \boldsymbol{F}_r$, and $\boldsymbol{P}$ are the total frequency-integrated radiation energy, flux, and pressure, respectively. Equations (4) and (5) express the conservation of mass and momentum, respectively, where the radiation applies a force proportional to $\kappa _P \rho \boldsymbol{F}_r/c$, where κ_P is the Planck mean opacity defined in Equation (9) below.

The hydrodynamic equations must also be closed via an equation of state. We use a gamma-law equation of state for simple ideal gases. The gas pressure P, density ρ, internal energy , and gas temperature T are related by the equations,

$$\begin{equation} P = (\gamma - 1)\rho \epsilon = \frac{N_a k_B}{\mu }\rho T, \end{equation} \tag{ 8 }$$

with adiabatic index γ = 5/3. N_a is Avogadro's number, k_B is the Boltzmann constant, and μ is the mean molecular weight.

The coupling between internal and radiation energy is expressed by Equations (6) and (7). Coupling comes through the emission and absorption of radiation energy. The assumption of LTE lets us write the source function as the Planck function, B, so that matter emits thermal radiation proportional to the rate,

$$\begin{equation*} 4 \pi \kappa _P \rho B(T) \equiv \kappa _P \rho c a T^4, \end{equation*}$$

where κ_P is the Planck mean opacity,

$$\begin{equation} \kappa _P = \frac{\int _0^\infty d \nu \kappa _{\nu }B_{\nu }(T)}{B(T)}, \end{equation} \tag{ 9 }$$

with

$$\begin{equation} B_{\nu }(T) = \frac{2 h \nu ^3 / c^2}{e^{h\nu /k_B T} - 1}. \end{equation} \tag{ 10 }$$

Meanwhile, matter is absorbing radiation out of the thermal field at a rate κ_PρcE_r. Energy lost by the radiation field is gained by the matter and vice versa. We will use the flux-limited diffusion approximation to relate the radiation flux $\boldsymbol{F}_r$ to the radiation energy E_r.

In practice, many grid codes implement the flux-limited diffusion approximation (Turner & Stone 2001; Krumholz et al. 2007b; Commerçon et al. 2011), which relates the radiation flux $\boldsymbol{F}_r$ to the radiation energy E_r. In the FLD approximation, the radiation and matter components may be thought of as two fluids, each with their own equation of state, internal energy, temperature, and pressure. The matter component may be further split into ion and electron components, which exchange energy via collisions. This is three-temperature, or "3T" approach is used in the modeling of laboratory plasmas. The flash code was extended in version 4 to include 3T radiation hydrodynamics (Lamb et al. 2010; Fatenejad et al. 2012; Kumar et al. 2011). Energy exchange under the 3T model is included for reference in Appendix A. In astrophysics, we are more concerned with modeling the gas/dust mixture of the ISM, with the gas sometimes ionized. By setting a unified matter temperature in flash, with T_ion = T_ele ≡ T_gas/dust, we return to a 2T model, with only matter and radiation temperatures varying and exchanging energy.

Radiation only behaves like a fluid when it is tightly coupled to matter, which occurs only in very optically thick regions. Hybrid radiation schemes more accurately transport the radiation energy through the computational domain.

Under the flux-limited diffusion approximation (Levermore & Pomraning 1981; Bodenheimer et al. 1990), the flux of radiation is proportional to the gradient in the radiation energy density,

$$\begin{equation} \boldsymbol{F}_r = - D \boldsymbol{\nabla }E_r, \end{equation} \tag{ 11 }$$

where

$$\begin{equation} D = \frac{\lambda c}{\kappa _R \rho } \end{equation} \tag{ 12 }$$

is the diffusion coefficient. κ_R is the Rosseland mean opacity and λ is the flux limiter that bridges the radiation diffusion rate between the optically thin and optically thick regimes. In the extreme optically thick limit, λ → 1/3, which is the diffusion limit. In the extreme optically thin limit, the flux of radiation becomes F_r = cE_r, the free-streaming limit.

We use the Levermore & Pomraning (1981) flux limiter, one of the most commonly used,

$$\begin{equation} \lambda = \frac{2+R}{6+3R+R^2}, \end{equation} \tag{ 13 }$$

with

$$\begin{equation} R = \frac{|\boldsymbol{\nabla }E_r|}{\kappa _R \rho E_r}. \end{equation} \tag{ 14 }$$

Other flux limiter are possible, such as Minerbo (1978), another popular choice. Different flux limiters result from different assumptions about the angular distribution of the specific intensity (Turner & Stone 2001).

flash solves the radiation diffusion Equation (7) using a general implicit diffusion solver. Using the method of operator splitting, terms such as advection and hydrodynamic work are handled by the code's hydro solver. The radiation diffusion solver then updates the radiation energy equation by solving

$$\begin{equation} \frac{\partial E_r}{\partial t} + \boldsymbol{\nabla }\boldsymbol{\cdot } \left(\frac{\lambda c}{\kappa _R \rho } \boldsymbol{\nabla }E_r\right) = \kappa _P \rho c (a T^4 - E_r), \end{equation} \tag{ 15 }$$

where T is the gas temperature.

To solve this equation, we use the method of the generalized minimum residual (GMRES) by Saad & Schultz (1986). It is included in flash via the hypre library (Falgout & Yang 2002). The GMRES method is the same one employed by Kuiper et al. (2010b) in their hybrid radiation-transport scheme. It belongs to the class of Krylov subspace methods for iteratively solving systems of linear equations of the form $A {x} = \boldsymbol{b}$, where the matrix A is large and must be inverted. Matrix inversion is, in general, a very computationally expensive operation, but because A is also a sparse matrix, methods exist for rapidly computing an approximate inverse with relatively high accuracy, outperforming methods such as conjugate gradient (CG) and successive over-relaxation (SOR). The hypre library that flash uses is a collection of sparse matrix solvers for massively parallel computers.

A difference in matter and radiation temperatures in Equation (15) results in an energy "excess" on the right-hand side. We must update the matter temperature due to the action of the combined radiation fields. Restating Equations (6) and (7), with operator-split hydrodynamic terms removed, and a stellar radiation source term added, we have

$$\begin{eqnarray} \partial _t \rho \epsilon &= - \kappa _P \rho c (a_R T^4 - E_r) - \boldsymbol{\nabla }\boldsymbol{\cdot } \boldsymbol{F}_* \end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray} \partial _t E_r + \boldsymbol{\nabla }\cdot \boldsymbol{F}_r &= + \kappa _P \rho c (a_R T^4 - E_r). \end{eqnarray} \tag{ 17 }$$

Discretizing Equation (16)) and (17), we have

$$\begin{eqnarray} &&\frac{\rho c_V T^{n+1} - \rho c_V T^{n}}{\Delta t} = -\kappa _P^n \rho c \left(a_R (T^{n+1})^4 - E_r^{n+1} \right) - \boldsymbol{\nabla }\boldsymbol{\cdot } \boldsymbol{F}_*, \nonumber\\ \end{eqnarray} \tag{ 18 }$$

and

$$\begin{equation} \hspace{-3pt}\frac{E_r^{n+1} - E_r^{n}}{\Delta t} - \nabla \cdot \left(D^{n} \nabla E_r^{n+1} \right) = + \kappa \rho c \left(a_R (T^{n+1})^4 - E_r^{n+1}\right), \end{equation} \tag{ 19 }$$

recalling that the specific internal energy = c_VT, with c_V being the specific heat capacity of the matter. Variables with superscript indices n and n + 1 take their values from before and after the implicit update, respectively. It is important to remember that the opacity κ_P is temperature-dependent. The source term for stellar radiation takes the form $\boldsymbol{\nabla }\boldsymbol{\cdot } \boldsymbol{F}_*$ and represents the amount of stellar radiation energy absorbed at a given location in the grid from all sources.

The presence of the nonlinear term (T^{n + 1})⁴ makes it difficult to solve for the temperature in (18). Assuming that the change in temperature is small, Commerçon et al. (2011) linearizes this term,

$$\begin{eqnarray} (T^{n+1})^4 &= (T^n)^4\left(1+\frac{T^{n+1} - T^n}{T^n}\right)^4 \nonumber \\ &\approx 4 (T^n)^3 T^{n+1} - 3 (T^n)^4. \end{eqnarray} \tag{ 20 }$$

This allows for Equation (20) to be substituted into Equation (18) and for us to write an expression for the updated temperature,

$$\begin{equation} T^{n+1} = \frac{3 a_R \alpha \left(T^n\right)^4 + \rho c_V T^n + \alpha E_r^{n+1} - \Delta t \, \boldsymbol{\nabla }\boldsymbol{\cdot } \boldsymbol{F}_*}{\rho c_V + 4 a_R \alpha \left(T^n\right)^3}, \end{equation} \tag{ 21 }$$

where we have used $\alpha \equiv \kappa _P^n \rho c \Delta t$.

The irradiation term, $\boldsymbol{\nabla }\boldsymbol{\cdot } \boldsymbol{F}_*$, is supplied by the raytracer. Its calculation is the subject of Section 3.4.

The raytrace is performed first in order to calculate $\boldsymbol{\nabla }\boldsymbol{\cdot } \boldsymbol{F}_*$. By substituting Equation (21) into (19) via the approximation of Equation (20), and making appropriate simplifications, the FLD equation becomes

$$\begin{eqnarray} && \frac{E_r^{n+1} - E_r^{n}}{\Delta t} - \nabla \cdot \left(D^{n} \nabla E_r^{n+1}\right)\nonumber\\ &&\quad + \kappa \rho c \left[\frac{ \rho c_V T^{n}}{\rho c_V T^{n} + 4 \alpha a_R \left(T^{n}\right)^4}\right] E_r^{n+1} \nonumber\\ &&\quad = \kappa \rho c a_R \left(T^{n}\right)^4 \left[\frac{\rho c_V T^{n} - 4 \Delta t \nabla \cdot F_*}{\rho c_V T^{n} + 4 \alpha a_R \left(T^{n}\right)^4} \right], \end{eqnarray} \tag{ 22 }$$

where we have gathered implicit terms on the left-hand side and explicit terms on the right-hand side.

We solve Equation (22) via the diffusion solver to calculate the updated radiation energy density and temperature, then update the gas/dust temperature via (21). This completes the radiation transfer update for the current timestep.

When the effects of scattering and emission are ignored, the equation for the specific intensity of radiation from a point source along a ray assumes a very simple form,

$$\begin{equation} I(r) = I_0 e^{-\tau (r)}, \end{equation} \tag{ 23 }$$

where I₀ is the specific intensity of the source and τ(r) is the optical depth from the source up to r:

$$\begin{equation} \tau (r) = \int _{R_*}^r \kappa _P(T_*) \rho (r^{\prime }) \, dr^{\prime }. \end{equation} \tag{ 24 }$$

The opacity κ_P is a function of temperature and the optical depth is calculated using the temperature of the radiation source, e.g., the effective temperature of a star. ρ is the usual gas density. Rays are traced from individual cells in the computational grid back to the point sources (usually stars) via a characteristics-based method described in Rijkhorst et al. (2006).

In characteristics-based raytracing, there are two approaches to computing the column density toward sources, each with its own advantages and disadvantages. Long characteristics involve tracing rays from the source to each cell in the computational domain. This method is very accurate, but involves many redundant calculations of the column density through cells intersected by similar rays near the source. Short characteristics try to mitigate the redundant calculations by tracing rays only across the length of each grid cell in the direction of the source, then interpolating upwind toward the source to calculate the total column density. The disadvantage of this approach is that the upwind values need to be known ahead of time, which imposes an order on the calculations and makes scaling the code to many processors very problematic. Some amount of numerical diffusivity is also introduced in this approach on account of interpolating column densities at each cell.

As a compromise between these two methods, Rijkhorst et al. (2006) developed a method they called hybrid characteristics that minimizes the number of redundant calculations and can be parallelized. The local contributions to the column density are computed for each block of n × n × n cells. Then the values on the block faces are computed and communicated to all processors. Finally, the total column density at each cell is computed by adding the local and interpolated face values.

Peters et al. (2010a) have improved the method to allow collapse simulations with arbitrary many refinement levels and sources of radiation. They added the propagation of (optically thin) nonionizing radiation and coupled the respective heating terms to a prescription of molecular and dust cooling (Banerjee et al. 2006). Furthermore, they have linked the radiation module to sink particles (Federrath et al. 2010) as sources of radiation and implemented a simple prestellar model to determine the value of the stellar and accretion luminosities. In our current implementation, we no longer use a separate cooling function because the thermal evolution is already implicit in the coupled equations of internal and radiation energy, as described in Section 3.3.

The stellar flux a distance r is given by

$$\begin{equation} F_*(r) = F_*(R_*) \left(\frac{R_*}{r}\right)^2 \exp \left(-\tau (r)\right), \end{equation} \tag{ 25 }$$

where R_* is the stellar radius, and F_*(R_*) is the flux at the stellar surface, given by

$$\begin{equation} F_*(R_*) = \sigma T_*^4. \end{equation} \tag{ 26 }$$

T_* is the effective temperature of the stellar surface and σ is the Stefan-Boltzmann constant.

At this time, our raytracing is purely gray; we use the Draine & Lee (1984) dust model and compute the temperature-dependent frequency-averaged opacities, assuming a 1% dust-to-gas ratio. Multifrequency raytracing is relatively straightforward to implement, but the computational costs increase linearly with the number of frequency bins. We leave the implementation of multifrequency effects to a future paper.

Using the stellar flux computed at each cell, we estimate an "irradiation" term, i.e., the amount of stellar flux absorbed by a given cell,

$$\begin{equation} \boldsymbol{\nabla }\cdot F_*(r) \approx - \frac{\left(1 - e^{\tau _{\mathrm{local}}}\right)}{\Delta r} F_*(r), \end{equation} \tag{ 27 }$$

where τ_local is the "local" contribution to the optical depth through the cell, i.e.,

$$\begin{equation} \tau _{\mathrm{local}} \approx \kappa _P(T_*) \rho \, \Delta r, \end{equation} \tag{ 28 }$$

and Δr is the distance traced by the ray through a grid cell. For numerical stability, when τ_local is very small, e.g., 10⁻⁵, we estimate the irradiation by Taylor-expanding the exponential,

$$\begin{equation} \boldsymbol{\nabla }\cdot F_*(r) \approx - \kappa _P \rho F_*(r), \end{equation} \tag{ 29 }$$

which is the optically thin limit.

The primary source of opacity in the interstellar medium is dust grains. We use the tabulated optical properties of graphite and silicate dust grains from Draine & Lee (1984). They evaluate absorption cross-sections for particles with sizes between 0.003 and 1.0 μm at wavelengths between 300 Å and 1000 μm. These tabulated dust properties are widely used both in simulations (e.g., Pascucci et al. 2004; Offner et al. 2009; Flock et al. 2013) and in observational work (e.g., Nielbock et al. 2012).

We assume a dust to gas ratio of 1% and integrate the dust opacity over a wavelength range of 1 Å–1000 μm. This is done as a function of temperature, resulting in "gray" (frequency-averaged) temperature-dependent Planck and Rosseland mean opacity tables. Figure 1 shows, e.g., the opacity κ over a range of temperatures from 0.1 K to 2000 K. We see that these opacities cover many orders of magnitude.

Zoom In Zoom Out Reset image size

Both the direct radiation field and the diffuse radiation field contribute to the radiation pressure. In the case of the diffuse (thermal) radiation field, flash makes the Eddington approximation,

$$\begin{equation} P_{\mathrm{rad}} = \frac{1}{3} E_r = \frac{aT_{\mathrm{rad}}^4}{3}, \end{equation} \tag{ 30 }$$

where E_r is the radiation energy density and T_rad is the corresponding temperature. Additional momentum is added to the gas from the direct component of the radiation field. We take the stellar fluxes computed by the raytracer (Equation (25)). These exert a body force

$$\begin{equation} f_{\mathrm{rad}} = \rho \kappa _P(T_*) \frac{F_*}{c} = - \frac{\boldsymbol{\nabla }\boldsymbol{\cdot } \boldsymbol{F}_*}{c}, \end{equation} \tag{ 31 }$$

where f_rad is the force density and F_* is the direct stellar radiation flux. The temperature of the direct radiation field is that of the source, T_* = T_eff, the effective stellar surface temperature.

From Equation (31) we see that the radiation force density is proportional to the absorbed radiation energy computed by the raytracer, i.e., the force felt by the gas and dust in a local region is proportional to its absorbed stellar radiation energy.

We do not expect our radiation force to change much if it were to be computed with a multifrequency raytracer. The radiation force is proportional to the absorbed flux. In the multifrequency case, infrared radiation can penetrate further into the gas, but also carries less energy. Most of the higher frequencies are all absorbed in the same (few) grid cells. Tests by Kuiper et al. (2010b) also showed very little deviation between the gray and multifrequency cases.

In summary, we treat radiation hydrodynamics by complementing the fluid dynamics equations to include the effects of radiation (Equations (4) through (7)). We have decomposed the radiation field into a direct stellar component and an indirect diffuse component. The method of raytracing is used to calculate the direct stellar radiation field. The method of flux-limited diffusion is used to transport radiation in the diffuse radiation field.

Matter and radiation are coupled by emission and absorption processes, while stars represent sources of radiation energy. This relationship is given in Equations (16) and (17).

We have discretized and linearized Equations (16) and (17) to derive an expression for evolving the matter temperature, given absorption and emission of radiation, the presence of a thermal radiation field, and the flux from discrete stellar sources. This is expressed in Equations (21) and (22).

The next sections describe tests of the radiation transport routines and the reliability of our method.

To test the accuracy of the matter/radiation coupling, we set up a unit simulation volume with a uniform gas density of ρ = 10⁻⁷ g cm⁻³ initially out of equilibrium with the radiation field (Turner & Stone 2001). If the radiation energy dominates the total energy, then any radiation energy absorbed or emitted is relatively small and the radiation field can be said to be unchanging. The matter energy evolves according to Equation (16), but without the source term, i.e.,

$$\begin{equation} \frac{\partial e}{\partial t} = \chi c (E_r - a_R T^4), \end{equation} \tag{ 32 }$$

where e = ρ is the volumetric matter energy density, and χ = κ_Pρ = 4 × 10⁻⁸ cm⁻¹ is the absorption coefficient.

We set the radiation energy density E_r = 10¹² erg cm³, as in Turner & Stone (2001). We solve Equation (32) assuming a constant E_r and plot the results in Figure 2 as the solid black lines. The initial matter energy density is e = 10¹⁰ erg cm³ in the case where the gas cools to equilibrium. In the gas warming case, the initial matter energy density is e = 10² erg cm³.

Zoom In Zoom Out Reset image size

The gas is assumed to be an ideal gas with γ = 5/3 and a mean mass per particle of μ = 0.6.

Figure 2 shows the results of numerical simulations with FLASH. The simulation has an adaptive timestep initially set to Δt = 10⁻¹⁴ s. The maximum timestep size for this simulation was set to Δt = 10⁻⁸ s. The total energy in the simulation volume is conserved to better than 1% over the duration of the simulation.

The treatment of radiative shocks is an important benchmark in many radiative transfer codes (Hayes & Norman 2003; Whitehouse & Bate 2006; González et al. 2007; Kuiper et al. 2010b; Commerçon et al. 2011). We follow the setup described in Ensman (1994) of a streaming fluid impinging on a wall, represented by a reflective boundary condition. The fluid is compressed and a shock wave travels in the upstream direction. The hot fluid radiates thermally, and the radiation field preheats the incoming fluid. By varying the speed of the incoming fluid, sub- or supercritical shocks can be formed. Criticality occurs when there is sufficient upstream radiation flux that the preshock temperature is equal to the postshock temperature. The fluid velocity at which this occurs is called the critical velocity. Numerical simulations can be compared to analytic arguments by Mihalas & Mihalas (1984) for the gas temperature in various parts of the shock to check how well these shock features are being reproduced by the code.

The initial conditions are as follows: an ideal fluid (γ = 5/3) has a uniform mass density of ρ₀ = 7.78 × 10⁻¹⁰ g cm⁻³, a mean molecular weight μ = 1, and is at a uniform temperature of T₀ = 10 K. The domain size is L = 7 × 10¹⁰ cm.

For the subcritical shock, the fluid moves to the left with a speed v = 6 km s⁻¹. For the supercritical shock, v = 20 km s⁻¹. The fluid is given a uniform absorption coefficient of σ = 3.1 × 10⁻¹⁰ cm⁻¹.

Upon compression of the gas at the leftmost boundary, the postshock temperature T₂ increases and a radiative flux of the order of $\sigma _B T_2^4$ is produced, where σ_B is the Stefan-Boltzmann constant. This radiation penetrates upstream to preheat the gas ahead of the shock front to a temperature T₋. This preheating can be clearly seen in Figure 3. The shock is considered subcritical so long as T₋ < T₂.

Zoom In Zoom Out Reset image size

If the speed of the incoming gas v is increased, there is greater preheating and the preshock temperature approaches the postshock temperature, T₋ ∼ T₂. These temperatures are equal for a critical shock. If the incoming gas speed is increased still further, the temperature of the preshock gas remains steady, while the radiative precursor is extended. This is a supercritical shock, as seen in Figure 4.

Zoom In Zoom Out Reset image size

For the subcritical case, if T₋ ≪ T₂, an approximate solution is given in Mihalas & Mihalas (1984) for the postshock temperature,

$$\begin{equation} T_2 \approx \frac{2\left(\gamma - 1\right)v^2}{R\left(\gamma + 1\right)^2}, \end{equation} \tag{ 33 }$$

where R = k_B/μm_H is the ideal gas constant.

For our subcritical case, with an incoming gas velocity of v = 6 km s⁻¹, this gives T₂ ≈ 812 K. In our test case, the postshock temperature at the far-left of the domain is 706 K, or about 13% less than the approximate value.

The preshock temperature can be approximated (Mihalas & Mihalas 1984) by

$$\begin{equation} T_- \approx \frac{\gamma - 1}{\rho v R} \frac{2 \sigma _B T_2^4}{\sqrt{3}} \sim 279\, {\rm K.} \end{equation} \tag{ 34 }$$

The temperature spike at the subcritical shock front is approximated by

$$\begin{equation} T_+ \approx T_2 + \frac{3 - \gamma }{\gamma + 1}T_- \sim 874\, {\rm K.} \end{equation} \tag{ 35 }$$

Our preshock temperature reaches 336 K, which is about 20% warmer than the analytic approximation. Similar calculations in Kuiper et al. (2010b) did not reproduce any preheating due to the one-temperature (1T) approach in FLD. Our numerical shock temperature T₊ ≈ 871 K agrees to within better than 1%.

In the supercritical case, the radiation spike has collapsed to a thickness of less than about a photon mean free path. The temperature of the spike is T₊ ≈ 5420 K, whereas the analytic approximation (Mihalas & Mihalas 1984) gives

$$\begin{equation} T_+ \approx (3 - \gamma)T_2 \sim 4612\, {\rm K,} \end{equation} \tag{ 36 }$$

with which our numerical calculations agree to within 18%. Comparing to other FLD implementations on AMR grids, Commerçon et al. (2011) matches the analytic estimates of the subcritical postshock temperature a little closer (to within 2%). Our scheme captures the preshock heating that the 1T scheme described in Kuiper et al. (2010b) could not (makemake has since been made into a 2T scheme), but not as closely as Commerçon et al. (2011), which captures it to within 1%. The subcritical shock temperature we capture comparably well to Commerçon et al. (2011). While flash has a different refinement criterion from the one used in Commerçon et al. (2011), this is not likely the cause of the differences. More likely is the different choice of flux limiter (Levermore & Pomraning 1981 versus Minerbo 1978), which, because of different assumptions about the angular dependence of the radiation field in a particular problem, can yield slightly different solutions, with the greatest differences seen in regions of intermediate to low optical depth (Turner & Stone 2001). It is not obvious which flux limiter is best for a given problem, but we have opted for one shared by Kuiper et al. (2010b).

We compute the radiation force density for the stellar radiation component through the midplane of our disk as in Equation (31). A similar test was done by Kuiper et al. (2010b), although their cut was not done through the midplane but along a polar angle of θ ≈ 27°. flash computes an isotropic radiation pressure via the Eddington approximation (Equation (30)), so we estimate the thermal radiation force density via

$$\begin{equation} f_{\mathrm{rad,thermal}} = - \frac{d p_{\mathrm{rad}}}{dx}. \end{equation} \tag{ 42 }$$

We use the makemake code (revised to be a 2T solver) to perform the same calculation through the midplane (θ ≈ 0°) and compare to flash. The results are shown in Figure 9. Most of the radiation pressure is due to the direct stellar component, with the diffuse, reprocessed field adding only a tiny contribution that is more apparent far from the source. The outer boundary is treated as optically thin outflow in makemake. In flash, we set the "vacuum" outer boundary condition, dF/dx = −2F.

Zoom In Zoom Out Reset image size

In the inner 40 AU, the two codes match to within about 20%, despite the radiation force density varying over six orders of magnitude. Because the radiation force is proportional to the absorbed radiation flux, the peak in the force density lies at the inner edge of the disk around 1 AU. Beyond 40 AU, the relative difference between the two calculations grows on account of boundary conditions, but the radiation force density is negligible beyond this point (six orders of magnitude below the peak value).

The flash code and updated makemake code produce comparable measurements of the radiation pressure. The advantage and key innovation with flash, however, is its AMR-capability that allows it to solve more general problems with high accuracy.

The first multi-source test involves two sources separated by about 658 AU. Each source has a stellar radiation field with T_eff = 5800 K. The medium is uniform and homogeneous, with a density of ρ = 8.321 × 10⁻¹⁸ g cm⁻³. We perform this test to contrast our hybrid radiation transfer method with M1 moment methods (Levermore 1984; González et al. 2007; Vaytet et al. 2010). M1 moment methods are similar to FLD methods in that they are both based on taking angular moments of the radiative transfer equation. FLD takes only the zeroth-order moment, and the conservation equations are closed using a diffusion relation based on the gradient of the radiation energy. M1 methods take the first moment of the radiative transfer equation, and the closure of the conservation equations takes a form that preserves the bulk directionality of photon flows. Problems arise, however, when multiple sources are present. With two sources, photons flow in opposing directions, canceling the opposing components of their flow and introducing spurious flows in the direction perpendicular to the line between the two sources (Rosdahl et al. 2013). Raytracers do not suffer from this artifact because the irradiation at each grid location is calculated along sightlines to each of the sources present in the simulation.

Figure 10 shows the simulation setup as well as the temperature of the gas after it has been allowed to equilibrate with the radiation field. We measure the temperature of the gas along the vertical dashed line and plot the result in Figure 11.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 11 shows the temperature along the dashed line from Figure 10, the locus of points equidistant from both radiation sources. Here we compare the simulated temperatures of our hybrid scheme in flash against a Monte Carlo calculation performed with RADMC-3D⁶ (Dullemond 2012). The Monte Carlo calculation was done with 10⁹ photons using the "modified random walk" method with scattering neglected, on an uneven grid (101 × 1 × 101 cells for the x × y × z dimensions) using the same Draine & Lee (1984) dust model as in the rest of this paper.

The flash temperatures differ from the Monte Carlo results by less than 5%, with flash calculating slightly warmer temperatures directly between the sources and slightly cooler temperatures at the edges of the simulation volume. These small differences are most likely due to multifrequency effects.

Monte Carlo calculations are considered the benchmark for accuracy, but are extremely computationally intensive and cannot, in general, be used in dynamic calculations. Simulations of star formation cannot be addressed with full Monte Carlo methods using present-day machines and reasonable time constraints. Moreover, they are highly nontrivial to parallelize. In these cases, approximate schemes such as flux-limited diffusion must be used.

We next create a test setup inspired by Rijkhorst et al. (2006), where a central concentration was irradiated by two sources. Since FLD methods are based on local gradients in the radiation energy, they cannot properly treat shadowing. In this test we place a dense clump of material at the center of our simulation volume and irradiate it from two angles, creating a shadowed region. We compute the irradiation and equilibrium temperatures.

Our setup contains a mix of optically thin and optically thick regions, creating strong gradients in the radiation energy density. We compare the equilibrium temperatures computed by the hybrid scheme to the temperature computed without diffusion.

The two radiation sources have a temperature T_eff = 8000 K and are situated in a low-density (ρ = 10⁻²⁰ g cm⁻³) medium to the side and below a central density concentration (ρ = 10⁻¹⁷ g cm⁻³) with radius R_c = 4 × 10¹⁵ cm ≈267 AU. The side length of the simulation box is about 2000 AU. Hydrodynamics are disabled.

Figure 12 shows the simulation setup and the specific irradiation of the gas in a slice through the midplane of the simulation volume. The overlaid grid shows the flash AMR block structure, with each block containing 8³ cells and refining on gradients in density, radiation energy, and matter temperature using the flash error estimator described in Section 3.1. The simulation was run with a maximum eight levels of refinement, achieving a maximum resolution of 1.96 AU. The black circle in the center marks the central density concentration. The central density concentration casts a shadow on the side opposite of the central concentration.

Zoom In Zoom Out Reset image size

Figure 13 shows the gas temperature after the simulation has been given time to reach equilibrium. Here we see warming of a region just interior to the central overdensity where strong density gradients exist.

Zoom In Zoom Out Reset image size

To visualize the radiation flux, we plot a vector field of the radiation flux, given by Equation (11). This is shown in Figure 14 overplotted on the total radiation energy density. The vector field represents how the radiation energy is being diffused by the FLD solver. The two radiation sources irradiate the central overdensity, heating it. It then drives radiation back into the surrounding medium via diffusion. There is also a radiation energy gradient across the central overdensity. The black circle marks the central overdensity, as in the previous two figures. In Figure 14 we see a halo of radiation energy just outside the circle on the side nearest to the two radiation sources also resulting from the diffusion of radiation.

Zoom In Zoom Out Reset image size

We now compare this to the zero-diffusion case by performing the same calculation with the diffusion term set to zero. The gas equilibrates with the radiation field imposed by the sources, but this energy is not allowed to diffuse. Figure 15 shows the difference in the resultant equilibrium temperature as a percentage of the relative difference. The result is that the central overdensity is approximately 50% hotter. It cannot cool by diffusing its radiation energy into the surrounding medium.

Zoom In Zoom Out Reset image size

Because the re-emission terms are often left out of raytracing-methods for the sake of simplicity and computational cost, they are unsuitable and inaccurate in simulations containing optically thick material, although re-emission is often treated using a cooling function to extract the energy. In optically thick envelopes surrounding young stellar objects, radiation is reprocessed as the envelope of gas absorbs and re-emits the radiation in the infrared.

The development of this hybrid radiation transfer code benefited from helpful discussions with Romain Teyssier, Christoph Federrath, Klaus Weide, and Petros Tzeferacos, and from two visits by M.K. to the flash Center at the University of Chicago. M.K. acknowledges financial support from the National Sciences and Engineering Research Council (NSERC) of Canada through a Postgraduate Scholarship. R.K. acknowledges funding from the Max Planck Research Group Star Formation throughout the Milky Way Galaxy at the Max Planck Institute for Astronomy. R.E.P. is supported by an NSERC Discovery Grant. T.P. acknowledges financial support through a Forschungskredit of the University of Zürich, grant No. FK-13-112. The flash code was in part developed by the DOE-supported Alliances Center for Astrophysical Thermonuclear Flashes (ASCI) at the University of Chicago. This work was made possible by the facilities of the Shared Hierarchical Academic Research Computing Network (SHARCNET: www.sharcnet.ca) and Compute/Calcul Canada.

We are grateful to KITP, Santa Barbara, for supporting M.K. as an affiliate (for 3 weeks), and R.E.P. as an invited participant (for 2.5 months) in the program "Gravity's Loyal Opposition: The Physics of Star Formation Feedback," where some of this research was supported by the National Science Foundation under grant No. NSF PHY11-25915.

Much of the analysis and data visualization was performed using the yt toolkit⁷ by Turk et al. (2011).