Radiation Feedback in ULIRGS: Are Photons Movers and Shakers?

# Radiation Feedback in ULIRGS: Are Photons Movers and Shakers?

Shane W. Davis11affiliation: Canadian Institute for Theoretical Astrophysics. Toronto, ON M5S3H4, Canada , Yan-Fei Jiang(å§çé£)22affiliation: Harvard-Smithsonian Center for Astrophysics, 60 Garden Street, Cambridge, MA 02138, USA 33affiliation: Einstein Fellow , James M. Stone44affiliation: Department of Astrophysical Sciences, Princeton University, Princeton, NJ 08544, USA , and Norman Murray11affiliation: Canadian Institute for Theoretical Astrophysics. Toronto, ON M5S3H4, Canada
###### Abstract

galaxies: ISM – ISM: jets and outflows – hydrodynamics – radiative transfer – methods: numerical
{CJK*}

UTF8gbsn

## 1. Introduction

Observations of star-formation in the Milky way and other galaxies provide consistent evidence that some feedback mechanism or mechanisms hamper the collapse of interstellar gas to form stars. For example, The Kenicutt-Schmidt law (Kennicutt, 1998) implies that, on average, only a few percent of the available gas actually collapses to form stars per dynamical time. Observations of molecular gas in rapidly star-forming ultraluminous infrared galaxies (ULIRGs) indicate that turbulent velocities of up to are present (e.g. Downes & Solomon, 1998). And most dramatically, galaxy scale outflows of cold, neutral gas are inferred in galaxies ranging from nearby dwarf starbursts to ULIRGs and rapidly star-forming galaxies at high redshift (e.g. Heckman et al., 1990; Pettini et al., 2001; Schwartz & Martin, 2004; Rupke et al., 2005).

Although a number of promising mechanisms have been proposed to explain these observations, we restrict our attention to the potential role of radiation pressure on dust grains in driving turbulence, hampering gravitational collapse, and launching outflows in such environments (Scoville et al., 2001; Murray et al., 2005; Thompson et al., 2005). This possibility has been explored extensively in recent years, with a number of studies considering how (in)effective radiation driving may be in various environments (e.g. Krumholz & Matzner, 2009; Andrews & Thompson, 2011; Hopkins et al., 2011, 2012; Wise et al., 2012; Krumholz & Thompson, 2012, 2013; Socrates & Sironi, 2013)

Here we focus specifically on the implications of the Rayleigh-Taylor instability (hereafter RTI; see e.g. Chandrasekhar, 1961). In recent work, Krumholz & Thompson (2012, hereafter KT12) and Krumholz & Thompson (2013) have argued that the RTI may play a significant role in limiting the exchange of momentum between radiation and dusty gas, ultimately reducing the role of radiation feedback in star-forming environments. A general, analytical calculation of the linear growth of the radiative RTI does not exist (see e.g. Mathews & Blumenthal, 1977; Krolik, 1977; Jacquet & Krumholz, 2011; Jiang et al., 2013), and non-linear evolution can generally only be explored via numerical simulations (KT12; Jiang et al., 2013).

In this paper, we attempt to replicate the results of KT12, who numerically solve the equations of radiation hydrodynamics using an implementation of the flux-limited diffusion algorithm in the ORION code (Krumholz et al., 2007). We utilize both a variable Eddington tensor method (Davis et al., 2012; Jiang et al., 2012) and our own version of the flux-limited diffusion method (Jiang et al., in preparation), each of which are implemented as part of the Athena astrophysical fluid dynamics code (Stone et al., 2008). Although we reproduce some aspects of the KT12 results, we find significant discrepancies that arise from innacuracies of the flux-limited diffusion treatment.

The plan of this work is as follows: We review the equations solved and summarize our numerical methods in section 2. We describe the setup of our numerical simulations, including formulations for the opacities, boundary conditions, and initial conditions in section 3. We summarize our key results from our numerical simulation in section 4 and discuss their implications in section 5. We provide our conclusions in section 6.

## 2. Equations and Numerical Methods

In this work we solve the equations of radiation hydrodynamics using the variable Eddington tensor (VET) implementation (Sekora & Stone, 2010; Davis et al., 2012; Jiang et al., 2012). The primary set of equations to be solved are the coupled systems of hydrodynamics and the radiation moment equations. The standard fluid equations are solved using a finite volume formulation as discussed in Stone et al. (2008). We solve the radiation system in the mixed frame approach (see e.g. Mihalas & Mihalas, 1984), where the radiation moments are Eulerian frame quantities, while the emissivities and opacities correspond to the comoving frame quantities (c.f. Mihalas & Klein, 1982). We follow the derivation of Lowrie et al. (1999), which retains all terms of order and some terms of order . The equations correspond to mass conservation

 (1)

momentum conservation

 ∂(ρv)∂t+∇⋅(ρvv+P)=ρg−Sr(P), (2)

and energy conservation

 ∂E∂t+∇⋅(Ev+P⋅v)=ρg⋅v−cSr(E). (3)

In the above, is the gas density, is the fluid velocity, is the gravitational acceleration, and is the total (fluid) energy density . The pressure tensor is defined as , where is the gas pressure and is the identity matrix.

The quantities and are the radiation momentum and energy source terms, respectively. They are given by

 Sr(P) = −σF[Fr−(vEr+v⋅Pr)]/c (4) +v(σParT4−σEEr)/c,

and

 Sr(E) = (σParT4−σEEr) (5) +σFvc2⋅[Fr−(vEr+v⋅Pr)].

Here is the radiation energy density, the radiation flux, is the radiation pressure tensor, is the radiation constant, is the gas temperature, is energy mean opacity, is the Planck mean opacity, and is the flux mean opacity. In this work we neglect any scattering contribution to the opacity. The radiation source terms on the right hand side of equations (2) and (3) are included via a modified Godunov method, as described in Jiang et al. (2012).

In order to compute these quantities, we solve the radiation moment equations representing conservation of radiation energy

 ∂Er∂t+∇⋅Fr=cSr(E), (6)

 1c2∂Fr∂t+∇⋅Pr=Sr(P). (7)

This radiation system is solved using a backward-Euler formulation (Jiang et al., 2012).

### 2.1. The Short-Characteristics Method

The above equations are incomplete because there is no evolution equation for the radiation pressure tensor. We address this by computing the eponymous variable Eddington tensor . We approximate by solving the time-independent equation of radiation transfer of the form

 ^n⋅∇I=σF(arc4πT4−I), (8)

where is the specific intensity for an angle defined by the unit vector . Hence, we have adopted a grey opacity treatment in which the characteristic opacity is assumed to correspond to the flux mean opacity. For the computation of the , we make no distinction between the comoving and Eulerian frames, since equation (8) drops all order terms or higher, including a term that would otherwise be present. The neglect of the term is typically a good approximation whenever but the neglect of velocity-dependent terms generally requires if is a characteristic optical depth for the system. In the results presented here and , so the neglect of these terms is a good approximation.

