Self-Consistent Solution of Cosmological Radiation-Hydrodynamics and Chemical Ionization

# Self-Consistent Solution of Cosmological Radiation-Hydrodynamics and Chemical Ionization

Daniel R. Reynolds John C. Hayes Pascal Paschos Michael L. Norman Mathematics, Southern Methodist University, Dallas, TX 75275-0156 Lawrence Livermore National Lab, PO Box 808, L-551, Livermore, CA 94551 Ctr. for Astrophysics and Space Sciences, U.C. San Diego, La Jolla, CA 92093 Physics Department, U.C. San Diego, La Jolla, CA 92093
###### Abstract

We consider a PDE system comprising compressible hydrodynamics, flux-limited diffusion radiation transport and chemical ionization kinetics in a cosmologically-expanding universe. Under an operator-split framework, the cosmological hydrodynamics equations are solved through the Piecewise Parabolic Method, as implemented in the Enzo community hydrodynamics code. The remainder of the model, including radiation transport, chemical ionization kinetics, and gas energy feedback, form a stiff coupled PDE system, which we solve using a fully-implicit inexact Newton approach, and which forms the crux of this paper. The inner linear Newton systems are solved using a Schur complement formulation, and employ a multigrid-preconditioned conjugate gradient solver for the inner Schur systems. We describe this approach and provide results on a suite of test problems, demonstrating its accuracy, robustness, and scalability to very large problems.

###### keywords:
Radiation, Hydrodynamics, Ionization, Cosmology, Numerical Methods, Implicit Methods
###### Pacs:
98.80.Bp, 02.60.Cb, 02.60.Lj
journal: Journal of Computational Physics

## 1 Introduction

A fundamental physics component in cosmology and astrophysics applications is the transport of ionizing radiation along with the interaction of that radiation with the hydrodynamic motion and ionization state of the surrounding gas. One example currently receiving a great deal of attention is cosmic reionization BarkanaLoeb2007 (), our motivation in this work. Observations indicate that an early population of UV emitting galaxies photoionized the intergalactic hydrogen and helium gas when the universe was about 800 million years old (redshift ). A computational challenge is to calculate this process self-consistently, coupling the radiative transfer of ionizing photons, the ionization kinetics and photo-heating of the gas, and the attendant hydrodynamic motions. This problem is challenging because the physics is numerically stiff and cosmic reionization is intrinsically three-dimensional, involving the growth, percolation, and overlap of ionization zones around an irregular and evolving distribution of galaxies with time-dependent luminosities. In addition, the problem inherits the large range of spatial scales () intrinsic in galaxy formation simulations, necessitating the use of spatially adaptive mesh or particle–tree methods and large-scale parallel computing NormanEtAl2007 (); Springel2005 ().

A variety of 3D radiative transfer methods are under development to tackle this problem IlievEtAl2006 (). These necessarily simplify the description of the radiation field in order to render the problem tractable. These methods include ray tracing using long and short characteristics MellemaEtAl2006 (); SusaUmemura2004 (); WhalenNorman2006 (); RijkhorstEtAl2006 (), Monte Carlo MaselliEtAl2003 (); SemelinEtAl2007 (), and moment methods GnedinAbel2001 (); HayesNorman2003 (); PaschosEtAl2007 (). However, only some of these codes allow solution of the coupled problem on spatially-adaptive grids, and very few allow distributed-memory parallelism. More importantly, all of these codes handle the interactions between hydrodynamic, radiative, and chemical processes in an explicit, operator-split fashion, which ignores stiff couplings that often arise between these components. When this happens, such codes must dramatically reduce allowable timesteps, or resort to subcycling, to ensure stability and accuracy of the coupled simulations.

Radiation transport and chemical ionization time scales are much faster than typical hydrodynamic time scales. This is particularly evident in dense gas bound to galaxies where recombination and light crossing times are short compared to the age of the universe (Hubble time). Moreover, these processes are very tightly coupled since radiation induces ionization that in turn affects opacities. While time-explicit, adaptive subcycling schemes have been developed that are capable of returning accurate solutions in all regimes of interest WhalenNorman2006 (), it is our view that for both computational efficiency and solution accuracy, tightly-coupled implicit methods require investigation. Here we present such a method.

We solve ionizing radiation transport, chemical ionization kinetics, and gas photo-heating using a fully implicit inexact Newton method. This algorithm is coupled to a cosmological hydrodynamics code through an explicit, operator-split formalism. The inner linear Newton systems are solved using a Schur complement formulation, which neatly decouples the local microphysics from the transport calculation. Radiation transport is modeled in the flux-limited diffusion approximation for simplicity, although our approach can be easily generalized to higher-order moment schemes. The use of Schur complements allows the application of optimal and scalable multigrid methods for the solution of the scalar radiation diffusion equation. We describe our algorithm in detail. We then illustrate our method’s accuracy, robustness, and parallel scalability against a suite of verification tests of increasing size and complexity. In its current implementation, we are restricted to calculations on uniform Cartesian grids. An extension of our algorithm on block structured adaptive meshes is under development, and requires only modifications to the inner multigrid linear solver.

The paper is organized as follows. In Sec. 2 the governing equations for cosmological radiation hydrodynamics are presented under two different assumptions about the radiation–matter coupling: a two temperature model assuming local thermodynamic equilibrium (LTE), and a non-LTE ionization kinetics multispecies model. In Sec. 3 we describe our solution procedures for splitting off the hydrodynamic calculation, and our coupled implicit radiation-ionization-gas energy system. Results from solution verification tests are presented in Sec. 4, as well as parallel scalability tests. Conclusions follow in Sec. 5.

