Radiation pressure driving of a dusty atmosphere
MMon. Not. R. Astron. Soc. , 000–000 (0000) Printed 14 November 2018 (MN L A TEX style file v2.2)
Radiation pressure driving of a dusty atmosphere
Benny T.-H. Tsang and Miloˇs Milosavljevi´c
Department of Astronomy, University of Texas at Austin, Austin, TX 78712, USA
14 November 2018
ABSTRACT
Radiation pressure can be dynamically important in star-forming environments such as ultra-luminous infrared and submillimeter galaxies. Whether and how radiation drives turbulenceand bulk outflows in star formation sites is still unclear. The uncertainty in part reflects thelimitations of direct numerical schemes that are currently used to simulate radiation trans-fer and radiation-gas coupling. An idealized setup in which radiation is introduced at thebase of a dusty atmosphere in a gravitational field has recently become the standard test forradiation-hydrodynamics methods in the context of star formation. To a series of treatmentsfeaturing the flux-limited-diffusion approximation as well as a short-characteristics tracingand M1 closure for the variable Eddington tensor approximation, we here add another treat-ment that is based on the Implicit Monte Carlo radiation transfer scheme. Consistent withall previous treatments, the atmosphere undergoes Rayleigh-Taylor instability and readjuststo a near-Eddington-limited state. We detect late-time net acceleration in which the turbulentvelocity dispersion matches that reported previously with the short-characteristics-based ra-diation transport closure, the most accurate of the three preceding treatments. Our technicalresult demonstrates the importance of accurate radiation transfer in simulations of radiativefeedback.
Key words: star: formation – ISM: kinematics and dynamics – galaxies: star formation –radiative transfer – hydrodynamics – methods: numerical
The forcing of gas by stellar and dust-reprocessed radiation hasbeen suggested to reduce star formation efficiency and drive super-sonic turbulence and large-scale outflows in galaxies (e.g., Thomp-son et al. 2005, 2015; Murray et al. 2010, 2011; Faucher-Gigu`ereet al. 2013; Kuiper et al. 2015). Generally the net effect of radiationpressure is to counter the gravitational force and modulate the rateof infall and accretion onto star-forming sites. In its most extremepresentation, radiation pressure accelerates gas against gravity sointensely that the gas becomes unbound. For example, Geach et al.(2014) have recently suggested that stellar radiation pressure drivesthe high-velocity, extended molecular outflow seen in a starburstgalaxy at z = 0 . . Theoretical and observational evidence thussuggests that radiation may profoundly influence the formation andevolution of star clusters and galaxies. While direct radiation pres-sure from massive young stars may itself be important in some,especially dust-poor environments (Wise et al. 2012), the trappingof the stellar radiation that has been reprocessed by dust grains intothe infrared (IR) should be the salient process enabling radiationpressure feedback in systems with the highest star formation ratedensities.Usually, the amplitude of radiative driving of the interstellarmedium (ISM) in star-forming galaxies is quantified with the av-erage Eddington ratio defined as the stellar UV (or, alternatively,emerging IR) luminosity divided by the Eddington-limited lumi- nosity computed with respect to the dust opacity. However, in re-ality, the ISM is turbulent and dust column densities vary widelybetween different directions in which radiation can escape. The lo-cal Eddington ratio along a particular, low-column-density direc-tion can exceed unity even when the average ratio is below unity.Thompson & Krumholz (2014) argue that this is sufficient for radi-ation pressure to accelerate gas to galactic escape velocities. An-drews & Thompson (2011) surveyed star forming systems on alarge range of luminosity scales, from star clusters to starbursts,and found that their dust Eddington ratios are consistent with theassumption that radiation pressure regulates star formation.Murray et al. (2005) highlighted the importance of radiationmomentum deposition on dust grains in starburst galaxies. Theyshowed that the Faber-Jackson relation L ∝ σ and the black holemass-stellar velocity dispersion relation M BH ∝ σ could both bemanifestations of self-regulation by radiation pressure. Thompsonet al. (2005) argued that radiation pressure on dust grains can pro-vide vertical support against gravity in disks of starburst galaxiesif the disks are optically thick to the reprocessed IR radiation. Ina study of giant molecular cloud (GMCs) disruption, Murray et al.(2010) found that radiation pressure actually dominates in rapidlystar-forming galaxies such as ULIRGs and submillimeter galaxies.Hopkins et al. (2010) identified a common maximum stellar masssurface density Σ max ∼ M (cid:12) kpc − in a variety of stellar sys-tems ranging from globular clusters to massive star clusters in star-burst galaxies and further to dwarf and giant ellipticals. These sys- c (cid:13) a r X i v : . [ a s t r o - ph . GA ] J u l B. T.-H. Tsang & M. Milosavljevi´c tems spanned ∼ orders of magnitudes in stellar mass and ∼ orders of magnitude in effective radius. The universality of maxi-mum stellar mass surface density can be interpreted as circumstan-tial evidence for the inhibition of gaseous gravitational collapse byradiation pressure.The preceding studies were based on one-dimensional or oth-erwise idealized models. To understand the dynamical effects ofradiation pressure in a dusty ISM, however, multi-dimensional ra-diation hydrodynamics (RHD) simulations are required. One spe-cific setup has emerged as the testbed for radiation hydrodynam-ics numerical methods used in simulating the dusty ISM, specifi-cally in the regime in which the gas (assumed to be thermally cou-pled to dust) is approximately isothermal and susceptible to com-pressive, high-mach-number perturbations. Krumholz & Thomp-son (2012, hereafter KT12) and Krumholz & Thompson (2013,hereafter KT13) designed a two-dimensional model setup to inves-tigate the efficiency of momentum transfer from trapped IR radi-ation to a dusty atmosphere in a vertical gravitational field. Usingthe flux-limited diffusion (FLD) approximation in the ORION code(Krumholz et al. 2007), they found that the optically thick gas layerquickly developed thin filaments via the radiative Rayleigh-Taylorinstability (RTI). The instability produced clumping that allowedradiation to escape through low-density channels. This significantlyreduced net momentum transfer from the escaping radiation to thegas, and the gas collapsed under gravity at the base of the compu-tational box where radiation was being injected.Davis et al. (2014, hereafter D14) then followed up by sim-ulating the same setup with the
ATHENA code (Davis et al. 2012)using the more accurate variable Eddington tensor (VET) approx-imation. They constructed the local Eddington tensor by solvingthe time-independent radiative transfer equation on a discrete setof short characteristics (Davis et al. 2012). Similar to the simula-tions of KT12, those of D14 developed filamentary structures thatreduced radiation-gas momentum coupling. However, in the long-term evolution of the radiation-pressure-forced atmosphere, D14detected significant differences, namely, the gas continued to accel-erate upward, whereas in KT12, it had settled in a turbulent steadystate confined near the base of the box. D14 interpreted this dif-ference of outcome by referring to an inaccurate modeling of theradiation flux in the optically thick-to-thin transition in FLD.Rosdahl & Teyssier (2015, hereafter RT15), simulated thesame setup with the new
RAMSES - RT RHD code using the com-putationally efficient M1 closure for the Eddington tensor. Thismethod separately transports radiation energy density and flux andassumes that the angular distribution of the radiation intensity is aLorentz-boosted Planck specific intensity. One expects this to pro-vide a significant improvement of accuracy over FLD, however stillnot approaching the superior accuracy of the short characteristicsclosure. The M1 results are qualitatively closer to those obtainedwith the FLD than with the short-characteristics VET. RT15 ar-gue that the differences between FLD and M1 on one hand and theshort characteristics VET on the other may be more subtle than sim-ply arising from incorrectly approximating the flux at the opticallythick-to-thin transition.In this paper, we revisit the problem of radiative forcing ofa dusty atmosphere and attempt to reproduce the simulations ofKT12, D14, and RT15, but now with an entirely different numericalscheme, the implicit Monte Carlo (IMC) method of Abdikamalovet al. (2012) originally introduced by Fleck & Cummings (1971).The paper is organized as follows. In Section 2 we review the equa-tions of radiation hydrodynamics and the IMC method. In Section 3we then assess the reliability of our approach in a suite of standard radiation hydrodynamics test problems. In Section 4 we describethe simulation setup and details of numerical implementation. Wepresent our results in Section 5 and provide concluding reflectionsin Section 6.
We start from the equations of non-relativistic radiation hydrody-namics. The hydrodynamic conservation laws are ∂ρ∂t + ∇ · ( ρ v ) = 0 , (1) ∂ρ v ∂t + ∇ · ( ρ vv ) + ∇ P = ρ g + S , (2) ∂ρE∂t + ∇ · [( ρE + P ) v ] = ρ v · g + cS , (3)where ρ , v , and P are the gas density, velocity, and pressure of thegas and E = e + 12 | v | (4)is the specific total gas energy defined as the sum of the specificinternal and kinetic energies of the gas. On the right hand side, g isthe gravitational acceleration and S and cS are the gas momentumand energy source terms arising from the coupling with radiation.The source terms, written here in the lab frame, generally de-pend on the gas density, thermodynamic state, and velocity. Wehere investigate a system that steers clear of the dynamic diffusionregime, namely, in our simulations τ v/c (cid:28) is satisfied at alltimes. Here, τ (cid:46) is the maximum optical depth across the boxand the velocity is non-relativistic v/c (cid:46) − . Therefore, we cansafely drop all O ( v/c ) terms contributing to the momentum sourcedensity and can approximate the lab-frame momentum source den-sity with a velocity-independent gas-frame expression S = 1 c (cid:90) ∞ d(cid:15) (cid:90) d Ω [ k ( (cid:15) ) I ( (cid:15), n ) − j ( (cid:15), n )] n . (5)Here, I ( (cid:15), n ) is the specific radiation intensity, (cid:15) is the photon en-ergy, n is the radiation propagation direction, c is the speed oflight, d Ω is the differential solid angle in direction n , and k ( (cid:15) ) and j ( (cid:15), n ) are the total radiation absorption and emission coefficients.The coefficients are also functions of the gas density ρ and temper-ature T but we omit these parameters for compactness of notation.Our energy source term includes the mechanical work per unitvolume and time v · S that is performed by radiation on gas. Sincethe gas-frame radiation force density is used to approximate thelab-frame value, the source term is correct only to O ( v/c ) cS = (cid:90) ∞ d(cid:15) (cid:90) d Ω [ k ( (cid:15) ) I ( (cid:15), n ) − j ( (cid:15), n )] + v · S . (6)Because in this scheme radiation exerts force on gas but gas doesnot on radiation, the scheme does not conserve energy and momen-tum exactly. However, it should be accurate in the non-relativistic,static-diffusion limit; we test this accuracy in Section 3. We splitthe absorption and emission coefficients by the nature of radiativeprocess k ( (cid:15) ) = k a ( (cid:15) ) + k s ( (cid:15) ) , (7) j ( (cid:15), n ) = j a ( (cid:15) ) + j s ( (cid:15), n ) . (8)The subscript ‘a’ refers to thermal absorption and emission, and ‘s’ c (cid:13) , 000–000 adiation pressure driving of a dusty atmosphere refers to physical scattering (to be distinguished from the effectivescattering that will be introduced in the implicit scheme).Equations (1–3) couple to the radiation subsystem via the radi-ation source terms defined in Equation (5) and (6). Assuming localthermodynamic equilibrium (LTE), the radiation transfer equationcan be written as c ∂I ( (cid:15), n ) ∂t + n · ∇ I ( (cid:15), n ) = k a ( (cid:15) ) B ( (cid:15) ) − k ( (cid:15) ) I ( (cid:15), n )+ j s ( (cid:15), n ) + j ext ( (cid:15), n ) (9)where B ( (cid:15) ) is the Planck function at temperature T and j ext ( (cid:15), n ) is the emissivity of external radiation sources. Note that the j s termdepends implicitly on the specific intensity I which makes Equa-tion (9) an integro-differential equation.Since the absorption and emission coefficients depend on thegas temperature, and the temperature in turn evolves with the ab-sorption and emission of radiation, the system is nonlinear. Wesolve the system by operator-splitting (Section 2.1), by replacinga portion of absorption and emission with effective scattering (thusmaking the solution implicit; Section 2.2), and by discretizing theradiation field with a Monte-Carlo (MC) scheme (Section 2.3). Our numerical method is based on the adaptive-mesh refinement(AMR) code
FLASH (Fryxell et al. 2000; Dubey et al. 2008), ver-sion 4.2.2. We use operator-splitting to solve Equations (1–3) and(9) in two steps:(i)
Hydrodynamic update : Equations (1–3) without the radiationsource terms ∂ρ∂t + ∇ · ( ρ v ) = 0 , (10) ∂ρ v ∂t + ∇ · ( ρ vv ) + ∇ P = ρ g , (11) ∂ρE∂t + ∇ · [( ρE + P ) v ] = ρ v · g (12)are solved using the HYDRO module in
FLASH .(ii)
Radiative transport and source deposition update : Equation(9) coupled to the radiative momentum ρ ∂ v ∂t = S (13)and energy ρ ∂E∂t = cS (14)deposition equations is solved with the implicit method that we pro-ceed to discuss. Under LTE conditions, the tight coupling between radiation and gasis stiff and prone to numerical instability. This limits the applicabil-ity of traditional, explicit methods unless very small time steps areadopted. The method of Fleck & Cummings (1971) for nonlinearradiation transport relaxes the limitation on the time step by treatingthe radiation-gas coupling semi-implicitly. Effectively, this methodreplaces a portion of absorption and immediate re-emission withelastic scattering, thus reducing the amount of zero-sum (quasi-equilibrium) energy exchange between gas and radiation. Numer-ous works have been devoted to investigating the semi-implicit scheme’s numerical properties. Wollaber (2008) provides a detaileddescription of the approximations made and presents a blueprint forimplementation and stability analysis. Cheatham (2010) providesan analysis of the truncation error. Recently, Roth & Kasen (2015)developed variance-reduction estimators for the radiation sourceterms in IMC simulations.In this section, we describe a choice of formalism for solvingthe coupled radiative transport and source term deposition equa-tions. The radiation transport equation is revised to reduce thermalcoupling between radiation and gas by replacing it with a pseudo-scattering process. The detailed derivation of the scheme can befound in Fleck & Cummings (1971) and Abdikamalov et al. (2012).Here we reproduce only the main steps. Our presentation followsclosely Abdikamalov et al. (2012) and the approximations are as inWollaber (2008).Given an initial specific intensity I ( (cid:15), n , t n ) and gas specificinternal energy e ( t n ) at the beginning of a hydrodynamic time step t n , our goal is to solve Equations (9) and (14) to compute the time-advanced values I ( (cid:15), n , t n +1 ) and e ( t n +1 ) at the end of the timestep t n +1 = t n + ∆ t , where ∆ t is the hydrodynamic time step.During this partial update, we assume that the gas density ρ andvelocity v remain constant. For mathematical convenience we in-troduce auxiliary parametrizations of the thermodynamic variables:the gas internal energy density u g = ρe, (15)the energy density that radiation would have if it were in thermo-dynamic equilibrium with gas u r = 4 πc (cid:90) ∞ B ( (cid:15) ) d(cid:15) (16)(for compactness of notation and at no risk of confusion, we do notexplicitly carry dependence on the gas temperature T ), the normal-ized Planck function b ( (cid:15) ) = B ( (cid:15) )4 π (cid:82) ∞ B ( (cid:15) ) d(cid:15) , (17)the Planck mean absorption coefficient k p = (cid:82) ∞ k a ( (cid:15) ) B ( (cid:15) ) d(cid:15) (cid:82) ∞ B ( (cid:15) ) d(cid:15) , (18)and a dimensionless factor, β , quantifying the nonlinearity of thegas-radiation coupling β = ∂u r ∂u g . (19)At a risk of repetition, we emphasize that the u r defined in Equation(16) is not the energy density of the radiation field; it is simply asan alternate parametrization of the gas internal energy density. Theabsorption coefficient and the gas-temperature Planck function, inparticular, can now be treated as functions of u r .Taking the physical scattering to be elastic, Equations (9) and(13), and (14) can be rewritten as c ∂I ( (cid:15), n ) ∂t + n · ∇ I ( (cid:15), n ) = k a ( (cid:15) ) b ( (cid:15) ) cu r − k a ( (cid:15) ) I ( (cid:15), n ) − k s ( (cid:15) ) I ( (cid:15), n ) + j s ( (cid:15), n )+ j ext ( (cid:15), n ) , (20) β ∂u r ∂t + ck p u r = (cid:90) ∞ d(cid:15) (cid:90) d Ω k a ( (cid:15) ) I ( (cid:15), n ) . (21) c (cid:13) , 000–000 B. T.-H. Tsang & M. Milosavljevi´c
We linearize the equations in u r and I ( (cid:15), n ) and denote the corre-sponding constant coefficients with tildes, c ∂I ( (cid:15), n ) ∂t = − n · ∇ I ( (cid:15), n ) + ˜ k a ( (cid:15) )˜ b ( (cid:15) ) cu r − ˜ k a ( (cid:15) ) I ( (cid:15), n ) − ˜ k s ( (cid:15) ) I ( (cid:15), n ) + j s ( (cid:15), n ) + j ext ( (cid:15), n ) , (22) β ∂u r ∂t = − c ˜ k p u r + (cid:90) ∞ d(cid:15) (cid:90) d Ω ˜ k a ( (cid:15) ) I ( (cid:15), n ) . (23)The scattering emission coefficient is j s ( (cid:15), n ) = (cid:90) d Ω (cid:48) ˜ k s ( (cid:15) ) Ξ( (cid:15), n , n (cid:48) ) I ( (cid:15), n (cid:48) ) , (24)where Ξ( (cid:15), n , n (cid:48) ) is the elastic scattering kernel. We evaluate theconstant coefficients explicitly at t n , the beginning of the time step.Next, for t n (cid:54) t (cid:54) t n +1 , we expand u r to the first order intime u r ( t ) (cid:39) u n r + ( t − t n ) u (cid:48) r n , (25)where u n r = u r ( t n ) and u (cid:48) r n = ∂u r /∂t ( t n ) . It is worth notingthat the implicitness in ‘IMC’ refers to that introduced in Equation(25). Substituting u r ( t ) from Equation (25) into (23), solving for u (cid:48) r n , and substituting the result back into Equation (25), we obtain u r = fu n r + 1 − fc ˜ k p (cid:90) ∞ d(cid:15) (cid:90) d Ω ˜ k a ( (cid:15) ) I ( (cid:15), n ) , (26)where f is a time-dependent factor f ( t ) = 11 + ( t − t n ) ˜ βc ˜ k p . (27)Since it is desirable to work with time-independent coefficients,we approximate f ( t ) with the so-called Fleck factor that remainsconstant during the time step f (cid:39)
11 + α ∆ t ˜ βc ˜ k p , (28)where (cid:54) α (cid:54) is a coefficient that interpolates between thefully-explicit ( α = 0 ) and fully-implicit ( α = 1 ) scheme for updat-ing u r . For intermediate values of α , the scheme is semi-implicit.The scheme is stable when . (cid:54) α (cid:54) (Wollaber 2008).We substitute u r from Equation (26) into Equation (20) to ob-tain an equation for I ( (cid:15), n ) in the form known as the implicit radi-ation transport equation c ∂I ( (cid:15), n ) ∂t + n · ∇ I ( (cid:15), n ) =˜ k ea ( (cid:15) )˜ b ( (cid:15) ) cu n r − ˜ k ea ( (cid:15) ) I ( (cid:15), n ) − ˜ k es ( (cid:15) ) I ( (cid:15), n )+ ˜ k a ( (cid:15) )˜ b ( (cid:15) )˜ k p (cid:90) ∞ d(cid:15) (cid:48) (cid:90) d Ω (cid:48) ˜ k es ( (cid:15) (cid:48) ) I ( (cid:15) (cid:48) , n (cid:48) ) − ˜ k s ( (cid:15) ) I ( (cid:15), n ) + j s ( (cid:15), n ) + j ext ( (cid:15), n ) , (29)where the effective absorption and scattering coefficients are ˜ k ea ( (cid:15) ) = f ˜ k a ( (cid:15) ) , (30) ˜ k es ( (cid:15) ) = (1 − f ) ˜ k a ( (cid:15) ) . (31)Equation (29) admits an instructive physical interpretation.The first two terms on the right-hand side represent the emissionand absorption of thermal radiation. Direct comparison with Equa-tion (20) shows that both terms are now a factor of f smaller.The following two terms containing ˜ k es are new; their functionalform mimics the absorption and immediate re-emission describing a scattering process. Meanwhile, the physical scattering and exter-nal source terms have remained unmodified.Since the effective absorption ˜ k ea and scattering ˜ k es coeffi-cients sum to the actual total absorption coefficient ˜ k a , we can in-terpret Equation (29) as replacing a fraction ( − f ) of absorp-tion and the corresponding, energy-conserving fraction of emissionby an elastic psedo-scattering process. The mathematical form ofthe Fleck factor can be rearranged to make this physical interpreta-tion manifest. Assuming an ideal gas equation of state, the radiativecooling time is t cool = 4 c ˜ β ˜ k p (32)and the Fleck factor is f = 11 + 4 α ∆ t/t cool . (33)When ∆ t/t cool (cid:29) so that f (cid:28) , the absorbed radiation is re-radiated within the same time step at zero net change in the gasenergy density; the only change is randomization of the radiationpropagation direction. The stability of the scheme rests precisely onthis reduction of the stiff thermal coupling. However, excessivelylarge time steps can still produce unphysical solutions (Wollaber2008).After the radiation transport equation has been solved usingthe IMC method (see Section 2.3), the net momentum and energyexchange collected during the radiative transport solve, which read S = 1 c ∆ t (cid:90) ∞ d(cid:15) ˜ k ( (cid:15) ) (cid:90) d Ω (cid:90) t n +1 t n dtI ( (cid:15), n ) n (34)and cS = − πcu n r (cid:90) ∞ d(cid:15) ˜ k ea ( (cid:15) ) ˜ b ( (cid:15) )+ 1∆ t (cid:90) ∞ d(cid:15) ˜ k ea ( (cid:15) ) (cid:90) d Ω (cid:90) t n +1 t n dtI ( (cid:15), n )+ v · S , (35)are deposited in the hydrodynamic variables. Therefore step (ii) inour operator splitting scheme has now been further split into twosub-steps:(ii (cid:48) ) Radiative transport and hydrodynamical source term col-lection : Solve Equation (29) with the IMC method while accumu-lating the contribution of radiative processes to gas source terms asin Equations (34) and (35).(ii (cid:48)(cid:48) ) Hydrodynamical source term deposition : Update gas mo-mentum and energy density using Equations (13) and (14).
The transition layer between the optically thick and thin regimesstrains the adequacy of numerical radiation transfer methods basedon low-order closures. In this transition layer, the MC radia-tive transfer method should perform better than computationally-efficient schemes that discretize low-order angular moments ofEquation (9). In the MC approach, one obtains solutions of theradiation transport equation by representing the radiation fieldwith photon packets and modeling absorption and emission withstochastic events localized in space and/or time. This permits accu-rate and straightforward handling of complicated geometries and,in greater generality than we need here, angle-dependent physi-cal processes such as anisotropic scattering. The specific intensity c (cid:13) , 000–000 adiation pressure driving of a dusty atmosphere I ( (cid:15), n ) is represented with an ensemble of a sufficiently large num-ber of MCPs. In the radiation transfer update, starting with the radiation fieldat an initial time t n , we wish to compute the coupled radiation-gassystem at the advanced time t n +1 = t n + ∆ t . Our MC schemefollows closely that of Abdikamalov et al. (2012) and Wollaber(2008). The radiation field is discretized using a large number ofMonte Carlo particles (MCPs), each representing a collection ofphotons. We adopt the grey approximation in which we track onlythe position and the collective momentum of the photons in eachMCP. MCPs are created, destroyed, or their properties are modifiedas needed to model emission, absorption, scattering, and propaga-tion of radiation.If a finite-volume method is used to solve the gas conserva-tion laws, the physical system is spatially decomposed into a finitenumber of cells. For the purpose of radiative transport, gas prop-erties are assumed to be constant within each cell. Particles arecreated using cell-specific emissivities. Each MCP is propagatedalong a piecewise linear trajectory on which the gas properties (ab-sorption and scattering coefficients) are evaluated locally. Our MCscheme computes an approximation to the solution of Equation (29)in two steps: by first creating MCPs based on the emissivities andboundary conditions, and then transporting MCPs through spaceand time. In Equation (29), the term ˜ k ea ˜ bcu n r on the right-hand side is thefrequency-dependent thermal emissivity. Assuming that thermalemission is isotropic ( ˜ k ea is angle-independent), the total thermalradiation energy emitted by a single cell of gas ∆ E can be calcu-lated as ∆ E = 4 π ∆ t ∆ V (cid:90) ∞ ˜ k ea ( (cid:15) ) B ( (cid:15) ) d(cid:15), (36)where ∆ t is the time step size, ∆ V is the cell volume, and B ( (cid:15) ) isthe Planck function at the gas temperature ˜ T ( t n ) . Since we furtherassume that the opacity is grey (independent of (cid:15) ) then ∆ E = c ∆ t ∆ V ˜ k ea u n r . (37)The net momentum exchange due to thermal emission is zero be-cause the thermal radiation source is isotropic.We specify that in thermal emission, N new MCPs are cre-ated in each cell in each time step. The energy carried by each newMCP is then ∆ E / N . The emission time of each such MCP is sam-pled uniformly within the interval [ t n , t n +1 ] . The spatial positionof the MCP is sampled uniformly within the cell volume and thepropagation direction is sampled uniformly on a unit sphere. Ev-ery MCP keeps track of its time t i , current position r i , momentum p i , and fraction of the energy remaining since initial emission ς i , One disadvantage of the MC scheme is low computational efficiency inthe optically thick regime where the photon mean free path is short. Effi-ciency in such regions can be improved by applying the diffusion approxi-mation (Fleck & Canfield 1984; Gentile 2001; Densmore et al. 2007). Re-cently, Abdikamalov et al. (2012) interfaced the IMC scheme at low opticaldepths with the Discrete Diffusion Monte Carlo (DDMC) method of Dens-more et al. (2007) at high optical depths. This hybrid algorithm has beenextended to Lagrangian meshes (Wollaeger et al. 2013). In the present ap-plication, the optical depths are relatively low and shortness of the meanfree path is not a limitation. where the index i ranges over all the MCPs active in a given hydro-dynamic time step. Newly created MCPs are added to the pool ofMCPs carried over from previous hydrodynamic time steps. To minimize noise, we treat absorption deterministically. This ‘con-tinuous absorption’ method is a variance reduction technique com-mon in practical implementations of IMC (Abdikamalov et al.2012; Hykes & Densmore 2009). Specifically, when an MCP trav-els a distance cδt i inside a cell with absorption coefficient ˜ k ea , itsmomentum is attenuated according to p i ( t i + δt i ) = p i ( t i ) e − ˜ k ea cδt i , (38)where we denote an arbitrary time interval with δt i to distinguishit from the hydrodynamic time step ∆ t . In a single hydrodynamic time step, the simulation transports theMCPs through multiple cells. The MCP-specific time remaininguntil the end of the hydrodynamic time step is t n +1 − t i . For eachMCP, we calculate or sample the following four distances:(i) The free streaming distance to the end of the hydrodynamictime step d t = c ( t n +1 − t i ) .(ii) Distance to the next scattering event assuming cell-localscattering coefficients d s = − ln ξk s + ˜ k es , (39)where ξ is a random deviate uniformly distributed on the interval (0 , .(iii) Distance to near-complete absorption d a defined as the dis-tance over which only a small fraction ς min = 10 − of the initialenergy remains.(iv) Distance to the current host cell boundary d b .We repeatedly update the four distances, select the shortestone, and carry out the corresponding operation until we reach theend of the hydrodynamic time step t n +1 . If in such a sub-cyclethe shortest distance is d t , we translate the MCP by this distance r i → r i + d t n i , where n i = p i /p i is the propagation direction.We also attenuate its momentum according to Equation (38) and ac-crue the momentum − ∆ p i, a and energy | ∆ p i, a | c transferred to thegas. If the shortest distance is d s , we do the same over this distance,but at the end of translation, we also randomize the MCP’s direc-tion n i → n (cid:48) i and accrue the corresponding additional momentum − ∆ p i, s = p i ( n (cid:48) i − n i ) and kinetic energy − v · ∆ p i, s transferredto gas. As a further variance-reduction tactic, given the statisticalisotropy of n (cid:48) i , we compute the momentum deposited in a scatter-ing event simply as − ∆ p i, s = − p i . If the shortest distance eitheris d a or d b , we translate the MCP while attenuating its momentumand accruing the deposited energy and momentum. Then we eitherremove the MCP while instantaneously depositing the remainingmomentum and energy to the gas (if the shortest distance is d a ), ortransfer the MCP to its new host cell (or removed the MCP if it hasreached a non-periodic boundary of the computational domain).As the MCPs are transported over a hydrodynamic time step ∆ t , the energy and momentum source terms for each cell are accu- c (cid:13) , 000–000 B. T.-H. Tsang & M. Milosavljevi´c
0 0.1 0.2 0.3 0.4 0.5 1 2 3 4 5 6 7 8
Radiation Energy Denstiy [erg cm -3 ] Refinement Level Radius [cm]
Figure 1.
Spherically-averaged radiation energy density profiles in thethree-dimensional radiative diffusion test (Section 3.1). The analytical so-lutions (solid lines) and the numerical solutions (crosses) are shown at fourdifferent times, (0 . , . , . , . × − s . The dashed line and theright axis show the refinement level of the AMR grid. mulated using cS = 1∆ t ∆ V (cid:88) ( | ∆ p i, a | c − v · ∆ p i, s ) , (40) S = − t ∆ V (cid:88) (∆ p i, a + ∆ p i, s ) , (41)where the sums are over all the absorption and scattering events thatoccurred in a specific computational cell during the hydrodynamictime step. The source terms are then substituted into Equations (13)and (14) to compute the gas momentum and energy at the end ofthe hydrodynamic time step. To assess the validity of our radiation hydrodynamics implementa-tion, we performed a series of standard tests: a test of radiative dif-fusion in a scattering medium (Section 3.1), a test of gas-radiationthermal equilibration (Section 3.2), a test of thermal wave propaga-tion (Marshak wave; Section 3.3), and a radiative shock test (Sec-tion 3.4).
Here we test the spatial transport of MCPs across the AMR gridstructure in the presence of scattering. In the optically thick limit,radiation transfer proceeds as a diffusion process. The setup is a cu-bical L = 1 cm three-dimensional AMR grid with no absorptionand a scattering coefficient of k s = 600 cm − . The scattering is as-sumed to be isotropic and elastic. We disable momentum exchangeto preclude gas back-reaction and focus on testing the evolution ofthe radiation field on a non-uniform grid.At t = 0 s, we deposit an initial radiative energy ( E init =3 . × erg) at the grid center in the form of 1,177,600 MCPs withisotropically sampled propagation directions. We lay an AMR grid -12 -10 -8 -6 -4 -2 Energy Density [erg cm -3 ] Time [s] gasradiation
Figure 2.
Evolution of the gas (diamonds) and radiation (squares) energydensity in the one-zone radiative equilibration test (Section 3.2). The exactsolution is shown with a solid (gas) and dashed (radiation) line. hierarchy such that the refinement level decreases with increasingradius as shown on the right axis of Figure 1. The grid spacing is ∆ x = 2 − (cid:96) − L , where (cid:96) is the local refinement level. A constanttime step of ∆ t = 2 × − s is used and the simulation is run for × − s.The diagnostic is the radiation energy density profile as a func-tion of distance from the grid center and time ρe rad ( r, t ) . The exactsolution in d spatial dimensions is given by ρe rad ( r, t ) = E init (4 πDt ) d/ exp (cid:18) − r Dt (cid:19) , (42)where D = c/ ( dk s ) is the diffusion coefficient.Figure 1 shows the spherically-averaged radiation energy den-sity profile at four times. Excellent agreement of our MC resultswith the analytical expectation shows that our algorithm accuratelycaptures radiation transport in a scattering medium. We have re-peated the test in one and two spatial dimensions and find the sameexcellent agreement. We have also checked that in multidimen-sional simulations, the radiation field as represented with MCPspreserves the initial rotational symmetry. Here, in a one-zone setup, we test the radiation-gas thermal cou-pling in LTE. We enabled IMC with an implicitness parameter of α = 1 . Defining u r = aT as in Section 2.2, where a is the ra-diation constant and T is the gas temperature, the stiff system ofequations governing the gas and radiation internal energy densityevolution is dedt = k a c (cid:18) e rad − u r ρ (cid:19) , (43) de rad dt = k a c (cid:18) u r ρ − e rad (cid:19) , (44)where k a is the absorption coefficient. We assume an ideal gas withadiabatic index γ = 5 / . c (cid:13) , 000–000 adiation pressure driving of a dusty atmosphere -5 -4 -3 -2 -1 u k a z Figure 3.
Radiation energy density profiles in the Marshak wave test prob-lem at times θ = (3 , , from left to right. Numerical integrations areshown by the data points. The solid lines are the reference solutions of Su& Olson (1996). We perform the one-zone test with parameters similar to thoseof Turner & Stone (2001) and Harries (2011), namely, the absorp-tion coefficient is k a = 4 . × − cm − and the initial energydensities are ρe = 10 erg cm − and ρe rad = 0 , respectively. Theresults and the corresponding exact solutions are shown in Figure2. The MC solution agrees with the exact solutions within (cid:46) throughout the simulation. It shows that the physics of radiation-gas thermal exchange is captured well by our scheme, and that instatic media, the scheme conserves energy exactly.As noted by Cheatham (2010), the order of accuracy associ-ated with the IMC method depends on the choice of α and on thespecifics of the model system. When the above test problem is re-peated with α = 0 . , the error is (cid:46) . because the O (∆ t ) -residuals cancel out and the method is O (∆ t ) -accurate. In this test, we simulate the propagation of a non-linear thermalwave, known as the Marshak wave, in a static medium in one spa-tial dimension (Su & Olson 1996; Gonz´alez et al. 2007; Krumholzet al. 2007; Zhang et al. 2011). The purpose of this standard testis to validate the code’s ability to treat nonlinear energy couplingbetween radiation and gas when the gas heat capacity is a functionof gas temperature. We employ IMC with α = 1 .Initially, a static, uniform slab with a temperature of 10 K oc-cupying the interval (cid:54) z (cid:54) cm is divided into 256 equalcells. An outflow boundary condition is used on the left and a re-flective one on the right. A constant incident flux F inc = σ SB T of k B T inc = 1 keV thermal radiation, where σ SB and k B are theStefan-Boltzmann and Boltzmann constants, is injected from theleft (at z = 0 ). The gas is endowed with a constant, grey absorp-tion coefficient of k a = 1 cm − and a temperature-dependent vol-umetric heat capacity of c v = αT . The constant α is related to theSu-Olson retardation parameter (cid:15) via α = 4 a/(cid:15) and we set (cid:15) = 1 .The diagnostics for this test problem are the spatial radia-tion and gas energy density profiles at different times. Su & Olson(1996) provided semi-analytical solutions in terms of the dimen-sionless position x = √ k a z (45) -5 -4 -3 -2 -1 v k a z Figure 4.
The same as Figure 3, but for the gas energy density. and time θ = 4 ack a tα . (46)The radiation and gas internal energy density are expressed in termsof the dimensionless variables u ( x, θ ) = c E r ( x, θ ) F inc , (47) v ( x, θ ) = c aT ( x, θ ) F inc , (48)where E r and T are the radiation energy density and gas temper-ature, respectively (note that the relation of v to the specific gasinternal energy e is nonlinear).Figures 3 and 4 show the dimensionless numerical solutionprofiles at three different times, over-plotting the correspondingsemi-analytical solutions of Su & Olson (1996). At early times,we observe relatively large deviations from the Su-Olson solutionsnear the thermal wavefront. This is in not surprising given that thesolutions were obtained assuming pure radiative diffusion, yet atearly times and near the wavefront, where the gas has not yet heatedup, the optical depth is only about unity and the transport is not dif-fusive. Gonz´alez et al. (2007) observed the same early deviationsin their computations based on the M1 closure. At later times whentransport is diffusive, both the thermal wave propagation speed andthe maximum energy density attained agree well with the Su-Olsonsolutions. To finally test the fully-coupled radiation hydrodynamics, we sim-ulate a radiative shock tube. As in the preceding tests, we useIMC with α = 1 , but now, the gas is allowed to dynamicallyrespond to the radiation. We adopt the setup and initial condi-tions of Ensman (1994) and Commerc¸on et al. (2011) and simu-late both subcritical and supercritical shocks. The setup consists ofa one-dimensional × cm-long domain containing an idealgas with a mean molecular weight of µ = 1 and adiabatic in-dex of γ = 7 / . The domain is initialized with a uniform massdensity ρ = 7 . × − g cm − and a uniform temperatureof T = 10 K. The gas has a constant absorption coefficient of k a = 3 . × − cm − and a vanishing physical scattering coef-ficient.Initially, the gas is moving with a uniform velocity toward theleft reflecting boundary. An outflow boundary condition is adopted c (cid:13) , 000–000 B. T.-H. Tsang & M. Milosavljevi´c
Table 1.
Simulation parameters for the radiation-driven gaseous atmosphere.Run Initial Perturbation Σ g h ∗ t ∗ τ ∗ f E , ∗ t max /t ∗ [ L x × L y ]/ h ∗ ∆ x/ h ∗ (cid:96) max (g cm − ) ( − dyne g − ) ( − pc) kyrT10F0.02 sin 4.7 37 0.25 0.045 10 0.02 80 512 ×
256 0.5 7T03F0.50 sin, χ × Temperature [K] z [10 cm] GasRadiation Figure 5.
Gas (solid curve) and radiation (dashed curve) temperature in thesubcritical radiative shock test at t = 3 . × s. The initial gas velocityis v = 6 km s − and the profiles are plotted as a function of z = x − vt . on the right to allow inflow of gas at fixed density ρ and fixed tem-perature T and also to allow the free escape of radiation MCPs. Asgas collides with the reflecting boundary a shock wave starts propa-gating to the right. The thermal radiation in the compressed hot gasdiffuses upstream and produces a warm radiative precursor. Theshock becomes critical when the flux of thermal radiation is highenough to pre-heat the pre-shock gas to the post-shock temperature(Zel’dovich & Raizer 1967). We choose the incoming speed to be v = 6 km s − and km s − in the subcritical and the supercriti-cal shock tests, respectively.Mihalas & Mihalas (1984) provide analytical estimates for thecharacteristic temperatures of the radiative shocks. For the subcrit-ical case, the post-shock temperature T is estimated to be T (cid:39) γ − v R ( γ + 1) . (49)Using the parameters for the subcritical setup, the analytical esti-mate gives T (cid:39) K. In our simulation, the post-shock tem-perature at t = 3 . × s is T (cid:39) K, which agrees withthe analytical solution. The immediate pre-shock temperature T − is estimated to be T − (cid:39) γ − √ Rρv σ SB T . (50)Our simulation gives T − ∼ K while T − is estimated to be T − = 270 K. Finally, the amplitude of the temperature spike can
Temperature [K] z [10 cm] GasRadiation Figure 6.
The same as Figure 5, but for the supercritical radiative shocktest at t = 7 . × s and with initial velocity v = 20 km s − . be estimated to be T + (cid:39) T + 3 − γγ + 1 T − , (51)which gives T + (cid:39) K. It also close to the value we find, T + (cid:39) K. In both cases, our simulations reproduce the ex-pected radiative precursors. Also, in the supercritical case, the pre-shock and the post-shock temperatures are identical, as expected.
We turn to the problem of how radiation drives an interstellargaseous atmosphere in a vertical gravitational field. The problemwas recently investigated by KT12 and KT13, by D14, and byRT15, using the FLD, VET, and M1 closure, respectively. Our aimis to attempt to reproduce these authors’ results, which are all basedon low-order closures, using an independent method that does notrely on such a closure. Critical for the hydrodynamic impact of ra-diation pressure is the extent of the trapping of IR radiation by dustygas. Therefore we specifically focus on the radiation transfer aspectof the problem and assume perfect thermal and dynamic couplingbetween gas and dust grains, T g = T d = T and v g = v d = v .We follow the setup of KT12 and D14 as closely as possible.Taking that UV radiation from massive stars has been reprocessedinto the IR at the source, we work in the grey approximation inwhich spectral averaging of the opacity is done only in the IR partof the spectrum. We set the Rosseland κ R and Planck κ P mean dust c (cid:13) , 000–000 adiation pressure driving of a dusty atmosphere opacities to κ R , P = (0 . , . (cid:18) T K (cid:19) cm g − . (52)This model approximates a dusty gas in LTE at T (cid:54) K (Se-menov et al. 2003). Diverging slightly from KT12 and D14, whoadopted the pure power-law scaling in Equation (52), to approxi-mate the physical turnover in opacity, we cap both mean opacitiesto their values at
K above this threshold temperature. Overall,our opacity model is reasonable below the dust grain sublimationtemperature ∼ L x × L y . The grid is adaptively refined using the standard FLASH second derivative criterion in the gas density. The dustygas is initialized as a stationary isothermal atmosphere. A time-independent, vertically incident radiation field is introduced at thebase of the domain ( y = 0 ) with a flux vector F ∗ ˆ y . The gravita-tional acceleration is − g ˆ y .For notational convenience, we define a reference temperature T ∗ = [ F ∗ / ( ca )] / , sound speed c ∗ = (cid:112) k B T ∗ / ( µm H ) , scaleheight h ∗ = c ∗ /g , density ρ ∗ = Σ /h ∗ (where Σ is the ini-tial average gas surface density at the base of the domain), andsound crossing time t ∗ = h ∗ /c ∗ . In the present setup F ∗ =2 . × L (cid:12) kpc − and the mean molecular weight is µ = 2 . as expected for molecular hydrogen with a 10% helium molar frac-tion. The characteristic temperature is T ∗ = 82 K and the corre-sponding Rosseland mean opacity is κ R , ∗ = 2 .
13 cm g − .Following KT12 and KT13, we adopt two dimensionless pa-rameters to characterize the system: the Eddington ratio f E , ∗ = κ R , ∗ F ∗ gc (53)and the optical depth τ ∗ = κ R , ∗ Σ . (54)The atmosphere is initialized at a uniform temperature T ∗ . The gasdensity is horizontally perturbed according to ρ ( x, y ) = (cid:20) χ (cid:18) πxλ x (cid:19)(cid:21) × (cid:40) ρ ∗ e − y/h ∗ , if e − y/h ∗ > − ,ρ ∗ − , if e − y/h ∗ (cid:54) − , (55)where λ x = 0 . L x . D14 introduced χ , a random variate uniformlydistributed on [ − . , . , to provide an additional perturbationon top of the sinusoidal profile. If χ = 0 , the initial density distri-bution reduces to that of KT12.KT12 found that at a given τ ∗ , the preceding setup has a hy-drostatic equilibrium solution when f E , ∗ is below a certain criticalvalue f E , crit . Note that KT12 defined f E , crit assuming the purepower-law opacity scaling in Equation (52). With our capping ofthe opacities above 150 K, which breaks the dimensionless natureof the KT12 setup, the exact KT12 values for f E , crit cannot bedirectly transferred to our model. Nevertheless, we use their defini-tion of f E , crit simply to normalize our values of τ ∗ and f E , ∗ . Weattempt to reproduce the two runs performed by KT12. The firstrun T10F0.02 with τ ∗ = 10 and f E , ∗ = 0 .
02 = 0 . f E , crit lies inthe regime in which such a hydrostatic equilibrium solution exists.The second run T03F0.5 with τ ∗ = 3 and f E , ∗ = 0 . . f E , crit corresponds to the run performed by both KT12 and D14 that hadthe smallest ratio f E , ∗ /f E , crit and was still unstable. The latter runprobes the lower limit for the occurrence of a dynamically unstablecoupling between radiation and gas. Figure 7.
Gas density snapshots at four different times in the stable runT10F0.02. The fully simulation domain is larger than shown, 512 × h ∗ ;here, we only show the bottom quarter. The stable outcome of this run isconsistent with the cited literature. The boundary conditions are periodic in the x direction, re-flecting at y = 0 (apart from the flux injection there), and out-flowing (vanishing perpendicular derivative) at y = L y . The re-flecting condition does not allow gas flow or escape of radiation.The outflow condition does allow free inflow or outflow of gasand escape of radiation. Unlike the cited treatments, we used non-uniform AMR. The AMR improves computational efficiency earlyin the simulation when dense gas occupies only a small portion ofthe simulated domain. In low-density cells, radiation streams al-most freely through gas; there, keeping mesh resolution low min-imizes the communication overhead associated with MCP han-dling while still preserving MCP kinematic accuracy. The appli-cation of AMR in conjunction with IMC is clearly not essential intwo-dimensional, low-dynamic-range setups like the one presentedhere, but should become critical in three-dimensional simulationsof massive star formation; thus, we are keen to begin validating iton simple test problems.To further economize computational resources, we require adensity (cid:62) − ρ ∗ for thermal emission, absorption, and scatteringcalculations; below this density, the gas is assumed to be adiabaticand transparent. We also apply a temperature floor of 10 K.As the simulation proceeds, the total number of MCPs in-creases. To improve load balance, we limit the maximum numberof MCPs allowed in a single computational block ( × cells) atthe end of the time step to 64, or on average ∼ MCP per cell. (Amuch larger number of MCPs can traverse the block in the courseof a time step.) If the number exceeds this specified maximum atthe end of the radiation transport update, we merge some of theMCPs in a momentum- and energy-conserving fashion. We, how-ever, do not properly preserve spatial and higher-angular-momentstatistical properties of the groups of MCPs subjected to merging.This deficiency is tolerable in the present simulation where merg-ing takes place only at the lowest level of refinement, where theradiation no longer affects the gas. In future applications, however,we will develop a manifestly more physical MCP merging strategy.The simulation parameters of the two runs are summarizedin Table 1. The quoted cell width ∆ x is that at the highest level ofmesh refinement (the cells are square). Gas with density (cid:38) − ρ ∗ always resides at the maximum refinement level (cid:96) max throughoutthe simulation of duration t max . c (cid:13) , 000–000 B. T.-H. Tsang & M. Milosavljevi´c
Figure 9.
Gas density in the unstable run T03F0.50 at four different times as indicated. The panels display the full simulation domain. The flow morphologyis consistent with that observed by D14 and RT15. Downward Rayleigh-Taylor plumes are visible in the second panel from the left.
Density snapshots at four different times in the simulation areshown in Figure 7. The simulation closely reproduces the quantita-tive results of both FLD (KT12) and VET (D14). Shortly after thebeginning of the simulation, the trapping of radiation at the bottomof the domain by the dense dusty gas produces a rise in radiationenergy density. As we assume LTE and perfect thermal couplingbetween gas and dust, the gas temperature increases accordingly.Specifically, after being heated up by the incoming radiation, theeffective opacity at the midplane rises by a factor of ∼ . Theopacity rise enhances radiation trapping and the temperature risesstill further to ∼ (3 − T ∗ ∼
300 K . The heating drives the atmo-sphere to expand upward, but radiation pressure is not high enoughto accelerate the slab against gravity. After the initial acceleration,the atmosphere deflates into an oscillatory, quasi-equilibrium state.This outcome is consistent with what has been found with FLD andVET. To better quantify how the dynamics and the degree of turbu-lence in the gas compare with the results of the preceding investi-gations, we compute the mass-weighted mean gas velocity (cid:104) v (cid:105) = 1 M (cid:90) L y (cid:90) L x ρ ( x, y ) v ( x, y ) dxdy (56)and linear velocity dispersion σ i = 1 M (cid:90) L y (cid:90) L x ρ ( x, y )( v i ( x, y ) − (cid:104) v i (cid:105) ) dxdy, (57)where M is the total mass of the atmosphere and i indexes thecoordinate direction. We also define the total velocity dispersion σ = (cid:112) σ x + σ y .The time evolution of (cid:104) v y (cid:105) and the linear dispersions is shownin Figure 8. All the velocity moments are expressed as fractions ofthe initial isothermal sound speed c ∗ = 0 . km s − . Both panelsclosely resemble those in Figure 2 of D14. Early on at ∼ t ∗ , ra-diation pressure accelerates the gas and drives growth in (cid:104) v y (cid:105) , σ y ,and σ x . After this transient acceleration, (cid:104) v y (cid:105) executes dampledoscillations about zero velocity (the damping is likely of a numer- c (cid:13) , 000–000 adiation pressure driving of a dusty atmosphere -0.5 0 0.5 1 1.5 < v y > / c * [ σ , σ x , σ y ] / c * t / t * σσ y σ x Figure 8.
Top panel:
Time evolution of the mass-weighted mean velocityin the vertical direction in the stable run T10F0.02.
Bottom panel:
The cor-responding time evolution of the mass-weighted mean velocity dispersion.The linear dispersions σ x and σ y are shown with the dotted and dashedlines, respectively. -5 0 5 10 15 < v y > / c * IMCM1 - trappingM1 - no trapping VETFLD 0 2 4 6 8 10 12 0 20 40 60 80 100 120 [ σ , σ x , σ y ] / c * t / t * σ σ y σ x Figure 10.
Top panel:
Time evolution of the mass-weighted mean velocityin the vertical direction in the unstable run T03F0.50 (black line). The col-ored lines are tracks from the cited references (see text and legend).
Bottompanel:
Mass-weighted mean velocity dispersions (see legend). In both pan-els, the late-time net acceleration and velocity dispersions are in agreementonly with the results obtained by D14 with their short characteristics-basedVET method. ical origin). The linear velocity dispersions also oscillate, but withsmaller amplitudes (cid:46) . c ∗ . The oscillation period in σ x is justslightly longer than that reported by D14. The agreement of our andVET results demonstrates the reliability of both radiation transfermethods. To facilitate direct comparison with D14 and RT15, we introducea random initial perturbation on top of the initial sinusoidal pertur-bation as in Equation (55). The grid spacing ∆ x is twice of thatadopted by KT12, D14, and RT15, but as we shall see, this coarserspacing is sufficient to reproduce the salient characteristics of theevolving system. MCP merging is activated at t = 36 t ∗ . Figure9 shows density snapshots at four different times. As in the stablecase, the incoming radiation heats the gas at the bottom of the do- f E , V IMCM1 - trappingM1 - no trapping VETFLD 5 6 7 8 9 10 11 τ V τ F / τ V t / t * Figure 11.
Top panel:
Time evolution of the volume-weighted Eddingtonratio in the unstable run T03F0.50. The colored lines are tracks from thecited references (see text and legend).
Middle panel:
The volume-weightedmean total vertical optical depth.
Bottom panel:
The flux-weighted meanoptical depth. main and opacity jumps. The flux soon becomes super-Eddingtonand a slab of gas is lifted upward. At t ∗ , fragmentation of theslab by the RTI is apparent; most the gas mass becomes concen-trated in dense clumps.We note that the slab lifting and the subsequent fragmenta-tion are consistently observed in all radiative transfer approaches;differences become apparent only in the long-term evolution. Asin the VET simulation (D14), coherent gaseous structures in oursimulation continue to be disrupted and accelerated. Qualitatively,radiation drives gas into dense, low-filling-factor filaments embed-ded in low density (10 − − − ) ρ ∗ gas. At 115 t ∗ , the bulk ofthe gas has a net upward velocity and has been raised to altitudes y ∼ h ∗ .Figure 10 compares the time evolution of the bulk velocity (cid:104) v y (cid:105) and velocity dispersions σ x,y with the corresponding tracksfrom the published FLD/VET and M1 simulations (respectively,D14 and RT15). Initially, (cid:104) v y (cid:105) rises steeply as the gas slab heats upand the incoming flux becomes super-Eddington. All simulationsexcept for the one performed with the M1 closure without radia-tion trapping (RT15) exhibit a similar initial rise. At ∼ t ∗ , theRTI sets in and the resulting filamentation reduces the degree ofradiation trapping. This in turn leads to a drop in radiation pres-sure and (cid:104) v y (cid:105) damps down under gravity. The transient rise anddrop in (cid:104) v y (cid:105) is observed with all the radiative transfer methods,although the specific times of the acceleration-to-deceleration tran-sition differ slightly. The bulk velocity peaks at (cid:104) v y (cid:105) (cid:39) c ∗ inIMC and at (cid:39) c ∗ in VET. The subsequent kinematics differs sig-nificantly between the methods. In IMC and VET, the gas filamentsrearrange in a way that enables resumption of upward accelerationafter ∼ (50 − t ∗ . At late times, the secondary rise in (cid:104) v y (cid:105) doesnot seem to saturate in IMC as it does in VET. Otherwise, the IMCand VET tracks are very similar to each other. In FLD and M1,however, gas is not re-accelerated after the initial transient accel-eration. Instead, it reaches a turbulent quasi-steady state in whichgas is gravitationally confined at the bottom of the domain and (cid:104) v y (cid:105) fluctuates around zero.The evolution of velocity dispersions in IMC is also in closeagreement with VET. Before the RTI onset, the dispersions riseslightly to σ y (cid:38) c ∗ . Once the RTI develops and the slab frag-ments, σ y increases rapidly and σ x somewhat more gradually. A c (cid:13) , 000–000 B. T.-H. Tsang & M. Milosavljevi´c drop in σ y is observed at ∼ t ∗ , but after that time, the verticaldispersion rises without hints of saturation. Velocity dispersions atthe end of our simulations are consistent with those in VET. In FLDand M1, on the other hand, the asymptotic turbulent quasi-steadystates have smaller velocity dispersions (cid:39) c ∗ .To further investigate the coupling of gas and radiation, wefollow KT12 and KT13 to define three volume-weighted quantities:the Eddington ratio f E , V = (cid:104) κ R ρF y (cid:105) V cgρ , (58)the mean total vertical optical depth τ V = L y (cid:104) κ R ρ (cid:105) V , (59)and the flux-weighted mean optical depth τ F = L y (cid:104) κ R ρF y (cid:105) V (cid:104) F y (cid:105) V , (60)where F y is the flux in the y direction and (cid:104)·(cid:105) V = L − x L − y (cid:82) L y (cid:82) L x · dxdy denotes volume avearge.Figure 11 compares the time evolution of the volume-weighted quantities in our IMC run with those in FLD, M1, andVET. The evolution of f E , V in IMC matches both qualitatively andquantitatively that in VET over the entire course of the run. Com-mon to all the simulations except the one carried out with the M1closure without radiative trapping, the mean Eddington ratio in-creases from its initial value of f E , V = 0 . to super-Eddingtonvalues soon after the simulation beginning. Then it immediatelydeclines toward f E , V (cid:46) . . After ∼ t ∗ , all methods becomesub-Eddington, with the M1 with radiative trapping and the FLDexhibiting the most significant decline. Then beyond ∼ t ∗ , allsimulations attain near-unity Eddington ratios. D14 pointed out thatthe time evolution of (cid:104) v y (cid:105) is in general sensitive to the value of f E , V , namely, (cid:104) v y (cid:105) increases when f E , V > and decreases other-wise. It is observed that IMC stays slightly super-Eddington at latetimes, similar to VET. The observed continuous acceleration of thegas with IMC suggests that gas dynamics can be very different be-tween simulations with similar volume-average Eddington rationsas long as the simulations are performed with different radiativetransfer methods.The middle panel of Figure 11 shows the evolution of thevolume-weighted mean total vertical optical depth τ V . Since thisquantity depends only on the gas state but not on the noisier radi-ation state, the IMC track is smooth. It is a global estimate of theoptical thickness of the gas layer and we expect its behavior to berelated to that of f E , V . The IMC track seems a flattened and down-scaled version of the others. This should be an artifact of the precisechoice of the opacity law. We cap the opacities κ R , P at their val-ues at T = 150 K, whereas the other authors allow the κ ∝ T scaling to extend at T >
K. Therefore, our choice of opac-ity underestimates the strength of radiation pressure compared tothe cited studies, but this discepancy does not appear to affect thehydrodynamic response of the gas.The bottom panel of Figure 11 shows the ratio τ F /τ V . Notethat τ F is the true effective optical depth felt by the radiation. There-fore, a small τ F /τ V implies a higher degree of flux-density anti-correlation. The evolution of this ratio is similar in all radiationtransfer methods. We applied the Implicit Monte Carlo radiative transfer method to astandard two-dimensional test problem modeling the radiation hy-drodynamics of a dusty atmosphere that is accelerated against grav-ity by an IR radiation field. The atmosphere is marginally capableof trapping the transiting radiation. We consider this idealized sim-ulation a necessary stepping stone toward characterizing the dy-namical impact of the radiation emitted by massive stars and activegalactic nuclei. We compare our IMC-derived results with those us-ing low-order closures of the radiative transfer hierarchy that havebeen published by other groups. Our particle-based approach en-ables independent validation of the hitherto tested methods.Sufficiently strong radiation fluxes universally render the at-mosphere turbulent, but its bulk kinematics differs between theVET and IMC methods on the one hand and the FLD and M1 meth-ods on the other. We find that the former continue to accelerate theatmosphere against gravity in the same setup in which the latter reg-ulate the atmosphere into a gravitationally-confined, quasi-steadystate. This exposes shortcomings of the local closures. Namely, incomplex geometries, the FLD seems to allow the radiation to moreeasily escape through optically thin channels. This can be under-stood in terms of a de facto artificial re-collimation of the radiationfield diffusing into narrow, optically-thin channels from their moreoptically thick channel walls. In the limit in which the radiationfreely streams in the channels, the flux in the channels becomesequal to what it would be for a radiation field in which the photonmomenta are aligned with the channel direction. Indeed, D14 arguethat in the optically thin regime, the FLD’s construction of the radi-ation flux is inaccurate in both its magnitude and direction, and hasthe tendency to reinforce the formation of such radiation-leakingchannels.Whether outflowing or gravitationally-confined, the turbulentatmosphere seems to reach a state approximately saturating the Ed-dington limit. The nonlinearity arising from the increase of dustopacity with temperature introduces the potential for bi-stability inthe global configuration. Subtle differences between numerical clo-sures can be sufficient to force the solution into degenerate, qual-itatively different configurations. Robust radiation-hydrodynamicmodeling seems to demand redundant treatment with distinct nu-merical methods including the IMC.Future work will of course turn to more realistic astrophysicalsystems. For example, the role of radiation trapping and pressure inmassive star forming regions remains a key open problem, both inthe context of the nearby (Krumholz & Matzner 2009; Krumholzet al. 2012, 2014; Coker et al. 2013; Lopez et al. 2014) and thedistant (Riechers & et al. 2013) universe. Radiative reprocessingby photoionization and dust requires a frequency-resolved treat-ment of the radiation field as well as a generalization the IMCmethod to nonthermal processes. The assumption of perfect gas-dust thermal coupling can be invalid and the respective temper-atures must be tracked separately. Numerical treatments may berequired to resolve dust sublimation fronts (Kuiper et al. 2010)and radiation pressure on metal lines (Tanaka & Nakamoto 2011;Kuiper & Yorke 2013). On the small scales of individual massive-star-forming cores, multifrequency radiative transfer may be ofessence for robust estimation of the final characteristic stellar massscale and the astronomically measurable accretion rate (Yorke &Sonnhalter 2002; Tan et al. 2014). Photoionization can set the fi-nal stellar masses through fragmentation-induced starvation (Pe-ters et al. 2010). The star formation phenomenon spans a huge dy-namic range that can be effectively treated with telescopic AMR c (cid:13) , 000–000 adiation pressure driving of a dusty atmosphere grids constructed to ensure that the local Jeans length is always ad-equately resolved. It will likely be necessary to invent new acceler-ation schemes for improving the IMC method’s efficiency in suchheterogeneous environments. One promising direction is the intro-duction of MCP splitting (see, e.g., Harries 2015, where MCP split-ting is applied in methods developed to simulate radiation transferin massive star forming systems). ACKNOWLEDGMENTS
We are grateful to the referee M. Krumholz for very helpful com-ments, to E. Abdikamalov for generously sharing details of his IMCradiative transfer implementation, to C. Ott for inspiring discus-sions, and to S. Davis and J. Rosdahl for consultation and sharingsimulation data with us. B. T.-H. T. is indebted to V. Bromm forencouragements throughout the course of this research. He also ac-knowledges generous support by The University of Hong Kong’sHui Pun Hing Endowed Scholarship for Postgraduate ResearchOverseas. The
FLASH code used in this work was developed in partby the DOE NNSA-ASC OASCR Flash Center at the University ofChicago. We acknowledge the Texas Advanced Computing Centerat The University of Texas at Austin for providing HPC resources,in part under XSEDE allocation TG-AST120024. This study wassupported by the NSF grants AST-1009928 and AST-1413501.
REFERENCES
Abdikamalov E., Burrows A., Ott C. D., L¨offler F., O’Connor E., DolenceJ. C., Schnetter E., 2012, ApJ, 755, 111Andrews B. H., Thompson T. A., 2011, ApJ, 727, 97Cheatham J. R., 2010, PhD thesis, The University of MichiganCoker C. T., Thompson T. A., Martini P., 2013, ApJ, 778, 79Commerc¸on B., Teyssier R., Audit E., Hennebelle P., Chabrier G., 2011,A&A, 529, A35Davis S. W., Jiang Y.-F., Stone J. M., Murray N., 2014, ApJ, 796, 107Davis S. W., Stone J. M., Jiang Y.-F., 2012, ApJS, 199, 9Densmore J. D., Urbatsch T. J., Evans T. M., Buksas M. W., 2007, Journalof Computational Physics, 222, 485Dubey A., Reid L. B., Fisher R., 2008, Physica Scripta Volume T, 132,014046Ensman L., 1994, ApJ, 424, 275Faucher-Gigu`ere C.-A., Quataert E., Hopkins P. F., 2013, MNRAS, 433,1970Fleck Jr. J. A., Canfield E. H., 1984, Journal of Computational Physics,54, 508Fleck Jr. J. A., Cummings J. D., 1971, J. Comp. Phys., 8, 313Fryxell B., Olson K., Ricker P., Timmes F. X., Zingale M., Lamb D. Q.,MacNeice P., Rosner R., Truran J. W., Tufo H., 2000, ApJS, 131, 273Geach J. E., Hickox R. C., Diamond-Stanic A. M., Krips M., RudnickG. H., Tremonti C. A., Sell P. H., Coil A. L., Moustakas J., 2014, Nature,516, 68Gentile N. A., 2001, Journal of Computational Physics, 172, 543Gonz´alez M., Audit E., Huynh P., 2007, A&A, 464, 429Harries T. J., 2011, MNRAS, 416, 1500Harries T. J., 2015, MNRAS, 448, 3156Hopkins P. F., Murray N., Quataert E., Thompson T. A., 2010, MNRAS,401, L19Hykes J. M., Densmore J. D., 2009, J. Quanti. Specstrosc. & Radiat.Transfer, 110, 1097Krumholz M. R., Bate M. R., Arce H. G., Dale J. E., Gutermuth R., KleinR. I., Li Z.-Y., Nakamura F., Zhang Q., 2014, Protostars and Planets VI,pp 243–266Krumholz M. R., Klein R. I., McKee C. F., 2012, ApJ, 754, 71Krumholz M. R., Klein R. I., McKee C. F., Bolstad J., 2007, ApJ, 667, 626 Krumholz M. R., Matzner C. D., 2009, ApJ, 703, 1352Krumholz M. R., Thompson T. A., 2012, ApJ, 760, 155Krumholz M. R., Thompson T. A., 2013, MNRAS, 434, 2329Kuiper R., Klahr H., Beuther H., Henning T., 2010, ApJ, 722, 1556Kuiper R., Yorke H. W., 2013, ApJ, 763, 104Kuiper R., Yorke H. W., Turner N. J., 2015, ApJ, 800, 86Lopez L. A., Krumholz M. R., Bolatto A. D., Prochaska J. X., Ramirez-Ruiz E., Castro D., 2014, ApJ, 795, 121Mihalas D., Mihalas B. W., 1984, Foundations of radiation hydrodynam-ics. Oxford University PressMurray N., M´enard B., Thompson T. A., 2011, ApJ, 735, 66Murray N., Quataert E., Thompson T. A., 2005, ApJ, 618, 569Murray N., Quataert E., Thompson T. A., 2010, ApJ, 709, 191Peters T., Banerjee R., Klessen R. S., Mac Low M.-M., Galv´an-MadridR., Keto E. R., 2010, ApJ, 711, 1017Riechers D. A., et al. 2013, Nature, 496, 329Rosdahl J., Teyssier R., 2015, MNRAS, 449, 4380Roth N., Kasen D., 2015, ApJS, 217, 9Semenov D., Henning T., Helling C., Ilgner M., Sedlmayr E., 2003, A&A,410, 611Su B., Olson G. L., 1996, J. Quanti. Specstrosc. & Radiat. Transfer, 56,337Tan J. C., Beltr´an M. T., Caselli P., Fontani F., Fuente A., Krumholz M. R.,McKee C. F., Stolte A., 2014, Protostars and Planets VI, pp 149–172Tanaka K. E. I., Nakamoto T., 2011, ApJ, 739, L50Thompson T. A., Fabian A. C., Quataert E., Murray N., 2015, MNRAS,449, 147Thompson T. A., Krumholz M. R., 2014, ArXiv e-printsThompson T. A., Quataert E., Murray N., 2005, ApJ, 630, 167Turner N. J., Stone J. M., 2001, ApJS, 135, 95Wise J. H., Abel T., Turk M. J., Norman M. L., Smith B. D., 2012, MN-RAS, 427, 311Wollaber A. B., 2008, PhD thesis, The University of MichiganWollaeger R. T., van Rossum D. R., Graziani C., Couch S. M., Jordan IVG. C., Lamb D. Q., Moses G. A., 2013, ApJS, 209, 36Yorke H. W., Sonnhalter C., 2002, ApJ, 569, 846Zel’dovich Y. B., Raizer Y. P., 1967, Physics of shock waves and high-temperature hydrodynamic phenomena. Academic PressZhang W., Howell L., Almgren A., Burrows A., Bell J., 2011, ApJS, 196,20c (cid:13)000