On each timestep, this equation is solved using the method of short characteristics (Kunasz & Auer, 1988), as described in Davis et al. (2012). As the radiation transfer calculation proceeds, the mean intensity , the first moment , and the second moment are tabulated in each grid zone. These are then used to compute for the subsequent timestep. We integrate the radiation field along 84 rays, distributed nearly uniformly over the half unit sphere, taking advantage of the reflection symmetry of two dimensional Cartesian domains. Further discussion of the angular grid and impact of angular resolution is provided in the appendix.

### 2.2. The Flux-limited Diffusion Method

We have also implemented a flux-limited diffusion (FLD) algorithm in Athena. Here, we only briefly describe the equations and implementation. A more detailed description of the algorithm and its implementation will be provided in a forthcoming work (Jiang et al., in preparation).

In FLD, one drops the radiation momentum equation from the evolution equations, replacing it with a diffusion equation for the flux

 Fr=−cλσF∇Er, (9)

where is the flux-limiter. The radiation energy equation is then solved by substituting this relation for in equations (4)-(6). Following Levermore & Pomraning (1981), we assume a flux limiter of the form

 λ = 1R(cothR−1R) R = |∇Er|+βσFEr. (10)

To specify we use the Eddington tensor

 f=12(1−f)I+12(3f−1)^f^f (11)

with and .

Following Shestakov & Offner (2008) we have added a parameter to the definition of , which acts as an effective floor on in optically thin regions. This addition was necessary to robustly obtain convergence in our backward-Euler scheme, in which one needs to solve a large matrix equation. We find poor convergence for our multigrid solver unless . These equations are implemented in a manner analogous to the VET method described above. The source terms are coupled to the hydro equations using our modified Godunov method and the radiation energy equation is solved using a backward-Euler algorithm. Finally, we note that is the Eulerian frame flux in the above equations. Strictly speaking, it is the comoving frame flux for which the diffusion approximation applies in high optical depth media. However, the low vales of in the simulations lead to very small differences between the Eulerian and comoving frame fluxes.

Due to the ad hoc aspects of this approximation, it is not as reliable as our VET method, which evolves the radiation momentum equation and uses a solution of the radiation transfer equation to estimate . Modulo the velocity dependent terms that we neglect in computing (and only in computing ), the accuracy of our approximation can always be improved by increasing the angular resolution. In contrast, FLD will always be inaccurate at some level in regions of the flow where the characteristic optical depths are near unity or less.111Although FLD is often claimed to be accurate in the optically thin limit, it assumes , which only applies to a radiation field with perfectly parallel rays. This is typically only obtained at large distances from a point source emitter and is not general. Well known deficiencies of this method include its inability to cast shadows (see e.g. Hayes & Norman, 2003) because the radiation field diffuses around opaque barriers, even in optically thin environments. We have implemented FLD in Athena only to facilitate comparison of our VET calculation with previous FLD-based results, since it aids in determining which aspects of the calculation depend on the transfer method.

## 3. Simulation Setup

Our simulation setup follows the computations performed in KT12. The goal is to simulate the evolution of the interstellar gas in ULIRGs, although the results may generalize to other rapidly star forming environments. We assume that the gas and dust are strongly coupled, both dynamically and thermally. This assumption should be quite reasonable for ULIRGs, as discussed in appendix A of Krumholz & Thompson (2013). Hence, the gas and dust are assumed to share a common temperature (hereafter simply referenced as the gas temperature ), and any momentum exchanged between the dust and radiation field is rapidly shared via collisions with the gas.

Following KT12, the Planck and Rosseland mean opacities are given by

 κP = 0.1(T10K)2cm2g−1 κR = 0.0316(T10K)2cm2g−1. (12)

This scaling approximately holds for dust in thermal equilibrium with a blackbody radiation field for (Semenov et al., 2003). This implicitly assumes the characteristic temperature of the radiation and are equal. Since we evolve the radiation energy density, we can define a characteristic radiation temperature and check a posteriori if . This is generally a good approximation although modest differences are present in some regions, which motivated us to perform alternative simulations with . The two sets of simulations agreed very closely with each other so we report only the here. These assumptions imply that the radiation field is characterized by blackbodies with K. Since the main source of radiation is massive stars emitting primarily in the UV, this analysis assumes the direct stellar radiation has already been reprocessed by the dust and converted to infrared radiation.

The simulations are performed on a two-dimensional, Cartesian grid. For simplicity, the gravitational acceleration is assumed to be constant and the radiation field is sourced by a constant flux , incident at the lower boundary. This flux defines a characteristic temperature , from which we can define a characteristic sound speed , scale height and sound crossing time . Following KT12, we assume and choose K, corresponding to and . The remaining free parameters are and the surface mass density , which can be specified in terms of the dimensionless Eddington ratio222Note that this constitutes only a local assessment of the force balance. Since varies with height and radius in real systems, a local balance between radiation and gravitational forces simply means the disk is radiation pressure supported. The global Eddington ratio of the system may be much less than unity.

 fE,∗=κR,∗F∗gc (13)

and optical depth

 τ∗=κR,∗Σ. (14)

### 3.1. Initial Conditions

As in KT12, we assume an initially static and isothermal atmosphere with . The density is initialized as an exponential profile with scale height , which would correspond to hydrostatic equilibrium if radiation pressure support was negligible. For both the initial condition and subsequent evolution we impose a density floor of , where is the maximum initial density at the lower boundary. We assume everywhere, which gives a constant . Hence, for even moderate values of and , this initial condition is neither in thermodynamic nor hydrostatic equilibrium.

On top of this initial density profile we include a perturbation of the general form

 δρρ=0.25(1±χ)sin(2πx/λ), (15)

where , and is the domain width. Here, is identically zero for sinusoidal perturbations or can be a random number uniformly distributed between -0.25 and 0.25. If , the perturbation is identical to KT12, but we found it useful to consider simulations with random perturbations on top of the sinusoidal variation.

Simulations parameters for our primary simulations are summarized in Table 1.

### 3.2. Boundary Conditions

For all simulations, both the radiation and hydrodynamic quantities obey periodic boundary conditions in the horizontal (-coordinate) direction. For the hydrodynamic quantities, we impose reflecting boundary conditions at the lower vertical (-coordinate) boundary and outflow boundary conditions on the upper vertical boundary. These boundary conditions match those of KT12, modulo implementation differences between Athena and ORION. The exception is that we simply continue and into the ghost zones at our upper boundary whereas KT12 fix and there.

For the lower boundary condition of the VET runs, we apply reflecting boundary conditions on and require

 Fry=F∗+vyEr+(v⋅Pr)y. (16)