We consider the coupled system of partial differential equations

 ∂tρb+1avb⋅∇ρb =−1aρb∇⋅vb, (1) ∂tvb+1a(vb⋅∇)vb =−˙aavb−1aρb∇p−1a∇ϕ, (2) ∂te+1avb⋅∇e =−2˙aae−1aρb∇⋅(pvb)−1avb⋅∇ϕ+G−Λ (3) ∂tni+1a∇⋅(nivb) =αi,jnenj−niΓphi,i=1,…,Ns (4) ∂tE+1a∇⋅(Evb) =∇⋅(D∇E)−m˙aaE+4πη−cκE. (5)

Here, the first three equations (1)-(3) correspond to the equations of ideal gas dynamics in a coordinate system that is comoving with the expanding universe BryanEtAl1995 (). These equations correspond to mass, momentum and energy conservation, respectively, in which is the proper peculiar baryonic velocity, is the proper pressure, and the total gas energy per unit mass is given by . The modified gravitational potential is given by , which satisfies the comoving form of Poisson’s equation,

 ∇2ϕ=4πga(ρb+ρdm−⟨ρ⟩), (6)

where provides the gravitational constant, and are the baryonic and dark matter densities, respectively, and is the cosmic mean density. The densities are comoving, relating to the proper densities through the relation . Here denotes the cosmological expansion parameter for a smooth homogeneous background, where the redshift is a function of time only. All spatial derivatives are taken with respect to the comoving position . The hydrodynamics equations are closed as usual with the ideal gas equation of state,

 e=pρb(γ−1)+12|vb|2, (7)

where is the ratio of specific heats, taken to be throughout this work.

The hydrodynamics equations are coupled with the elemental species rate equations (4) and an equation describing the flux-limited diffusion (FLD) approximation of radiation transport in a cosmological medium (5) HayesNorman2003 (); Paschos2005 (). In these equations denotes the chemical species (out of total), is the comoving number density, is the electron number density, corresponds to other ions that react with the species , and are the rate coefficients defining these interactions AbelEtAl1997 (); HuiGnedin1997 (); corresponds to the comoving radiation energy density. The parameter controls whether is monochromatic at a specified frequency (), or an integrated grey radiation energy density ().

The baryonic gas is coupled to collisionless dark matter solely through their self–consistent gravitational field via (6). The dark matter density is evolved using the Particle-Mesh method described in HockneyEastwood1988 (); NormanBryan1999 (); OSheaEtAl2004 (). As the N-body method is standard and not the focus of this paper, we do not elaborate on it here.

### 2.1 Model Coupling

In addition to the advective components of the chemistry and radiation equations, coupling between these equations arise through a number of spatially-local reaction terms. The radiation energy density and chemical number densities affect the gas energy through the heating and cooling rates and , respectively. The ionization and recombination rates and emissivity depend on the gas temperature,

 T=(γ−1)pμmpρbkb, (8)

where corresponds to the mass of a proton, corresponds to the local molecular weight, and is Boltzmann’s constant. Finally, the photoionization rate depends on the radiation energy density, and the opacity depends on the state of chemical ionization.

In determining these coupling terms we distinguish between two cases: those in local thermodynamic equilibrium (LTE) and those that are not (nLTE). In the LTE case the chemical species are assumed to be in equilibrium, and hence their equations (4) may be omitted from the time-dependent system (1)-(5). For problems in this regime, the coupling terms resemble those typically encountered in radiation-hydrodynamics simulations HayesNorman2003 (); HayesEtAl2006 (); HowellGreenough2003 (),

 ηLTE(T)=κPB=κPσSBπT4,GLTE(ρb,E)=cκρbE,ΛLTE(ρb,T)=4πρbηLTE(T), (9)

where is the speed of light, is the Stefan-Boltzmann constant, and and correspond to the problem-dependent Planck mean and total opacities for the gas.

For simulations that may not be approximated as being in local thermodynamic equilibrium, these coupling terms involve the dynamically-changing chemical ionization states. Here, the combined opacity depends on the local ionization states , the emissivity depends on both and , the gas heating rate depends on both and , and the gas cooling rate depends on and , with the corresponding formulas given in the references AbelEtAl1997 (); Black1981 (); Osterbrock1989 (); Cen1992 (); RazoumovEtAl2002 ().

### 2.2 Cosmological Flux-Limited Radiation Diffusion Model Details

We derive the cosmological flux-limited radiation diffusion equation (5) from the general multi-frequency version Paschos2005 (),

 ∂tEν+1a∇⋅(Eνvb)=∇⋅(D∇Eν)+ν˙aa∂νEν+4πην−cκνEν. (10)

Through assumption of a given radiation frequency spectrum, , the frequency-dependent radiation energy density may be written in the form . With this, we define the single “grey” radiation energy density used in the model (1)-(5) as

 E(x,t)=∫∞ν0Eν(x,t,ν)dν=~E(x,t)∫∞ν0χE(ν)dν. (11)

The radiation equation (5) may then be derived through integration of the equation (10) over frequencies ranging from the ionization threshold of Hydrogen ( eV) to infinity; integration of the term gives rise to the term in (5). We note that this approximation (11) is valid only if the assumed spectrum is defined such that the indefinite integral exists, as is the case for quasar and stellar type spectra where it scales with frequency as where . However, through this formulation we may also consider problems involving a monochromatic radiation energy density, since such energy densities may also be expanded as , where for radiation at the monochromatic frequency , the assumed spectrum is given through the Dirac-delta function . In such cases, the term vanishes, giving rise to the parameter in (5). For standard grey radiation approximations, we assume a radiation spectrum of the form of either a power law, , or as a blackbody spectrum, .

