Smoothed Particle Radiation Hydrodynamics: Two-Moment method with Local Eddington Tensor Closure
MMon. Not. R. Astron. Soc. , 000–000 (0000) Printed 18 February 2021 (MN L A TEX style file v2.2)
Smoothed Particle Radiation Hydrodynamics: Two-Momentmethod with Local Eddington Tensor Closure
T. K. Chan (cid:63) , Tom Theuns , Richard Bower and Carlos Frenk
Institute for Computational Cosmology, Department of Physics, Durham University, South Road, Durham DH1 3LE, UK
18 February 2021
ABSTRACT
We present a new radiative transfer method (
SPH-M1RT ) that is coupled dynamicallywith smoothed particle hydrodynamics (
SPH ). We implement it in the (tasked-basedparallel)
SWIFT galaxy simulation code but it can be straightforwardly implementedin other
SPH codes. Our moment-based method simultaneously solves the radiationenergy and flux equations in
SPH , making it adaptive in space and time. We mod-ify the M1 closure relation to stabilize radiation fronts in the optically thin limitwhich performs well even in the case of head-on beam collisions. We also introduceanisotropic artificial viscosity and high-order artificial diffusion schemes, which allowthe code to handle radiation transport accurately in both the optically thin and opti-cally thick regimes. Non-equilibrium thermochemistry is solved using a semi-implicitsubcycling technique. The computational cost of our method is independent of thenumber of sources and can be lowered using the reduced speed of light approximation.We demonstrate the robustness of our method by applying it to a set of standard testsfrom the cosmological radiative transfer comparison project of Iliev et al. The SPH-M1RT scheme is well-suited for modelling situations in which numerous sources emitionising radiation, such as cosmological simulations of galaxy formation or simulationsof the interstellar medium.
Key words:
Physical Data and Processes: radiative transfer— software: development— galaxies: ultraviolet — radiation: dynamics — Interstellar Medium: H II regions
Almost everything we know about galaxies and most of whatwe know about stars comes from studying their radiation.However, radiation is not just a messenger informing usabout the sources and sinks of radiation, but may impactgas directly, for example through photoheating or the sup-pression of cooling, or by affecting its chemistry. Radiationpressure on gas and dust can also affect the dynamics ofthe gas directly. Unfortunately, including the effects of radi-ation in numerical models is challenging: the equation thataccounts for the change of intensity of a light ray result-ing from emission and absorption is 7 dimensional. To makematters worse, radiation travels at the speed of light, re-quiring dramatically shorter timesteps than those requiredto solve the associated hydrodynamics equations.Progress has been made by concentrating on particu-lar aspects of the impact of radiation. We briefly mentionsome of these aspects and the codes in which they are im- (cid:63)
Email: (TKC)[email protected] plemented, without aiming to be exhaustive. The cloudy code, last described by Ferland et al. (2017), implements ingreat detail the interaction between radiation and matter insimple geometries assuming equilibrium conditions. cloudy has been instrumental in interpreting the spectra of galax-ies. Accounting for absorption and re-emission of light bydust in more complex geometries has been implemented us-ing Monte-Carlo ray-tracing in for example the skirt (Baes& Camps 2015), sunrise (Jonsson 2006), and
CMacIonize (Vandenbroucke & Wood 2018) codes. The resonant scat-tering of Lyman- α has been implemented by, for exampleZheng & Miralda-Escudé (2002); Cantalupo et al. (2005);Verhamme et al. (2006) and others. Radiation can also regu-late star formation through radiative feedback. The infraredradiation on the interstellar medium (ISM) is modelled in,e.g., Turner & Stone (2001); Davis et al. (2012). Radiativefeedback is also important in the formation of the first starsand galaxies (e.g Bromm et al. 2009; Kim et al. 2019).In this paper we concentrate on the propagation of(hydrogen) ionizing photons. Radiative transfer (henceforthRT) of ionizing photons is important in the context of galax- © 0000 RAS a r X i v : . [ a s t r o - ph . I M ] F e b T. K. Chan et al. ies, governing the evolution of HII regions in the interstel-lar medium (ISM), and in the physics of the intergalacticmedium (IGM), which is highly ionized (Gunn & Peterson1965) by radiation from active galactic nuclei (AGN, Sargentet al. 1980) and massive stars in galaxies (Shapiro & Giroux1987; Madau & Meiksin 1994). In both situations, the follow-ing considerations are relevant to the design of a successfulRT implementation: (1) there is no useful symmetry to beexploited; (2) radiation is emitted by numerous sources; and(3) gas and radiation interact under non-equilibrium condi-tions. In addition, even without including RT, simulating theISM and the IGM is computationally demanding requiringthe inclusion of many other physical processes. These con-siderations motivate us to build RT on top of an existinghydrodynamics code, and implement a method that is inde-pendent of the number of sources.Smooth Particle Hydrodynamics (SPH) (Lucy 1977;Gingold & Monaghan 1977) is a Lagrangian hydrodynamicsscheme that has been applied to a large variety of astro-physical problems (from planet to star to galaxy formationsimulations) as well as non-astrophysical problems. In thisscheme, the hydrodynamic properties of a fluid are carriedby a set of discrete particles that move with the fluid andare used to interpolate physical quantities such as densityusing a smooth function called the ‘kernel’. The method iscomputationally efficient, highly adaptive in space and time,and can easily be coupled to gravity. Many current state-of-the art astrophysical hydrodynamics codes are SPH based(e.g. Springel 2005; Hopkins et al. 2014; Schaye et al. 2015;Wadsley et al. 2017; Price et al. 2018; Springel et al. 2020).We briefly discuss available options for including RT inhydrodynamical codes, especially for transporting ionizingphotons. Conceptually most intuitive is direct ray tracing(also called ‘long characteristics’), where each source castsa number of rays and the equation for RT is solved alongall rays simultaneously. With a computational cost scalingas N source × N , where N source and N sink are the numberof sources and sinks, this method may be accurate but itis also computationally extremely demanding. Approxima-tions to full ray-tracing are possible though, for example us-ing short characteristic (Mihalas et al. 1978; Mellema et al.2006), hybrid characteristic (Rijkhorst et al. 2006), or adap-tive ray-tracing (Abel & Wandelt 2002; Wise & Abel 2011).While the short characteristic method is faster than the longcharacteristic methods, its angular resolution is lower (Fin-lator et al. 2009); it is difficult to handle bright point sources(e.g Davis et al. 2012); and it has not yet been implementeddirectly on top of irregular meshes or particle-based codessuch as SPH (though see Finlator et al. 2009). Adaptiveray-tracing is fast and can be applied to irregular meshes(and particle methods), so it remains a viable option forRT in SPH, although its computational cost still scales withthe number of sources. In cases where the radiation field islargely known, reverse ray-tracing (Kessel-Deynet & Burk-ert 2000; Altay & Theuns 2013) has been used to calculatethe attenuation of ionizing radiation in high density regions.A variation of reverse ray-tracing has proved to be efficientlyparallelizable in SPH (e.g. Susa 2006; Hasegawa & Umemura2010).An alternative to ray-tracing is to discretize radiationdirections in a finite number of cones (Pawlik & Schaye 2008;Petkova & Springel 2011). The scaling of this implementa- tion is independent of the number of sources provided a ‘conemerging scheme’ is implemented. The method has been ap-plied in reionization simulations (Pawlik et al. 2017). How-ever, the method is still relatively expensive, given that ahigh number of cones is required to avoid excessive loss ofangular resolution. It also requires substantial modificationsto the hydrodynamics code, e.g. virtual particles and ro-tation of cones to improve the angular sampling and avoidartificial radiation spikes. Another strategy, the Monte Carlomethod (e.g. Altay et al. 2008; Baek et al. 2009; Grazianiet al. 2013), is even more expensive requiring a large numberof photon packets to reduce shot noise to acceptable levels.A different starting point for an RT algorithm is to com-pute angular and spectral moments of the radiative trans-fer equation and integrate the resulting ‘moment’ equationsnumerically. It is the radiation equivalent of integrating thefluid equations rather than the full Boltzmann equation. Inboth cases, doing so leads to a dramatic reduction in the di-mensionality of the problem. Just as in the case of the fluidequations, there is an infinite hierarchy of moment equa-tions which needs to be truncated by a ‘closure relation’.The closure relation is not unique and obtaining a good clo-sure relation is challenging, because it needs to be able tohandle the very different nature of the transport of opticallythick and optically thin radiation.RT moment methods vary in terms of the order of themoments used and in the choice of closure relation. Ideally,the closure relation uses only local properties of the gas andthe radiation: this makes the computational cost indepen-dent of the number of sources and makes the implementa-tion easily parellelisable. Moment methods do not requirefine angular discretization - unlike cone-based or short char-acteristic methods - so the computational cost per cell orparticle can be lower.The ‘Flux Limited Diffusion’ method (FLD, Levermore& Pomraning 1981) solves only for the zeroth moment ofthe RT equation, which is a diffusion equation providedthe time derivative of the first moment is neglected. Thespeed with which a radiation front propagates is not lim-ited by the speed of light but can be infinite; however, itis possible to impose a ‘flux limiter’ to enforce causality.FLD has been used in many astrophysical simulations (e.g.Turner & Stone 2001; Reynolds et al. 2009; Commerçonet al. 2011; Krumholz & Thompson 2012), some of whichuse SPH (Whitehouse & Bate 2004). Gnedin & Abel (2001)developed the otvet method which also only uses the zerothmoment, but with a closure relation applicable to opticallythin radiation. FLD and otvet are fast with a computetime that is largely independent of the number of sources.However, the relatively diffusive nature of transport inFLD makes it hard to preserve the propagation directionof radiation accurately. As a consequence, neither standardFLD nor otvet cast sharp shadows behind optically thickregions, albeit for slightly different reasons (Hayes & Nor-man 2003; Gnedin & Abel 2001). The timestep for prop-agating radiation in these methods is very restrictive as aconsequence of the infinite propagation speed of information;thus, it may be more efficient to use an implicit integrationscheme. Unfortunately, an implicit method is computation-ally inefficient in a scheme like SPH for a large neighbournumber (Whitehouse & Bate e.g 2004; Petkova & Springel © 0000 RAS, MNRAS , 000–000 moothed Particle Radiation Hydrodynamics—SPH-M1RT e.g 2009; see also the discussion about the efficiency of theimplicit FLD scheme in Skinner & Ostriker 2013).The ‘Two Moment’ method solves the zeroth and firstorder moments of the RT equation simultaneously. A popu-lar closure relation for this method was introduced by Lev-ermore (1984) to which we will refer as the ‘M1’ closurerelation . The M1 method was first used in astrophysicsby González et al. (2007), and has also been implementedin other hydrodynamics codes, e.g. grid-based (Aubert &Teyssier 2008, Skinner & Ostriker 2013, Rosdahl et al. 2013,and Kannan et al. 2019) and hybrid schemes (Hopkins et al.2020).This computational scheme is accurate up to order v/c (the fluid velocity divided by the speed of light; Buchler1983) for a single source in the optically thin or thick lim-its (Levermore 1984). In the optically thick case, it capturesthe minimum entropy (production) principle in the pres-ence of one preferred direction (Levermore 1996; Dubroca& Feugeas 1999). In the optically thin case, it preserves theradiation’s direction - and hence it can cast shadows - withradiation fronts moving at the speed of light. It may be sur-prising at first, but this second-order method is generally faster than FLD or otvet if solved explicitly. This is a con-sequence of the hyperbolic nature of the equations whichresult in a much less restrictive timestep (Thomas 1998).The speed and accuracy of the method makes it a promis-ing scheme for including RT in astrophysical hydrodynamicscalculations.While the M1 method works well on structured and un-structured meshes, to date, it has not been implemented inSPH. One reason is that SPH has zeroth-order errors un-der irregular particle distributions meaning that the SPHestimate does not converge to the true value in the limit ofvanishing smoothing length (Lucy 1977; Raviart 1985; Lan-son & Vila 2008; Dehnen & Aly 2012). Secondly, devisinga good artificial dissipation scheme in SPH is not trivial.However, such artificial dissipation is necessary in order tosuppress numerical oscillations around discontinuities. As wewill demonstrate, the usual scheme, e.g. Price (2008), for im-plementing artificial dissipation fails when applied to the M1scheme. Finally, the original M1 closure relation artificiallyamplifies noise in the optically-thin regime which requireschanges to the closure relation.Despite these difficulties, implementing the M1 RTmethod in SPH would be highly desirable: SPH is highlyadaptive and ideal for problems that are characterised by avery large dynamic range, whereas the M1 method is effi-cient and accurate in both the optically thick and thin lim-its . Furthermore, the M1 method can be straightforwardlyimplemented on top of any SPH code, since the structures ofthe hydrodynamics equations and radiation moment equa-tions are quite similar. This results in an accurate and fastcode that can handle a very large number of sources ina computationally efficient way. As such, the method de-scribed in this paper goes some way towards enabling the The ‘M’ in ‘M1’ refers to Minerbo, who introduced the maxi-mum entropy closure in Minerbo (1978). But it has issues in handling multiple sources in the opticallythin region; see §4. inclusion of RT in simulations of galaxy formation as a mat-ter of course.In this paper, we describe how to incorporate the M1method into SPH and examine its performance throughstandard RT problems. We begin in §2 by briefly illuminat-ing the analogy between taking moments of the Boltzmannequation to derive the fluid equations, and taking momentsof the RT equation to derive the two-moment method. Wethen discuss closure relations and discuss our modification tothe M1 closure. Next, we show how the SPH equations canbe dicretized to yield the more accurate gradients requiredfor implementing the two-moment method and discuss waysof capturing discontinuities in the radiation field. We finish§2 by discussing the coupling of radiation to the thermo-dynamics and chemistry of the gas, explain and discuss theadvantages and drawbacks of the ‘reduced speed-of-light’ ap-proximation, and discuss how we inject radiation. In §3, wepresent the results of tests with a known solution and com-pare more realistic tests without a known solution to thosein the RT code comparison project (Iliev et al. 2006, 2009).In §4, we comment on the strengths and weaknesses of ourscheme and compare with other radiative transfer methods.In §5, we briefly summarise our findings and foresee possibleimprovements in the future.
In the following, we first describe the two moment methodincluding modifications to the M1 closure relation. We con-tinue by discussing the implementation in SPH as well asthe thermochemistry solver we employ. The equations con-tain numerous variables which we have collated for easy ref-erence in Table 1.
The radiative transfer equation expresses the constancy ofthe specific intensity ( I ; in erg cm − s − Hz − sr − ) of a beamof light in the absence of sources or sinks (e.g. Pomraning1973) c ∂∂t I + ˆ n · ∂∂ x I = (cid:18) DIDt (cid:19) SS . (1)Here, I is a function of position ( x ), direction ( ˆ n ), frequency( ν ) and time ( t ); the right hand side represents sources andsinks of light. In the astronomical literature I is usuallycalled the ‘surface brightness’, and just as surface bright-ness, I also suffers from cosmological dimming or Dopplereffects due to the redshifting of photons, which is not cap-tured by the approximate Eq. (1) (see Buchler 1983 for thederivation of the more complete RT equation). Moment methods drastically simplify the solution of thisequation by multiplying Eq. (1) with some function of direc-tion and integrating the resulting equation over solid angle.This yields an infinite number of moment equations, with thehierarchy closed after a finite number of moments by a ‘clo-sure relation’. Solving the resulting RT moment equations isthe radiation equivalent of solving the fluid equations ratherthan the collisional Boltzmann equation. We point the inter-ested reader to a sketch of the derivation of these momentequations and the relation to fluid equations in Appendix © 0000 RAS, MNRAS000
The radiative transfer equation expresses the constancy ofthe specific intensity ( I ; in erg cm − s − Hz − sr − ) of a beamof light in the absence of sources or sinks (e.g. Pomraning1973) c ∂∂t I + ˆ n · ∂∂ x I = (cid:18) DIDt (cid:19) SS . (1)Here, I is a function of position ( x ), direction ( ˆ n ), frequency( ν ) and time ( t ); the right hand side represents sources andsinks of light. In the astronomical literature I is usuallycalled the ‘surface brightness’, and just as surface bright-ness, I also suffers from cosmological dimming or Dopplereffects due to the redshifting of photons, which is not cap-tured by the approximate Eq. (1) (see Buchler 1983 for thederivation of the more complete RT equation). Moment methods drastically simplify the solution of thisequation by multiplying Eq. (1) with some function of direc-tion and integrating the resulting equation over solid angle.This yields an infinite number of moment equations, with thehierarchy closed after a finite number of moments by a ‘clo-sure relation’. Solving the resulting RT moment equations isthe radiation equivalent of solving the fluid equations ratherthan the collisional Boltzmann equation. We point the inter-ested reader to a sketch of the derivation of these momentequations and the relation to fluid equations in Appendix © 0000 RAS, MNRAS000 , 000–000
T. K. Chan et al.
D/Dt
LagrangianDerivative 2 ρ gas density 2 v gas velocity 2 p gas pressure 3 χ opacity 3 ˜ c reduced speed oflight 6 f radiation flux /gas density 5 F radiation flux 5 u gas internal en-ergy 4 ξ radiation energydensity / gasdensity 5 P radiation stresstensor 5 F Eddington ten-sor 7 e gas thermal en-ergy 4 E radiation energydensity 5 f Edd
Eddington fac-tor 7 ˆ n radiation direc-tion 13 ν photon fre-quency 9 Ω solid angle 8 I specific inten-sity 1 (cid:126) reduced Planckconstant 9 π (cid:126) ¯ ν spectral meanphoton energy 9 ε optical thick-ness estimator 18 h SPH particlesize 24 m SPH particlemass 24 D diffusion coeffi-cient 27 W SPH kernel 24 v sig SPH signalspeed 30 φ slope limiter 31 n γ photon numberdensity 9 f γ photon numberflux 9 x , x HI neutral hydro-gen fraction 38 α f , α ξ artificial dissi-pation factor 36 Λ combined heat-ing and coolingrate 4 S injection sourcerate 4 c physical speedof light (cid:15) γ photon-ionizationheating perionization A5 &41 σ γ photon-ionizationcross-section 39 &A2 Ω i correctionfor variablesmoothinglength 23 Table 1.
A list of variables, including their symbols, names, and the equations they are first used.
D. It is worth recalling that fluid equations, being differen-tial equations, do not properly describe the behaviour of aset of particles in case of discontinuities such as shocks orcontact discontinuities. Their numerical integration requiresthe addition of extra terms (such as ‘artificial viscosity’ or‘artificial conduction’). The same is true for moments of theRT equation, and we describe the discontinuity capturingscheme below, after we introduce the full moment equationsfor fluid and radiation combined in the next section.
The moment equations describing the interaction of a fluidwith radiation are (e.g Buchler 1983; Mihalas & Mihalas1984):
DρDt + ρ ∇· v = 0 , (2) D v Dt = − ∇ pρ − ∇ φ + χρ ˜ c f + S v , (3) DuDt = − pρ ∇ · v + Λ u + S u , (4) DξDt = − ρ ∇ · ( ρ f ) − ∇ v : P ρ + Λ ξ + S ξ , (5) c DDt f = − ∇ · P ρ − χρ ˜ c f + S f , (6) P = F E = F ρξ . (7)These equations are series expansions including all terms upto v/c ; they include Doppler shifting but not cosmologicalredshifting. The fluid variables are mass density ( ρ ), velocity( v ), pressure ( p ), and thermal energy per unit mass ( u ); ∇ φ is the gravitational acceleration. The radiation variables arethe energy density per unit volume ( E ), flux ( F ) and the‘radiation stress tensor’ ( P ); D/Dt is the Lagrangian timederivative.The first few angular moments of I are E ν = 1˜ c (cid:90) I dΩ , F iν = (cid:90) ˆ n i I dΩ , P ijν = (cid:90) ˆ n i ˆ n j I dΩ , (8)When additionally integrated over frequency (in the ‘grey’approximation, e.g. Turner & Stone 2001), those angularmoments become E , F and P , respectively. Integrating oversmall frequency intervals, in order to mimic multi-frequencyRT, is challenging when Doppler shifts or redshifts are large.This case is not considered here. © 0000 RAS, MNRAS , 000–000 moothed Particle Radiation Hydrodynamics—SPH-M1RT We further defined the ratio of the radiation energy den-sity over the fluid’s density, ξ ≡ E/ρ , and the ratio of radia-tion flux over the fluid’s density as f ≡ F /ρ . The tensor F isthe Eddington tensor, which we will discuss in §2.3; ∇ v : P is short hand for the contraction P ij v i,j .There are some further terms on the right hand sideof Eqs. 5 and 6 that represent sources and sinks. Λ u and Λ ξ are the combined heating and cooling rates for gas andradiation, respectively; χ is the opacity related to the opticaldepth per unit length as dτ /dr = χ ρ . We will describe theseterms in more detail in §2.7. S X is the source term thatcorresponds to quantity X: the injection of radiation will bedescribed in more detail in §2.8.The ( χρ/ ˜ c ) f term in Eq. (3) represents the radiationpressure acting on the gas. Currently, we apply radiationpressure inferred from the quantities averaged over the vol-ume of each particle, as described below. However, Hopkins& Grudić (2019) demonstrated that it is more accurate toapply radiation pressure to the interface between particles,an improvement we intend to implement in future. In thecase of ionizing radiation propagating through a low resolu-tion simulation - for example when simulating cosmic reion-ization - the resulting differences are expected to be smallbecause the radiation imparts little momentum. However, amore accurate treatment of radiation pressure may be re-quired in high-resolution simulations to capture radiationpressure from massive stars or active galactic nuclei.In this paper we propagate radiation at the speed ˜ c < c ,which is a ‘reduced’ speed of light. The motivation, validityand limitations of this approximation are discussed in §2.10.In the two-moment method, the time derivatives of theradiation density and radiation flux are kept, unlike in thecase of flux limited diffusion (FLD, Levermore & Pomran-ing 1981). There are two advantages in keeping this term.Firstly, Buchler (1983) showed that the time derivative of f may be significant in the optically thin (free streaming)regime, making the two-moment method more accurate thanFLD. Secondly, including the time derivative yields hyper-bolic equations rather than the parabolic equation of FLD.Solving a parabolic differential equation explicitly requiresa more restrictive timestep, ∆ t ∝ (∆ x ) / ˜ c , compared to thehyperbolic case where ∆ t ∝ ∆ x/ ˜ c ; where ∆ x is the spatialresolution. Combined with using a reduced speed of lightapproximation ( ˜ c rather than c ) improves the efficiency ofthe RT implementation compared to FLD .The relations between the photon number density, n γ ,and the radiation energy density, E , and between the photonflux, F γ and the radiation flux, F , are n γ = E π (cid:126) ¯ ν ; F γ ≈ F π (cid:126) ¯ ν , (9)where the mean photon energies π (cid:126) ¯ ν are (cid:126) ¯ ν = (cid:20)(cid:90) I dΩd ν (cid:21) (cid:20)(cid:90) ( I/ (cid:126) ν )dΩd ν (cid:21) − ; (cid:126) ¯ ν ≈ (cid:20)(cid:12)(cid:12)(cid:12)(cid:12)(cid:90) ˆ n I dΩd ν (cid:12)(cid:12)(cid:12)(cid:12)(cid:21) (cid:20)(cid:12)(cid:12)(cid:12)(cid:12)(cid:90) ˆ n ( I/ (cid:126) ν )dΩd ν (cid:12)(cid:12)(cid:12)(cid:12)(cid:21) − . (10) The M1 method is faster even if FLD is solved implicitly be-cause the inversion step in the implicit solver is expensive (see,e.g. Skinner & Ostriker 2013).
The second relation is a good approximation when the radi-ation is either isotropic or optically thin. For reference, themean photon energy π (cid:126) ¯ ν of ionizing radiation is 29.6eV fora black-body spectrum at T = 10 K ; (cid:126) is Planck’s constantdivided by π (see Appendix A for details). Taking successive angular moments of the RT relation leadsto an infinite set of coupled moment equations (Mihalas &Mihalas 1984). A ‘closure relation’ which relates higher or-der moments to lower-order ones, is required to break thishierarchy. Unfortunately, the closure relation is not uniqueand depends on the problem at hand. Levermore (1984) de-rived a closure relation as follows: consider the RT equationsin the fluid frame by setting v = 0 , in which case Eqs. (5)and (6) simply to the following two moment equations, ∂E∂t = −∇ · F , (11) c ∂∂t F = −∇ · P − χρ ˜ c F , (12)in the absence of sources and sinks. Provided that the radi-ation field is symmetric around a given direction ˆ n , Lever-more (1984) demonstrated that the second moment can bewritten as P = E − f Edd ) I + E f Edd − ˆnˆn , (13)where f Edd is called the ‘Eddington factor’.When the radiation field is almost isotropic, P ij ≈ ( E/ δ ij which corresponds to f Edd = 1 / . Combining thetwo moment equations with this relation yields ∂E∂t = ∇ · (cid:18) ˜ c χρ ∇ E (cid:19) − ∇ · (cid:18) cχρ ∂∂t F (cid:19) . (14)This describes isotropic diffusion of the energy density, E ,in case the rate of change of the flux (the last term on theright hand side) can be neglected. Of course, if the radi-ation field were exactly isotropic everywhere it has to beuniform as well - but this diffusion approximation can beused, provided E varies sufficiently slowly in space and time(Levermore 1984). This case corresponds to the classical ‘Ed-dington’ approximation for the propagation of radiation inthe isotropic case, and we will refer to as the ‘optically-thick’solution.In contrast, the value f Edd = 1 leads to anisotropicradiation propagation with ∇ · P = ˆ n (ˆ n · ∇ E ) . (15)In this ‘optically-thin’ case, ∂ E∂t + ˜ c (ˆ n · ∇ ) E = − ˜ cχρ ∂E∂t , (16)and radiation ‘streams’ in direction ˆ n with speed ˜ c , with itsintensity decreasing due to absorption as quantified by theright hand side of the equation.The closure relation of Eq. (13) therefore captures thepropagation correctly in the two limiting cases of (1) high-optical depth, with the solution describing isotropic diffu-sion, and (2) the optically-thin regime of negligible opticaldepth, where the solution describes free propagation at the © 0000 RAS, MNRAS000
The second relation is a good approximation when the radi-ation is either isotropic or optically thin. For reference, themean photon energy π (cid:126) ¯ ν of ionizing radiation is 29.6eV fora black-body spectrum at T = 10 K ; (cid:126) is Planck’s constantdivided by π (see Appendix A for details). Taking successive angular moments of the RT relation leadsto an infinite set of coupled moment equations (Mihalas &Mihalas 1984). A ‘closure relation’ which relates higher or-der moments to lower-order ones, is required to break thishierarchy. Unfortunately, the closure relation is not uniqueand depends on the problem at hand. Levermore (1984) de-rived a closure relation as follows: consider the RT equationsin the fluid frame by setting v = 0 , in which case Eqs. (5)and (6) simply to the following two moment equations, ∂E∂t = −∇ · F , (11) c ∂∂t F = −∇ · P − χρ ˜ c F , (12)in the absence of sources and sinks. Provided that the radi-ation field is symmetric around a given direction ˆ n , Lever-more (1984) demonstrated that the second moment can bewritten as P = E − f Edd ) I + E f Edd − ˆnˆn , (13)where f Edd is called the ‘Eddington factor’.When the radiation field is almost isotropic, P ij ≈ ( E/ δ ij which corresponds to f Edd = 1 / . Combining thetwo moment equations with this relation yields ∂E∂t = ∇ · (cid:18) ˜ c χρ ∇ E (cid:19) − ∇ · (cid:18) cχρ ∂∂t F (cid:19) . (14)This describes isotropic diffusion of the energy density, E ,in case the rate of change of the flux (the last term on theright hand side) can be neglected. Of course, if the radi-ation field were exactly isotropic everywhere it has to beuniform as well - but this diffusion approximation can beused, provided E varies sufficiently slowly in space and time(Levermore 1984). This case corresponds to the classical ‘Ed-dington’ approximation for the propagation of radiation inthe isotropic case, and we will refer to as the ‘optically-thick’solution.In contrast, the value f Edd = 1 leads to anisotropicradiation propagation with ∇ · P = ˆ n (ˆ n · ∇ E ) . (15)In this ‘optically-thin’ case, ∂ E∂t + ˜ c (ˆ n · ∇ ) E = − ˜ cχρ ∂E∂t , (16)and radiation ‘streams’ in direction ˆ n with speed ˜ c , with itsintensity decreasing due to absorption as quantified by theright hand side of the equation.The closure relation of Eq. (13) therefore captures thepropagation correctly in the two limiting cases of (1) high-optical depth, with the solution describing isotropic diffu-sion, and (2) the optically-thin regime of negligible opticaldepth, where the solution describes free propagation at the © 0000 RAS, MNRAS000 , 000–000 T. K. Chan et al. speed ˜ c in the characteristic direction ˆn . The expectation isthen that Eq. (13) also provides a good approximation toany intermediate case (Levermore & Pomraning 1981).One disadvantage of the scheme is that radiation be-haves as a ‘collisional’ fluid: beams of light with differentpropagation directions ˆ n that intersect will collide. This isbecause the local Eddington tensor closure relation of Eq.(13) can only handle one direction ˆ n at a time (in addition toan isotropic component). We will discuss this issue in moredetails in §4.1.1. Next we turn to the choice of f Edd . As shown by Levermore(1984), given that F / (˜ cE ) and P /E are the first and secondmoments of a non-negative unit density variable requiresthat ε = (cid:12)(cid:12)(cid:12)(cid:12) F ˜ cE (cid:12)(cid:12)(cid:12)(cid:12) = (cid:12)(cid:12)(cid:12)(cid:12) f ˜ cξ (cid:12)(cid:12)(cid:12)(cid:12) (cid:54) f Edd (cid:54) , (17)which we term the ‘original’ closure.Of course, even if we demanded that the Eddington fac-tor should only depend on the local values of F and E , thenthis would not specify f Edd uniquely (see Levermore 1984for a summary of reasonable choices). One particular choiceis the ‘M1’ closure, which Levermore (1984) derived by as-suming that there exist inertial frames in which the radia-tion density is isotropic (not necessarily isotropic in the labor fluid frame). This original M1 relation is f Edd = 3 + 4 ε √ − ε . (18)Dubroca & Feugeas (1999) showed that this corresponds tothe simplest moment closure that maximizes the entropyand is anisotropic .The evaluation of this expression for M1 is computa-tional efficient as well as highly parallelisable, as comparedto e.g. ray-tracing or Monte Carlo methods, because f Edd de-pends only on local quantities. Because of this, several astro-physical RT implementations use this M1 closure relation,e.g. González et al. (2007); Aubert & Teyssier (2008); Skin-ner & Ostriker (2013); Rosdahl et al. (2013); Kannan et al.(2019); Skinner et al. (2019). However, this choice is notwithout its problems (as are other variants based on localvariables). Firstly, consider the case of two otherwise iden-tical beams of radiation propagating in opposite directions.Where the beams hit the net flux is zero, | f / (˜ cξ ) | ∼ so that (cid:15) = 0 and f Edd = 1 / : this corresponds to the optically-thicksolution, even in the system were optically thin. It is as if thebeams of radiation collide with each other (see also Rosdahlet al. 2013). Clearly, this behaviour is incorrect.This choice of closure relation also results in artificialdispersion, since radiation does not move at the same speedwhen | f / (˜ cξ ) | varies: radiation propagates with speed be-tween ˜ c/ and ˜ c , when Eq. (14) or Eq. (16) applies, respec-tively.An improved closure relation can be derived from the Note that in the mathematics community, the entropy has theopposite sign compared to that in the physics community. following considerations. In the optically thick limit, we de-sire that f Edd = 1 / , since the corresponding isotropic dif-fusion captures the random walks of photons through themedium as a consequence of numerous independent scat-tering events. In the opposite limit of an optically thinmedium, we desire that f Edd = 1 , since that correctly de-scribes streaming of radiation at the speed of light. Finally,we require that f Edd (cid:54) .Our proposed modification to the closure relation is ε = max [exp( − χρh ) , | f / (˜ cξ ) | ] , (19)where τ ≡ χρh is the local optical depth across the extent h of a resolution element. This choice satisfies ε (cid:54) , andhas the correct limiting behaviour: ε → in the opticallythin limit of τ → and ε → | f / (˜ cξ ) | → in the opticallythick case of τ → ∞ when the flux | f | is small according toEq. (6). In case | f | is small due to the ‘collision’ of two beamsof radiation, ε can still be of order 1 and describe radiationstreaming rather than diffusion provided the optical depthis small. We will demonstrate in §4.1.1 that our modifiedM1 closure (Eq.19) can handle head-on beam collisions andmore generally, 1D RT problems.We choose to modify M1 by the factor exp( − χρh ) tomimic the diffusion of radiation when the optical depth islarge. The choice is also motivated by a desire to help numer-ical convergence: the combined contributions of two resolu-tion elements, for example two SPH particles with extents h i and h j , is approximately exp( − χρh i ) × exp( − χρh j ) =exp[ − χρ ( h i + h j )] , which corresponds to the approximateeffect of a lower resolution SPH particle with size ( h i + h j ) .We will also show in Fig.1 that our scheme is more sta-ble than the original M1 closure in optically thin regionswhen simulated with SPH. However, the scheme does notsolve the problem of the artificial collision of radiation beamsin case they are not head on (§4.1.1), since the radiation di-rections will still merge locally according to Eq. (13). Fortu-nately, even in this case, our closure (Eq.19) will still preventthe numerical diffusion in the optically thin limit.Finally, we justify the use of physical quantities otherthan radiation energy and flux in the Eddington factor. TheEddington tensor should be derived from the RT equation,which contains information about the gas, e.g. density, veloc-ity, and the opacity (through the collisional term). Thus, theEddington tensor should be also a function of these gas prop-erties. In fact, in the absence of the collision term ( χ = 0 ),the radiation should always be free streaming at the thespeed of light regardless of the value of ξ and f . In the standard formulation of SPH (e.g Monaghan 2002),the density, ρ i , at the location of particle i is calculatedthrough interpolating over ‘neighbouring’ particles in agather approach, ρ i = (cid:88) j m j W ij ( h i ) , (20)where the kernel W ij ( h i ) = W ( | r i − r j | , h i ) is a functionwith compact support (by default the M cubic B-splinefunction), h i is the smoothing length and m i the mass of © 0000 RAS, MNRAS , 000–000 moothed Particle Radiation Hydrodynamics—SPH-M1RT particle i . We follow the variable smoothing length treat-ment similar to that in Springel & Hernquist (2002) suchthat the number of neighbour particles that contribute tothe sum is N ngb (=48 in 3D) (see the GADGET-2
SPH sec-tion in Schaller et al. 2015 for more details, including theSPH formulation of the hydrodynamics in
SWIFT ).In the radiation hydrodynamics tests presented below,we do not use the standard SPH formulation but ratherthe modifications introduced by Borrow et al. (2020) called
SPHENIX , which uses the density and energy hydrodynamicvariables, rather than pressure and energy.
SPHENIX appliesthe Cullen & Dehnen (2010) shock detector to minimize ar-tificial viscosity away from shocks and the artificial diffusionterm to capture fluid mixing described by Price (2008).One of the main hurdles to over come for implementinga moment method in SPH is that such a higher-order methodrequires the calculation of derivatives, and these tend to benoisy when the particle distribution is irregular. For exam-ple, there are several ways to estimate the divergence of avector field X , which include (e.g. Tricco & Price 2012) the‘ symmetric ’ estimate: ( ∇ · X ) i (cid:12)(cid:12) sym = ρ i (cid:88) j m j (cid:20) X i Ω i ρ i · ∇ i W ij ( h i )+ X j Ω j ρ j · ∇ i W ij ( h j ) (cid:21) , (21)and the ‘ difference ’ estimate: ( ∇ · X ) i (cid:12)(cid:12) diff = − (cid:88) j m j Ω i ρ i ( X i − X j ) · ∇ i W ij ( h i ) , (22)where X is an arbitrary vector or tensor associated witheach particle, and Ω i = 1 + h i ρ i (cid:88) j m j ∂W ij ( h i ) ∂h i , (23)is a correction factor introduced by Springel & Hernquist(2002) to account for spatial variations in the value of thesmoothing length, h .We use the difference form to evaluate Eqs.(5) and (6): (cid:18) Dξ i Dt (cid:19)(cid:12)(cid:12)(cid:12)(cid:12) diff = (cid:20) − ρ ∇ · ( ρ f ) (cid:21)(cid:12)(cid:12)(cid:12)(cid:12) i = − (cid:88) j m j Ω i ρ i ( ρ i f i − ρ j f j ) · ∇ i W ij ( h i ) , (24)and (cid:18) c D f i Dt (cid:19)(cid:12)(cid:12)(cid:12)(cid:12) diff = (cid:20) − ∇ · P ρ − χρ ˜ c f (cid:21)(cid:12)(cid:12)(cid:12)(cid:12) i = − (cid:88) j m j Ω i ρ i ( ρ i ξ i F i − ρ j ξ j F j ) · ∇ i W ij ( h i ) − χ i ρ i ˜ c f i . (25)The difference formulation subtracts the zeroth-order errorsthat occur in SPH explicitly, yielding first-order accuracyregardless of the underlying particle distribution. This re-sults in superior numerical estimates of the divergence par-ticularly near steep gradients. However, the difference esti-mate does not manifestly conserve flux, unlike the ‘symmet-ric’ estimate . Fortunately, we find that the level of non- The ‘symmetric’ SPH form can also help to regularize the par- conservation of flux is small in our experiments (Typicallyless than one percent, better accuracy could be reached byincreasing the order of the scheme, if required.). There is noknown formulation that simultaneously avoids zeroth-ordererrors and is manifestly conservative in SPH (see the discus-sion in Price 2012).We add the term − χ i ρ i f i / ˜ c to Eq. (6) using operatorsplitting, f i ( t + ∆ t ) = exp( − χ i ρ i ˜ c ∆ t ) × f i ( t ) . (26)Though this scheme in unconditionally stable, it neverthe-less yields the wrong answer when the timestep, ∆ t , is toolong. This could be avoided by limiting the timestep to ∆ t (cid:54) / ( χρ ˜ c ) , but that would result in unacceptably shorttimesteps in regions of high optical depth. Since in such re-gions the impact of radiation may be small anyway, we willlimit the timestep by the usual Courant–Friedrichs–Lewy(CFL) condition, ∆ t (cid:54) . h/ ˜ c . To ensure that our results arephysically meaningful, numerically stable and satisfy causal-ity, we apply the following additional limiters at the begin-ning of each timestep: ( i ) | f | (cid:54) ˜ cξ , ( ii ) ξ (cid:62) , and ( iii ) wezero unused components of f in 1D or 2D simulations. Thelatter limiter corrects for any numerical scatter of radiationinto unused dimensions, as may happen if the Eddingtontensor is non-zero but | f | is small. In general, we set the propagation direction, ˆ n , to be thatof the local flux, ˆ n = ˆ f . However, radiation propagates in aconstant direction in the optically thin case. Since the fluxis computed numerically, round-off errors or numerical noisemay rotate the flux vector so that imposing ˆ n = ˆ f does notguarantee that radiation travels in a straight line, even inthe optically thin case.In some special cases, for example of light emanatingfrom a single point source or the propagation of a planeparallel radiation front, the direction ˆ n is known a priori, andwe can therefore choose to simply impose the propagationdirection, and use that direction to compute the opticallythin Eddington tensor, Eq. (13).Not surprisingly, imposing the direction of radiationpropagation yields spherical ionization regions around apoint source (Fig.8) and casts sharp shadows behind op-tically thick absorbers even at low resolution (Fig.11). Ofcourse in general, the propagation direction ˆ n is not gen-erally known, for example, there may be several sources oran additional isotropic background, which require improve-ments of our scheme. The fluid equations encoded in SPH are differential equa-tions and hence need to be supplemented with extra terms ticle distribution in hydrodynamics calculations, albeit by intro-ducing purely numerical forces (Price 2012). This is less importantwhen propagating ionising radiation which does not usually exertstrong forces on the gas particles. In the timestep determination, we will use the smallest h of allneighbouring SPH particles and of the particle itself, to ensurestability and conservation.© 0000 RAS, MNRAS000
SPHENIX appliesthe Cullen & Dehnen (2010) shock detector to minimize ar-tificial viscosity away from shocks and the artificial diffusionterm to capture fluid mixing described by Price (2008).One of the main hurdles to over come for implementinga moment method in SPH is that such a higher-order methodrequires the calculation of derivatives, and these tend to benoisy when the particle distribution is irregular. For exam-ple, there are several ways to estimate the divergence of avector field X , which include (e.g. Tricco & Price 2012) the‘ symmetric ’ estimate: ( ∇ · X ) i (cid:12)(cid:12) sym = ρ i (cid:88) j m j (cid:20) X i Ω i ρ i · ∇ i W ij ( h i )+ X j Ω j ρ j · ∇ i W ij ( h j ) (cid:21) , (21)and the ‘ difference ’ estimate: ( ∇ · X ) i (cid:12)(cid:12) diff = − (cid:88) j m j Ω i ρ i ( X i − X j ) · ∇ i W ij ( h i ) , (22)where X is an arbitrary vector or tensor associated witheach particle, and Ω i = 1 + h i ρ i (cid:88) j m j ∂W ij ( h i ) ∂h i , (23)is a correction factor introduced by Springel & Hernquist(2002) to account for spatial variations in the value of thesmoothing length, h .We use the difference form to evaluate Eqs.(5) and (6): (cid:18) Dξ i Dt (cid:19)(cid:12)(cid:12)(cid:12)(cid:12) diff = (cid:20) − ρ ∇ · ( ρ f ) (cid:21)(cid:12)(cid:12)(cid:12)(cid:12) i = − (cid:88) j m j Ω i ρ i ( ρ i f i − ρ j f j ) · ∇ i W ij ( h i ) , (24)and (cid:18) c D f i Dt (cid:19)(cid:12)(cid:12)(cid:12)(cid:12) diff = (cid:20) − ∇ · P ρ − χρ ˜ c f (cid:21)(cid:12)(cid:12)(cid:12)(cid:12) i = − (cid:88) j m j Ω i ρ i ( ρ i ξ i F i − ρ j ξ j F j ) · ∇ i W ij ( h i ) − χ i ρ i ˜ c f i . (25)The difference formulation subtracts the zeroth-order errorsthat occur in SPH explicitly, yielding first-order accuracyregardless of the underlying particle distribution. This re-sults in superior numerical estimates of the divergence par-ticularly near steep gradients. However, the difference esti-mate does not manifestly conserve flux, unlike the ‘symmet-ric’ estimate . Fortunately, we find that the level of non- The ‘symmetric’ SPH form can also help to regularize the par- conservation of flux is small in our experiments (Typicallyless than one percent, better accuracy could be reached byincreasing the order of the scheme, if required.). There is noknown formulation that simultaneously avoids zeroth-ordererrors and is manifestly conservative in SPH (see the discus-sion in Price 2012).We add the term − χ i ρ i f i / ˜ c to Eq. (6) using operatorsplitting, f i ( t + ∆ t ) = exp( − χ i ρ i ˜ c ∆ t ) × f i ( t ) . (26)Though this scheme in unconditionally stable, it neverthe-less yields the wrong answer when the timestep, ∆ t , is toolong. This could be avoided by limiting the timestep to ∆ t (cid:54) / ( χρ ˜ c ) , but that would result in unacceptably shorttimesteps in regions of high optical depth. Since in such re-gions the impact of radiation may be small anyway, we willlimit the timestep by the usual Courant–Friedrichs–Lewy(CFL) condition, ∆ t (cid:54) . h/ ˜ c . To ensure that our results arephysically meaningful, numerically stable and satisfy causal-ity, we apply the following additional limiters at the begin-ning of each timestep: ( i ) | f | (cid:54) ˜ cξ , ( ii ) ξ (cid:62) , and ( iii ) wezero unused components of f in 1D or 2D simulations. Thelatter limiter corrects for any numerical scatter of radiationinto unused dimensions, as may happen if the Eddingtontensor is non-zero but | f | is small. In general, we set the propagation direction, ˆ n , to be thatof the local flux, ˆ n = ˆ f . However, radiation propagates in aconstant direction in the optically thin case. Since the fluxis computed numerically, round-off errors or numerical noisemay rotate the flux vector so that imposing ˆ n = ˆ f does notguarantee that radiation travels in a straight line, even inthe optically thin case.In some special cases, for example of light emanatingfrom a single point source or the propagation of a planeparallel radiation front, the direction ˆ n is known a priori, andwe can therefore choose to simply impose the propagationdirection, and use that direction to compute the opticallythin Eddington tensor, Eq. (13).Not surprisingly, imposing the direction of radiationpropagation yields spherical ionization regions around apoint source (Fig.8) and casts sharp shadows behind op-tically thick absorbers even at low resolution (Fig.11). Ofcourse in general, the propagation direction ˆ n is not gen-erally known, for example, there may be several sources oran additional isotropic background, which require improve-ments of our scheme. The fluid equations encoded in SPH are differential equa-tions and hence need to be supplemented with extra terms ticle distribution in hydrodynamics calculations, albeit by intro-ducing purely numerical forces (Price 2012). This is less importantwhen propagating ionising radiation which does not usually exertstrong forces on the gas particles. In the timestep determination, we will use the smallest h of allneighbouring SPH particles and of the particle itself, to ensurestability and conservation.© 0000 RAS, MNRAS000 , 000–000 T. K. Chan et al. in order to correctly capture discontinuities such as shocksand contact discontinuities. These terms broaden disconti-nuities by introducing numerical dissipation so that they canbe resolved by the interpolation scheme (see e.g Monaghan1997; Agertz et al. 2007; Price 2008).The SPH implementation of the moment method needsto be extended with similar dissipation terms to handle dis-continuities in the radiation field, and we base these on theartificial diffusion and artificial viscosity terms of the SPHfluid equations. The energy dissipation term is similar to the artificial diffusion term in SPH hydrodynamics, (cid:18) D ξ i D t (cid:19)(cid:12)(cid:12)(cid:12)(cid:12) diss = N (cid:88) j =1 D ξ,ij m j ¯ ρ ( ˜ ρ i ξ i − ˜ ρ j ξ j ) ˆr ij · ∇ i W ij r ij , (27)where ¯ ρ = √ ρ i ρ j is the geometric mean of the densitiesof the pair of interacting particles i and j . If the densitycontrast is larger than 10, we found the scheme to be morestable with the choice ¯ ρ = min( ρ i , ρ j ) , but this choice is notused in the tests in this paper.Flux dissipation can modelled on the artificial viscosity term in SPH hydrodynamics , (cid:18) D f i D t (cid:19)(cid:12)(cid:12)(cid:12)(cid:12) diss = (cid:80) Nj =1 D f ,ij m j ¯ ρ ( ρ i f i − ρ j f j ) · ˆr ij ∇ i W ij r ij , if ( ρ i f i − ρ j f j ) · ˆr ij < , , otherwise . (28)In these expression, D ξ,ij and D f ,ij are the ‘artificial dissipa-tion’ coefficients, they have the units of a diffusion constantand we write them as D ξ,ij = α ξ,i v sig ,i h i + α ξ,j v sig ,j h j ,D f ,ij = α f ,i v sig ,i h i + α f ,j v sig ,j h j . (29)Here, v sig is the signal speed, v sig = | ˆf · ˆr ij | ˜ c, , (30)and ( α f , α ξ ) (cid:54) are dimensionless numbers that quan-tify the strength of the numerical dissipation. The formsof Eqs. (27) and (28) are consistent with the Riemannsolver across the boundary of two SPH particles (Mon-aghan 1997). The kernel averaged over smoothing length is ∇ i W ij = 0 . ∇ i W ij ( h j ) + ∇ i W ij ( h j )] .Eqs. (27) and (28) are diffusion equations (see e.g.Jubelgas et al. 2004; Price 2008). The maximum value ofthe diffusion speed is h ˜ c , where h is the particle size and ˜ c is the propagation speed of the radiation. A numerical dif-fusion coefficient larger than this maximum value will resultin numerical instabilities if the timestep is set by the CFLcondition, ∆ t (cid:46) h/ ˜ c .Price (2008) set ˜ ρ i ξ i = ρ i ξ i , the values associated withindividual SPH particles, and minimized the amount of nu-merical dissipation by choosing how the signal speed de-pends on local quantities. However, another way to minimizeartificial dissipation is by reconstructing fluid quantities atthe interface between particles (see e.g. Frontiere et al. 2017; This is not our default choice, see Eq. (32) below.
Rosswog 2020a). To do so, we reconstruct the radiation en-ergy density at the interface using a Taylor series expansion, ˜ ρ i ξ i − ˜ ρ j ξ j = | ˆ f · ˆ r ij |{ ( ρ i ξ i − ρ j ξ j )+ φ [ h i h i + h j r ji · ∇ ( ρ i ξ i ) − h j h i + h j r ij · ∇ ( ρ j ξ j )] } , (31)where φ is the slope limiter (implemented using the minmod function, minmod ( x )= max (0, min ( x ,1)) to minimize spuri-ous oscillations. The term ˆ f · ˆ r ij limits unwanted dissipa-tion perpendicular to direction of the flux. We find that for α ξ = 1 , discontinuity-capturing is good while dissipation issmall in smooth regimes. A discontinuity ‘detector’for artificial diffusion is therefore not required.Chow & Monaghan (1997) suggested turning off artifi-cial dissipation when ( ρ i f i − ρ j f j ) · ˆr ij > , in order to re-duce unnecessary diffusion, e.g. behind a discontinuity. How-ever in our experiments we found that such a switch makesthe scheme unstable, in particular in cases where radiationbeams collide in the optically thin regime. The instabilityresults in significant non-conservation of energy. We there-fore do not use the Chow & Monaghan (1997) switch, butinstead apply a discontinuity detector to minimize the arti-ficial viscosity as described in §2.5.1.The artificial flux dissipation term of Eq. (28) causesnumerical dissipation of the flux in directions perpendicularto the flux. This is problematic, in particular in the opti-cally thin case where it leads to the destruction of a packetof radiation as shown in Fig.1. To avoid this requires thatany artificial flux dissipation should be in the direction ofthe flux itself. Simply multiplying the right hand side ofEq. (28) by ˆ f · ˆ r ij does not work: any component of numeri-cal flux perpendicular to the actual flux, for example due tonumerical noise, will still lead to the artificial destruction ofan optically thin radiation packet .A better solution is to implement the dissipation schemeas as anisotropic diffusion .Here we outline our default choice of artificial flux dis-sipation . We begin by rewriting Eq. (28) in the form of ananisotropic diffusion equation, (cid:18) D f D t (cid:19)(cid:12)(cid:12)(cid:12)(cid:12) diss = 1 ρ ∇ · (cid:104) D f ∇ · ( ρ f ) (cid:105) , ≡ ρ ∇ · (cid:104) ρ D f ψ (cid:105) , (32)where the tensor D f and the scalar ψ are given by D f = α f v sig h ˆ n ˆ n ψ ≡ ρ − ∇ · ( ρ f ) . (33) In the optically thick regime, χ already provides the necessarydissipation. This is one of the reason why the flux-limited diffusiondoes not require artificial dissipation. Another reason is that thereare no artificial oscillations when solving a diffusion equation). It is possible to use Eq. 28 without disrupting radiation direc-tions if the optically thin direction is imposed, as in §2.4.1. Inthis case, we will only consider the dissipation component (in Eq.28) along the optically thin direction, and only consider the fluxdifference in that direction. ‘Anisotropic artificial viscosity’ was also considered by Owen2004, but our SPH form is simpler and different from theirs.© 0000 RAS, MNRAS , 000–000 moothed Particle Radiation Hydrodynamics—SPH-M1RT We implement the diffusion equation in SPH as ψ i = − (cid:88) j m j Ω i ρ i ( ρ i f i − ρ j f j ) · ∇ i W ij ( h i ) , (34) (cid:18) D f D t (cid:19)(cid:12)(cid:12)(cid:12)(cid:12) i, diss = − (cid:88) j m j Ω i ρ i (cid:16) ρ i ψ i D f i − ρ j ψ j D f j (cid:17) · ∇ i W ij ( h i ) . (35)This formulation of anisotropic viscosity in SPH is noveland we suggest that it may be applicable to other situ-ations as well, for example when implementing magneto-hydrodynamics or cosmic ray propagation. Clearly it would be advantageous to activate flux dissipationonly near discontinuities in the radiation, which requires effi-cient detection of such discontinuities. Such switches are alsoregularly used to activate dissipation in the SPH equationsfor hydrodynamics itself.Morris & Monaghan (1997) proposed to use the diver-gence of the velocity as a measure of how discontinuous thefluid flow evolves, but this cannot distinguish compression -which conserves entropy - from true discontinuities. In addi-tion, flux dissipation may be activated unnecessarily in thecase of wave-like disturbances. Rosswog (2020b) suggestedto use changes in entropy as a discontinuity detector, how-ever it is not clear how to apply this in the case of radiation.(Cullen & Dehnen 2010) suggested to track the time deriva-tive of the velocity divergence, ∇ · v , so that the diffusioncoefficient is of the form h | ˙ ∇ v | /v , where v sig is the sig-nal velocity. This effectively corresponds to a switch that isbased on the second time-derivative of the density and hencecan distinguish between gas in the pre- and post-shock re-gions. Such a switch is implemented in SPHENIX .Inspired by the Cullen & Dehnen (2010) switch andafter experimenting with various forms of how their expres-sion can be applied to the case of radiation, we settled onthe following target value for the diffusion coefficient, α f , aim = min (cid:18) max (cid:26) − Ah ρξ ˜ c D[ ∇ · ( ρ f )]D t , (cid:27) , (cid:19) . (36)The denominator is ρξ ˜ c rather than ρ | f | ˜ c because ρcξ (cid:29)| ρ f | in the optically thick limit where artificial dissipation isnot needed. A (= 200) is a constant multiplication factor tocompensate for the large ˜ c in the equation.Upstream from a discontinuity, we require that the dif-fusion coefficient be large enough so that the discontinuitycan be captured by the interpolation scheme. Downstreamfrom the discontinuity, a smaller level of diffusion is stillrequired to suppress any numerical oscillations. We followMorris & Monaghan (1997) and implement this by makingthe diffusion coefficient time-dependent, as follows: (1) when α (cid:54) α aim , we set α = α aim ; (2) when α (cid:62) α aim , we evolve α back to α aim by solving DD t ( α − α aim ) = − τ relax = − (˜ c/h + ˜ cχρ ) , (37)where τ relax is the relaxation time scale. The term ˜ cχρ en-sures that a large value for the α quickly relaxes back tothe target value in the optically thick yet smooth regime.Finally, for gas particles in which we inject radiation we set Figure 1.
A packet of radiation propagating upwards (from x = 0 to x = 2 ) over a 2D glass-like distribution of × particlesshown at time t = 0 . ; the reduced speed of light is ˜ c = 10 andthe extent of the simulation volume is ∆ x = 2 in the verticaldirection and ∆ y = 0 . in the horizontal direction. Colours rep-resent the radiation energy density ξ and small red-arrows theradiation flux density, f . Panels from left to right illustrate the default , isomax , none and original choices for the artificial dis-sipation and closure scheme (see text). The optical depth of themedium is zero and the packet should propagate freely at thespeed of light while retaining its square form. With isotropic arti-ficial dissipation ( isomax ) or with the original closure scheme, theradiation packet incorrectly dissolves quickly. The case withoutartificial dissipation ( none ) results in strong oscillations whichleads to significant non-conservation of energy. The default choicecorrectly maintains the morphology of the beam while suppress-ing artificial numerical oscillations. α = 1 to better capture any discontinuities associated withradiation injection.Before ending this section, we comment on the requirednumber of SPH neighbour loops associated with our RTscheme. If the radiation moment and hydrodynamics equa-tions (Eqs. 2-6) are solved simultaneously, then at least threeneighbour loops are required to compute the right-hand-sizeof the anisotropic diffusion equation, Eq. (32): the variable ψ in Eq. (33) requires (1) a loop to compute ρ and (2) a sec-ond loop to compute the gradient; and finally the scheme re-quires (3) a third loop to compute the gradient of ψ Eq. (32).Similarly, the dissipation switch of Eq. (36) requires threeloops. The scheme may be optimized by solving the radia-tion transport equation on a shorter time-scale than usedto update the hydrodynamics. During such sub-cycling, thedensity is kept a constant, in which case the radiation trans-port only requires two SPH neighbour loops. We will reporton such improvements elsewhere. © 0000 RAS, MNRAS000
A packet of radiation propagating upwards (from x = 0 to x = 2 ) over a 2D glass-like distribution of × particlesshown at time t = 0 . ; the reduced speed of light is ˜ c = 10 andthe extent of the simulation volume is ∆ x = 2 in the verticaldirection and ∆ y = 0 . in the horizontal direction. Colours rep-resent the radiation energy density ξ and small red-arrows theradiation flux density, f . Panels from left to right illustrate the default , isomax , none and original choices for the artificial dis-sipation and closure scheme (see text). The optical depth of themedium is zero and the packet should propagate freely at thespeed of light while retaining its square form. With isotropic arti-ficial dissipation ( isomax ) or with the original closure scheme, theradiation packet incorrectly dissolves quickly. The case withoutartificial dissipation ( none ) results in strong oscillations whichleads to significant non-conservation of energy. The default choicecorrectly maintains the morphology of the beam while suppress-ing artificial numerical oscillations. α = 1 to better capture any discontinuities associated withradiation injection.Before ending this section, we comment on the requirednumber of SPH neighbour loops associated with our RTscheme. If the radiation moment and hydrodynamics equa-tions (Eqs. 2-6) are solved simultaneously, then at least threeneighbour loops are required to compute the right-hand-sizeof the anisotropic diffusion equation, Eq. (32): the variable ψ in Eq. (33) requires (1) a loop to compute ρ and (2) a sec-ond loop to compute the gradient; and finally the scheme re-quires (3) a third loop to compute the gradient of ψ Eq. (32).Similarly, the dissipation switch of Eq. (36) requires threeloops. The scheme may be optimized by solving the radia-tion transport equation on a shorter time-scale than usedto update the hydrodynamics. During such sub-cycling, thedensity is kept a constant, in which case the radiation trans-port only requires two SPH neighbour loops. We will reporton such improvements elsewhere. © 0000 RAS, MNRAS000 , 000–000 T. K. Chan et al.
A test of the artificial dissipation scheme and the choiceof the Eddington tensor closure by propagating a singleshort beam of light is illustrated in Fig. 1. The underlying × particle distribution is glass-like. The figure showsfour choice of the artificial dissipation and closure relationlabelled default , isomax , none , and original : • default : uses the default artificial dissipation of Eqs. 27,31, and 32 described in §2.5.1, and the default modified clo-sure M1 scheme of Eq. (19); • isomax : uses isotropic artificial dissipation of Eqs. 27and 28, setting α f = α ξ = 1 , see Eq. (29), and the default modified closure scheme; • none : does not use any artificial dissipation (i.e. α f = α ξ = 0 in Eq. 29) and the default closure scheme; • original : uses the default artificial dissipation but the‘ original ’ closure Eq. (17).Figure 1 demonstrates that artificial diffusion is needed tosuppress the artificial oscillations seen in panel none : os-cillations lead to non-physical negative values of ξ , which,if zeroed, lead to a catastrophic artificial increase in radia-tion energy. However, such diffusion should be anisotropic toavoid that the beam artificially diffuses perpendicular to thepropagation direction as seen in panel isomax . The original scheme fails to preserve the beam’s shape: it does not handlewell non-uniform particle distributions. Fortunately, our de-fault scheme preserves the beam shape, suppresses artificialoscillations and conserves energy. In this section we briefly describe how we implement theinteraction between matter and radiation, limiting the dis-cussion to the particular case of ionizing radiation in a purehydrogen gas.
The processes of collisional ionization and photoionization,photoheating, and collisional and radiative cooling in a hy-drogen gas are (e.g. Aubert & Teyssier 2008): ∂n γ ∂t = ρ (Λ ξ + S ξ )2 π (cid:126) ¯ ν = − n HI ˜ cσ γ n γ + n e n HII ( α A − α B ) + S γ , (38) ∂ f γ ∂t = − χρ ˜ c f γ + ρ S f π (cid:126) ¯ ν = − n HI ˜ cσ γ f γ + S γ , (39) ∂n HI ∂t = − n HI ˜ cσ γ n γ + n e n HII α A − n e n HI β , (40) ∂e tot ∂t = ρ ∂u∂t = ρ ( S u + Λ u )= (cid:15) γ n HI ˜ cσ γ n γ − n HI n e Γ e HI − n HII n e Γ e HII . (41)In addition to the variables described in section 2.2, n X isthe number density of element species X ( i.e neutral andionised hydrogen, and electrons), e tot ≡ ρu = n tot k B T isthe thermal energy density of the fluid, n tot = ρ/ ( µm H ) = n HI + 2 × n HII is the total number density of all parti-cles, µ = 1 / (2 − x ) is the mean molecular weight and x HI ≡ n HI /n H is the neutral fraction (sometimes we justrefer the neutral fraction as x for simplicity when there isno confusion). Furthermore, α A and α B are the ‘case A’ and‘case B’ recombination coefficients, respectively, β is the col-lisional ionization coefficient, σ γ is the photoionization crosssection, Γ e HI is the collisional cooling rate per unit volume(combining collisional ionization and collisional excitation), Γ e HII represents the sum of recombination cooling and ther-mal Bremsstrahlung, and (cid:15) γ is the photoionization heatingenergy per ionization. The values of these coefficients, to-gether with their dependence on photon frequency, ν , and/orgas temperature, T , are listed in Appendix A. Finally, S ξ and S γ are external sources of photons.The above set of differential equations is in general nu-merically stiff, meaning that the numerical solution is un-stable unless the equations are integrated in time using avery short timestep, ∆ t . The reason is that the coefficients inthese equations are large in some situations, e.g. n HI ˜ cσ γ n γ islarge in the neutral region near radiation sources. The usualremedy is to use an implicit scheme because this is stable,however its solution may not be sufficiently accurate. Ourstrategy described below is to combine explicit and implicitmethods. To illustrate the solution method we make the ‘on-the-spot’approximation by assuming that recombinations directly tothe ground state produce an ionizing photon that is ab-sorbed close to where it was emitted ( i.e. ‘on the spot’).In this approximation we set α A = α B , resulting in the fol-lowing set of three coupled differential equations, ∂n γ ∂t = − xn H ˜ cσ γ n γ , (42) ρ ∂u∂t = (cid:15) γ xn H ˜ cσ γ n γ − n x (1 − x )Γ e HI − n (1 − x ) Γ e HII , (43) ∂x∂t = − x ˜ cσ γ n γ + (1 − x ) n H α B − x (1 − x ) n H β, (44)where x = n HI /n H (= x HI ) is the neutral hydrogen fraction.The partial time derivatives refer to changes due to in-teraction between radiation and gas only. There may be ad-ditional terms, for example, due to other photon sources orsinks, and heating and cooling due to adiabatic processes orshocks. Here, we restrict ourselves to solving these radiativeequations, treating any other source/sink terms in operatorsplit fashion.We integrate these equations following the approach ofPetkova & Springel (2009): solve the first two equations ex-plicitly and use that solution to solve the third equation(the chemistry equation) implicitly. However, we addition-ally perform subcycling , requiring that n γ and u do notchange by more than 10% in each subcycle.We do so by requiring that ∆ t (cid:54) . / ( xn H cσ γ ) , u/ | ∂u/∂t | ) . While the implicit solverfor the neutral fraction x is unconditionally stable, it can © 0000 RAS, MNRAS , 000–000 moothed Particle Radiation Hydrodynamics—SPH-M1RT Figure 2.
A schematic diagram illustrating the subcyclingmethod for solving non-equilibrium thermochemistry equations.In the main loop, we advance the radiation, hydrodynamics, grav-ity, and radiation injection equations, and any other equations notrelated to radiation ( e.g. subgrid schemes). In between we sub-cycle the thermochemistry equations (Eqs.42, 43 & 44), with asmaller timestep. Note that the whole set of the thermochemistryequations is solved sequentially within one subcycle, accuratelyaccounting for any rapid changes of ionization state or photondensity. be inaccurate if the time-step is too large. Therefore, wefurther limit the subcycle timestep to ∆ t (cid:54) C x/ | ∂x/∂t | ,where C is a parameter we choose to be 0.1 (but can belarger depending on the tolerance). So in summary, we takethe subcycling step to be ∆ t = 0 . n γ | ∂n γ /∂t | , u | ∂u/∂t | , x | ∂x/∂t | ) . (45)Fig. 2 illustrates the semi-implicit subcycle scheme.We update n γ using the analytic solution of Eq. (42) incase x and n H are held constant, n γ ( t n +1 ) = n γ ( t n ) exp( − ∆ tσ γ cn H x ) . (46)Doing so guarantees that n γ is always positive, as it shouldbe, and that the solution is asymptotically correct when re-combinations are negligible. We update u using the corre-sponding analytical solution of Eq.43.The implicit solution to the chemistry equation uses theupdated values for the radiation and internal energies, x ( t n +1 ) − x ( t n )∆ t = − ˜ cσ γ n γ ( t n +1 ) x ( t n +1 )+ n H α B ( t n +1 )[1 − x ( t n +1 )] − n H β ( t n +1 ) x ( t n +1 )[1 − x ( t n +1 )] . (47)which is a quadratic equation for x ( t n +1 ) . This method can be generalised to the case of moreelements by adopting the approach of Anninos et al. (1997), i.e. updating each element implicitly one by one in orderof increasing timescale. A possible alternative scheme usesthe CVODE library (Cohen et al. 1996) to solve these stiffequation, e.g.
Kannan et al. (2019).For now, we have restricted this discussion to the on-the-spot approximation. This limitation can be relaxed byadding recombination radiation as a source terms to each gasparticle. Given that the computing time of our RT scheme isindependent of the number of sources, this could be feasiblyimplemented in the future.Our semi-implicit subcycling scheme is accurate, as weshow below, as well as computationally efficient, as shown inAppendix B. Subcycles are initiated only when the systemis out of equilibrium, for example when a source of pho-tons is suddenly switched on. But even in such situations,we find that only a few dozen subcycles occur. Subcyclingshould also help with load balancing the computation. With-out subcycling, most of the computing time will be spent inthe vicinity of ionization fronts, where the thermochemistryis highly out of equilibrium. With subcycling enabled, thetime-step of the main loop is instead limited by the over-all CFL timestep. We proceed by showing some tests of thethermochemistry implementation.
This test is a variation of Test 0 in Iliev et al. (2006): aninitially neutral parcel of pure hydrogen gas at low temper-ature is suddenly ionized and heated by a source of ionizingradiation with a specified spectrum of ionizing photons. Af-ter a specified time, the ionizing source is switched off. Thetotal density of the gas parcel is kept constant. The radia-tion is assumed to be optically thin at all times. The testinvolves following the evolution of the neutral fraction, x ,and of the temperature of the gas, T . To enable a fair com-parison between codes, it is, of course, important to makesure that the physical constants used - such as, for example,the frequency dependence of the ionization cross section -are the same.As the source is switched on and the hydrogen gas getsionized, the temperature increases to a value that dependson the shape of the ionizing spectrum. The gas is then inionization equilibrium. However, it takes longer for the gasto be in thermal equilibrium - where photoheating balancesradiative cooling. When the source is switched off, the gasstarts to recombine and cool. We do not include molecule for-mation in the calculation, and hence the cooling rate dropswith decreasing T . • Analytical description
The evolution can be understoodby writing Eq. (44) as a Ricatti equation, dxdt = − xτ i + (1 − x ) τ r − x (1 − x ) τ e ≡ τ ( x − x )( x − x ) , (48)where we defined three characteristic time-scales, τ i ≡ cσ γ n γ ; τ e ≡ n H β ; τ r ≡ n H α B , (49) © 0000 RAS, MNRAS000
The evolution can be understoodby writing Eq. (44) as a Ricatti equation, dxdt = − xτ i + (1 − x ) τ r − x (1 − x ) τ e ≡ τ ( x − x )( x − x ) , (48)where we defined three characteristic time-scales, τ i ≡ cσ γ n γ ; τ e ≡ n H β ; τ r ≡ n H α B , (49) © 0000 RAS, MNRAS000 , 000–000 T. K. Chan et al. and τ = 1 τ e + 1 τ r ; x + x = 2 + τ r τ i ; x x = ττ r . (50)Choosing initial condition x = x at t = 0 , the generalsolution in case all τ ’s are constant, is x ( t ) = x ( x − x ) − x ( x − x ) f ( t )( x − x ) − ( x − x ) f ( t ) f ( t ) = exp[ − ( x − x ) t/τ ] . (51)In the special case where collisional ionizations are ne-glected, τ e → ∞ and when τ i (cid:28) τ r , this solution simplifiesto approximately x ( t ) = x exp( − t/τ i ) + τ i /τ r : the neutralfraction approaches its ionization equilibrium exponentiallyon the ionization timescale, τ i .This approximate description assumes that τ r - andhence the temperature of the gas - remains a constant. Anestimate of the change in temperature following rapid ion-ization, τ i (cid:28) τ r , follows from neglecting radiative cooling inthe short time it takes to ionize the gas, so that ρ dudt ≈ − (cid:15) γ dn HI dt = − (cid:15) γ n H dxdt , (52)where (cid:15) γ is the mean energy injected into the gas perphotoionization, see §A. Writing the thermal energy perunit mass, u , in terms of the neutral fraction, x , as u = k B T / ( µm H ) with µ = (2 − x ) − , yields the following rela-tion between the initial temperature T = T at t = 0 , andthe temperature T when the neutral fraction is x : k B T = 2 − x − x k B T + 2 (cid:15) γ x − x − x . (53)In the test described below, the gas is initially neutral, x = 1 , and at low temperature, k B T (cid:28) (cid:15) γ , and the pho-toionization rate is high so that in photoionization equilib-rium, x (cid:28) . In this case, the photoionization equilibriumtemperature is T ≈ (cid:15) γ / (3 k B ) , which depends only on theionization energy per photon, regardless of radiation inten-sity or gas density.On a longer timescale, the parcel of gas will reach ther-mal equilibrium (temperature T ), where photoheating bal-ances radiative cooling. Provided T > T , the timescale toreach this equilibrium can be estimated by simply neglectingcooling and noting that the rate at which the gas is heatedis approximately the product of the ionization rate, x/τ i ,times the energy injected per photoionization per hydrogenatom, (cid:15) γ /m H , hence d u d t = (cid:15) γ m H xτ i ≈ (cid:15) γ m H τ r . (54)Therefore it takes approximately a recombination time toreach thermal equilibrium (see also Pawlik & Schaye 2011).When the source is suddenly switched off, gas will startto recombine. This is still described by a Ricatti equation ofthe form of Eq. (51), except that now x ( t ) = x ( x eq − x ) − x ( x eq − x ) f ( t )( x eq − x ) − ( x eq − x ) f ( t ) f ( t ) = exp[ − ( x − x )( t − t eq )) /τ ] , (55)where x eq is the neutral fraction in thermal equilibrium, t eq is the time that the ionizing source is switched off, and x + 1 /x = 1 + τ /τ r with x = 1 /x . If the gas is no longer heated, it will of course simply keep on cooling and there isno further equilibrium state. • Numerical solution
For the numerical values for thistest, we take ˜ c = c , n H = 1 cm − , x = 1 and T = 100 K .From time t = 0 , the parcel of gas is being irradiated witha black body spectrum of temperature T = 10 K (for which (cid:15) γ = 6 .
33 eV in the optically thin limit) and photon flux F photon = 10 photons s − cm − . The source is switchedoff after a time t s = 5 × yr and we follow the evolutionuntil time t e = 10 yr .For a reference temperature of T = 10 K , the threetimescales are τ i ≈ − . yr , τ r ≈ . yr and τ e ≈ . yr ,so that τ i (cid:28) τ r ≈ τ (cid:28) t e < τ e . The temperature in pho-toionization equilibrium is T = 6 .
33 eV / (3 k B ) ≈ . K ,and the temperature in thermal equilibrium is about twicethat. The timescale for the gas to reach thermal equilibriumcan be estimated as follows. The heating rate of the gas,when in photoionization equilibrium, is du/dt ≈ (cid:15) γ /τ r ≈ . × − erg Myr − . Therefore the timescale to reach ther-mal equilibrium is approximately t eq ≈ ∆ udu/dt = u eq − u i du/dt ≈ u i du/dt ≈ . yr , (56)where we have used the fact that for the parameters of thistest, u eq ≈ u i . In summary: the gas should reach its pho-toionization equilibrium temperature by a time τ i , reach itsthermal equilibrium temperature by the time t eq , and startto cool and recombine after time t s .We want to verify that the combination of explicit sub-cycling and implicitly solving the chemistry equations yieldsthe correct solution, independently of a globally imposedtimestep. To demonstrate the accuracy of the integrationscheme, we also want to compare to a run in which we inte-gration the equations with a short, fixed timestep. However,the ionization timescale τ i is much smaller than the evolu-tion timescale t e , and it is impractical to simulate the wholetime evolution with a timestep much shorter than τ i . Herewe follow Pawlik & Schaye (2011) and perform a dozen sim-ulations with different (fixed) timesteps, from ∆ t (cid:28) τ i to ∆ t (cid:29) τ i .Results are shown in Fig. 3, where differently colouredcurves show the evolution for different values of the globaltimestep, ∆ t . All curves follow the analytical expectation:gas heats and gets almost fully ionized on a timescale τ i (vertical dashed red line), reaching its photoionization equi-librium temperature (horizontal dashed line), continues tobe heated on a timescale τ r (vertical dashed black line) toreach thermal equilibrium (horizontal dashed black line),and finally starts to recombine and cool when the sourceis switched off. All simulation runs fall on top of each other,demonstrating that the numerical solution is independent ofthe global fixed value of ∆ t .After the radiation is switched off, gas recombines andcools rapidly to T ∼ K , below which the cooling ratedrops rapidly as we only include cooling by neutral hydro-gen. There are two major takeaways from Fig. 3. Firstly,the simulated evolution follows the analytical expectation Note that the value of ˜ c is irrelevant to the solution since thepure thermochemistry equation is independent of ˜ c if F photon isgiven. © 0000 RAS, MNRAS , 000–000 moothed Particle Radiation Hydrodynamics—SPH-M1RT Figure 3.
A variation of Iliev et al. (2006) Test 0: photoheat-ing of a single gas parcel irradiated with black body radiationof temperature T = 10 K . The source is switched off aftertime t s = 5 × yr . Vertical dashed lines indicate the ionizationtimescale τ i ( red ) and the time to reach thermal equilibrium, t eq ( black ), horizontal dashed lines indicate the expected photoion-ization equilibrium temperature (or neutral fraction) ( red ) andthe thermal equilibrium temperature ( black ). Points with differ-ent colours show the evolution computed with different, fixed,global timesteps, as per the legend. Simulations with different ∆ t step sizes yield the same evolution, which also matches the an-alytical estimate a well as the curves in Iliev et al. (2006) andPawlik & Schaye (2011) (not shown here). as well as the results from other simulation codes in Ilievet al. (2006) (see their Fig. 5), therefore the scheme isaccurate. Secondly, the evolution is independent of theglobal timestep, demonstrating the convergence of the sub-cycle thermochemistry solver. With this solver, the (main)timestep of the simulation is only limited by the Courantcondition in solving the moment equation (§2.2). Radiation is injected by star particles. In our SPH imple-mentation, star particles have a smoothing length h , whichis calculated in the same way as that of gas particles, byrequiring that each star particle overlaps with a specifiednumber of weighted gas neighbours. A star particle i , whichhas a given energy injection rate, ˙ e i, rad , transfers during atime-step ∆ t i an amount of radiation ∆ e ij into its gas neigh-bour j , where ∆ e ij = m j ∆ ξ j = m j N nor ρ j r ij ˙ e i, rad ∆ t i . (57)This kernel-weighted transfer is normalized by N nor , suchthat (cid:80) j ∆ e ij = ˙ e i, rad ∆ t i , summing over all neighbor parti-cles j .We simultaneously inject isotropic radiation flux, m j ∆ f j , as if the surrounding medium were optically thin: ∆ f j = ˜ c ˆ r ji ∆ ξ j . (58) Because the distribution of gas neighbours around asource is generally relatively disordered, the resulting radia-tion field may not be isotropic unless energy is injected overa sufficient number of gas particles. Therefore we increasethe smoothing lengths of star particles to be a few timesthe smoothing lengths of gas particles, e.g. h star = 2 h gas (see Fig. 7). An alternative way of ensuring isotropic radia-tion around sources is to impose the radiation direction inthe optically thin limit (§2.4.1). In some tests, e.g. tests ofStrömgren spheres, we calculate the total radiation energywithin the injection region, and then reset the radiation dis-tribution according to the optically thin expectation. Our RT scheme is implemented in the public version ofthe
SWIFT (SPH with interdependent fine-grained tasking)(Schaller et al. 2016) code, which has been applied ingalaxy formation and planetary giant impacts (Kegerreiset al. 2019). The target application of SWIFT is large-volumecosmological simulations with subgrid physics modules sim-ilar to eagle (Schaye et al. 2015).
SWIFT is an SPH code that solves cosmologi-cal or non-cosmological hydrodynamic equations, includ-ing self-gravity, and is designed to work on hybridshared/distributed memory computer architectures. Loadbalance is optimised using task-based parallelism, with tasksassigned by a graph-based domain decomposition, and usingdynamic, asynchronous communication. For hydrodynam-ics, Borrow et al. (2018) found
SWIFT to have good weakscaling from 1 to 4096 codes (losing only 25% performance)in low redshift cosmological galaxy simulations (with
EA-GLE physics from Schaye et al. 2015).
When radiation travels at the speed of light, the timestepto advance a radiation front correctly is of order ∆ t c ∼ h/c ,for a smoothing length h of an SPH particle. This is, ofcourse, much shorter than the CFL step, which is of order ∆ t s ∼ h/v s , where v s is the sound speed. However, ionisingradiation with flux F moves at the speed of the ionizationfront, v I ∼ F/ (2 πn H (cid:126) ν ) , through neutral gas with density n H . When v I (cid:28) c , the code can be speeded up by a largefactor by ‘reducing’ the speed of light, from c to ˜ c . As longas ˜ c > v I , the speed of an ionization front can still be correctfor a given F (see e.g. the discussion in Rosdahl et al. 2013).This ‘reduced speed of light’ (RSL) approximation wasintroduced by Gnedin & Abel (2001) to simulate radiativetransfer efficiently and has been applied to other radiativetransfer simulations, e.g. by Aubert & Teyssier (2008). Theydemonstrate that RSL performs well in problems involvingionization, photoheating, and expansion of HII regions.However, there is no unique way to implement RSL.Skinner & Ostriker (2013) implemented the RSL approx-imation in simulating RT in the interstellar medium, e.g.modeling radiation reprocessed by dust. However, their ap-proach does not conserve total radiation plus matter energy http://swift.dur.ac.uk/ © 0000 RAS, MNRAS000
When radiation travels at the speed of light, the timestepto advance a radiation front correctly is of order ∆ t c ∼ h/c ,for a smoothing length h of an SPH particle. This is, ofcourse, much shorter than the CFL step, which is of order ∆ t s ∼ h/v s , where v s is the sound speed. However, ionisingradiation with flux F moves at the speed of the ionizationfront, v I ∼ F/ (2 πn H (cid:126) ν ) , through neutral gas with density n H . When v I (cid:28) c , the code can be speeded up by a largefactor by ‘reducing’ the speed of light, from c to ˜ c . As longas ˜ c > v I , the speed of an ionization front can still be correctfor a given F (see e.g. the discussion in Rosdahl et al. 2013).This ‘reduced speed of light’ (RSL) approximation wasintroduced by Gnedin & Abel (2001) to simulate radiativetransfer efficiently and has been applied to other radiativetransfer simulations, e.g. by Aubert & Teyssier (2008). Theydemonstrate that RSL performs well in problems involvingionization, photoheating, and expansion of HII regions.However, there is no unique way to implement RSL.Skinner & Ostriker (2013) implemented the RSL approx-imation in simulating RT in the interstellar medium, e.g.modeling radiation reprocessed by dust. However, their ap-proach does not conserve total radiation plus matter energy http://swift.dur.ac.uk/ © 0000 RAS, MNRAS000 , 000–000 T. K. Chan et al. and momentum, and the non-equilibrium solution might notbe correct.Ocvirk et al. (2019) examined the ‘dual speed of light’(DSL) approximation, where c → ˜ c in the propagation equa-tion but not in the thermo-chemistry equations. Unfortu-nately, DSL fails to reproduce the correct equilibrium gasproperties. This is because when c is reduced to ˜ c in thepropagation equation, the photon-matter interaction ratedoes not change accordingly in DSL. This can be seen byconsidering the analytical solution of the Strömgren sphere(Appendix C), n HI cσ γ n γ = (cid:16) c ˜ c (cid:17) n HI σ γ ˙ N γ πr exp (cid:18) − (cid:90) r n HI σ γ d r (cid:19) = n e n HII α B , (59)where the factor ˜ c arises from the propagation equation inthe optically thin limit, and the factor c comes from thethermo-chemistry equation. The equilibrium neutral frac-tion will deviate from the correct solution due to the c/ ˜ c factor.Thus, our default treatment is to replace c → ˜ c in allequations (Aubert & Teyssier 2008; Rosdahl et al. 2013),including the propagation and thermo-chemistry equations.For a fixed photon flux, F , (or photon injection rate), thischoice reduces the interaction strength between light andmatter (e.g. σ γ ˜ c ) to compensate for higher photon density(due to the slower photon propagation speed). As a result,the photoionization rate will be independent of ˜ c (as long as ˜ c is larger than other speeds). Furthermore, the choice of ˜ c will not affect the equilibrium gas properties, as demonstrated inEq. 59 (and see the tests in next sections).However, there are limitations to RSL. First, ˜ c shouldexceed the speed v I of any ionization front. For example,Bauer et al. (2015) showed that using ˜ c = c/ affects thetiming of reionization. Another issue of RSL is that using ˜ c < c increases the momentum term ∇ ( F E ) , and if this isnot corrected for then the radiation pressure will be too large(see also Jiang et al. 2012; Jiang & Oh 2018). This may beproblematic in cases where radiation pressure is crucial, forexample when modelling radiation pressure from AGN. Inthe case of reionization simulations, the photon-density islow and radiation pressure is mostly neglected anyway. This section contains an extensive series of tests to vali-date the numerical scheme and its implementation in the swift code. The tests combine the default scheme for radia-tion (§2.6) with the
SPHENIX
SPH formulation for hydrody-namics (§2.4), unless explicitly stated otherwise. Some testimpose the optically thin direction in the Eddington ten-sor (§2.4.1). Otherwise, the flux propagates in the direction ˆ n = ˆ f , as computed for each gas particle. We will make theon-the-spot approximation in all of the tests (§2.7.2). We donot use the reduced speed of light (RSL) approximation in§3.1, since we are interested in the radiation distribution.But we use RSL in §3.2-3.4, where we are concerned withgas properties (see §2.10). Figure 4.
Propagation in 1D: radiation energy density, ξ ( x ) of apackage of radiation propagating to the right at the speed of light( c = 1 in code unit). The black dashed line is the initial profilein units of the initial value of E , the blue line is the edge of theinitial profile shifted to the right by ∆ x = 5 , the red points are thesimulated values of E for individual SPH particles at time t = 50 .The simulation uses 400 SPH particles located on the vertices ofa regular grid. The test shows that the front moves at the correctspeed. The setup propagates a finite radiation packet in an opticallythin medium. We test whether the radiation front travels atthe correct speed, c = 1 , without excessive smoothing of thefront and without generating artificial oscillations in the ra-diation density, E ( x ) , behind the radiation packet. Initially,the radiation energy density and flux are uniform and non-zero only for x < , with the radiation flux is pointing inthe + x direction initially. Fig. 4 shows the initial conditionand the configuration at t = 5 , when the radiation front haspropagated 100 times the mean interparticle spacing. Whilethere is small broadening of the radiation front caused by theartificial dissipation, numerical oscillations are suppressedsignificantly and the scheme is stable. The front propagatesat the correct speed ( c = 1 ) and the radiation energy density E remains constant inside the radiation package, unaffectedby the artificial dissipation. We repeat the previous test but now in two-dimensions: arectangular radiation package propagates in empty space in2D. This tests the extent to which radiation leaks out ofthe package artificially, either perpendicular or parallel tothe propagation direction, as a consequence of the artificialdissipation. The propagation direction of the radiation onindividual particles is not imposed, but computed from ˆ n =ˆ f . Results are shown in Fig. 5 where the radiation packagemoves upwards from the initial state at time t = 0 (leftpanel) to time t = 0 . (right panel).The top panel of the figure demonstrates the ability of © 0000 RAS, MNRAS , 000–000 moothed Particle Radiation Hydrodynamics—SPH-M1RT Figure 5.
Propagation in 2D: propagation of a packet of radia-tion in two dimension at the speed of light, c = 1 . Upper panels :Radiation energy density of the packet in units of its initial value.The left panel shows the initial state, where the packet has length δy = 1 . Panels to the right show the system at different times. TheSPH particle distribution is glass-like, with × particles fill-ing the computational volume of horizontal extent ∆ x = 0 . andvertical extent ∆ y = 2 . . Colors represent the radiation energydensity and red arrows show the radiation fluxes. Lower panel : ra-diation energy density, ξ ( x, y ) , of SPH particles with . (cid:54) y (cid:54) . at time t = 1 . ( red line ), the dashed-black line is the initial shapeof the packet. the implementation to maintain the direction of the radia-tion, with little artificial leakage of radiation perpendicularto the beam. In the right panel, the packet has propagatedupwards over 128 mean particle spacings. Two small ‘tails’of radiation trail the package, where radiation leaked out ofthe beam. Smoothing perpendicular to the beam is quan-tified in more in detail in the lower panel, which is a cutthrough the middle of the beam at time t = 1 . This profilehas approximately Gaussian-shaped edges, as expected from artificial diffusion (§2.5); the diffusion coefficient is propor-tional to h c . The dependence on the smoothing length, h ,means that the beam can propagate further without distor-tion at higher resolution. Due to the finite resolution and, ingeneral, non-uniform underlying SPH particle distribution,it is not possible to completely eliminate radiation leakageperpendicular to the propagation direction without causinginstabilities; higher-order shock-capturing schemes (e.g. Liuet al. 1994) might suppress such artificial leakage more effi-ciently.The next test is that of radiation propagating isotropi-cally away from a source in two dimensions, see Fig. 6. Thefigure confirms that the radiation front preserves rotationalsymmetry as it propagates out at the speed of light. The ra-diation energy density is smooth in the radial shell, with noappreciable noise even though the underlying particle distri-bution is non-uniform. The energy density is small behindthe shell. The absence of significant artificial ‘left over’ ra-diation results from the artificial dissipation switch. The first test is Test 1 in Iliev et al. (2006). This tests theradiative transfer scheme and the thermo-chemistry solveragainst an analytical solution: uniform density, neutral gasis photoionized by a source that emits ionizing photons ata constant rate. We keep the density and temperature ofthe gas constant, i.e. the gas is not allowed to move, heator cool. The ionization front propagates into the gas clouduntil it reaches its Strömgren radius. The analytical solu-tion (assuming grey opacity) is derived and summarised inAppendix C.The numerical parameters are taken to be identicalto those used by Iliev et al. (2006) to allow for a directcomparison: the gas cloud consists of pure hydrogen gaswith density n H = 10 − cm − , the collisional ionization co-efficient β = 3 . × − cm s − and the recombinationcoefficient is α = 2 . × − cm s − (with the on-the-spot approximation); the photoionization cross section is σ γ HI = 8 . × − cm . The source emits ionizing radiationat a constant rate of ˙ N γ = 5 × photons s − . The compu-tational volume has linear extent of 20 kpc; the SPH particledistribution is glass-like with approximately particles.In this problem we also test the RSL approximation, using ˜ c = c/ .We first compare several implementations of the injec-tion of radiation energy by the source (§2.8) and of the op-tically thin closure relation (§2.4.1); results are shown inFig. 7. Without requiring that the optically thin directionbe radial, the ionization front is not spherical when radiationis injected over one smoothing length (left panel), a conse-quence of the fact that the SPH particles are not exactlyuniformly distributed around the source. This can be reme-died by either injecting radiation into gas particles up to twosmoothing lengths away from the source (middle panel) orby requiring that the radiation should move radially awayfrom the source, i.e. ˆ n = ˆ r (right panel).For the actual test, we inject radiation in all gas parti-cles within two smoothing lengths from the source but with-out imposing a propagation direction, as in the middle panelof Fig. 7. The simulation results at time t = 500 Myr are © 0000 RAS, MNRAS000
Propagation in 2D: propagation of a packet of radia-tion in two dimension at the speed of light, c = 1 . Upper panels :Radiation energy density of the packet in units of its initial value.The left panel shows the initial state, where the packet has length δy = 1 . Panels to the right show the system at different times. TheSPH particle distribution is glass-like, with × particles fill-ing the computational volume of horizontal extent ∆ x = 0 . andvertical extent ∆ y = 2 . . Colors represent the radiation energydensity and red arrows show the radiation fluxes. Lower panel : ra-diation energy density, ξ ( x, y ) , of SPH particles with . (cid:54) y (cid:54) . at time t = 1 . ( red line ), the dashed-black line is the initial shapeof the packet. the implementation to maintain the direction of the radia-tion, with little artificial leakage of radiation perpendicularto the beam. In the right panel, the packet has propagatedupwards over 128 mean particle spacings. Two small ‘tails’of radiation trail the package, where radiation leaked out ofthe beam. Smoothing perpendicular to the beam is quan-tified in more in detail in the lower panel, which is a cutthrough the middle of the beam at time t = 1 . This profilehas approximately Gaussian-shaped edges, as expected from artificial diffusion (§2.5); the diffusion coefficient is propor-tional to h c . The dependence on the smoothing length, h ,means that the beam can propagate further without distor-tion at higher resolution. Due to the finite resolution and, ingeneral, non-uniform underlying SPH particle distribution,it is not possible to completely eliminate radiation leakageperpendicular to the propagation direction without causinginstabilities; higher-order shock-capturing schemes (e.g. Liuet al. 1994) might suppress such artificial leakage more effi-ciently.The next test is that of radiation propagating isotropi-cally away from a source in two dimensions, see Fig. 6. Thefigure confirms that the radiation front preserves rotationalsymmetry as it propagates out at the speed of light. The ra-diation energy density is smooth in the radial shell, with noappreciable noise even though the underlying particle distri-bution is non-uniform. The energy density is small behindthe shell. The absence of significant artificial ‘left over’ ra-diation results from the artificial dissipation switch. The first test is Test 1 in Iliev et al. (2006). This tests theradiative transfer scheme and the thermo-chemistry solveragainst an analytical solution: uniform density, neutral gasis photoionized by a source that emits ionizing photons ata constant rate. We keep the density and temperature ofthe gas constant, i.e. the gas is not allowed to move, heator cool. The ionization front propagates into the gas clouduntil it reaches its Strömgren radius. The analytical solu-tion (assuming grey opacity) is derived and summarised inAppendix C.The numerical parameters are taken to be identicalto those used by Iliev et al. (2006) to allow for a directcomparison: the gas cloud consists of pure hydrogen gaswith density n H = 10 − cm − , the collisional ionization co-efficient β = 3 . × − cm s − and the recombinationcoefficient is α = 2 . × − cm s − (with the on-the-spot approximation); the photoionization cross section is σ γ HI = 8 . × − cm . The source emits ionizing radiationat a constant rate of ˙ N γ = 5 × photons s − . The compu-tational volume has linear extent of 20 kpc; the SPH particledistribution is glass-like with approximately particles.In this problem we also test the RSL approximation, using ˜ c = c/ .We first compare several implementations of the injec-tion of radiation energy by the source (§2.8) and of the op-tically thin closure relation (§2.4.1); results are shown inFig. 7. Without requiring that the optically thin directionbe radial, the ionization front is not spherical when radiationis injected over one smoothing length (left panel), a conse-quence of the fact that the SPH particles are not exactlyuniformly distributed around the source. This can be reme-died by either injecting radiation into gas particles up to twosmoothing lengths away from the source (middle panel) orby requiring that the radiation should move radially awayfrom the source, i.e. ˆ n = ˆ r (right panel).For the actual test, we inject radiation in all gas parti-cles within two smoothing lengths from the source but with-out imposing a propagation direction, as in the middle panelof Fig. 7. The simulation results at time t = 500 Myr are © 0000 RAS, MNRAS000 , 000–000 T. K. Chan et al.
Figure 6.
Propagation in 2D: propagation of a shell of light away from a source at the speed of light, c = 1 . Upper panel : radiationenergy density of the shell of light in units of its initial value. The left panel shows the initial state, where radiation is filled uniformlywithin a circle with radius r = 0 . and points out radially. Panels to the right show the system at times t = 0 . and t = 0 . . The SPHparticle distribution is glass-like, with × particles filling the computational volume of horizontal extent ∆ x = 2 and verticalextent ∆ y = 2 . Colours represent the radiation energy density and red arrows show the radiation fluxes; the colour scale is not the samein each panel to bring out the smoothness of the radiation density as the shell moves out. Lower panel : radiation flux ( left panel ) andradiation density ( right panel ) as a function of radius at time t = 0 . . The red points represent values of all individual SPH particles, blue points show binned values. The blue vertical lines at r = 0 . show the location of the radiation front at time t = 0 . . compared to the analytical solution derived in Appendix Cin Fig. 8. By this time, the system is in a steady-state whereionizations balance recombinations and the ionization frontis at the location of the Strömgren radius. We can comparethe run of the neutral fraction in the simulation to the exactanalytical solution.The mean value of the neutral fraction as a functionof radius, x ( r ) = n HI ( r ) /n H , follows the analytical resultclosely with relatively small scatter. There are small sys-tematic deviations from the analytical solution near r = 0 ,where radiation is injected, and at r (cid:62) , where the analyt-ical value of x drops faster than the simulated value. Thelatter is due to radiation ‘leaking’ beyond the Strömgren ra-dius in the simulation due to the artificial dissipation. Theoverall performance of the scheme is relative good: (1) the neutral fraction is approximately spherically symmetric; (2)the scheme is photon-conserving and therefore the locationof the Strömgren radius agrees well with the analytical so-lution; (3) the scheme is accurate both in the optically thinregion near the source as well as in the optically thick regionoutside the Strömgren radius, as well as in the intermediateregion. We note that some cone-based (e.g. Pawlik & Schaye2008) or short-characteristic (e.g. Finlator et al. 2009) RTschemes display artificial ‘ray-like’ features in the neutralfraction if the angular resolution is insufficient; no such fea-tures appear in the present scheme.The effect of using the RSL approximation on the timeevolution of the I-front is illustrated in Fig. 9. We use asimilar setup as above, but to capture the early phase of theexpansion of the ionization front (1) we inject radiation over © 0000 RAS, MNRAS , 000–000 moothed Particle Radiation Hydrodynamics—SPH-M1RT Figure 7.
Isothermal Strömgren sphere from Iliev et al. (2006) Test 1: a source of radiation, located at the centre of the panels,photoionizes hydrogen gas, kept at constant density and constant temperature. Panels show the neutral fraction, n HI /n H , at time t = 500 Myr in a slice through the centre of the three-dimensional volume. Left panel : radiation is injected in all gas particles within onesmoothing length from the source, see §2.8 for the method of injection.
Central panel: radiation is injected in all gas particles within twosmoothing lengths from the source.
Right panel: as in left panel, but the direction in which radiation travels is set to ˆ n = ˆ r , see §2.4.1 fordetails. The ionization front becomes more spherical if the injection region is enlarged, or if the radiation is forced to propagate radially. Figure 8.
Isothermal Strömgren sphere from Iliev et al. (2006)Test 1: a source of radiation, located at radius r = 0 , photoion-izes hydrogen gas, kept at constant density and constant tem-perature. The system is shown a time t = 500 Myrs after thesource is switched on. Red points are the neutral hydrogen frac-tion, x HI = n HI /n H , of individual gas particles, with the thickblue points showing binned values, with horizontal error bars in-dicating the bin width, and vertical error bars the standard devi-ation; the black line is the analytical solution derived in AppendixC. Small blue points , and thick black points with black error barsshow the corresponding ionized fraction, − x HI , with the yellowline the analytical solution. The vertical blue line is the approxi-mate location of the Strömgren radius from Eq. (C5). The meaninterparticle separation of the SPH particles is 0.625 kpc. Herewe inject radiation over two smoothing lengths. only one smoothing length but impose that the radiationpropagates radially outwards, i.e. ˆ n = ˆ r ; (2) we use higherresolution, particles.The traditional analytical solution for the time-dependent location of the ionization front, r I ( t ) , of Eq. (C5) Figure 9.
Test 1 (Iliev et al. 2006): pure hydrogren isothermalHII region expansion from time t = 0 to
500 Myr . Top : the evo-lution of the radius of the ionization front, defined by where 50%hydrogen in the spherical shell is ionized, divided by the Ström-gren radius (Eq. C6).
Bottom : the evolution of the velocity ofthe ionization front, divided by the Strömgren radius over the re-combination time t rec = ( n H α B ) − . The c/ ˜ c lines represent oursimulation results with different reduced speeds of light and thefixed optically thin direction (§2.4.1). The black lines show theanalytic solutions (Eq. C5 and its derivative; note that the ana-lytic solution is only approximately correct). Finally, the ‘C2-ray’lines represent the C2-ray (Mellema et al. 2006) result in Ilievet al. (2006). assumes that the front is infinitely thin and that the down-stream gas is fully ionized. In reality, the down-stream gasis not completely ionized and the analytical solution of © 0000 RAS, MNRAS000
Bottom : the evolution of the velocity ofthe ionization front, divided by the Strömgren radius over the re-combination time t rec = ( n H α B ) − . The c/ ˜ c lines represent oursimulation results with different reduced speeds of light and thefixed optically thin direction (§2.4.1). The black lines show theanalytic solutions (Eq. C5 and its derivative; note that the ana-lytic solution is only approximately correct). Finally, the ‘C2-ray’lines represent the C2-ray (Mellema et al. 2006) result in Ilievet al. (2006). assumes that the front is infinitely thin and that the down-stream gas is fully ionized. In reality, the down-stream gasis not completely ionized and the analytical solution of © 0000 RAS, MNRAS000 , 000–000 T. K. Chan et al.
Figure 10.
Iliev et al. (2006) Test2: HII region expansion in a uniform gas with varying temperature and static gas particles at time t = 100 Myr . Top left : a slice of the hydrogen neutral fraction through the center;
Top right : a slice of the gas temperature through thecenter;
Bottom left : hydrogen neutral fraction as a function of radius;
Bottom right : temperature as a function of radius. The red pointsrepresent the values of individual particles. The vertical error bars show the mean and standard deviation. The horizontal error barsshow the smoothing length. The solid black line is the ‘TT1D thin’ result (Pawlik & Schaye 2011) and the shaded region is the upperand lower bounds of the Iliev et al. (2006) results. Here we inject radiation over three smoothing lengths.
Eq. (C5) is only approximately correct . Because of this,we compare our simulation results to another simulationcode, namely C2-ray (Mellema et al. 2006) (see also Ilievet al. 2006) and, follow them by defining the position of theI-front as the radius at which x HI = 0 . .Fig. 9 demonstrates that our results converge for ˜ c → c ,and even when ˜ c = c/ are close to those obtained with C2-ray . Using ˜ c = c/ , we notice deviations of a few tensof percent at times less than a recombination time, and muchsmaller than 10 per cent later on. This matches our expec- The mean free path of ionizing photons is λ = n H σ γ ≈ .
04 kpc and not resolved in the simulation setup. tation discussed in §2.10 (see also Rosdahl et al. 2013). Thescheme works well at low resolution even when using non-uniform particle distributions. This is important because intypical applications ( e.g. reionization simulations or simula-tions of the interstellar medium) the gas distribution aroundthe ionizing sources is often at best marginally resolved.
We repeat the previous test, but now allow photoheating ofthe ionized gas, testing the interaction between the radiativetransfer and the photochemistry solver. The parameters ofthe test are identical to the previous case; the photoheating © 0000 RAS, MNRAS , 000–000 moothed Particle Radiation Hydrodynamics—SPH-M1RT Figure 11.
Iliev et al. (2006) Test 3: I-front trapping in a dense clump.
Left : a slice of mass-weighted hydrogen neutral fraction throughthe midplane of the simulation volume at t = 15 Myr ; Right : same but for temperature.
Upper panels : for the default optically thindirection, ˆ n = ˆ f ; Lower panels : for an imposed fixed optically thin direction in the x direction (see §2.4.1), as shown by the arrows. and cooling processes are detailed in Appendix A (the testmakes the on-the-spot approximation). We use the opticallythin value for the photoheating energy rate per ionization, (cid:15) γ (see Appendix A). The underlying particle distributionis glass-like with 32 particles on a side; the injection radiusis three smoothing lengths (see §2.8); and ˜ c = c/ . Fig. 10summarises our simulation results.We compare our results to results for the same setuppublished by Iliev et al. (2006), since this test has no knownanalytical solution. Unfortunately the comparison is notstraightforward because some codes in that paper use differ-ent values for the thermo-chemistry coefficients (see Fig. 2in Iliev et al. 2006), and some codes use multi-frequency RTto account for spectral hardening. Spectral hardening leadsto pre-heating of gas upstream from the ionization front. Tomake the comparison appropriate, we also compare to TT1D (TestTraphic1D), which is a 1D radiative transfer code de-veloped by Pawlik & Schaye (2008, 2011) for testing the
TRAPHIC code. We will compare to their ‘
TT1D thin ’ result,which uses the same assumptions as ours, i.e. one frequencybin (grey approximation) and photoheating in the opticallythin limit. Fig. 10 demonstrates that our scheme producesa roughly spherical morphology at this resolution with themean value at a given radius matching those obtained by
TT1Dthin . In addition, our result falls within the grey-banddefined by the range of simulation results published by Ilievet al. (2006).Next, we illustrate our scheme’s ability to trap an ion-ization front and cast a shadow, by repeating ‘Test 3’ inIliev et al. 2006. The setup is as follows: a cubic volumeof linear extent . is filled with hydrogen gas of den-sity n out = 2 × − cm − and temperature T out = 8000 K . © 0000 RAS, MNRAS000
TT1Dthin . In addition, our result falls within the grey-banddefined by the range of simulation results published by Ilievet al. (2006).Next, we illustrate our scheme’s ability to trap an ion-ization front and cast a shadow, by repeating ‘Test 3’ inIliev et al. 2006. The setup is as follows: a cubic volumeof linear extent . is filled with hydrogen gas of den-sity n out = 2 × − cm − and temperature T out = 8000 K . © 0000 RAS, MNRAS000 , 000–000 T. K. Chan et al.
Figure 12.
Iliev et al. (2006) Test 3: I-front trapping in dense clump at t = 15 Myr . Left : the neutral fraction along the axis of symmetrythrough the centre of the clump, where x = 0 ; Right : the corresponding plot for temperature. The red points represent the values ofindividual particles (particles with distance < 0.1 kpc from the axis). The vertical errorbars show the mean and standard deviation andthe horizontal errorbars show the smoothing length. The shaded region encompasses the upper and lower bounds of the Iliev et al. (2006)results. Vertical dashed lines show the initial clump front position, whereas horizontal dashed lines show the initial neutral fraction andtemperature. We only show the result with the fixed optically thin direction (§2.4.1).
At spherical cloud with radius . , hydrogen density n clump = 0 .
04 cm − and temperature T clump = 40 K isplaced at the centre of the volume. The system is irradiatedfrom the top with a black-body spectrum with temperature T BB = 10 K at a constant photon flux of s − cm − . TheSPH particle distribution is glass-like with approximately particles (which is lower than the tests published in Ilievet al. (2006) but enough to demonstrate the performance ofour scheme). We model the overdensity of the gas in theclump by increasing the hydrogen fraction of the SPH par-ticles in the clump. We use the RSL approximation, setting ˜ c = c/ .We show the simulation results in Fig. 11, with (lowerpanels) and without (upper panels) imposing the opticallythin propagation direction for the radiation, ˆ n = ˆ x , ( §2.4.1).As expected, the gas in front of the high-density region ishighly ionized; inside the high-density region but upstreamfrom the ionization front it is ionised to a level n HI /n H ∼ − ; the ionization front is trapped inside the high-densityregion at x ∼ . , and finally the gas is mostly neutralbehind the ionization front in the shadow behind the high-density region.We follow the reasoning of §2.7.3 to estimate the tem-perature immediately after the ionization front has passed, T = 2 − x − x T + 23 x − x − x (cid:15) γ k B , (60)where x and T are the initial neutral fraction and temper-ature, x and T are the neutral fraction and temperaturewhen the ionization front has passed, and (cid:15) γ ≈ .
33 eV isthe photoheating per ionization. In front of the high-densityregion, we take x = 1 , T = 8000 K and x ≈ to find T = 10 . . K ; upstream from the ionization front inside thehigh-density region, we take x = 1 , T = 40K , x = 10 − tofind T = 10 . K . Once the front has passed, the gas will continue to be photoheated while also cooling radiatively.In front of the high-density region, the recombination time, τ r ≡ / ( α B n H ) ≈
600 Myr , is long compared to the sim-ulation time, and the gas is still heating slowly to its newequilibrium temperature of . K . Inside the high-densityregion, τ r ≈ is short compared to the simulation timeand the gas should be in thermal equilibrium. Scaling theflux so that the neutral fraction is − . , close to whatis found in the simulation at the front of the high-densityregion, yields an equilibrium temperature of ≈ . K ,consistent with what we find in the simulation (Fig. 12).When the propagation direction is imposed, the over-dense gas traps the I-front at time t = 15 Myr , and therun of the neutral fraction and temperature from the topof the volume in the x -direction follows our analytical es-timates. It also falls within the grey-region, defined by thethe locus of the simulation results for the different codesfor Test 3 in Iliev et al. (2006), as shown in Fig. 12. Theshadow is relatively sharp, with a small level of ionizationat its boundary due to the numerical/artificial diffusion. Asexpected, is also cooler in the shadow with the temperaturethere agreeing with the results of some of the codes in Ilievet al. (2006). As in the previous test case, our results arenot directly comparable to some of the codes in Iliev et al.(2006), in particular codes that perform multi-frequency RTcapture the pre-heating ahead of the ionization front due tothe smaller optical depth of higher-energy photons.When the direction of propagation is not imposed (toppanels in Fig. 11), the ionization front is still trapped in thehigh-density region and there is still a significant fractionof self-shielded gas. However, the numerical and/or artifi-cial diffusion wipes-out the shadow. Some high-density gasdown-stream from the ionization front also gets ionized as aconsequence of artificial diffusion. © 0000 RAS, MNRAS , 000–000 moothed Particle Radiation Hydrodynamics—SPH-M1RT Figure 13.
Slice plots of neutral fraction ( upper left ), temperature ( upper right ), gas density ( lower left ), and Mach number ( lowerright ) of the HII region expansion in an initially uniform medium (Iliev et al. 2009 Test 5) at t = 200 Myr through the midplane of thesimulation box. Hydrodynamics is turned on and we also fix the optically thin direction (§2.4.1) in this test. If we choose to inject overtwo/three smoothing lengths, the profiles look similar but the front will propagate slightly faster than the reference solution. In this section we present tests with heating, cooling, andhydrodynamics with radiation. The first test is a repeat ofthe HII region, but now allowing the gas to expand as it isheated; this is Test 5 in Iliev et al. (2009) . The problemsetup and conditions are exactly identical to the variabletemperature Strömgren sphere test in §3.3, but we simulatealso the hydrodynamics response of the gas. We use the de-fault SPHENIX
SPH formulation for hydrodynamics (§2.4)with particles in the computational volume. We ignore the Doppler shift and radiation pressure terms (see§2.2), since they are not important in this test problem.
If we do not impose that radiation propagates in the op-tically thin direction and we limit the injection region to onesmoothing length, we find that at the relatively low resolu-tion of particles in the computational volume, the ionizationfront shows strong deviations from sphericity due to parti-cle noise (see also Fig. 7). If instead we inject radiation overa few smoothing lengths, then we find that the ionizationregion will be almost spherical as desired, yet the locationof the front exceeds slightly that found in the reference pro-file computed with
ZEUS-MP , since photons move out toofast initially. To remedy these limitations, we require thatphotons propagate radially, and inject radiation in particlesless than one smoothing length away from the source. Withthis injection scheme, Fig. 13 shows that the gas profiles are © 0000 RAS, MNRAS000
ZEUS-MP , since photons move out toofast initially. To remedy these limitations, we require thatphotons propagate radially, and inject radiation in particlesless than one smoothing length away from the source. Withthis injection scheme, Fig. 13 shows that the gas profiles are © 0000 RAS, MNRAS000 , 000–000 T. K. Chan et al.
Figure 14.
Gas density ( upper left ), gas velocity ( upper right ), neutral fraction ( lower left ), and pressure ( lower right ) as a function ofradius of the HII region expansion in an initially uniform medium (Iliev et al. 2009 Test 5) at t = 200 Myr . The red points representthe values of individual particles. The vertical errorbars show the mean and standard deviation. The horizontal errorbars show thesmoothing length. The lines are the results from the ZEUS-MP code (Whalen & Norman 2006) in Iliev et al. (2009). Black solid representsmonochromatic light, i.e. a single frequency bin, similar to our implementation, whereas green dashed lines represent multifrequencytransfer. Hydrodynamics is turned on and we also fix the optically thin direction (§2.4.1) in this test. spherical and there are no numerical instabilities or irregu-larities.In Fig. 14 we compare our simulation results (red pointsshow values for individual SPH particles, blue points showthe mean and scatter of values in radial bins) to those of ZEUS-MP single-frequency bin result (black solid lines), aspublished by Iliev et al. (2009). For the gas neutral fractionand pressure, the single-frequency bin result is not avail-able, so we compare to the multi-frequency
ZEUS-MP resultinstead.As the gas is ionized and heated, the surrounding gasis swept-up in a dense shell. The location of the shell agreeswell with that found by
ZEUS-MP , as do gas neutral fraction,density and velocity profiles. In our simulation, the shell issomewhat wider than that found by
ZEUS-MP , partly be- cause of our much lower resolution and possibly due to theapplied artificial viscosity in the hydrodynamics solver.Our pressure profile agrees with
ZEUS-MP-multi until r = 8 kpc , where ZEUS-MP-multi shows a slower falloff inpressure. This is because
ZEUS-MP-multi includes also spec-tral hardening , where high frequency photons penetrate fur-ther into the neutral medium. Our current single-bin methodcannot model this effect.
In the previous sections, we have demonstrated that our RTimplementation accurately propagates radiation in the opti-cally thin limit, preserving the direction of propagation and © 0000 RAS, MNRAS , 000–000 moothed Particle Radiation Hydrodynamics—SPH-M1RT advancing the radiation front at the correct speed (§2 and§3.1). In three dimensions, it accurately reproduces the ini-tial expansion of an ionization front around a source, andits asymptotic slow down to the Strömgren radius. The im-plementation also handles ionization front trapping (§3), re-producing results accurately also when the hydrodynamicsof the gas is accounted for, even at moderate numerical res-olution (§3.4).Importantly, the method has favourable computationalscaling, which is proportional only to the number of gasparticles, and independent of the number of sources. By im-plementing a ‘reduced speed of light’ method, the timestepassociated with radiation propagation can be dramaticallyincreased (see §2.10). The thermochemistry uses subcyclingin order to decouple the short ionization timescales from themuch longer radiation propagation timescale (see §2.7).The RT implementation inherits the full spatial andtemporal adaptivity of the underlying gas dynamics scheme(e.g. §3.4). The method as described can be combined withany SPH code, without any need for extra structures (e.g.grids, rays, photon packets, or angular discretization).However, the method also has limitations, e.g. thoseassociated with the closure relation and numerical noise. Inthis section, we discuss these limitations in more detail. As discussed in §2, the two-moment method results fromtruncating an infinite hierarchy of moment equations bypostulating a closure relation for the Eddington tensor. Thechoice of closure relation affects the accuracy of the method.The ‘M1 closure relation’ is not exact in the regime in-termediate between optically thin and optically thick, evenin the case of a single source. As shown by Levermore (1984),this is because in such a situation the closure relation cannotbe uniquely determined by the first two moments. Secondly,this closure relation cannot handle situations where parti-cles receive radiation from two or more directions, even inthe optically thin regime (see §2.3 and §4.1.1). In such cases,beams of radiation ‘collide’ with each other rather than sim-ply pass through one another as they should do. The reasonfor this is twofold. Firstly, the M1 closure relation of Eq. (17)erroneously implies that the radiation is optically thick whentwo beams collide (see also Fig. 1 in Rosdahl et al. 2013).Secondly, the form of the Eddington tensor in Eq. (13) im-plicitly assumes that radiation is moving in a single direction(single stream, plus an additional isotropic component).We illustrate the ‘collision of radiation’ in Fig. 15. Thesetup is as follows: two sources emit a burst of radiationisotropically into an optically thin medium. This resultsin two spherical shells of radiation, propagating away fromeach source. When these shells overlap, they should simplypass through each other. When this situation is simulatedwith the M1 closure relation, the shells ‘collide’ instead,erroneously producing two spikes of radiation. Such colli-sions could significantly distort the radiation morphology e.g. when sources are associated with multiple star clustersor multiple galaxies.There are several ways to improve the method to re-duce the impact of such collisions. Firstly, we introduced the modified M1 closure relation in Eq. (19). This new closurerelation does not lead to radiation collision in one dimen-sion, as shown in Fig. 16. The reason is that this modifiedclosure relation correctly identifies that the radiation is op-tically thin even where the two beams collide (unlike theoriginal M1 close relation).Unfortunately, the modified closure scheme does not re-solve the inability of the closure relation to identify correctly,and propagate several overlapping optically thin beams, andhence it does not resolve the problem illustrated in Fig. 15.A possible way forward would be to resort to higher or-der methods, e.g.
Vikas et al. (2013) and Levermore (1996).However, a stable and efficient high-order method has notbeen discussed in the astrophysics literature, as far as weare aware. Another avenue might be to calculate the closurerelation itself more accurately, for example using short char-acteristic (Finlator et al. 2009; Jiang et al. 2012) or usinga Monte Carlo scheme (Foucart 2018). Unfortunately, boththese schemes are computationally more expensive in thecase of multiple bright sources.
Another limitation of the moment method, which is not re-stricted to the two-moment or M1 methods, is numericalnoise. Such noise can destroy the coherence of the radiationeven in the optically thin regime (see e.g.
Fig. 1), or causeradiation to propagate into a shadow (see e.g.
Fig. 11).There are two major sources of noise. Firstly, the dis-cretization of the density field into a disordered set of SPHparticles, which is especially problematic when the symme-try of the particle distribution differs from that of the ra-diation field and/or when the numerical resolution is low.In addition, SPH suffers from ‘zeroth-order errors’ that arealso particularly severe when the particle distribution is ir-regular (see §1, Dehnen & Aly 2012); our SPH formulationis designed to minimize these (see §2.4). Higher-order SPHschemes could reduce such SPH noise further ( e.g.
Vila 1999;Gaburov & Nitadori 2011; Rosswog 2015; Hopkins 2015;Frontiere et al. 2017).The second major source of noise is the discontinuity-capturing artificial dissipation, which can introduce unphysi-cal diffusion and damping. We have tried to limit the severityof such artefacts by introducing anisotropic viscosity as wellas a higher-order artificial diffusion scheme. Even so, our ar-tificial diffusion is not completely anisotropic, so that a smallamount of radiation still diffuses artificially perpendicular tothe propagation direction of a beam of radiation (see Figs. 5and 11). We suggest that higher order artificial anisotropicdiffusion, similar to that used in the HLL Riemann solver(Harten et al. 1983), might reduce these artefacts.Our method can further reduce unphysical diffusion incases where it is possible to impose the direction of propa-gation of the radiation.
Another limitation of the M1 method is that the computa-tional cost may be high in some physically interesting situ-ations, e.g. when capturing the final stages of reionizationwhen the reduced speed of light needs to be close to c to © 0000 RAS, MNRAS000
Another limitation of the M1 method is that the computa-tional cost may be high in some physically interesting situ-ations, e.g. when capturing the final stages of reionizationwhen the reduced speed of light needs to be close to c to © 0000 RAS, MNRAS000 , 000–000 T. K. Chan et al.
Figure 15.
Radiation collision in three dimensions: spherical radiation fronts emanate from two sources into an optically thin medium.The setup is the same as Fig. 6, but with two spherical fronts instead of one. The panels from left to right show the initial conditionsand the state at later times of t = 0 . and . when the spherical fronts overlap. The simulation uses the standard M1 closure relation(Eq. 13). Colors represent the radiation energy density, red arrows show the radiation fluxes. Whereas the shells should pass througheach other, this closure relation results in a collision of the two fronts, resulting in the formation of horizontal beams of light. Figure 16.
A short beam of radiation propagating to the top,colliding with a short beam of radiation propagating to the bot-tom, simulated with the default two-moment method (see §2.6)in an optically thin medium (with the same condition as in Fig.1,except there is an extra radiation beam pointing downward). Col-ors represent the radiation energy density, whereas the red arrowsshow the radiation fluxes. This shows our modified M1 closure(Eq.19) can handle the head-on beam collision problems. capture the speed of ionization fronts correctly (e.g.Baueret al. 2015). Possible improvements include using a ‘variablespeed of light approximation’ (see §2.10, Katz et al. 2017) orimplementing the radiation transfer on graphics processingunits (Ocvirk et al. 2016). Subcycling the radiative transfermodule can further improve the performance of the overallcode (Rosdahl et al. 2013; Kannan et al. 2019).However, we emphasise that the computational cost of our method scales with the number of gas particles, N gas ,more favourably than ray-tracing methods (which addition-ally scale with the number of sources) and OTVET (whichscales with the N log N scaling of its Poisson solver). The two-moment M1 method is a popular scheme in thefield of galaxy formation simulations and is implemented in e.g.
ATON (Aubert & Teyssier 2008),
RAMSE-RT (Rosdahlet al. 2013), and
AREPO-RT (Kannan et al. 2019). The firsttwo codes are fixed grid AMR schemes, whereas
AREPO-RT is a moving mesh scheme (Springel 2010).Our implementation in the swift code takes advan-tages of the Lagrangian nature of SPH, which is impor-tant particularly when gas flows at high speeds. This is agreat advantage compared to AMR codes, especially at rel-atively low resolution (Robertson et al. 2010) and high red-shift (Pontzen et al. 2020). While fixed grid codes can gainadaptivity through the adaptive refinement, the refinementand derefinement are not trivial and can be noisy.
AREPO-RT and our code can follow fluid motions andhighly adaptive, so both of them are advantageous in galaxyformation simulations, where the large dynamic range andhigh speeds of the gas are both numerically challenging.We think that the main advantages of our scheme areits efficiency and its parallelizability. The moving mesh coderequires mesh reconstruction, which can be computationallyexpensive; such reconstruction is not necessary in SPH. SPHcodes can also be parallelized efficiently through task-basedparallelism, e.g.
ChaNGa (Jetley et al. 2008) and
SWIFT (Schaller et al. 2016), in which our method has been imple-mented. © 0000 RAS, MNRAS , 000–000 moothed Particle Radiation Hydrodynamics—SPH-M1RT We have developed a numerical radiation hydrodynamicsscheme based on the two-moment method and using and im-proved closure relation. The two-moment method is basedon the first two moments of the radiation moment hier-archy, and the hierarchy is closed with an Eddington clo-sure relation. The M1 Eddington tensor closure is the sim-plest anisotropic closure that maximizes the entropy (see§2), our modified closure relation improves the stability ofthe method, in particular in the optically thin regime. Thenumerical scheme is implemented in the smoothed particlehydrodynamics (SPH) code
SWIFT (Schaller et al. 2016).The interaction of radiation with a pure hydrogen gas is im-plemented in non-equilibrium, using subcycling to improvecomputational speed.Key aspects of the method and its implementation in-clude: • the first stable and accurate SPH implementation of thetwo moment method (Eqs. 24 and 25); • improvements to the M1 closure relation (Eq. 19), mak-ing the method (1) less prone to noise in the optically thinregime (enabling the implementation of the two momentmethod in SPH) and (2) correct in the optically thin limit • anisotropic artificial viscosity and high-order artificialdiffusion schemes (§2.5) to capture discontinuities in theradiation. These are essential to propagate radiation accu-rately in the optically thin regime. • an efficient non-equilibrium thermochemistry solverwith subcycling (§2.7) • implementation in SWIFT , a task-based parallel SPHgalaxy simulation code (Schaller et al. 2016).The accuracy and stability of the scheme is demon-strated in §3. The scheme: • preserves the directions and speed of propagation in theoptically thin regime; • is able to cast shadows if the optically thin direction isgiven; • can simulate radiation hydrodynamics in dynamicalmultiscale problems; • has accuracy comparable to other radiative transfercodes in the cosmological code comparison papers (Ilievet al. 2006, 2009).Note that the method yields robust results for sphericallysymmetric problems, e.g. Strömgren sphere, even withoutimposing that radiation propagates radially, provided thatthe injection region is sufficiently well sampled; see §2.8 &3.2. The main advantage of our scheme (and of momentmethods in general) as compared to other radiative transferschemes, is that it is computationally efficient in large-scalesimulations with numerous sources, since this cost is inde-pendent of the number of sources (which has been demon-strated in many other studies; see §1). Our scheme is alsohighly adaptive in both space and time, since the radiationfield is directly sampled by individual SPH particles.Our scheme can be also implemented in other SPHcodes without requiring substantial structural changes. Thisis because the two-moment equations resemble the set of hy-drodynamic equations themselves, e.g. the radiation energy density equation resembles the gas density equation, and theradiation flux equation resembles the momentum equationof hydrodynamics. The artificial dissipation terms for thepropagation of radiation are also very similar to the artifi-cial viscosity and diffusion terms in SPH.This method was developed to enable tracking of thepropagation of ionizing radiation through the intergalac-tic medium or interstellar medium in galaxy simulations. Itshould be sufficiently accurate to correctly model the prop-agation of ionization fronts and cast shadows behind opti-cally thick absorbers, while at the same time being suffi-ciently efficient to be able to handle thousands of sourceswithout overly slowing down the calculations. These designrequirements may make the implementation suited for tack-ling other problems in computational astrophysics. For ex-ample, with a computation time independent of the num-ber of sources, the code would be suited to study the ef-fect of diffuse radiation or account for multiple scatteringevents. Such problems arise, for example, when studying thetransport of IR radiation in the interstellar medium, wheremultiply-scattered IR photons are thought to be importantin transferring momentum to gas (e.g. Skinner & Ostriker2013; Rosdahl & Teyssier 2015; Kannan et al. 2019).We are aiming to release this code to the community inthe near future, after we have completed the following twoimprovements. First, we are implementing multifrequencyradiative transfer so that we can include Helium in the cal-culations. Multi-frequency radiation also enables the mod-elling of spectral hardening. Secondly, we are improving theefficiency of the implementation by decoupling the update ofthe hydrodynamics variables and the radiation variables, i.e.sub-cycling the radiation transport module. Taking advan-tage of the very different time-scales for the motion of gasand of radiation should lead to a large speed-up in comput-ing time. We will report on these improvements and demon-strate the efficiency of our code elsewhere. ACKNOWLEDGEMENTS
We would like to thank J. Borrow, M. Schaller, M. Ivkovic,L. Hausammann, and C. Correa for their assistance withthe implementation in
SWIFT
SWIFT © 0000 RAS, MNRAS000
SWIFT © 0000 RAS, MNRAS000 , 000–000 T. K. Chan et al. (Jones et al. 2001), swiftsimio (Borrow & Borrisov 2020),and NASA’s Astrophysics Data System.
DATA AVAILABILITY
The data underlying this article will be shared on reasonablerequest to the corresponding author (TKC).
REFERENCES
Abel T., Haehnelt M. G., 1999, ApJL, 520, L13Abel T., Wandelt B. D., 2002, MNRAS, 330, L53Agertz O., et al., 2007, MNRAS, 380, 963Altay G., Theuns T., 2013, MNRAS, 434, 748Altay G., Croft R. A. C., Pelupessy I., 2008, MNRAS, 386,1931Anninos P., Zhang Y., Abel T., Norman M. L., 1997, NewAstron, 2, 209Aubert D., Teyssier R., 2008, MNRAS, 387, 295Baek S., Di Matteo P., Semelin B., Combes F., Revaz Y.,2009, A& A, 495, 389Baes M., Camps P., 2015, Astronomy and Computing, 12,33Bauer A., Springel V., Vogelsberger M., Genel S., TorreyP., Sijacki D., Nelson D., Hernquist L., 2015, MNRAS,453, 3593Borrow J., Borrisov A., 2020, Journal of Open Source Soft-ware, 5, 2430Borrow J., Bower R. G., Draper P. W., Gonnet P., SchallerM., 2018, Proceedings of the 13th SPHERIC InternationalWorkshop, Galway, Ireland, June 26-28 2018, pp 44–51Borrow J., Schaller M., Bower R. G., Schaye J., 2020, arXive-prints, p. arXiv:2012.03974Bromm V., Yoshida N., Hernquist L., McKee C. F., 2009,Nature, 459, 49Buchler J. R., 1983, Journal of Quantitative Spectroscopyand Radiative Transfer, 30, 395Cantalupo S., Porciani C., 2011, MNRAS, 411, 1678Cantalupo S., Porciani C., Lilly S. J., Miniati F., 2005,ApJ, 628, 61Cen R., 1992, ApJS, 78, 341Chow E., Monaghan J. J., 1997, Journal of ComputationalPhysics, 134, 296Cohen S. D., Hindmarsh A. C., Dubois P. F., 1996, Com-puters in Physics, 10, 138Commerçon B., Teyssier R., Audit E., Hennebelle P.,Chabrier G., 2011, A& A, 529, A35Cullen L., Dehnen W., 2010, MNRAS, 408, 669Davis S. W., Stone J. M., Jiang Y.-F., 2012, ApJS, 199, 9Dehnen W., Aly H., 2012, MNRAS, 425, 1068Dubroca B., Feugeas J., 1999, Academie des Sciences ParisComptes Rendus Serie Sciences Mathematiques, 329, 915Ferland G. J., Peterson B. M., Horne K., Welsh W. F.,Nahar S. N., 1992, ApJ, 387, 95Ferland G. J., et al., 2017, RMxAA, 53, 385Finlator K., Özel F., Davé R., 2009, MNRAS, 393, 1090Foucart F., 2018, MNRAS, 475, 4186Frontiere N., Raskin C. D., Owen J. M., 2017, Journal ofComputational Physics, 332, 160Gaburov E., Nitadori K., 2011, MNRAS, 414, 129 Gingold R. A., Monaghan J. J., 1977, MNRAS, 181, 375Gnedin N. Y., Abel T., 2001, New Astron, 6, 437Gnedin N. Y., Ostriker J. P., 1997, ApJ, 486, 581González M., Audit E., Huynh P., 2007, A& A, 464, 429Graziani L., Maselli A., Ciardi B., 2013, MNRAS, 431, 722Gunn J. E., Peterson B. A., 1965, ApJ, 142, 1633Harten A., Lax P. D., Leer B. v., 1983, SIAM Review, 25,35Hasegawa K., Umemura M., 2010, MNRAS, 407, 2632Hayes J. C., Norman M. L., 2003, ApJS, 147, 197Hopkins P. F., 2015, MNRAS, 450, 53Hopkins P. F., Grudić M. Y., 2019, MNRAS, 483, 4187Hopkins P. F., Kereš D., Oñorbe J., Faucher-Giguère C.-A., Quataert E., Murray N., Bullock J. S., 2014, MNRAS,445, 581Hopkins P. F., Grudić M. Y., Wetzel A., Kereš D., Faucher-Giguère C.-A., Ma X., Murray N., Butcher N., 2020, MN-RAS, 491, 3702Hui L., Gnedin N. Y., 1997, MNRAS, 292, 27Hunter J. D., 2007, Computing in Science and Engineering,9, 90Iliev I. T., et al., 2006, MNRAS, 371, 1057Iliev I. T., et al., 2009, MNRAS, 400, 1283Jetley P., Gioachin F., Mendes C., Kale L. V., Quinn T.,2008, in 2008 IEEE International Symposium on Paralleland Distributed Processing. pp 1–12Jiang Y.-F., Oh S. P., 2018, ApJ, 854, 5Jiang Y.-F., Stone J. M., Davis S. W., 2012, ApJS, 199, 14Jones E., Oliphant T., Peterson P., et al., 2001, SciPy:Open source scientific tools for Python,
Jonsson P., 2006, MNRAS, 372, 2Jubelgas M., Springel V., Dolag K., 2004, MNRAS, 351,423Kannan R., Vogelsberger M., Marinacci F., McKinnon R.,Pakmor R., Springel V., 2019, MNRAS, 485, 117Katz H., Kimm T., Sijacki D., Haehnelt M. G., 2017, MN-RAS, 468, 4831Kegerreis J. A., Eke V. R., Gonnet P., Korycansky D. G.,Massey R. J., Schaller M., Teodoro L. F. A., 2019, MN-RAS, 487, 5029Kessel-Deynet O., Burkert A., 2000, MNRAS, 315, 713Kim J.-h., Wise J. H., Abel T., Jo Y., Primack J. R., Hop-kins P. F., 2019, ApJ, 887, 120Krumholz M. R., Thompson T. A., 2012, ApJ, 760, 155Lanson N., Vila J.-P., 2008, SIAM Journal on NumericalAnalysis, 46, 1912Levermore C. D., 1984, J. Quant. Spec. Radiat. Transf.,31, 149Levermore C. D., 1996, Journal of Statistical Physics, 83,1021Levermore C. D., Pomraning G. C., 1981, ApJ, 248, 321Liu X.-D., Osher S., Chan T., 1994, Journal of Computa-tional Physics, 115, 200Lucy L. B., 1977, AJ, 82, 1013Madau P., Meiksin A., 1994, ApJL, 433, L53Mellema G., Iliev I. T., Alvarez M. A., Shapiro P. R., 2006,New Astron, 11, 374Mihalas D., Mihalas B. W., 1984, Foundations of radiationhydrodynamicsMihalas D., Auer L. H., Mihalas B. R., 1978, ApJ, 220,1001 © 0000 RAS, MNRAS , 000–000 moothed Particle Radiation Hydrodynamics—SPH-M1RT Minerbo G. N., 1978, Journal of Quantitative Spectroscopyand Radiative Transfer, 20, 541Monaghan J. J., 1997, Journal of Computational Physics,136, 298Monaghan J. J., 2002, MNRAS, 335, 843Morris J. P., Monaghan J. J., 1997, Journal of Computa-tional Physics, 136, 41Ocvirk P., et al., 2016, MNRAS, 463, 1462Ocvirk P., Aubert D., Chardin J., Deparis N., Lewis J.,2019, A& A, 626, A77Osterbrock D. E., 1989, Astrophysics of Gaseous NebulaeandActive Galactic Nuclei. Palgrave Macmillan.Owen J. M., 2004, Journal of Computational Physics, 201,601Pawlik A. H., Schaye J., 2008, MNRAS, 389, 651Pawlik A. H., Schaye J., 2011, MNRAS, 412, 1943Pawlik A. H., Rahmati A., Schaye J., Jeon M., Dalla Vec-chia C., 2017, MNRAS, 466, 960Petkova M., Springel V., 2009, MNRAS, 396, 1383Petkova M., Springel V., 2011, MNRAS, 415, 3731Pomraning G. C., 1973, The equations of radiation hydro-dynamicsPontzen A., Rey M. P., Cadiou C., Agertz O., TeyssierR., Read J. I., Orkney M. D. A., 2020, arXiv e-prints, p.arXiv:2009.03313Price D. J., 2008, Journal of Computational Physics, 227,10040Price D. J., 2012, Journal of Computational Physics, 231,759Price D. J., et al., 2018, PASA, 35, e031Raviart P. A., 1985, in Brezzi F., ed., Numerical Methods inFluid Dynamics. Springer Berlin Heidelberg, Berlin, Hei-delberg, pp 243–324Reynolds D. R., Hayes J. C., Paschos P., Norman M. L.,2009, Journal of Computational Physics, 228, 6833Rijkhorst E. J., Plewa T., Dubey A., Mellema G., 2006,A& A, 452, 907Robertson B. E., Kravtsov A. V., Gnedin N. Y., Abel T.,Rudd D. H., 2010, MNRAS, 401, 2463Rosdahl J., Teyssier R., 2015, MNRAS, 449, 4380Rosdahl J., Blaizot J., Aubert D., Stranex T., Teyssier R.,2013, MNRAS, 436, 2188Rosswog S., 2015, MNRAS, 448, 3628Rosswog S., 2020a, MNRAS, 498, 4230Rosswog S., 2020b, ApJ, 898, 60Sargent W. L. W., Young P. J., Boksenberg A., Tytler D.,1980, ApJS, 42, 41Schaller M. e. a., 2018, SWIFT: SPH With Inter-dependentFine-grained Tasking, Astrophysics Source Code Library(ascl:1805.020)Schaller M., Dalla Vecchia C., Schaye J., Bower R. G., The-uns T., Crain R. A., Furlong M., McCarthy I. G., 2015,MNRAS, 454, 2277Schaller M., Gonnet P., Chalk A. B. G., Draper P. W.,2016, arXiv e-prints, p. arXiv:1606.02738Schaye J., et al., 2015, MNRAS, 446, 521Shapiro P. R., Giroux M. L., 1987, ApJL, 321, L107Skinner M. A., Ostriker E. C., 2013, ApJS, 206, 21Skinner M. A., Dolence J. C., Burrows A., Radice D., Var-tanyan D., 2019, ApJS, 241, 7Spitzer L., 1978, Physical processes in the interstellarmedium, doi:10.1002/9783527617722. Springel V., 2005, MNRAS, 364, 1105Springel V., 2010, MNRAS, 401, 791Springel V., Hernquist L., 2002, MNRAS, 333, 649Springel V., Pakmor R., Zier O., Reinecke M., 2020, arXive-prints, p. arXiv:2010.03567Strömgren B., 1939, ApJ, 89, 526Susa H., 2006, PASJ, 58, 445Theuns T., Leonard A., Efstathiou G., Pearce F. R.,Thomas P. A., 1998, MNRAS, 301, 478Thomas J., 1998, Numerical Partial Differential Equations:Finite Difference Methods. Vol. 22, Springer Science &Business MediaTricco T. S., Price D. J., 2012, Journal of ComputationalPhysics, 231, 7214Turner N. J., Stone J. M., 2001, ApJS, 135, 95Vandenbroucke B., Wood K., 2018, Astronomy and Com-puting, 23, 40Verhamme A., Schaerer D., Maselli A., 2006, A& A, 460,397Verner D. A., Ferland G. J., Korista K. T., Yakovlev D. G.,1996, ApJ, 465, 487Vikas V., Hauck C. D., Wang Z. J., Fox R. O., 2013, Journalof Computational Physics, 246, 221Vila J. P., 1999, Mathematical Models and Methods in Ap-plied Sciences, 09, 161Wadsley J. W., Keller B. W., Quinn T. R., 2017, MNRAS,471, 2357Whalen D., Norman M. L., 2006, ApJS, 162, 281Whitehouse S. C., Bate M. R., 2004, MNRAS, 353, 1078Wise J. H., Abel T., 2011, MNRAS, 414, 3458Zheng Z., Miralda-Escudé J., 2002, ApJ, 578, 33van der Walt S., Colbert S. C., Varoquaux G., 2011, Com-puting in Science Engineering, 13, 22
APPENDIX A: THERMO-CHEMISTRY RATECOEFFICIENTS
The hydrogen photo-ionization rate, Γ γ HI in units s − , is(e.g. Osterbrock 1989) Γ γ, HI = (cid:90) ∞ ν HI πJ ν π (cid:126) ν σ γ ( ν ) d ν , (A1)where hν HI ≈ . is the hydrogen binding energy, J ν is the angular-averaged specific intensity I (in unit of erg cm − s − Hz − ), and σ γ is the photo-ionization cross-section of hydrogen as a function of frequency, ν , where weused the fit from Verner et al. (1996). Defining the frequency-averaged photo-ionization cross-section (cid:104) σ γ (cid:105) ≡ (cid:20)(cid:90) ∞ ν HI πJ ν π (cid:126) ν σ γ ( ν ) d ν (cid:21) (cid:20)(cid:90) ∞ ν HI πJ ν π (cid:126) ν d ν (cid:21) − , (A2)the photo-ionization rate of Eq. (A1) is Γ γ, HI = (cid:104) σ γ (cid:105) (cid:90) ∞ ν HI πJ ν π (cid:126) ν d ν = (cid:104) σ γ (cid:105) ˜ cn γ , (A3)where the frequency-averaged photon flux is: ˜ cn γ ≡ (cid:90) ∞ ν HI πJ ν π (cid:126) ν d ν . (A4)For reference, the spectrum of a black-body (BB) with tem-perature T = 10 K has (cid:104) σ γ (cid:105) = 1 . × − cm . © 0000 RAS, MNRAS000
The hydrogen photo-ionization rate, Γ γ HI in units s − , is(e.g. Osterbrock 1989) Γ γ, HI = (cid:90) ∞ ν HI πJ ν π (cid:126) ν σ γ ( ν ) d ν , (A1)where hν HI ≈ . is the hydrogen binding energy, J ν is the angular-averaged specific intensity I (in unit of erg cm − s − Hz − ), and σ γ is the photo-ionization cross-section of hydrogen as a function of frequency, ν , where weused the fit from Verner et al. (1996). Defining the frequency-averaged photo-ionization cross-section (cid:104) σ γ (cid:105) ≡ (cid:20)(cid:90) ∞ ν HI πJ ν π (cid:126) ν σ γ ( ν ) d ν (cid:21) (cid:20)(cid:90) ∞ ν HI πJ ν π (cid:126) ν d ν (cid:21) − , (A2)the photo-ionization rate of Eq. (A1) is Γ γ, HI = (cid:104) σ γ (cid:105) (cid:90) ∞ ν HI πJ ν π (cid:126) ν d ν = (cid:104) σ γ (cid:105) ˜ cn γ , (A3)where the frequency-averaged photon flux is: ˜ cn γ ≡ (cid:90) ∞ ν HI πJ ν π (cid:126) ν d ν . (A4)For reference, the spectrum of a black-body (BB) with tem-perature T = 10 K has (cid:104) σ γ (cid:105) = 1 . × − cm . © 0000 RAS, MNRAS000 , 000–000 T. K. Chan et al.
Recombination rate (Hui & Gnedin 1997) by fitting Ferland et al. (1992) α A = . × − cm s − λ . [1 . λ/ . . ] − . α B = . × − cm s − λ . [1 . λ/ . . ] − . Collisional ionization rate Theuns et al. (1998) modified from Cen (1992) β = . × − cm s − T / exp( − . /T )(1 + T / ) − Collisional ionization cooling rates Theuns et al. (1998) modified from Cen (1992) Γ ion , eHI = 2 . × − erg cm s − T / exp( − . /T )(1 + T / ) − Collisional excitation cooling rates Theuns et al. (1998) modified from Cen (1992) Γ line , eHI = 7 . × − erg cm s − exp( − /T )(1 + T / ) − Recombination cooling rates taken from Hui & Gnedin (1997) (fitted from Ferland et al. (1992)) Γ A,e
HII = 1 . × − erg cm s − K − T λ . [1 . λ/ . . ] − . Γ B,e
HII = 3 . × − erg cm s − K − T λ . [1 . λ/ . . ] − . Bremsstrahlung cooling rate Theuns et al. (1998) modified from Cen (1992) and Spitzer (1978) Γ ff , eHII = 1 . × − erg cm s − T / { . .
34 exp( − [5 . − log( T )] / } Table A1.
Coefficients for heating and cooling of hydrogen. T n = T/ (10 n K) and temperature T is in K. λ = 315614K /T . If theon-the-spot approximation is applied, we used the case B recombination cooling (in Eq. 41). Otherwise, the case A recombinationcooling is used. Figure A1.
Injected energy, (cid:15) , in eV per photo-ionisation as afunction of the temperature, T , of the irradiating black-body spec-trum. Solid blue line is the optically thin case, dashed red line isthe optically thick case.
The energy injected into the gas per photo-ionization is (cid:15) γ = (cid:20)(cid:90) ∞ ν HI πJ ν (cid:126) ν σ γ ( ν )( (cid:126) ν − (cid:126) ν HI ) d ν (cid:21) × (cid:20)(cid:90) ∞ ν HI πJ ν π (cid:126) ν σ γ ( ν ) d ν (cid:21) − , (A5)in the ‘optically thin’ limit where the probability that a pho-ton of frequency ν is responsible for the ionization is set bythe photo-ionization cross-section.In the ‘optically thick’ limit, we simply assume that ev-ery photon with hν (cid:62) hν HI causes an ionization, and replace σ γ → . This increases the value of (cid:15) γ as higher-energy pho-tons contribute relatively more to the ionizations (e.g. Abel& Haehnelt 1999). This is only an approximation: in theoptically thick-limit, hard photons with hν (cid:29) hν HI tendto partially ionize gas upstream from the ionization front, which is not describe accurately by simply increasing thevalue of (cid:15) γ .In the special case of a BB spectrum, J ν = B ν , with B ν = 4 π (cid:126) ν c π (cid:126) ν/k B T ) − , (A6)the photo-heating per photo-ionization is (cid:15) γ kT = (cid:82) ∞ ζ T ζ (exp( ζ ) − − σ γ HI ( ζ ) dζ (cid:82) ∞ ζ T ζ (exp( ζ ) − − σ γ HI ( ζ ) dζ − ζ T , (A7)where ζ T ≡ π (cid:126) ν th / ( k B T ) ≈ . × K /T is a dimen-sionless fraction. The value of (cid:15) γ as a function of the tem-perature T of the BB is plotted in Fig. A1. For reference, (cid:15) γ ≈ .
33 eV ( . ) in the optically-thin (thick) case,when T = 10 K .Table A1 lists the interpolation formula for the ioniza-tion, recombination, heating, and cooling coefficients for hy-drogen, as used in the thermo-chemistry described in section§2.7. For reference, the cooling rate due to collisional line ex-citation, collisional ionization, thermal Bremsstrahlung andrecombination radiation (in the on-the-spot approximation),in units of erg cm − s − , are respectively ρ dudt (cid:12)(cid:12)(cid:12)(cid:12) line = − x (1 − x )Γ line , eHI n (A8) ρ dudt (cid:12)(cid:12)(cid:12)(cid:12) Ion = − x (1 − x )Γ ion , eHI n (A9) ρ dudt (cid:12)(cid:12)(cid:12)(cid:12) Bremss = − (1 − x ) Γ ff , eHII n (A10) ρ dudt (cid:12)(cid:12)(cid:12)(cid:12) Recomb = − (1 − x ) Γ B ,e HII n , (A11)whereas the photo-heating rate is ρ dudt (cid:12)(cid:12)(cid:12)(cid:12) Heat = x(cid:15) γ Γ γ, HI n H . (A12)These rates are plotted in Fig.A2, together with the ther-mal equilibrium values of the temperature and neutral frac-tion. For a gas temperature of K (cid:54) T (cid:54) K , line cool-ing typically dominates the cooling rate, whereas at lowertemperatures recombination cooling takes over. For refer-ence, the case-A and case-B recombination coefficients are © 0000 RAS, MNRAS , 000–000 moothed Particle Radiation Hydrodynamics—SPH-M1RT Figure A2. Left panel:
Heating and cooling rates from Table A1 as a function of temperature under thermal equilibrium: collisionalline cooling rate, x (1 − x )Γ line , eHI ( T ) ( orange solid ), collisional ionization cooling rate, x (1 − x )Γ ion , eHI ( T ) ( blue dashed ), Bremsstrahlungcooling rate, (1 − x ) Γ ff , eHII ( T ) ( red dotted ), recombination cooling rate (1 − x ) Γ B,e
HII ( T ) ( green long-dashed ) and photo-heating rate, (cid:15) γ x Γ γ, HI /n H ( purple dot-dashed ). Each rate is shown twice: for a hydrogen number density of n H = 1 cm − (lines with symbols),and for n H = 10 − cm − (lines without symbols). The assumed photo-ionisation rate is Γ = 10 − s − and the energy injected perphoto-ionization is (cid:15) γ = 6 .
33 eV , appropriate for a black-body spectrum of temperature K Right panel: corresponding equilibriumtemperature ( T eq , solid black line) and neutral fraction ( x eq , blue dashed line). Figure B1.
The number of sub-cycles per global time-step in thesingle particle photo-ionization heating test, described in §2.7.3,as a function of time. The global time-step ∆ t main = 0 . .This ratio can be up-to a factor of ten initially, reducing to tensof percents after the gas reaches equilibrium. α A = 4 . × − cm s − and α B = 2 . × − cm s − respectively, and the collisional ionization coefficient β =1 . × − cm s − at a gas temperature of T = 10 K . APPENDIX B: THERMO-CHEMISTRYSUB-CYCLING
Here we analyze the computational efficiency of using sub-cycling in the thermo-chemistry solver, as described in §2.7.Note that during sub-cycling, we only need to solve a fewequations (at most quadratic) for each active gas particle,making the computational cost small compared to calculat-ing hydrodynamical or gravitational forces, or performingthe radiative transfer, all of which require loops over neigh-bouring particles and potentially communication betweencompute nodes. Therefore, the number of sub-cycling steps
Figure B2.
The number of sub-cycles per global time-step forthe variable temperature Strömgren sphere test with static gasparticles described in §3.2. This ratio is plotted for every SPHparticle ( red dots ) as a function of its distance to the centre, attime t = 100 Myr ; blue points with error bars show the mean andvariance in radial bins. The speed-up due to sub-cycling is arounda factor of two, mainly near the location of the ionization front. per global time-step is a measure of the efficiency of sub-cycling.We plot this ratio for the single gas particle test (see§2.7.3) in Fig. B1, and for the static 3D Strömgren test(see §3.2) in Fig. B2. The number of sub-cycling stepsper global time-step is typically highest in nearly neutralregimes, where the ratio can be up-to an order of magnitude.This is because the gas in this regime is far from equilibrium,and the sub-cyle time-step is limited by the photo-ionizationtime scale, τ i = 1 / Γ γ, HI , which can be much shorter thanthe global time-step set by the cfl condition. Once the gasis highly ionized, sub-cyle time-step is usually set by the re-combination time-scale, which is typically much longer than τ i . In addition, the gas may be in ionization equilibrium oreven thermal equilibrium, so that the chemistry time-step islong. In the kind of astrophysical application that we have © 0000 RAS, MNRAS000
The number of sub-cycles per global time-step forthe variable temperature Strömgren sphere test with static gasparticles described in §3.2. This ratio is plotted for every SPHparticle ( red dots ) as a function of its distance to the centre, attime t = 100 Myr ; blue points with error bars show the mean andvariance in radial bins. The speed-up due to sub-cycling is arounda factor of two, mainly near the location of the ionization front. per global time-step is a measure of the efficiency of sub-cycling.We plot this ratio for the single gas particle test (see§2.7.3) in Fig. B1, and for the static 3D Strömgren test(see §3.2) in Fig. B2. The number of sub-cycling stepsper global time-step is typically highest in nearly neutralregimes, where the ratio can be up-to an order of magnitude.This is because the gas in this regime is far from equilibrium,and the sub-cyle time-step is limited by the photo-ionizationtime scale, τ i = 1 / Γ γ, HI , which can be much shorter thanthe global time-step set by the cfl condition. Once the gasis highly ionized, sub-cyle time-step is usually set by the re-combination time-scale, which is typically much longer than τ i . In addition, the gas may be in ionization equilibrium oreven thermal equilibrium, so that the chemistry time-step islong. In the kind of astrophysical application that we have © 0000 RAS, MNRAS000 , 000–000 T. K. Chan et al. in mind, for example reionization simulations, sub-cycling isessential since otherwise the short time-step required in gasbeing over run by an ionization front will grind the code to ahalt. Finally, our sub-cycling scheme parallelizes well, sinceit does not require any communication.
APPENDIX C: ANALYTIC SOLUTION OF THESTRÖMGREN SPHERE
In the classical Strömgren sphere problem (Strömgren 1939),a source emitting ionizing photons at a constant rate ˙ N γ is embedded in a spherical cloud, initially filled with com-pletely neutral hydrogen atoms with density n HI = n H . Asthe source switches on, an ionization front expands aroundthe source, and the gas inside the ionization front, radius R I , will be mostly ionized, x ≡ n HI /n H (cid:28) , and outside R I will be mostly neutral.The equations describing the evolution of the neutralfraction and photon density of such an idealized system are ∂n HI ∂t = − n HI cσ γ n γ + n e n HII α A − n e n HI β , (C1) ∂n γ ∂t = − n HI cσ γ n γ + n e n HII ( α A − α B ) + S γ , = ∂n HI ∂t − α B n e n HII + S γ . (C2)As a first approximation to describe such a system, weassume that the gas inside R I is fully ionized, x = 0 , andoutside R I is fully neutral, x = 1 , and that the ionizationfront is infinitely sharp. We further neglect collisional ion-ization, setting β = 0 . In this case, n HI ( r ) = n HI Θ( R I − r ) ,where Θ( x ) is the step function. Integrating each term ofEq. (C2) over the volume centered at the source, we find (cid:90) ∂n HI ∂t d V = (cid:90) ∂∂t [ n HI Θ( R I − r )] d V = n HI (cid:90) δ ( R I − r ) ∂R I ∂t d V = − πR I n HI ∂R I ∂t (cid:90) α B n e n HII d V ≈ α B n π R (cid:90) S γ d V = ˙ N γ . (C3)Combined, these yield the well-known equation, ˙ N γ = 4 πR I n H ∂R I ∂t + α B n π R I , (C4)with solution R I ( t ) = R S [1 − exp( − t/τ r )] / . (C5)Here, the recombination time τ r = ( α B n H ) − and theStrömgren radius R S ≡ (cid:18) N γ πα B n (cid:19) . (C6)For times t (cid:29) τ r , and greater than the ionization time scale τ i , τ i ≡ (4 π/ n H R s ˙ N γ , (C7) the ionization front reaches an equilibrium location whereionizations balance recombinations.To derive the profile of the hydrogen neutral fraction inequilibrium analytically, we work with the time independentequation in the on-the-spot approximation, α A → α B , ∂n γ ∂t = −∇ · f γ − n HI cσ γ n γ = 0 , (C8)In spherical symmetry (and neglecting the scattered radia-tion), f γ = n γ c ˆ r so that ∇ · ( cn γ ˆ r ) + n HI cσ γ n γ = 0 . (C9)Writing this is spherical coordinates yields r ∂∂r (cid:0) r n γ (cid:1) + n HI σ γ n γ = 0 , (C10)with formal solution n γ = ˙ N γ πc r exp (cid:18) − (cid:90) r n HI ( r (cid:48) ) σ γ d r (cid:48) (cid:19) . (C11)This derivation implicitly assumes that the mean free-pathof photons is much less than R S , which should be a goodapproximation for typical HII regions.The steady-state neutral hydrogen profile follows bybalancing ionizations and recombinations, i.e. substitutingthe previous relation into Eq. (C1) and setting ∂n HI /∂t = 0 , xσ γ ˙ N γ πr exp ( − τ ( r )) = (1 − x ) n H α B − x (1 − x ) β n , (C12)where the optical depth is given by τ ( r ) = n H σ γ (cid:90) r x ( r (cid:48) )d r (cid:48) . (C13)The nature of the solution is brought out better by cast-ing these equations in dimensionless form, τ S xq exp( − τ ) = 3(1 − x ) − βα B x (1 − x ) τ = τ S (cid:90) q x ( q (cid:48) ) dq (cid:48) , (C14)where the dimensionless radius q ≡ r/R s , and the ‘Ström-gren optical depth’ τ S ≡ n H σ γ R S . This is an integral equa-tion for the neutral fraction, x ( q ) . We follow Altay & Theuns(2013) to convert this into an easier to integrate differentialequation: take the logarithm of both sides and differentiatewith respect to q , which yield (cid:20) x + 2(1 − x ) − (1 − x ) β/α B (1 − x ) − x (1 − x ) β/α B (cid:21) dxdq = τ S x + 2 q , (C15)with boundary condition for q → x → τ S q . (C16)Provided that collisional ionizations can be neglected, thedifferential equation simplifies to dxdq = x (1 − x )1 + x (cid:18) τ S x + 2 q (cid:19) , (C17)which shows that there is a one-parameter family of solu-tions that are characterised by the value of τ S . The numeri-cal integration of Eq. C15 with its associated boundary con-ditions is plotted as the line labelled ‘Analytic’ in Fig. 8. © 0000 RAS, MNRAS , 000–000 moothed Particle Radiation Hydrodynamics—SPH-M1RT APPENDIX D: MOMENT DERIVATIONS
This short Appendix aims to elucidate the analogy betweentaking moments of the Boltzmann equation to derive thefluid equations, and taking moments of the RT equation toderive the moment equations for radiation. The collisionalBoltzmann equation is ∂∂t f + v · ∂∂ x f + a · ∂∂v f = (cid:18) DfDt (cid:19) coll , (D1)where the right hand side (R.H.S.) is the collision term, andthe distribution function f is a function of position, x , ve-locity v , and time, t . We will suppress this dependency toavoid clutter. Moments of the equation are derived by mul-tiplying Eq. (D1) with some function Q ( v ) and integratingover velocities, ∂∂t (cid:90) Q ( v ) f d v + ∂∂ x · (cid:90) v Q ( v ) f d v + a · (cid:90) (cid:26) ∂∂ v Q ( v ) f − f ∂∂ v Q ( v ) (cid:27) d v = (cid:90) Q ( v ) (cid:18) DfDt (cid:19) coll d v . (D2)The first term in curly brackets is the flux Q f evaluated atthe integration limits of the velocity. Provided we integrateover all velocities, we can assume that f goes to zero suffi-ciently fast that this term vanishes. Writing the velocity as v = V + w , where V is the mean and w is the random com-ponent of v , the density, momentum, pressure and viscousstress tensor, are defined as ρ ≡ m (cid:90) f d v P ≡ (cid:90) w f d v Π ij ≡ P δ ij − (cid:90) w i w j f d v . (D3)The fluid equations then follow by realising that integralsover collision term on the R.H.S. of Eq. (D2) vanishes forfunctions Q that are conserved during collisions. This is thecase for Q = m (the particle’s mass), and Q = m v (the par-ticle’s momentum), which then yield the first two momentequations (the continuity and Euler equations), ∂∂t ρ + ∂∂ x · ρ V = 0 ∂∂t ρ V i + ∂∂ x j (cid:16) ρ V j + P δ ij − Π ij (cid:17) − ρ a i = 0 . (D4)Contrast this derivation with taking moments of the radia-tive transfer equation (e.g Levermore & Pomraning 1981) , c ∂∂t I + n · ∂∂ x I = (cid:18) DIDt (cid:19) SS , (D5) Here we only briefly illustrate the concept, so we suppress theacceleration/Doppler shift term of Eq. (D5) for simplicity. Thismissing term is included in §2.2 which is based on the deriva-tion by Buchler (1983). Gnedin & Ostriker (1997) (see also e.g.Petkova & Springel 2009; Cantalupo & Porciani 2011) did includethe rate of change of frequency in Eq. (D5) but only to includeredshifting of photons while neglecting Doppler shifting (from gasvelocity). where the specific intensity I is a function of position, x , di-rection, θ, φ in spherical coordinates, frequency, ν , and time, t . The Cartesian coordinates of the unit vector in direction n are n = (sin( θ ) cos( φ ) , sin( θ ) sin( φ ) , cos( θ )) . The R.H.S.now represents photon sources and sinks.We proceed as before, by multiplying with some func-tion Q ( θ, φ ) , which can be a scalar or a tensor, and inte-grating the RT equation over solid angle d Ω , but not overfrequency,Taking Q = 1 and Q = n yield the first two momentequations, ∂∂t E + ∂∂x i F i = (cid:90) (cid:18) DIDt (cid:19) SS d Ω1 c ∂∂t F i + ∂∂x j P ij = (cid:90) n i (cid:18) DIDt (cid:19) SS d Ω E = 1 c (cid:90) I d
Ω = 3 P P ij ≡ P δ ij − Π ij F i = (cid:90) n i I d ΩΠ ij = P δ ij − (cid:90) n i n j Id Ω . (D6)Here, E is the energy density, P the radiation pressure, F theflux, and the trace-less tensor Π is the radiation equivalentof the viscous stress tensor. These are the moment equationsof Eq.(5) and Eq. (6) in the fluid frame, v = 0 .In the special case where the sources plus sinks termhave the form of isotropic absorption, ( DI/Dt ) SS = − κ I ,where κ is the isotropic absorption coefficient, the sink termsin Eq. (D6) are − κE and − κ F for the first and second equa-tion, respectively. Provided that the term (1 /c ) ∂ F /∂t canbe neglected, the moment equations then combine to ∂∂t E = ∇ · (cid:20) κ ∇ P (cid:21) − κE , (D7)which is a diffusion equation. In the isotropic case, Π ij = 0 , and a Gaussian package of the form E ( x , t ) = (cid:0) πσ (cid:1) − / exp[ − x / (2 σ )] exp( − κt ) is a solution, spread-ing out as σ ( t ) = σ ( t = 0) + 2 t/ (3 κ ) while dimming, ∝ exp( − κt ) . This paper has been typeset from a TEX/L A TEX file prepared by the author. © 0000 RAS, MNRAS000