Hence, the comoving flux in the ghost zones is equal to . We specify by integrating the time-independent vertical momentum equation, assuming the Eddington tensor at the lower boundary of the computational domain applies throughout the ghost zones. At the upper boundary, both components of are continued in to the ghost zones and we specify using

 Er=JHyFry, (17)

where and are the first moment and mean intensity (respectively) returned by the short-characteristics module. In other words, the ratios of the vertical component of the first moment to the “zeroth” moment for the two calculations are forced to match identically in the ghost zones.

For the lower boundary condition of the FLD computations, we compute the flux limiter using equation (10) and solve for using equation (9), assuming and no horizontal gradient in energy density. At the upper boundary we assume a fixed in the ghost zones for consistency with KT12.

For the lower boundary in the short-characteristics module, we integrate the outgoing (downward) radiation intensity over angle . We then assume isotropic incoming (upward) radiation. The intensity of the incoming radiation is normalized so that the corresponding angle integrated upward intensity obeys . At the upper boundary we assume that the intensity of incoming radiation is zero. Periodic boundary conditions in the horizontal direction are imposed by iterating the short-characteristics solution to convergence, as described in Davis et al. (2012).

## 4. Results

KT12 computed the equilibrium profiles for one-dimensional atmospheres obeying the assumptions of section 3, assuming equation (9) holds. They find that for a given value of , there is a maximum value for () for which their iterative method converges. Above this value, no equilibrium solution exists because the radiation and gravitational forces cannot be forced to match exactly due to the temperature dependence of the opacities. Their results for are summarized in their figure 2. Although their results are derived using the FLD approximation, we expect that their boundary curve should approximately demarcate the locus of stable equilibria solutions for our VET calculations, as well.

Here we consider two sets of and , as summarized in Table 1. The first (, ) falls in the range where an equilibrium exists, while the second (, ) falls in the unstable regime. We choose these parameters to match simulations performed by KT12.

### 4.1. Stable Case

We first consider the T10F0.02VET run, which we run using our VET module. The parameters for this run place it in the stable regime, according to the one-dimensional equilibrium solutions and two-dimensional simulation results of KT12. We initialized the simulation with a sinusoidal perturbation only (i.e no random fluctuations), as described in section 3.1 and ran it for .

Figure 1 shows four snapshots of the density from this run. At the beginning of the run, the large optical depth and constant incident flux require to increase. Since the radiation and dusty gas are well-coupled, the increase in drives a corresponding rise in . This leads to an increase in the opacity and a corresponding rise in the optical depth and radiation force. As a result, the atmosphere briefly expands upward, relaxes, and begins to oscillate around a quasi-equilibrium state.

This transient evolution and subsequent oscillatory behavior are most clearly seen in the mass-weighted mean velocity

 ⟨v⟩=1M∫ρvdV, (18)

and mass-weighted velocity dispersion

 σ2i=1M∫ρ(vi−⟨vi⟩)2dV, (19)

where is the total mass in the atmosphere. The top panel of figure 2 shows and , along with the total velocity dispersion () and the bottom panel shows .

The initial transient acceleration lead to growth in both the vertical and horizontal velocity dispersion, although the vertical component initially grows faster. After , both and have already reached their maximum for the simulation. After another , both quantities settle down into oscillations, with varying about zero. The horizontal velocity dispersion is also oscillatory, albeit on a longer timescale with a period greater than . Comparison with the top panel of figure 7 in KT12 shows close agreement between the results of the two different codes.

At the end of the run, the simulation is still oscillating, albeit with a slow decay in . Reductions in the kinetic energy due to “numerical viscosity” and diffusive damping of compressive motions by the radiation field will eventually damp the oscillatory behavior and the simulation will presumably settle into a steady state. However, this would seem to require a significantly longer run time and we are not interested in the detailed properties of equilibrium. We stop the run at , which is long enough to confirm that we reproduce the KT12 results with our VET formalism for this stable regime.

### 4.2. Unstable Case

We now focus on the unstable regime, considering simulations with and . This run has the lowest ratio of considered in KT12 for the unstable regime. Furthermore, it had the weakest transient acceleration at the beginning of the simulation, reached the lowest maximum vertical extent, and had the lowest maximum velocity dispersion. In this sense, the run represented the worst case for driving turbulence and outflows among the runs considered by KT12.

In our work, we consider two runs: one using our FLD module (T3F0.5FLD) and one using the VET module (T3F0.5VET), in order to assess the impact of the radiation transfer method. In principle, we could simply compare our VET method directly against the ORION FLD results, but we perform Athena FLD runs to control for differences in the hydrodynamics algorithms between Athena and ORION. The simulations were initialized with the goal of nearly reproducing the KT12 setup. One notable exception is that both runs have sinusoidal and random perturbations, as described in equation (15). We made this choice after running a simulation with the VET module that only included sinusoidal perturbations. In that case, the development of non-linear structure (due to the RTI, as discussed below) was slower than seen in KT12 and most of the mass in the shell was able to reach the top of the domain before this structure produced a significant feedback on the radiative acceleration of the shell. Since the sinusoidal perturbation is a artificial construction and because we were interested in the evolution of an atmosphere with significant non-linear structure, we seeded the RTI with additional random perturbations. This allowed for faster development of the RTI in the VET runs. We also initialize the FLD run with random perturbations to facilitate direct comparison between our FLD and VET runs.

We first consider the results from T3F0.5FLD, which we run for . This is run in a taller box than the T10F0.02VET simulation, with the initial condition and boundary conditions described in section 3. Note that we began the simulations with in equation (10) to place a floor on and avoid convergence problems in our backward Euler scheme. This was sufficient for the first , but we had to increase at to ensure convergence thereafter.

As in the T10F0.02VET run, there is an initial increase in at the base of the domain due to the finite optical depth, constant incoming flux, and close thermal coupling between gas, dust, and radiation. This increase in drives an increase in , resulting in an increase of the radiation force above the gravitational force. Hence, the lower part of the atmosphere quickly becomes super-Eddington, and the vast majority of the mass is driven upward in a thin shell.

Figure 3 shows snapshots of and from four different times. By , the shell has already begun to break up, with some material running behind the shell (or even falling back) in a number of plumes. This behavior is consistent with expectations that the shell should be subject to the RTI when accelerated against gravity by the radiation forces under these conditions. (See section 5.3 for further discussion.) Our results are also qualitatively consistent with those of KT12, who also attribute the non-linear structure to RTI. The growth of non-linear structure in our run is somewhat faster and the structures are less uniform, consistent with our additional random perturbations, which seed the growth of smaller scale structure.