As is standard with FLD approximations to radiation transfer, one must take special care in construction of the diffusion coefficient function . In its simplest form, may be written as , where is the total extinction coefficient, corresponds to total absorption (here taken to be the opacity ) and corresponds to scattering HayesEtAl2006 (). Use of this form for the diffusion coefficient, however, results in an infinite signal speed for the radiative flux . To preserve causality, the analytic form of is modified with a dimensionless flux-limiter whose particular form may be tuned to the specific problem of interest, but whose overriding purpose is to guarantee that the radiation transfer equation (5) gives the correct numerical behavior in the limiting cases of (nearly) isotropic and free-streaming radiation. Several choices for flux-limited forms of have been proposed in the literature Levermore1984 (); LevermorePomraning1981 (). We consider the diffusion coefficient to be of the form

 D(E)=⎡⎢⎣D1(E)000D2(E)000D3(E)⎤⎥⎦,whereDi(E)=c(2κT+Ri)6κ2T+3κTRi+R2i, (12)

with . We note that this function has been reformulated from its original version HayesNorman2003 () to allow increased numerical precision for scattering-free simulations involving extremely small opacities (i.e. ), as is typical in cosmology applications.

## 3 Solving the Coupled System

### 3.1 Operator-Split Hydrodynamics with Radiative Feedback

Since typical astrophysical and cosmological simulations involve the hydrodynamic motion of gases encountering shocks, whereas radiation diffusion and chemical kinetics processes are more of reaction-diffusion type, we choose to solve the coupled system (1)-(5) in an operator-split fashion. In this approach, a time step to is taken using the general steps

• Deposit the dark matter particles onto the mesh to calculate the .

• Solve for the gravitational potential resulting from the densities and using equation (6).

• Evolve the dark matter particles using the Particle Mesh Method HockneyEastwood1988 (); NormanBryan1999 (); OSheaEtAl2004 ().

• Evolve the hydrodynamics equations (1)-(3) with a high-order, explicit-time upwind method. In this step, use the velocity field to advect both the chemical number densities and radiation energy density .

• Using a high-order implicit-time method, solve a coupled reaction-diffusion system to obtain the time-evolved number densities , radiation energy density and gas energy .

In order to allow us to split the equations (1)-(5) into the two steps (iv) and (v) above, we consider the gas energy as consisting of two components, , where is the fluid energy arising from the hydrodynamic evolution of the system, and is the gas energy correction arising from the couplings with radiation and chemistry. Under this decomposition, the energy conservation equation (3) may be equivalently written as

 ∂t(eh+ec)+1avb⋅∇(eh+ec)= (13) −2˙aa(eh+ec)−1aρb∇⋅(pvb)−1avb⋅∇ϕ+G−Λ.

Under this splitting, the hydrodynamic solver used in step (iv) of the operator-split algorithm solves the system of equations

 ∂tρb+1avb⋅∇ρb =−1aρb∇⋅vb, (14) ∂tvb+1a(vb⋅∇)vb =−˙aavb−1aρb∇p−1a∇ϕ, (15) ∂teh+1avb⋅∇eh =−2˙aaeh−1aρb∇⋅(pvb)−1avb⋅∇ϕ (16) ∂tni+1a∇⋅(nivb) =0, (17) ∂tE+1a∇⋅(Evb) =0, (18)

to evolve the solution at , , to the time-updated variables at , , and the advected variables . For this step, we employ the Piecewise Parabolic Method (PPM) ColellaWoodward1984 (), on a regular finite-volume spatial grid, implemented in the community astrophysics code Enzo NormanEtAl2007 (); OSheaEtAl2004 (); enzo-site ().

The remainder of the coupled system,

 ∂tec =−2˙aaec+G−Λ, (19) ∂tni =αi,jnenj−niΓphi, (20) ∂tE =∇⋅(D∇E)−m˙aaE+4πη−cκE, (21)

is then solved using a fully implicit nonlinear solution approach to evolve the advected variables to the time-evolved quantities . Here we may assume that since the hydrodynamic solver uses the full energy at in its evolution, i.e. . Once this step is finished, we compute the time-evolved total energy as the sum of the hydrodynamic portion and the adjustments due to radiative feedback , i.e. . The treatment of the implicit radiation, chemical ionization, and gas energy feedback system (19)-(21) serves as our focus for the remainder of this section.

### 3.2 Solving the Radiation, Ionization and Energy Feedback System

Under a method-of-lines approach, we consider a two level, up-to-second order accurate theta-scheme for implicit integration of our system (19)-(21),

 en+1c+ΔtθLn+1e =enc+Δt(θ−1)Lne, (22) nn+1i+ΔtθLn+1ni =nni+Δt(θ−1)Lnni, (23) En+1+Δtθ[Dn+1E+Ln+1E] =En+Δt(θ−1)[DnE+LnE]. (24)

Here, the parameter defines the implicit integration method: =1 corresponds to a first-order implicit Euler method, =0.5 gives a second-order time-centered approach (i.e. Crank-Nicolson). We note that in the ensuing computational results from section 4, we have typically taken =0.51 to provide a nearly-second-order time integration while avoiding the “ringing” traditionally associated with fully time-centered approaches SimoArmero1994 (); Osterby2003 (). For the above equations, we have defined the diffusive operator

 DE=DE(E,ni) ≡−∇⋅(D∇E), (25)