Subsequent evolution is also consistent with KT12. The RTI plumes grow into high density filaments interspersed with lower density channels. Since is largest in the low density channels, there is an anti-correlation of and that reduces the momentum exchange between the radiation field and gas. The radiation forces become sub-Eddington in the optically-thick, high-density filaments and the majority of the matter falls back by . This behavior is reinforced by a drop in the volume averaged optical depth, which leads to a drop in (and therefore ) at the base of the atmosphere, further reducing the radiation force. By , almost all of the mass is contained in the bottom of the domain and remains there for the duration of the run, which we stop at .

As expected from the lack of a one-dimensional equilibrium, the matter at the base of the domain does not settle into a hydrostatic equilibrium, but instead remains in a turbulent state, with a moderate velocity dispersion and very little average vertical motion. Again, this quasi-steady state turbulent flow is qualitatively and quantitatively consistent with the results in KT12.

We now turn our attention to the VET run labeled T3F0.5VET. This simulation began in a domain with the same dimensions, resolution, and initial condition as in the T3F0.5FLD run. Figures 4 and 5 show four snapshots of and (respectively) from T3F0.5VET. Initially, the dynamics of the simulation are qualitatively consistent with the T3F0.5FLD run. Again, and rises rapidly at the base of the domain and the bulk of the mass is launched upward as a shell. The RTI seems to grow slightly slower than in the FLD run, but is still fairly rapid and significant growth of non-linear structure is apparent by .

The difference between the FLD and VET runs become more apparent in the subsequent evolution. As in the T3F0.5FLD run, there is a flux-density anti-correlation that allows high density filaments to become sub-Eddington even as the volume averaged structure remains super-Eddington. Some of these high-density filaments do fall back to the base of the domain, where they are reflected by the boundaries. Others are disrupted and reaccelerated by the radiation field as they disperse. The overall evolution leads to a gradual filling of the volume with gas, although most of the mass remains in several dense filaments.

At , the dense gas reaches the upper boundary of the initial domain and the assumptions made for the radiation boundary condition (near vacuum with no incoming radiation) cease to be valid. This leads to an abrupt, unphysical increase in at the upper boundary. Therefore, we restart the simulation at in a domain with the and (i.e. we double the height at fixed resolution). All variables are copied into the bottom half of the new domain. Velocities in the upper half of the domain are set to zero. All other material variables, , and are averaged on the upper most grid zones and the upper half of the domain is uniformly initialized with these values. The Eddington tensor is recomputed with the short characteristics module. Subsequent evolution matches the original evolution until when the breakdown of the vertical boundary condition alters the evolution in the smaller domain. We continue our integration in the larger domain until , when the dense gas again reaches the upper boundary of the larger domain. Again, most of the mass is concentrated in a few of the densest filaments, but the remainder is nearly volume filling with , except for transient periods when most of the volume near the bottom of the domain can have .

Figure 5 shows that there is relatively little horizontal variation in , consistent with T3F0.5FLD and KT12. In contrast to T3F0.5FLD, the vertical distribution spreads out vertically with time as matter fills the domain, but temperature at the base of the domain remains relatively constant, with . Modest, but short-lived increases in are seen, and correspond to fall back and reflection of dense filaments at the lower boundary (see e.g. the second panel at ). These values are modestly higher than the mean values of at the base in T3F0.5FLD, although there is greater fluctuation in in the T3F0.5FLD run.

Figure 6 shows a comparison of volume averaged quantities and their ratios from the T3F0.5FLD and T3F0.5VET simulations. The top panel compares the characteristic Eddington ratio

 fE,V=⟨κRρFry⟩cg⟨ρ⟩ (20)

where denote volume averages (e.g. ). Both runs have an initial super-Eddington period that is followed by a gradual decline to sub-Eddington values, although this decline is more rapid and deeper in the FLD case. Both simulations show a return to an Eddington ratio near unity, although T3F0.5FLD tends to have fluctuating near unity, while T3F0.5VET tends to remain moderately super-Eddington for most of the run. Although the differences are rather modest, the upshot is that T3F0.5VET receives a nearly continual, modest upward acceleration while T3F0.5FLD settles into a quasi equilibrium state. This evolution is apparent in the evolution of in figure 7. The trend in suggests that T3F0.5VET might settle into a similar equilibrium, but with a much larger scale height and higher velocity dispersion.

In order to explore the evolution and impact of optical depth variations, it is useful to define two characteristic averages for the optical depth at the base of the domain. The first is the volume-weighted mean optical depth

 τV=Ly⟨κRρ⟩, (21)

and the second is the flux-weighted mean optical depth

 τF=Ly⟨κRρFry⟩⟨Fry⟩. (22)

The middle panel of figure 6 shows that correlates closely with as both evolve over the course of the simulations. Large values of correspond to large values of , because and is roughly proportional to throughout much of the domain. Note, however, that it is not quite true that determines the Eddington ratio as can be larger in T3F0.5VET, even when is greater for T3F0.5FLD. This occurs primarily because the nature of anti-correlation between and differs in the two runs.

The bottom panel of figure 6 shows the ratio . Note that multiplies to give the momentum per unit area transfered from the radiation to the gas while would be the characteristic value for the total infrared optical depth in a uniform medium. Hence, this ratio gives an estimate of how much the flux – density anticorrelation reduces the momentum coupling between radiation and gas. Our , where is trapping factor discussed in KT12. Our estimate of agrees with KT12’s for , .

Overall, figure 6 leaves the impression that the FLD and VET simulations yield very similar results, particularly at late times, in contrast with the impression provided by comparison of figure 3 to figures 4 and 5. The key point is that even modest variations can lead to rather large differences in the outcome when systems are near an Eddington ratio of unity.

The differences are somewhat more apparent in figure 7, which shows the evolution of and . For both simulations, the evolution of is largely determined by the value of in figure 6. When , increases and when , decreases. In both simulations, there is an early period of transient growth, followed by a decline after the RTI sets in. For T3F0.5FLD, becomes negative as most of the mass falls back, but eventually settles into a quasi-steady state, with simply fluctuating near zero. For T3F0.5VET, grows for the majority of the run, but appears to be flattening out at late times when .

For both runs, the and components of show a modest initial increase, followed by a faster rises in as the RTI sets in and the shells start to break apart. The subsequent evolution for T3F0.5FLD involves a weak decline in compensated by a slow increase in . These trends roughly cancel so that fluctuates about for the majority of the simulation, consistent with the results of KT12. The T3F0.5VET evolution displays a more continuous rise in , with only a brief drop after , when increases quickly during a period of moderately high . We find a slow, continuous rise in , similar to the T3F0.5FLD. When we need to stop the run shortly after , is still increasing and already a factor of larger than in T3F0.5FLD.

## 5. Discussion

### 5.1. The Dependence on the Radiation Transfer Method