and we have defined the local “reaction” operators as

 Le=Le(ec,E,ni) ≡2˙aaec−G+Λ (26) Lni=Lni(ni,ec,E) ≡niΓphi−αi,jnenj (27) LE=LE(E,ec,ni) ≡m˙aaE−4πη+ckE. (28)

The equations (22), (23) and (24) form a coupled nonlinear system of reaction-diffusion equations for evolution of the fluid energy correction , the elemental number densities , and the radiation energy density . Denoting the vector of unknowns , we first define the nonlinear residual function for the time step , as

 (29)

where the vectors are formed using the previous time-level information from (22)-(24). In order to evolve the coupled implicit system, we solve the nonlinear problem for the updated vector of unknowns . For this nonlinear solve, we use a globalized Inexact Newton’s Method Kelley1995 (); KnollKeyes2004 (), in which we apply an iterative process for convergence toward the solution in the following manner.

Given an initial guess , we iterate toward a solution satisfying (we typically choose ):

1. Approximately solve the linearized Newton system, , to tolerance for the correction vector . Here, , and we typically choose the tolerance as .

2. Update the vector of unknowns as , where is the line-search parameter DennisSchnabel1996 (); OrtegaRheinboldt2000 ().

We measure convergence of the Newton iteration with the RMS norm

 ∥v∥=(∥v∥22N(Ns+2))1/2, (30)

where is the number of unknowns in ( spatial cells, variables), since such a norm does not grow artificially larger with mesh refinement. The key to efficiency of the inexact Newton algorithm lies in a fast and robust solver for the linear systems . Once such a solver has been provided, the algorithm exhibits very fast convergence – superlinear for this choice of Kelley1995 (); EisenstatWalker1996 (). Moreover, for diffusive PDE systems similar to the one solved here, the Newton convergence rate has been shown to be independent of spatial resolution WeiserEtAl2005 (), suggesting that this entire implicit algorithm should allow scalability to the limits of the inner linear solver.

#### 3.2.1 Linear Solver

In solving the system (29), we make one approximation within the Newton system matrices , wherein we lag the dependence of in (25) to the previous Newton iterate. Mathematically, this results in a full Newton step for all but the limiter’s dependence on the chemical opacities, which are instead converged through a fixed-point iteration. The resulting solution retains the accuracy and stability of the full Newton iteration, albeit with theoretically slower convergence. However, in practice we have not noticed any increase in nonlinear iterations due to this approximation, and most importantly it results in inexact Newton matrices with the form

 J(U)=I+Δtθ⎡⎢⎣Je,eJe,nJe,EJn,eJn,nJn,EJE,eJE,nJE,E⎤⎥⎦, (31)

where nearly all of the blocks are given by the spatially-local components,

 (32)

and the only block containing spatial couplings is . Thus, although the Jacobian matrix contains couplings both within and between variables, it has a very desirable structure: all inter-variable couplings occur locally in space, and the only nonlocal couplings are within the block , consisting of a scalar-valued reaction-diffusion operator.

In keeping with a block-structured view of the Jacobian (31), we rewrite the Newton system in the form

 [MULD](sMsE)=−(fMfE), (33)

where

 M=I+Δtθ[Je,eJe,nJn,eJn,n],U=Δtθ[Je,EJn,E],L=Δtθ[JE,eJE,n],D=I+Δtθ[JE,E], (34)

and . We note that the only matrix containing spatial dependencies is , so under an appropriate variable ordering the other sub-matrices are block diagonal. Hence, we may efficiently invert to obtain as a function of :

 MsM+UsE=−fM⇒sM=−M−1(fM+UsE).

Inserting this into the second row, we have the single equation for ,

 (D−LM−1U)sE=−fE+LM−1fM.

Therefore, this Schur complement formulation BrownWoodward2001 () for solution of the linear Newton system (33) proceeds with the following steps:

1. Set .

2. Solve the system for .

3. Recover the remaining solution pieces, .

We examine each of these steps below.

The step (i) corresponds to solving the linear system

 [I+ΔtθJe,eΔtθJe,nΔtθJn,eI+ΔtθJn,n](~fe~fn)=(fefn). (35)

Due to the spatial locality of each component in , we order the equations and unknowns in this system so that application of , and more notably , may be performed independently in every spatial cell. Such solves consist of dense matrix algebra on linear systems (for chemical densities). In addition, at this step we compute the matrix through one additional solve with , which will be used in the following steps.

The step (ii) corresponds to solving the system . This is denoted the Schur complement system, with the matrix . We note that due to the spatially-local nature of the matrices and , we may form by constructing the diffusive sub-matrix , followed by updates to the diagonal entries corresponding to the entries of . Similarly, construction of the right-hand side may occur independently at each spatial location. Once this system has been computed, we use a multigrid-preconditioned conjugate gradient parallel linear solver from the HYPRE library FalgoutYang2002 (); hypre-site () to perform the scalar-valued solve, . We note that this is the only step in the solution of the Jacobian systems (33) that requires communication between processors. Moreover, we point out that in recent tests the HYPRE library has demonstrated ideal weak scaling up to over 100,000 processors for diffusion problems similar to the one encountered in this work BakerFalgoutYang2006 (). As this solver comprises the majority of the non-local components within our nonlinear solver, we therefore expect similar scalability for the overall implicit solution approach described here.

The final step (iii) in solution of the system (33) is to recover the solution components via the system . Again, since we have already computed the spatially-local matrix and the vector in step (i), we may trivially obtain the remaining solution components through cell-local matrix-vector products and vector operations, .