In this work, we have revisited the calculations of KT12, but have been unable to reproduce all of their results with our more accurate VET radiation transfer algorithm. In contrast, our own implementation of the FLD algorithm seems to agree well with their calculations, indicating that discrepancies arise from the use of the FLD approximation, rather than implementation differences between the Athena and ORION codes. Although we reproduce several aspects of KT12’s results (see section 4), in this section we focus primarily on the discrepancies, and describe the deficiencies of the FLD algorithm that we believe are primarily responsible.

First, we note that our VET calculations do closely reproduce the KT12 results in the stable regime. Therefore, qualitative differences in the evolution arise from the development of non-linear structure, driven by the RTI. Although the development of the RTI is slightly slower in the VET case, some general aspects of the evolution are similar to the FLD simulations: a thin shell develops, becomes unstable to RTI, breaks up into high density filaments interspersed with low-density channels, and eventually settles to state with at late times.

In spite of these qualitative similarities, the differences in transfer methods have a fairly striking impact on the evolution of the velocity and spatial distributions of the gas. For the majority of the VET simulation, the radiation force exceeds gravity and the net vertical velocity is always positive. In contrast, the FLD results only see a transient acceleration phase followed by fallback of most of the gas, and settle into quasi-steady state turbulence with a scale height and velocity dispersion which are much too low to explain observed systems.

We would like to understand what aspects of the FLD method lead to these discrepancies in evolution. We note that since both simulations reach Eddington ratios near unity, only modest differences are needed to produce divergent outcomes. Figure 8 shows a snapshot of for the T3F0.5VET run at . We focus on a small fraction of the domain where the shell is being disrupted. The white vectors denote the values, scaled by magnitude, from the VET calculation. For comparison, we also plot green arrows showing what would be using equation (9) with the values in the VET run.

In the densest regions, our VET calculations are nearly in the diffusion limit where the FLD method is reliable, and the two sets of vectors are in approximate agreement. In lower density regions, where the diffusion approximation is poor, the magnitude of the FLD flux is much larger and the directions may be significantly different, even anti-parallel. For example, below and adjacent to the large plume in the center of the figure, the FLD fluxes point downward or horizontal while the VET fluxes point nearly upward almost everywhere in the domain. Since the radiation forces are proportional to , the FLD fluxes would not be opposing the downward motion of the plume to the same extent as the VET flux and may even be reinforcing it. In contrast, the underdense regions near and with have implied FLD fluxes that significantly exceed the VET fluxes. This means that (relative to VET) the FLD radiation forces would be more effective at reinforcing the development of the low density channels that are forming. Overall, the tendency is for FLD to reinforce the development of non-linear structure to a greater extent than the VET algorithm, which is consistent with the faster non-linear development of the RTI in the FLD runs.

One objection to the comparison in figure 8 is that the values of in the FLD run will not generally be the same as those in the VET run. Figure 9 avoids this issue by comparing and in the T3F0.5VET run to the same quantities in the T3F0.5FLD run. In both cases, we focus on a subsets of the simulation domains. The upper boundaries of these subsets are chosen to so that less than an optical depth of intervening matter lies between the top of the subset and the top boundary of the simulation domain, where vacuum boundary conditions are imposed.

Since the mean-free-path of photons is proportional to , will generally be larger when is smaller and vice-versa. Such an anti-correlation is apparent in both runs and is particularly clear when we scale with . Such anti-correlations are expected to be particularly strong if the gradient in is relatively uniform and the diffusion limit applies (equation [9] with ). However, the optical depths across most of the filaments seen in figure 9 are of order unity or less, so the diffusion limit is not valid and the radiation flux is expected to be less sensitive to local values of . Instead it is determined by the geometry of sources of strong emission and the integrated optical depth along different lines-of-site to these sources. The local value of the radiation field is determined by non-local properties of the flow. The VET algorithm accurately captures the variation of in this regime, because the Eddington tensor follows from the calculation of the angle dependent radiation transfer equation. In contrast, the FLD radiation field is highly constrained by the ad hoc assumptions that underlie the approximation. There is a single preferred direction (parallel to ) and the ratio is determined by (equation [10]). Since and are relatively uniform, the sharp variation of dominates the variation of and, therefore, . This stronge dependence of on leads to an unphysically high level of anti-correlation between and in low-to-moderate optical depth regions.

The upshot is that FLD tends to overestimate the contrast in between the high density filaments and the low density channels, where most of the radiation escapes. In effect, the FLD radiation field is more effective at “punching holes” through the distribution and then reinforces the resulting channels to a greater degree than in the more accurate VET calculations. This allows more of the radiation to escape through the low density channels, leading to a lower average radiation force for the majority of the gas.

Finally, we note that the FLD and VET methods agree well in regions where structures are moderately optically thick () and the diffusion approximation applies. This suggests that the FLD results may still be relevant for systems where the photospheric Eddington ratio is low. In these systems (i.e. systems with low and large ), the radiation force can only exceed the gravitational force at high optical depth, where , and therefore , is sufficiently large.

### 5.2. Implications for Observed Systems and Subgrid Models of Radiation Feedback

Since our VET solutions differ in important respects from earlier results using FLD, it is important to reexamine the implications of RTI for radiation feedback. Although the current simulation setup is conducive to exploring the role of RTI, the choice to start with a hydrostatic initial condition limits the relevance to observations of real star forming galaxies. This is because the hydrostatic length scale in the problem is (for and ) and . Hence the domain is nearly two orders of magnitude too small to contain a realistic ULIRG disk. Other simplifications, such as the assumption of a constant , approximate grey opacity, and the constant incident infrared flux may also have a strong effect on the evolution.

Nevertheless, we can reach some tentative conclusions. One quantity of interest is the mass-weighted velocity dispersion . Since is growing throughout most of the simulation and at the simulations’s end, we view our maximum as a lower bound rather than a characteristic estimate. We find . This is about an order-of-magnitude lower than the values typical inferred (e.g. Downes & Solomon, 1998), although this is unsurprising given the small size of the simulation domain. It is difficult to assess what implications the trend towards will have for the subsequent evolution of and . The system is very sensitive to the precise value of when is constant. If (on average) continues to be slightly greater than one, then acceleration will continue and we would expect the scale height and velocities to grow with time indefinitely. In a real system, the vertical variation of would presumably play a significant role in determining an equilibrium scale height and velocity dispersion.

A second question of interest is what is the correct efficiency (or “trapping”) factor to use in assessing the momentum imparted by radiation pressure. In feedback models for optically thick environments, the rate of momentum injection has been parameterized as (see e.g Hopkins et al., 2011, and references therein)

 dpdt=(1+ητIR)Lc, (23)