#### 3.2.2 Multiphysics/Cosmology Units

As with any multi-physics system, special care must be taken when solving such systems computationally due to disparate scales between variables and equations. This problem is especially evident in cosmology applications, where in CGS units one may typically enounter specific gas energies on the order of , number densities on the order of , and radiation energy densities on the order of , with all proper density values decreasing in time due to cosmological expansion. To this end, we define the scaled variables

 ~ec =ec/ue,~E=E/uE,~ni=ni/un, (36) ~x =x/ux,~t=t/ut,

where the constants , , , and correspond to the typical magnitudes of gas energy, radiation energy density, chemical number density, length and time at the start of the simulation. We note that due to our use of comoving values for , and , these constants are all redshift-independent, with the proper values of these quantities given by

 Eproper =E/a3(t)=~EuEa3(t), ni,proper =ni/a3(t)=~niuna3(t), (37) xproper =xa(t)=~xuxa(t).

The constants are supplied on a problem-dependent basis, to allow for adaptable, on-the-fly rescaling of simulations ranging from normalized test problems to cosmological reionization. With these rescaled variables, we rewrite our equations (22)-(24) as the normalized system

 ~en+1c+Δ~tθ~Ln+1e =~enc+Δ~t(θ−1)~Lne, (38) ~nn+1i+Δ~tθ~Ln+1n =~nni+Δ~t(θ−1)~Lnn, (39) ~En+1+Δ~tθ[~Dn+1E+~Ln+1E] =~En+Δ~t(θ−1)[~DnE+~LnE]. (40)

Here the operators , , and have correspondingly absorbed the renormalization constants . These equations, along with the normalized solution vector are then used within the solution strategy described in section 3.2.

#### 3.2.3 Adaptive Time Step Selection

A strong appeal of using implicit methods is their stability with respect to time step size; however such freedom gives rise to the questionof what time step should be used. At one extreme, we may choose a large step to achieve overall efficiency of the simulation, with little to no knowledge of the resulting temporal accuracy. At the other extreme, we may choose a very small time step for an accurate solution, resulting in inefficient simulations due to the increased cost of solving the nonlinear systems at each step. As the approach described here is operator-split, in which the hydrodynamics is solved using an explicit approach, we are therefore bound by the hydrodynamic CFL stability limit; however for most problems involving radiation and chemical ionization, the dynamic time scales of interest remain significantly faster than the hydrodynamic time scale. Thus the question of how to adaptively choose the time step size remains.

To that end, we adaptively choose the time steps as the largest possible that additionally satisfy a prescribed accuracy requirement. We estimate this accuracy through comparison of the updated solution with an explicit predictor for that solution . Defining the weighting vector in a spatial cell for the variable as

 ωi,v=√|Un+1i,vUpredi,v|+1, (41)

(which assumes normalized values of ), we estimate the local accuracy of the current time step as

 εloc=(1N(Ns+2)∥∥∥Un+1−Upredω∥∥∥pp)1/p, (42)

where we have used the standard -norm (including as the ‘max’ norm, in which case we do not divide by ), and where the quotient inside the norm is taken pointwise. With this estimate, we set the new time step to

 Δtn+1=τtolΔtnεloc, (43)

which should provide the maximal value that still satisfies the desired integration accuracy tolerance, , assuming that approximates the time-evolved solution to .

Here, the vector is designed so that estimates the average relative change in each solution component, and includes the harmonic mean of the predicted and new states to allow increased robustness in the case of cosmology-type problems where variables change by orders of magnitude across cells and time steps. The value of is typically taken to be in the ensuing test problems; however such a choice may limit parallel scaling since such a measure is sensitive to pointwise changes, of which there are many more as dynamics propagate throughout an increasingly refined domain. Lastly, we use the explicit predictor as the initial guess for the Newton method, , which we describe in the following section 3.2.4.

#### 3.2.4 Explicit Predictor

A well known property of Newton’s method is that its robustness and efficiency benefit greatly from an accurate initial guess. To this end, we provide the predicted initial Newton iterate

 U0=⎛⎜⎝ec,0ni,0E0⎞⎟⎠=⎛⎜⎝encnniEn⎞⎟⎠+Δt⎛⎜⎝LneLnniDnE+LnE⎞⎟⎠, (44)

i.e. we use an initial guess given by the -accurate explicit Euler update to the coupled system (19)-(21). As this provides only an initial guess to the solution, its instability at larger will not affect the temporal stability of the overall method, since the solution to each step must satisfy the implicit system (22)-(24). However, as we use an adaptive time-stepping strategy, for very fast dynamics (that give rise to very small ), such an initial guess may already satisfy the nonlinear tolerance and the solver will not require any Newton iterations, effectively allowing an adaptive explicit/implicit simulation of the coupled system (19)-(21).

The final detail that we describe in this solution approach relates to the choice of assumed radiation spectrum . As noted in section 2.2, we may choose either a monochromatic or an integrated “grey” radiation equation, based on the choice of this assumed spectrum. This choice affects all terms involving the radiation energy density in the general radiation energy equation (10) and in the coupling terms and . As each of these terms involve a product of the form , integration over converts these to

 ∫∞ν0f(ν)Eν(x,t,ν)dν=~E(x,t)∫∞ν0f(ν)χE(ν)dν. (45)

We therefore allow a user-defined functional form for , which we then numerically integrate to high accuracy upon initialization of the simulation, providing the relevant constants necessary to convert the -dependent equation (10) to the monochromatic or grey integrated equation (5).

## 4 Numerical Results

We present test problems designed to verify the accuracy of the radiation diffusion and chemical ionization modules in conjunction with hydrodynamical fluid motions. Since a number of distinct processes may compete for importance in a full simulation, we begin with simple tests that isolate various single components, and subsequently build upon those results with more sophisticated problems that couple additional physics. We begin with a radiation diffusion problem (§4.1) that exercises the diffusion term of the radiation equation (5) in the absence of energy, chemical, or hydrodynamic coupling; this test is followed by an examination of the matter-radiation coupling terms (§4.2) in an infinite uniform medium, a diffusion wave with material coupling (§4.3), and a non-equilibrium radiating shock problem (§4.4), all of which assume chemical equilibrium. Our attention then turns to problems including ionization, beginning with an HII ionization front propagating through a static, isothermal medium (§4.5), followed by a cosmological I-front propagation problem that exercises the cosmology terms and units (§4.6). We then consider a fully-coupled radiation-hydrodynamics-ionization calculation (§4.7). This section concludes with additional calculations (§4.8) demonstrating the parallel scalability of the radiation diffusion module, which of all components places the highest communication demands upon a domain-decomposed parallel calculation.

We note that for all test problems except (§4.6) and (§4.8), we use a non-cosmological problem (i.e.  and ). In problems (§4.6) and (§4.8), the cosmological parameters are described therein.

Because the standard diffusion equation is parabolic, the associated signal speed of the diffusion variable is formally infinite. However in reality radiation fronts propagate at speeds bounded by the speed of light in vacuum, so we modify the diffusion coefficient in our radiation energy equation (5) with a flux-limiter, as discussed in section 2.2. Our first test problem verifies the correct action of this limiter by examining the propagation of a planar radiation front through a transparent medium. Radiation is assumed to propagate along the -axis of our computational mesh; a Dirichlet boundary condition is imposed on the left boundary specifying an incident radiation energy density of 1.0 erg cm. Physically, the expectation is that with a sufficiently small (but nonzero, due to numerical constraints) opacity, a sharp radiation front will move through the domain at the speed of light. The Planck and Rosseland mean opacities are assigned a constant value of cm, ensuring an essentially transparent medium. The spatially uniform initial value of the radiation energy density is 10 erg cm.

The computational mesh has a domain length of 1.0 cm along the propagation direction of the light wave. We have run the problem for 8.3391 picoseconds, which is one quarter of the light-crossing time for this length. Figure 1 shows a series of curves resulting from calculations at mesh sizes of 128, 256, 512, 1024, and 2048 zones along the -axis. The dashed line indicates the expected location of the radiation front, , where is the speed of light and the evolution time. In the absence of the flux limiter, the numerical curves would give the formal solution for the diffusion equation, which for our problem parameters would be a nearly horizontal profile for throughout the domain (given the nearly zero opacity). That the curves capture the correct location of the radiation front is due entirely to the action of the limiter. The sharpening of this front with increased resolution is evident.

Also apparent is a slight lag (about 0.01 cm) between the analytical location of the radiation front and the numerical location, taken as the common intersection point of the numerical curves. The size of this lag depends upon the choice of the adaptive time step. In the language of (§3.2.3), we compute using and we vary , which here corresponds to the maximum allowed fractional change in the radiation energy density per timestep.

Figure 2 illustrates this effect by showing 3 curve pairs, each of which has been computed with a different choice of : 0.01 (green), 0.05 (blue), and 0.25 (red). In each colored pair, the solid curve was obtained using 512 zones, and the dashed curve shows the 128-zone result. Two curves at different mesh resolutions are provided to identify the “consensus” value of the light front location via their point of intersection. This location is seen to converge as .

### 4.2 Matter-Radiation Equilibration in a Homogeneous Medium

We now consider a problem in which and are spatially uniform but are initialized to values far away from equilibrium. This case thus isolates the matter-radiation coupling terms in the gas and radiation energy equations. The parameters for this test were published by Turner and Stone TurnerStone2001 (), who assumed an isotropic medium characterized by a single opacity of cm, a gas density of 10 g cm, and an average particle mass of 0.6 m, where m is the mass of a hydrogen atom. Coupled to this medium is a radiation field with a uniform value of 10 erg cm. From this value we compute a “radiation temperature”, , of about K. Here we have defined as the radiation constant, erg cm K. Two cases are considered: one in which the initial gas energy density is erg cm, and one in which the initial value is erg cm. For the stated parameters, these energies correspond to gas temperatures of roughly and 4.8 K, respectively, which therefore bracket the radiation temperature. In both cases, however, the initial radiation temperature is sufficiently high that the radiation energy density should remain constant to good approximation as the gas evolves to thermal equilibrium. To see this clearly, consider the effective heat capacity of a unit volume of the radiation “gas” as compared to that for the material. For radiation, this number is simply 1.0 cm , which evaluates to roughly erg K. In contrast, the gamma value and mean particle mass translate to a specific heat (erg g K) of roughly , which yields, for our assumed density, a material heat capacity of 20 erg K for a unit control volume. The physical result is that the radiation field has an effectively infinite thermal reservoir when compared to the material.

If the radiation energy density is formally assumed to be constant, the gas energy equation may be written as a simple ODE:

 ˙e=cκE−4πκB(e), (46)

where is the temperature-dependent Planck function

 B(T)=car4πT4. (47)

Using the ideal-gas law (7), we write as a function of and solve the simplified gas energy equation for the equilibrium value of such that ,

 eeq=32(ρkB0.6mH)(Ear)1/4. (48)

Notice that this expression is nothing more than the ideal gas formula for evaluated at the fixed radiation temperature, the expected result.

The results of our two test calculations are shown in figure 3. Both tests were run in a small box domain (4 zones) with triply-periodic boundaries. The case of is indicated by the solid curve; the low- case is shown by the dashed line. The horizontal dotted line has been placed at . Both energy curves converge to the correct result. Note that in this test the opacity serves only to control the timescale to reach thermal equilibrium; neither the value of the equilibrium energy nor the validity of our assumption of constant are dependent on the value of . While figure 3 demonstrates convergence to the correct asymptotic value of the gas energy, it provides no information as to whether the rate at which it approaches this value is correct. A quantitative assessment of this latter metric is provided by our next test problem.

We further use this test to examine the conservation properties of the coupled radiation and gas energy solver. In Table 1 we show the value of

 ∫|Etotal(t)−Etotal(0)|dx∫Etotal(0)dx (49)

for both tests at the final time =2.5e-7 sec., run using a variety of nonlinear and linear solver tolerances, and , respectively. We note that in all cases, the total energy is conserved to more than 10 digits of accuracy. Moreover, while the conservation is weakly dependent on the nonlinear solver tolerance, it is entirely independent of the linear solver tolerance. This behavior is most likely due to use of the Schur complement formulation (§3.2.1), that exactly solves for coupling between variables to floating point roundoff, leaving the iterative linear solver to handle only the radiation equation.

We further note that this is an ideal problem to test conservation of the coupled solver, since it is the only test considered that uses a closed system. We further comment that since the PPM finite-volume method is constructed to satisfy conservation, and the implicit subsolver achieves conservation to high accuracy, overall conservation of the coupled solver follows. However, we note that for problems utilizing non-periodic boundary conditions, chemical ionization cooling, gravitational heating, or cosmological expansion, the model no longer represents a closed system and therefore will not conserve energy.

### 4.3 Non-Equilibrium Marshak Waves

This test exercises both radiation diffusion and the physics of matter-radiation coupling. Non-equilibrium Marshak waves characterize the evolution of the radiation field in an initially cold, uniform halfspace on which a radiation source is imposed. The particular form of the Marshak problem described here is originally due to Pomraning Pomraning1979 (). The problem was re-examined by Su and Olson SuOlson1996 (), who derived semi-analytic exact solutions for the radiation and gas energy densities and tabulated select values of them on a grid of space and time values. The problem considers a formally 1-D semi-infinite domain in which denote dimensional space and time coordinates. Ignoring hydrodynamic motions, Su and Olson write simplified forms for the radiation and material energy equations as

 ∂tE(z,t)−∂z[c3κ∂zE(z,t)] =cκ[arT4(z,t)−E(z,t)], (50) cv(T)∂tT(z,t) =cκ[E(z,t)−arT4(z,t)], (51)

in which is the constant opacity, the radiation constant as defined previously, and the specific heat of the material. Note that the flux-divergence term in (50) assumes pure diffusion with no flux limiter. The matter temperature, , is assumed to be related to the gas energy, , via

 e = ρcvT. (52)

As written, (50) and (51) are coupled nonlinear PDEs in the dependent variables and . Su and Olson linearized the equations by choosing the following form for the specific heat:

 cv = αT3, (53)

where is an arbitrary constant. The dependence of on temperature has two effects: it allows equations (50) and (51) to be written as linear ODEs in and , and it gives the heat capacity of the material the same temperature dependence as the effective heat capacity of the radiation field. With an appropriate choice of , a problem can therefore be designed in which the material and radiation will both evolve significantly in space and time.

The description of the problem is completed with a specification of the boundary conditions. A Marshak boundary condition is applied to at :

 E(0,t) − 23κ∂zE(0,t) = 4cFinc, (54)

where is the incident flux at . The boundary condition at , and initial conditions at are

 E(∞,t) = 0,E(z,0) = T(z,0) = 0. (55)

Su and Olson construct the linearized equations by defining dimensionless independent and dependent variables, and such that

 X ≡ zκ√3, u(X,τ) ≡ (c4Finc)E(z,t), (56) τ ≡ (4arcκα)t, v(X,τ) ≡ (car4Finc)T4(z,t).

With these definitions, and letting , equations (50), (51), and (54)-(55) become

 ϵ∂τu(X,τ) − ∂2X2u(X,τ) = v(X,τ) − u(X,τ), (57) ∂τv(X,τ) = u(X,τ) − v(X,τ), (58) u(0,τ) − 2√3∂Xu(0,τ) = 1, (59) u(∞,τ) = 0, (60) u(X,0) = v(X,0) = 0. (61)

The Marshak boundary condition represented by (59) enforces the constraint of constant flux on the left boundary. This is an example of a “mixed” or Robin boundary condition, and as such requires special treatment in Enzo. For the purposes of this verification test, we implement this boundary condition by imposing a Dirichlet condition with a time-varying value of computed from SuOlson1996 ()’s equation 36, evaluated at . Because the integrands in their equation are highly oscillatory for , we substitute the asymptotic expression given by their equation 51 when .

Figure 4 shows results from a high-resolution simulation with 2048 zones along the coordinate. The exact solution values tabulated by SuOlson1996 () span the range . Since the right boundary condition is specified at , we choose our domain such that is sufficiently large (about 35) for the evolution time of interest that the Dirichlet condition may be reasonably applied. We choose opacity and coupling parameters 1.0 cm. The curves indicate profiles of (dashed lines) and (solid lines) for values of 1, 3, 10, 30, 100. The squares and circles indicate exact values of and , respectively. We have indicated these values on corresponding curves at evolution times sufficiently early that the material and radiation have not yet had time to equilibrate. Figure 5 shows a resolution study for the curves computed at . Curves at mesh sizes of 128 (red), 256 (orange), 512 (blue), 1024 (green), and 2048 (black) are shown. Each calculation is performed with the timestep restriction . Because this treatment allows for adaptive timesteps, the evident first-order rate of convergence measures the combined effect of time and space discretization methods.

### 4.4 Subcritical Radiating Shock Waves

We now add hydrodynamic motions to our mix of physics by examining the propagation of shock waves for which the radiation energy plays a significant role in the shock structure and evolution. Radiating shock waves represent a broad class of phenomena figuring prominently in both astrophysical and terrestrial applications. The particular formulation of the problem we present is due to Lowrie and Edwards LowrieEdwards2008 (), who considered the propagation of planar, steady shock waves in the grey nonequilibrium diffusion limit. Under the assumption of steady flow, LowrieEdwards2008 () transform the coupled gas and radiation energy equations into a set of nonlinear ODEs in dimensionless gas and radiation temperature variables, which must be integrated numerically to achieve semi-analytic solutions. Nonetheless, their radiation diffusion model corresponds identically to that implemented in Enzo in the grey LTE limit, and the unique structure of the post-shock material temperature profile for a given Mach number makes this problem an excellent verification test for computer codes.

We have run the Mach-2 test case described in LowrieEdwards2008 (). The computational domain has a length of 0.1 cm. The material has a uniform initial density of 1.0 g cm, a constant specific heat of erg g eV, and a uniform initial velocity of cm s. The material and radiation are assumed to be in thermal equilibrium at at a temperature of 121.6 eV. Outflow and reflecting boundary conditions are imposed upon the left and right boundaries, respectively, resulting in a shock wave that forms near the right boundary and propagates to the left. The total evolution time is 1.73325 nanoseconds.

Figure 6 shows the result of a high-resolution simulation (4096 zones along the propagation axis). The curves represent the dimensionless gas (solid curve) and radiation (dashed curve) temperatures, and . The circles and squares are taken from exact solution data kindly provided by R. Lowrie for this parameter set. Both the gas and radiation have dimensionless far-field temperatures of 1.0 in the pre-shock state. Examining the gas temperature curve, there are three features of particular interest: the precursor, in which the material is preheated ahead of the shock front by the radiation wave which travels ahead of the shock; the Zel’dovich spike, shown by the overshoot in temperature at the shock front, and the radiation relaxation region, delineated by the decline in the material temperature to its eventual far-field postshock value. Letting denote the maximum preshock value of the gas temperature in the precursor, and the asymptotic postshock value, we note that the property identifies this calculation as an example of a subcritical radiating shock. In the limit of high Mach number, can become equal to (but never exceed) , such a shock wave is referred to as supercritical.

As vividly demonstrated by LowrieEdwards2008 (), the strength of the precursor, the height of the Zel’dovich spike, and the precise temperature structure in the relaxtion region are extremely sensitive to the Mach number. While the case we have shown is subcritical, it lies near the limit for which a multidimensional code can reasonably capture this structure without resorting to adaptive mesh refinement. The degree to which we resolve this structure as a function of resolution is shown in figure 7, in which we magnify the region near the shock and show gas temperature curves for mesh sizes of 128, 256, 512, 1024, 2048 and 4096 zones. As shown in LowrieEdwards2008 (), raising the Mach number results in a dramatic increase in the height of the spike and narrowing of the relaxation region; a proper representation of the postshock structure in a supercritical shock with Enzo must await the implementation of adaptive mesh refinement in our radiation module.

Since this problem considers coupled radiation and hydrodynamics, we also examine how the adaptive time step selection strategy from Section 3.2.3 compares with the hydrodynamic CFL-limited time step. For this problem, the average radiation time step ranged from 1.1e-5 down to 9.4e-5 for the coarsest (128-cell) to finest (4096-cell) grids, exhibiting a near-constant time step selection that tracks evolution of the radiation field. For these same problems, the hydrodynamic CFL limits on the time steps were 3.3e-4, 1.7e-4, 4.2e-5, 2.1e-5, 1.0e-5 and 5.0e-6. Hence, for most problems the stiff radiation time scale limits the overall time step size, until very fine grids where the mesh-dependent hydrodynamic CFL stability condition becomes more restrictive.

### 4.5 Isothermal Ionization of a Static Neutral Hydrogen Region

Our first test problem incorporating ionization chemistry is due to Iliev et al. IlievEtAl2006 (). This problem combines radiative transfer and hydrogen ionization in a static astrophysical region. The physical situation of interest is the expansion of an ionized hydrogen (HII) region in a uniform gas around a single monochromatic ionizing source emitting photons per second at the ionization frequency of hydrogen ( eV). We enforce a fixed gas temperature of K, and a static hydrodynamic state (i.e. ). In such a problem, the radiation source should rapidly ionize the surrounding hydrogen, and then should develop a spherically-propagating ionization front (I-front) that propages quickly at first, slows, and then eventually stagnates at an equilibrium position referred to as the Strömgren radius, where ionization (HIHII) and recombinations (HIIHI) balance. For this scenario, the analytically-provided I-front radius is given by

 rI =rS[1−exp(−t/trec)]1/3,where (62) rS =[3˙