where is an estimate of the momentum transferred to the interstellar gas, is an estimate of the integrated optical depth through the dusty gas at infrared wavelengths, is luminosity of a star (or star cluster), and is an ad hoc reduction factor () included to account for inhomogeneities in the surrounding gas. Since corresponds to a volume-weighted average total optical depth, one can estimate , while represents the effective optical depth for momentum exchange. Therefore, we estimate the in equation (23) as the ratio , where we have time averaged the VET run for .

A value of is in the same range that was considered by Hopkins et al. (2011) in their calculations. As noted in section 4.1, our value of agrees with KT12’s (since ). However, in their discussion KT12 argue that Hopkins et al. (2011) and others likely overestimate . The origin of this apparent discrepancy lies partly in how one estimates and partly in the dependence of and on and . First, KT12 estimate , where is the temperature at the base of the simulation domain. For our simulation, this quantity is a factor of larger than for T3F0.5VET.

KT12 also perform a number of calculations with , which are associated with lower values of . Although we consider only a single parameter set with , we can derive approximate scalings for with and by assuming that all simulations will approach a characteristic horizontally average vertical profile. At late times our T3F0.5VET run can be fit approximately by333Note that equation (24) implies , which is roughly the ratio returned by our short-characteristics calculation at the top of the domain. For an approximately isotropically emitting, semi-infinite atmosphere, which is essentially what our periodic horizontal boundary conditions provide, this agrees well with standard grey treatments (Mihalas, 1978). This is greater than the assumed in FLD calculations, a limit that is only obtained at large distances from a point source.

 ¯¯¯¯T4=T4∗3(23+¯¯¯τF), (24)

where quantities with an overbar correspond to -dependent, horizontally average quantities (e.g. ). The exceptions are

 ¯¯¯τF ≡ ∫¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯κRρFry¯¯¯¯¯¯¯Frydz, ¯¯¯τ ≡ ∫¯¯¯¯¯¯¯¯¯κRρdz. (25)

To a reasonable approximation, we find that

 ¯¯¯τF≃{¯¯¯τ;¯¯¯τ<1η¯¯¯τ¯¯¯τ≳1, (26)

with a constant . Using the fact that

 ∂¯¯¯τ∂m=κR,∗(¯¯¯¯TT∗)2, (27)

with , we find

 ¯¯¯τ=34η(κR,∗m)2+√2κR,∗m. (28)

At the base of the domain, and for , which agrees with our measurement of in T3F0.5VET.

If we assume that all of our simulations will approach at late times we can infer that

 τF≈τ∗fE,∗=cgF∗Σ, (29)

which is equivalent to KT12’s equation (43) with . Thus, for a given input flux and constant gravitational acceleration , we can define a characteristic opacity , for which radiative acceleration and gravitational acceleration balance. Then the optical depth can be used in place of in equation (23). In our calculation near the photosphere where . This is in approximate agreement with conclusion of Krumholz & Thompson (2013) that the photospheric opacity is the relevant opacity for evaluating momentum feedback, although it is not clear how precisely this would generalize to other parameters.

Alternatively, we can provide an estimate for by combing equations (28) and (29) to obtain

 η≈23τ∗(√2+3τ∗/fE,∗−√2), (30)

which implies for and , in good agreement with the simulation results. Hence, would be a function of and , with for , although it should be kept in mind that approximation becomes poor for , so there is a limited range where this relation could be physically relevant. Finally, we can use equations (24) and (29) to estimate the optical depth using the temperature at the base of the domain

 κR,0Σ=τ∗(2+3τ∗fE,∗)1/2. (31)

This implies for and , nearly matching the measured value of from T3F0.5VET.

Another question of interest is whether these simulations are consistent with the radiation driving of large scale outflows. Our results show that the simulations approach an Eddington ratio near unity for a constant value of . However, if we assume that the gas in real ULIRGs is distributed in a disk-like geometry the vertical component of gravity will increase from near zero at the midplane (where ) to a maximum value set by the overall potential of the galaxy. It seems plausible that an Eddington ratio near unity will also be reached in this case, but it is not clear what would pick out the specific value of where the radiation forces and gravity balance. However, if this hypothetical equilibrium behaves like our simulations, we would expect a fraction of the gas to be accelerated well beyond this characteristic , because the opacity typically increases while decreases as we approach the midplane. Therefore, gas in low density channels can be efficiently accelerated to high velocities, conceivably approaching the escape velocity from the galaxy. Such acceleration of low density gas to high velocities is seen in our simulations. For example, when we stop the T3F0.5VET run, 4.1% and 1.2% of the gas had greater than and (respectively). In typical ULIRGs, only a few percent of the gas is observed in the cold neutral outflows, so even a small tail of high velocity gas could explain the observed outflow rates and velocities if can be obtained in more realistic simulations.

In this respect, a prominent role of RTI in the dynamics could help solve a potential problem with radiation driving – that it cannot simultaneously explain both the presence of a radiation supported, quasi-equilibrium gas disk and the launching of outflows. If all the gas experienced the same radiative acceleration, it would all be in state of hydrostatic equilibrium or it would all be accelerated in an outflow. In an RTI dominated picture, the bulk of the gas can be in a quasi-equilibrium turbulent state with , but the presence of low density channels with larger than average could still allow a modest fraction of the gas to be accelerated out the gas disk. In principle, this could allow radiation to launch outflows even in galaxies that are well below their global Eddington limit (c.f. Socrates & Sironi, 2013), although more realistic calculations are required to assess the effectiveness of this mechanism.

Finally, we note that all of our analysis and discussion assumes that only radiation and gravitational forces play a role. In real systems, mechanical feedback from stellar winds, supernovae, cosmic rays, magnetic fields, and other sources may be present. Even if they are not the dominant mechanisms of momentum transfer to the gas, they may still have a measurable impact on the gas velocity and density distributions. The implication that RTI will generically produce density inhomogeneities that drive towards a value may break down if any of these other mechanisms are strong enough to reorient the low density channels and dense filaments. It seems much more likely that such reorientation will interfere with, rather than enhance, the escape of photons, leading to increased coupling and conceivably to generically super-Eddington configurations. This hypothesis can be tested with future calculations.

### 5.3. The Role of the Rayleigh-Taylor Instability

We attribute the growth of structure in our simulations primarily to the RTI, but other instability mechanisms may play a role, in principle. Section 5.1 of KT12 provides a fairly comprehensive and persuasive discussion why the RTI is dominant and other instability mechanisms, such as radiative driving of unstable acoustic modes (Shaviv, 2001; Blaes & Socrates, 2003), are likely to be absent or unimportant. Here we simply summarize some of the most salient results and refer the reader to KT12 for further details.

Although a general solution for linear growth of the radiative RTI has never been formulated (see e.g. Mathews & Blumenthal, 1977; Krolik, 1977; Jacquet & Krumholz, 2011; Jiang et al., 2013), numerical experiments (Jiang et al., 2013) suggest that the instability will grow at a rate that can be approximated by , if is a characteristic horizontal wavelength. This is the analytic prediction in the optically thick limit of an incompressible flow if the density contrast at the base of the shell is large (i.e. Atwood number of unity). Radiative diffusion effects can slow the rate of growth, but typically only by factors of order unity unless radiation pressure is much greater that gas pressure. For our simulations, this timescale evaluates to and we have for . The growth of the instability on small scales is even faster and there is ample time for growth of the RTI to explain the large scale non-linear structure observed at (and earlier).

Further evidence that the evolution is due to RTI is provided by calculations we performed with constant, temperature-independent opacities in the limits of pure scattering or pure absorption. The simulation setup for these runs were identical to our T3F0.5VET calculations except that . In the standard runs, with , radiation forces are largest near the base of the domain and decrease with height. Therefore, gas closer to the base of the domain quickly reaches higher velocities than the overlying material and leads to the formation of dense shells with sharp density inversions. In contrast, the gas in the constant opacity simulations receives an acceleration that is significantly more uniform, and no significant density inversions form until very late in the run. As a result, there is very little RTI, and the small scale random perturbations are smoothed out by diffusion before they can grow appreciably. Only the largest scales show non-linear development from the initial sinusoidal perturbation, and the timescale for this growth is much longer. Thus, there is strong evidence that density inversions driven by the opacity law are essential to the dynamics, providing a clear indication that the RTI is the dominant instability mechanism.

## 6. Summary and Conclusions

We have considered the role of the Rayleigh-Taylor instability in the interaction of infrared radiation fields, dust, and gas in rapidly star-forming environments. We have focused on the regime of radiation supported, dense gas that may be present in some systems, such as ULIRGs. Our primary results stem from the numerical simulation of such environments, which are studied in a simplified problem setup with a constant gravitational acceleration, a constant incident infrared flux on the base of the domain, and initialized with a perturbed isothermal atmosphere to match previous calculations in KT12.

In the stable regime, we find that the atmosphere settles down into an equilibrium solution in agreement with the previous results. In the unstable regime, we confirm that the RTI develops and has a significant impact on subsequent evolution. However, we find that after the growth of the RTI, the evolution depends significantly on the choice of algorithm for modeling radiation transfer. Our VET simulations show a stronger coupling between radiation and dusty gas, leading to continuous net upward acceleration of the gas. No steady state is reached before the end of the calculation, when high density material had reached the top of the domain. The mean velocities and velocity dispersion are both increasing at the end of the run. In contrast, our FLD calculations broadly reproduce the FLD results of KT12, finding weaker coupling between gas and radiation. This leads to a short initial burst of acceleration, followed by a period of fallback, finally settling into quasi-steady state turbulence with zero mean velocity and low velocity dispersion. As a result, our VET calculations imply a much larger scale height and higher velocity dispersion than the FLD-based calculations.

We argue that these discrepancies result from limitations in the diffusion-based FLD algorithm, which lead to inaccuracies in modeling how the radiation field responds to structure in the gas distribution that is a few optical depths or smaller in size. These errors are related to the FLD radiation field’s tendency to diffuse around denser, optically-thicker structures even when the diffusion limit does not apply. Relative to our VET calculations, this leads the FLD radiation forces to be more effective at opening and reinforcing low density channels but less effective at disrupting high density filaments, ultimately reducing the coupling between radiation and dusty gas.

Despite the discrepancies in the scale height and velocity distributions, it appears that both the VET and FLD simulations trend towards an approximate balance between the volume-averaged radiation and gravitational forces at late times. If this behavior is general, it confirms one of the key results of KT12 (see also Krumholz & Thompson, 2013), suggesting the rate of momentum transfer between radiation and dusty gas may scale approximately as , where is the luminosity, , is the mass surface density, and is the opacity for which the radiation and gravitational accelerations balance. Since is generally larger nearer to the midplane, will be lower than estimates of that assume volume-average or mid-plane opacities. For example, in the simulation parameters considered here (), we infer a total optical depth using a volume-averaged opacity, which must be reduced by an efficiency factor to obtain the correct rate of momentum transfer. However, if the Eddington ratio is always near unity, this efficiency factor may be lower for gas disks or clouds with larger surface densities.

Although the simulation setup considered here is useful for studying the development and saturation of the RTI, it is not optimally suited for making observational predictions. Future work with more realistic assumptions, such as a vertically varying gravitational acceleration, non-grey opacity, distributed radiation sources, and physically relevant simulation volumes will be necessary to provide robust predictions and facilitate direct comparison with observations.

We thank C.-A. Faucher-Giguer, P. Hopkins, C. Matzner, E. Ostriker, and T. Thompson for useful discussions. N.M. is supported in part by the Canada Research Chair program. Y.F.J. is supported by NASA through Einstein Postdoctoral Fellowship grant number PF-140109 awarded by the Chandra X-ray Center, which is operated by the Smithsonian Astrophysical Observatory for NASA under contract NAS8-03060. Computations were performed on the GPC supercomputer at the SciNet HPC Consortium. SciNet is funded by the Canada Foundation for Innovation under the auspices of Compute Canada, the Government of Ontario, Ontario Research Fund-Research Excellence, and the University of Toronto.

### .1. Impact of Angular Resolution

The discussion in section 5.1 assumes that the VET radiative transfer algorithm provides a more accurate realization of the radiation field than the FLD algorithm. Since we solve the momentum equation (7), this should be true, as long as our method for computing the Eddington tensor is valid. As discussed in section 2.1, the main approximation made in our computation of the Eddington tensor is the neglect of order terms (which are retained in the radiation moment equations). We confirm a posteriori that this is a very good approximation.

The only significant potential problem for our VET calculation would be if our angular grid of rays under-resolves the radiation field. Our experience has been that angular resolution used here is more than adequate for previous problems, but angular resolution requirements can vary from problem to problem. In the short characteristics method, the most apparent features in under-resolved systems are usually referred to as “ray-effects”, which correspond to unphysical anisotropies in the radiation field correlated with the angular grid. Ray-effects can be particularly problematic when bright point sources dominate the radiation field and scattering opacity is negligible (see e.g. Larsen & Wollaber, 2008; Finlator et al., 2009). The intensity can be overestimated in grid zones that lie along ray directions that point directly back to the point source and underestimated in regions between rays, with a relative error that increases in magnitude with the number of grid zones traversed. This leads to the appearance of ‘rays’ or ‘spokes’ in radiation variables emanating radially outward from the point source.

Scattering introduces angular diffusion that can mitigate ray-effects (see e.g. Larsen & Wollaber, 2008), but is absent in the current calculations. Although bright point sources are also absent, dense filaments can lead to rather sharp variations in the emissivity. Ray-effects are most clearly apparent in the low density regions well above the photosphere, reflecting variations in the emissivity near the photosphere that propagate along rays in the angular grid into the regions above. However, there is very little mass in this region and radiation is less important because the radiation forces are sub-Eddington here.

Ray-effects are less apparent below the photosphere, but must be present at some level. In order to assess their impact (if any), we reran the T3F0.5VET run with an alternative ray distribution that has effectively higher resolution in the plane. Our default ray distribution is shown in the left panel of figure 10. It attempts to cover the unit sphere uniformly, treating all three spatial directions, including the implied third dimension () on equal footing. This is formally required even though the domain is two dimensional, because rays traveling at different angles relative to will travel different total path lengths for a given displacement in the plane. This means that many of the rays nearly overlap when projected onto the plane, and the effective angular resolution in the plane is closer to rays per radians, even though we have 84 rays total. The right panel shows our alternative distribution, in which all rays have the same projection onto the -axis ( so that ). We distribute 80 total rays uniformly in azimuthal angle (). This amounts to a factor increase in the effective angular resolution in the plane even though the total number of angles is nearly identical.

The higher effective azimuthal angular resolution significantly reduced ray-effects above the photosphere, but had relatively little effect on the Eddington tensor below the photosphere, suggesting ray effects were already minor in the original run. Overall, the global properties and evolution of the two simulation were nearly identical. A sense of how similar the runs are is provide by figure 11, which compares snapshots of from the two simulation at . We therefore conclude that our results are well-converged in terms of the effective angular resolution of our VET calculations.

## References

• Andrews & Thompson (2011) Andrews, B. H., & Thompson, T. A. 2011, ApJ, 727, 97
• Blaes & Socrates (2003) Blaes, O., & Socrates, A. 2003, ApJ, 596, 509
• Chandrasekhar (1961) Chandrasekhar, S. 1961, Hydrodynamic and hydromagnetic stability
• Davis et al. (2012) Davis, S. W., Stone, J. M., & Jiang, Y.-F. 2012, ApJ
• Downes & Solomon (1998) Downes, D., & Solomon, P. M. 1998, ApJ, 507, 615
• Finlator et al. (2009) Finlator, K., Özel, F., & Davé, R. 2009, MNRAS, 393, 1090
• Hayes & Norman (2003) Hayes, J. C., & Norman, M. L. 2003, ApJS, 147, 197
• Heckman et al. (1990) Heckman, T. M., Armus, L., & Miley, G. K. 1990, ApJS, 74, 833
• Hopkins et al. (2011) Hopkins, P. F., Quataert, E., & Murray, N. 2011, MNRAS, 417, 950
• Hopkins et al. (2012) —. 2012, MNRAS, 421, 3522
• Jacquet & Krumholz (2011) Jacquet, E., & Krumholz, M. R. 2011, ApJ, 730, 116
• Jiang et al. (2013) Jiang, Y.-F., Davis, S. W., & Stone, J. M. 2013, ApJ, 763, 102
• Jiang et al. (2012) Jiang, Y.-F., Stone, J. M., & Davis, S. W. 2012, submitted to ApJS
• Kennicutt (1998) Kennicutt, Jr., R. C. 1998, ApJ, 498, 541
• Krolik (1977) Krolik, J. H. 1977, Physics of Fluids, 20, 364
• Krumholz et al. (2007) Krumholz, M. R., Klein, R. I., McKee, C. F., & Bolstad, J. 2007, ApJ, 667, 626
• Krumholz & Matzner (2009) Krumholz, M. R., & Matzner, C. D. 2009, ApJ, 703, 1352
• Krumholz & Thompson (2012) Krumholz, M. R., & Thompson, T. A. 2012, ApJ, 760, 155
• Krumholz & Thompson (2013) —. 2013, MNRAS, 434, 2329
• Kunasz & Auer (1988) Kunasz, P., & Auer, L. H. 1988, J. Quant. Spec. Radiat. Transf., 39, 67
• Larsen & Wollaber (2008) Larsen, E. W., & Wollaber, A. B. 2008, Nuclear Science and Engineering, 160, 267
• Levermore & Pomraning (1981) Levermore, C. D., & Pomraning, G. C. 1981, ApJ, 248, 321
• Lowrie et al. (1999) Lowrie, R. B., Morel, J. E., & Hittinger, J. A. 1999, ApJ, 521, 432
• Mathews & Blumenthal (1977) Mathews, W. G., & Blumenthal, G. R. 1977, ApJ, 214, 10
• Mihalas (1978) Mihalas, D. 1978, Stellar atmospheres /2nd edition/
• Mihalas & Klein (1982) Mihalas, D., & Klein, R. I. 1982, Journal of Computational Physics, 46, 97
• Mihalas & Mihalas (1984) Mihalas, D., & Mihalas, B. W. 1984, Foundations of radiation hydrodynamics, ed. Mihalas, D. & Mihalas, B. W.
• Murray et al. (2005) Murray, N., Quataert, E., & Thompson, T. A. 2005, ApJ, 618, 569
• Pettini et al. (2001) Pettini, M., Shapley, A. E., Steidel, C. C., Cuby, J.-G., Dickinson, M., Moorwood, A. F. M., Adelberger, K. L., & Giavalisco, M. 2001, ApJ, 554, 981
• Rupke et al. (2005) Rupke, D. S., Veilleux, S., & Sanders, D. B. 2005, ApJS, 160, 115
• Schwartz & Martin (2004) Schwartz, C. M., & Martin, C. L. 2004, ApJ, 610, 201
• Scoville et al. (2001) Scoville, N. Z., Polletta, M., Ewald, S., Stolovy, S. R., Thompson, R., & Rieke, M. 2001, AJ, 122, 3017
• Sekora & Stone (2010) Sekora, M. D., & Stone, J. M. 2010, Journal of Computational Physics, 229, 6819
• Semenov et al. (2003) Semenov, D., Henning, T., Helling, C., Ilgner, M., & Sedlmayr, E. 2003, A&A, 410, 611
• Shaviv (2001) Shaviv, N. J. 2001, ApJ, 549, 1093
• Shestakov & Offner (2008) Shestakov, A. I., & Offner, S. S. R. 2008, Journal of Computational Physics, 227, 2154
• Socrates & Sironi (2013) Socrates, A., & Sironi, L. 2013, ApJ, 772, L21
• Stone et al. (2008) Stone, J. M., Gardiner, T. A., Teuben, P., Hawley, J. F., & Simon, J. B. 2008, ApJS, 178, 137
• Thompson et al. (2005) Thompson, T. A., Quataert, E., & Murray, N. 2005, ApJ, 630, 167
• Wise et al. (2012) Wise, J. H., Abel, T., Turk, M. J., Norman, M. L., & Smith, B. D. 2012, MNRAS, 427, 311
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters