RADIATION FEEDBACK IN ULIRGs: ARE PHOTONS MOVERS AND SHAKERS?

The above equations are incomplete because there is no evolution equation for the radiation pressure tensor. We address this by computing the eponymous VET ${\sf f} \equiv {\sf P}_r/E_r$. We approximate ${\sf f}$ by solving the time-independent equation of radiation transfer of the form

$$\begin{equation} \hat{n} \cdot \nabla I = \sigma _{F} \left(\frac{a_r c}{4\pi } T^4 - I\right), \end{equation} \tag{ 8 }$$

where I is the specific intensity for an angle defined by the unit vector $\hat{n}$. Hence, we have adopted a gray opacity treatment in which the characteristic opacity is assumed to correspond to the flux mean opacity. For the computation of the $\sf f$, we make no distinction between the comoving and Eulerian frames, since Equation (8) drops all order v/c terms or higher, including a c⁻¹∂I/∂t term that would otherwise be present. The neglect of the c⁻¹∂I/∂t term is typically a good approximation whenever v/c ≪ 1, but the neglect of velocity-dependent terms generally requires τv/c ≪ 1 if τ is a characteristic optical depth for the system. In the results presented here v_max/c  ∼  10⁻⁴ and τ_max  ∼  15, 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 J, the first moment H, and the second moment $\sf K$ are tabulated in each grid zone. These are then used to compute ${\sf f} ={\sf K}/J$ for the subsequent timestep. For our two-dimensional (2D) simulations, we integrate the radiation field along 84 rays, distributed nearly uniformly over the half unit sphere, taking advantage of the reflection symmetry of 2D Cartesian domains. Our single three-dimensional (3D) run uses 168 rays. Further discussion of the choice of angular grid and impact of angular resolution is provided in the Appendix.

We have also implemented an 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 (Y.-F. 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

$$\begin{equation} \mathbf {F}_r=-\frac{c \lambda }{\sigma _{F}} \nabla E_r, \end{equation} \tag{ 9 }$$

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

$$\begin{eqnarray} \lambda &= &\frac{1}{R}\left(\coth {R} - \frac{1}{R}\right)\nonumber \\ R &=&\frac{| \nabla E_r |+\beta }{\sigma _{F} E_r}. \end{eqnarray} \tag{ 10 }$$

To specify ${\sf P}_r$, we use the Eddington tensor

$$\begin{eqnarray} &&{\sf f} = \frac{1}{2}(1-f){\sf I}+\frac{1}{2}(3 f-1)\hat{f}\hat{f} \end{eqnarray} \tag{ 11 }$$

with f = λ + λ²R² and $\hat{f} = \nabla E_r/|\nabla E_r|$.

Following Shestakov & Offner (2008), we have added a β parameter to the definition of R, which acts as an effective floor on R 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 β ≳ 10⁻⁴. 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 F_r 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 v/c 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 $\sf f$. Modulo the velocity-dependent terms that we neglect in computing $\sf f$ (and only in computing $\sf f$), 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.⁶ For realistic problems, we can enhance the accuracy of our VET solution by increasing the angular and spatial resolution of our grid (Y.-F. Jiang et al., in preparation). There is no comparable mechanism for improving accuracy in FLD. Simply increasing the spatial resolution will not make FLD converge to the correct solution. 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.

As in KT12, we assume an initially static and isothermal atmosphere with T = T_*. The density is initialized as an exponential profile with scale height h_*, 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 10⁻¹⁰ρ_*, where ρ_* = Σ_*/h_* is the maximum initial density at the lower boundary. We assume T_r = T everywhere, which gives a constant $E_r=a_r T_*^4=c F_*$. Hence, for even moderate values of τ_* and f_{E, *}, 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

$$\begin{equation} \frac{\delta \rho }{\rho } = \left\lbrace \begin{array}{@{}ll}0.25 (1 \pm \chi) \sin {(2 \pi x/\lambda _x)} & {\rm 2D}, \\ 0.25 (1 \pm \chi) \sin {(2 \pi x/\lambda _x)}\sin {(2 \pi y/\lambda _y)} & {\rm 3D}, \end{array}\right. \end{equation} \tag{ 15 }$$

where λ_{(x, y)} = 0.5L_{(x, y)} and L_{(x, y)} are the domain widths. Here, χ is identically zero for sinusoidal perturbations or can be a random number uniformly distributed between −0.25 and 0.25. If χ = 0, the perturbation is identical to KT12, but we found it useful to consider simulations with random perturbations on top of the sinusoidal variation.

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

For all simulations, both the radiation and hydrodynamic quantities obey periodic boundary conditions in the horizontal direction (x-coordinate in 2D; x- and y-coordinates in 3D). For the hydrodynamic quantities, we impose reflecting boundary conditions at the lower vertical (y-coordinate in 2D, z-coordinate in 3D) 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 E into the ghost zones at our upper boundary, whereas KT12 fix ρ = 10⁻¹³ρ_* and T = 10³T_* there.

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

$$\begin{equation} F_{ry}=F_*+v_y E_r+(\mathbf {v} \cdot {\sf P}_r)_y. \end{equation} \tag{ 16 }$$

Hence, the comoving flux in the ghost zones is equal to F_*. We specify E_r 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 F_r are continued in to the ghost zones and we specify E_r using

$$\begin{equation} E_r=\frac{J}{H_y} F_{ry}, \end{equation} \tag{ 17 }$$

where H and J 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 E_r using Equation (9), assuming F_ry = F_* and no horizontal gradient in energy density. At the upper boundary, we assume a fixed $E_r=a T_*^4$ 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 J₋. We then assume isotropic incoming (upward) radiation. The intensity of the incoming radiation is normalized so that the corresponding angle integrated upward intensity J₊ obeys H_y = J₊ − J₋ = F_*. 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).

We first consider the T10_F0.02 run, which we run using our VET module. The parameters for this run place it in the stable regime, according to the 1D equilibrium solutions and the 2D 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 80t_*.

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 E_r to increase. Since the radiation and dusty gas are well-coupled, the increase in T_rad drives a corresponding rise in T. 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.

Zoom In Zoom Out Reset image size

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

$$\begin{equation} \langle \mathbf {v} \rangle = \frac{1}{M}\int \rho \mathbf {v} dV, \end{equation} \tag{ 18 }$$

and mass-weighted velocity dispersion

$$\begin{equation} \sigma ^2_{i} = \frac{1}{M} \int \rho (v_i-\langle v_i\rangle)^2 dV, \end{equation} \tag{ 19 }$$

where M = ∫ρdV is the total mass in the atmosphere. The top panel of Figure 2 shows σ_x and σ_y, along with the total velocity dispersion ($\sigma =(\sigma _x^2+\sigma _y^2)^{1/2}$), and the bottom panel shows 〈v_y〉.

Zoom In Zoom Out Reset image size

The initial transient acceleration lead to growth in both the vertical and horizontal velocity dispersion, although the vertical component initially grows faster. After ∼3t_*, both 〈v_y〉 and σ_y have already reached their maximum for the simulation. After another ∼10t_*, both quantities settle down into oscillations, with 〈v_y〉 varying about zero. The horizontal velocity dispersion is also oscillatory, albeit on a longer timescale with a period greater than 50t_*. 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 〈v_y〉. 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 80t_*, which is long enough to confirm that we reproduce the KT12 results with our VET formalism for this stable regime.

We now focus on the unstable regime, considering simulations with τ_* = 3 and f_{E, *} = 0.5. This run has the lowest ratio of f_{E, *}/f_{E, crit} 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 (T3_F0.5_FLD) and one using the VET module (T3_F0.5) 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 an 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. As we shall see, the VET run shows stronger coupling between gas and radiation than found in the FLD run. The inferred coupling would have been even stronger if we had not added random perturbations to seed the RTI.

We first consider the results from T3_F0.5_FLD, which we run for 200t_*. This is run in a taller box than the T10_F0.02 simulation, with the initial condition and boundary conditions described in Section 3. Note that we began the simulations with β = 10⁻⁴ in Equation (10) to place a floor on R and avoid convergence problems in our backward Euler scheme. This was sufficient for the first ∼100t_*, but we had to increase β = 4 × 10⁻⁴ at t = 100t_* to ensure convergence thereafter.

As in the T10_F0.02 run, there is an initial increase in T_rad 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 T drives an increase in κ_R, 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 T_rad from four different times. By 25t_*, 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 structures.

Zoom In Zoom Out Reset image size

Subsequent evolution is also consistent with KT12. The RTI plumes grow into high-density filaments interspersed with lower density channels. Since F_ry is largest in the low-density channels, there is an anti-correlation of F_ry 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 50t_*. This behavior is reinforced by a drop in the volume averaged optical depth, which leads to a drop in T_rad (and therefore T) at the base of the atmosphere, further reducing the radiation force. By 75t_*, almost all of the mass is contained in the bottom ∼50h_* of the domain and remains there for the duration of the run, which we stop at 200t_*.

As expected from the lack of a 1D 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 T3_F0.5. This simulation began in a domain with the same dimensions, resolution, and initial condition as in the T3_F0.5_FLD run. Figures 4 and 5 show four snapshots of ρ and T_rad (respectively) from T3_F0.5. Initially, the dynamics of the simulation are qualitatively consistent with the T3_F0.5_FLD run. Again, T_rad and T 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 39t_*.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The difference between the FLD and VET runs become more apparent in the subsequent evolution. As in the T3_F0.5_FLD 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 ρ ≳ 10⁻⁴ρ_* gas, although most of the mass remains in several dense filaments.

At ∼83t_*, 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 T_rad at the upper boundary. Therefore, we restart the simulation at t = 80t_* in a domain with the L_y = 2048h_* and N_y = 4096 (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, F_r and E_r, 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 ∼83t_* when the breakdown of the vertical boundary condition alters the evolution in the smaller domain. We continue our integration in the larger domain until 158t_*, 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 ρ ≳ 10⁻⁴ρ_*, except for transient periods when most of the volume near the bottom of the domain can have ρ ≲ 10⁻⁶ρ_*.

Figure 5 shows that there is relatively little horizontal variation in T_rad, consistent with T3_F0.5_FLD and KT12. In contrast to T3_F0.5_FLD, 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 T ≃ T_rad  ∼  2.2T_*. Modest but short-lived increases in T_rad are seen and correspond to fall back and reflection of dense filaments at the lower boundary (see, e.g., the second panel at 78t_*). These values are modestly higher than the mean values of T_rad at the base in T3_F0.5_FLD, although there is greater fluctuation in T_rad in the T3_F0.5_FLD run.

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

$$\begin{equation} f_{\rm E, V}=\frac{\langle \kappa _{\rm R} \rho F_{ry} \rangle }{c g \langle \rho \rangle }, \end{equation} \tag{ 20 }$$

where 〈 〉 denote volume averages (e.g., $\langle \rho \rangle = L_x^{-1} L_y^{-1} \int \int \rho dx dy$). 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. The Eddington ratio is initially slightly greater in the VET run than FLD, but some of this is due to issues with our VET boundary condition and deviation of the initial condition from thermal equilibrium. The value of F_ry initially exceeds the target F_* (up to ∼37%), but declines to with 1% of F_* by 5t_*, accounting for some of the difference. Both simulations show a return to an Eddington ratio near unity, although T3_F0.5_FLD tends to have f_{E, V} fluctuating near but typically just below unity, while T3_F0.5 tends to remain moderately super-Eddington for most of the run. Although the differences are rather modest, the upshot is that T3_F0.5 receives a nearly continual, modest upward acceleration while T3_F0.5_FLD settles into a quasi equilibrium state. This evolution is apparent in the evolution of 〈v_y〉 in Figure 7. The trend in f_{E, V} suggests that T3_F0.5 might settle into a similar equilibrium, but with a much larger-scale height and higher velocity dispersion.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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

$$\begin{equation} \tau _{\rm V} = L_y \langle \kappa _{\rm R} \rho \rangle, \end{equation} \tag{ 21 }$$

and the second is the flux-weighted mean optical depth

$$\begin{equation} \tau _{\rm F} = L_y \frac{\langle \kappa _{\rm R} \rho F_{ry}\rangle }{\langle F_{ry} \rangle }. \end{equation} \tag{ 22 }$$

The middle panel of Figure 6 shows that τ_V correlates closely with f_{E, V} as both evolve over the course of the simulations. Large values of τ_V correspond to large values of κ_R∝T² because T ≃ T_rad and T_rad is roughly proportional to $\tau _{\rm V}^{1/4}$ throughout much of the domain. Note, however, that it is not quite true that τ_V determines the Eddington ratio as f_{E, V} can be larger in T3_F0.5, even when τ_V is greater for T3_F0.5_FLD. This occurs primarily because the nature of anti-correlation between F_ry and ρ differs in the two runs.

The bottom panel of Figure 6 shows the ratio τ_F/τ_V. Note that τ_F multiplies 〈F_ry〉 = F_* to give the momentum per unit area transferred from the radiation to the gas while τ_V would be the characteristic value for the total infrared optical depth τ_IR 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 τ_F = f_trap + 1, where f_trap is trapping factor discussed in KT12. Our estimate of τ_F ≃ 6 agrees with KT12's f_trap = 5 for τ_* = 3, f_{E, *} = 0.5.

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 F_r 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 〈v_y〉 and σ. For both simulations, the evolution of 〈v_y〉 is largely determined by the value of f_{E, V} in Figure 6. When f_{E, V} > 1, 〈v_y〉 increases and when f_{E, V} < 1, 〈v_y〉 decreases. In both simulations, there is an early period of transient growth, followed by a decline after the RTI sets in. For T3_F0.5_FLD, 〈v_y〉 becomes negative as most of the mass falls back, but eventually settles into a quasi-steady state, with 〈v_y〉 simply fluctuating near zero. For T3_F0.5, 〈v_y〉 grows for the majority of the run.

For both runs, the x and y components of σ show a modest initial increase, followed by a faster rises in σ_y as the RTI sets in and the shells start to break apart. The subsequent evolution for T3_F0.5_FLD involves a weak decline in σ_y compensated by a slow increase in σ_x. These trends roughly cancel so that σ fluctuates about ∼5c_* for the majority of the simulation, consistent with the results of KT12. The T3_F0.5 evolution displays a more continuous rise in σ, with declines at ∼75t_* and 170t_*, when 〈v_y〉 increases quickly during periods of higher f_{E, V}. When we stop the run near 250t_*, σ is increasing, but still below the maximum, which was a factor of ∼3 larger than in T3_F0.5_FLD.

Our initial focus on τ_* = 3 and f_{E, *} = 0.5 is motivated by a desire to compare with KT12, but also by the fact that this gives parameters that are broadly consistent with observational constraints on luminous ULIRGs, like Arp 220, which has one of the highest inferred surface densities and star formation rates in the local universe. Many local star forming galaxies have lower surface densities and, therefore, lower dust optical depths. Higher optical depths may also be obtained in local regions of star forming galaxies and may be common in galaxies at high redshift.

These considerations motivate runs with τ_* = 1 (T1_F0.5) and τ_* = 10 (T10_F0.5), both assuming f_{E, *} = 0.5. Figure 8 compares f_{E, V}, τ_V, and τ_F/τ_V for these runs, along with T3_F0.5. The T10_F0.5 (solid black) shows a similar profile for f_{E, V} as in T3_F0.5 (red dotted), but with a notably larger initial transient. This follows from the larger τ_V ≳ 50 in this run and the non-linear dependence of τ_V on τ_* results from the T² dependence of the opacity. However, we still have f_{E, V}  ∼  1 due to the regulating effects of the RTI, which show up as a lower ratio of τ_F/τ_V. Figure 9 shows the evolution of 〈v_y〉 and σ. Due to the large initial transient, 〈v_y〉 grows quickly for T10_F0.5 and a significant fraction of the mass reaches the top of simulation (L_y = 2048h_*) at t = 65t_*, considerably faster than the t  ∼  150t_* that it takes for the T3_F0.5 to reach the same height. Although σ increases faster in the T10_F0.5 run, σ  ∼  15c_* at the end of both simulations.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The evolution of the T1_F0.5 (dashed blue) is qualitatively different. There is a very brief initial transient after which f_{E, V} ≲ 1 for the remainder of the run. We find τ_V  ∼  2 and τ_F/τ_V is always slightly below unity. As a result, there is never any significant net upward acceleration. The mass that is initially accelerated upward falls back and 〈v_y〉 oscillates about zero. The mass settles down into a turbulent state reminiscent of the T3_F0.5_FLD run, but with even smaller velocity dispersion.

Since the fiducial Eddington ratio f_{E, *} = 0.5, the mass-weighted average opacity would have to increase by a factor of two in order to obtain f_{E, V}  ∼  1, which requires the mass-weighted average of T² to be $2 T_*^2$. If the flow were homogeneous this would be possible, but only a very small reduction in the trapping factor will prevent this. As we will see in the discussion section (Equation (32)), we need τ_F  ∼  0.94τ_V to maintain f_{E, V}  ∼  1. This value is approached in Figure 8, but not sustained.

We have also considered the impact of spatial resolution on our simulations. Most of our simulations assume L_{(x, y)}/N_{(x, y)} = 0.5h_*, leading to a large number of grid zones in both spatial dimensions. These high resolutions are attainable in 2D simulations but would be very computationally demanding for 3D calculations. Therefore, we have performed simulations at a medium resolution L_{(x, y)}/N_{(x, y)} = h_* (T3_F0.5_MR) and low resolution with L_{(x, y)}/N_{(x, y)} = 2h_* (T3_F0.5_LR) to study the impact of resolution.

Figure 10 compares T3_F0.5 (black solid), T3_F0.5_MR (red dotted), and T3_F0.5_LR (blue dashed). The lower-resolution runs evolve with a time dependence similar to the high-resolution runs, but spend most of their time at lower Eddington ratios. At early times, the RTI develops more slowly as resolution decreases. However, the development of inhomogeneities in the non-linear phase has a larger affect. The Eddington ratio falls more quickly and takes longer to recover to f_{E, V} > 1. Both τ_V and τ_F/τ_V are lower, consistent with a stronger flux-density anticorrelation in these lower-resolution runs. The mass is more concentrated in narrower and denser filaments in the lower-resolution runs. This can be seen by comparing the density snapshot from T3_F0.5_LR is shown at t = 100t_* in the right-hand panel of Figure 11 with t = 78t_* snapshot of T3_F0.5 in Figure 4.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In spite of these differences, the lower-resolution runs still have a significant amount of moderately dense (ρ  ∼  10⁻³ρ_*) material reaching the top of the box and we need to stop both the T3_F0.5_MR and T3_F0.5_LR runs at t  ∼  103t_*, when enough mass has reached the upper boundary to invalidate the vacuum boundary condition on the incoming radiation. In comparison, the T3_F0.5 simulation reached the same point at t = 83t_*. Furthermore, the bulk of the mass has received much less acceleration in the low-resolution runs and has a lower 〈v_y〉 and σ. Most of the mass remains in the lower half of the box in both T3_F0.5_MR and T3_F0.5_LR, whereas ρ is more uniformly distributed with height in the T3_F0.5 run.

The upshot is that higher resolution leads to a more diffuse density distribution and enhances the volume-averaged radiation force. There is a slight trend of increasing f_{E, V} from T3_F0.5_LR to T3_F0.5_MR, but the two low-resolution runs are rather similar. This may imply that even higher resolution would further enhance or feedback or could simply indicate that one needs to resolve the h_* scale in the simulation. We have not performed higher-resolution runs to test for convergence due to the computational expense, but this should be explored further in future applications.

It is well known that the non-linear development of hydrodynamics RTI differs in 2D and 3D simulations with and without magnetic fields (e.g., Young et al. 2001; Stone & Gardiner 2007). However, we know of no published work that considers the development of radiative RTI in 3D. Since our results thus far have suggested that modest changes in the angular distribution of the radiation field can have a large affect on the dynamics, it is important to see how well our conclusions from 2D simulations carry over the 3D case. Due to the substantial increase in computational cost, we must perform a single 3D calculation (T3_F0.5_3D) at lower resolution, equivalent to the T3_F0.5_LR run in 2D.

Figure 12 compares f_{E, V}, τ_V, and τ_F/τ_V for T3_F0.5_3D (solid black) and T3_F0.5_LR (dotted red). The evolution of all three parameters overlaps for the first t  ∼  25t_*, when the flow is in the linear or early non-linear phases of RTI growth. Subsequent growth shows a significant enhancement for all three quantities in the 3D run. In the 3D run, f_{E, V} remains above unity from t = 48t_* until the end of the run while the 2D run fluctuates about unity. This leads to a monotonic increase in 〈v_y〉. Both simulations must be stopped when ρ  ∼  10⁻³ρ_* material reaches the upper boundary. Both simulations finish with σ  ∼  7c_*, but 〈v_y〉 = 8c_* in T3_F0.5_3D, more than twice the final value of 3.5c_* in T3_F0.5_LR.

Zoom In Zoom Out Reset image size

These results suggest that the Eddington ratio is modestly enhanced in 3D. Figure 11, which shows a 2D slice through the ρ distribution of T3_F0.5_3D, provides some indication of why this enhancement occurs. While the 2D density distribution is dominated by a number of narrow filaments, ρ is more diffusely distributed in a larger number of lower density filaments. The overall density field is less structured inhibiting the escape of radiation through low-density channels and enhancing the coupling between matter and radiation. These results suggest that the dimensionality of the simulation will matter for properly assessing the role of radiative driving in dust environments.

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 f_{E, V}  ∼  1 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. The top panel of figure 13 shows a snapshot of ρ for the T3_F0.5 run at t = 25t_*. We focus on a small fraction of the domain where the shell is being disrupted. The bottom panel compares F_r with the FLD prediction for the flux using Equation (9) with the values E_r in the VET run. The color map shows the flux ratio on a logarithmic scale while arrows denote the directions of the VET (white) and FLD (green) fluxes.

Zoom In Zoom Out Reset image size

In the densest regions, our VET calculations are nearly in the diffusion limit where the FLD method is reliable, and the fluxes are in approximate agreement. In lower-density regions, where the diffusion approximation is poor, the magnitude of the FLD flux is much larger (up to factors ∼30) than VET 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 F_r, 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 x  ∼  0, − 115h_* and with 160 ⩽ y/h_* ⩽ 190 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 13 is that the values of E_r in the FLD run will not generally be the same as those in the VET run. Figure 14 avoids this issue by comparing ρ and F_ry/E_rc in the T3_F0.5 run to the same quantities in the T3_F0.5_FLD run. In both cases, we focus on a 256h_* × 300h_* 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.

Zoom In Zoom Out Reset image size

Since the mean-free path of photons is proportional to 1/(κρ), F_ry 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 F_ry with E_r. Such anti-correlations are expected to be particularly strong if the gradient in E_r is relatively uniform and the diffusion limit applies (Equation (9) with λ → 1/3). However, the optical depths across most of the filaments seen in Figure 14 are of the order of 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 sight 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 F_ry/E_rc 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 ∇E_r) and the ratio F_ry/E_rc is determined by R (Equation (10)). Since |∇E_r|/E_r and T are relatively uniform, the sharp variation of ρ dominates the variation of R and therefore λ(R). This strong dependence of λ(R) on ρ leads to an unphysically high level of anti-correlation between F_ry/E_rc and ρ in low-to-moderate optical depth regions.

The upshot is that FLD tends to overestimate the contrast in F_r 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 (τ ≳ 5) 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 f_{E, *} and large τ_*), the radiation force can only exceed the gravitational force at high optical depth, where T, and therefore κ_R, is sufficiently large.

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 h_* = 6.4 × 10⁻⁴ pc (for T_* = 82 K and f_{E, *} = 0.5) and L_y = 2048h_* = 1.3 pc. 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 g, the approximate gray 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 σ. For T3_F0.3, σ reaches a maximum ≃ 14c_* ≃ 7.5 km s⁻¹. However, since no steady-state is reached and the gas appears to be unbound in this simple setup, there is no evidence that this value of σ would be characteristic of realistic self-gravitating system for either T3_F0.3 or T10_F0.3. It is difficult to assess what implications the trend toward f_{E, V}  ∼  1 will have for the subsequent evolution of 〈v_y〉 and σ. The system is very sensitive to the precise value of f_{E, V} when g is constant. If, on average, f_{E, V} 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 g would presumably play a significant role in determining an equilibrium scale height and velocity dispersion.

In contrast, the T1_F0.3 run with τ_* = 1 seems to settle into a steady state with σ ≲ 3c_*. This is perhaps not surprising because for τ_*  ∼  1 the trapping factor is not expected to be particularly large to begin with. Hence, for lower optical depth we conclude that radiative driving is likely ineffective due to the RTI, although this conclusion may be sensitive to the simplifying assumptions discussed above. This could mean that radiative feedback does not play a dominant role in less extreme starbursts and lower surface density ULIRGs (see, e.g., Krumholz & Thompson 2013) unless the radiation approaches the Eddington limit. Or, it may be that feedback in systems with lower mean optical depths might be dominated by regions with significantly higher than average optical depth. Assessing the role of radiative feedback in such systems will be an important goal for future work.

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 g. However, if we assume that the gas in real ULIRGs is distributed in a disk-like geometry, the vertical component of gravity g(y) will increase from near zero at the midplane (where y = 0) 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 select the specific value of g₀ 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 g₀ because the opacity typically increases while g(y) 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 T3_F0.5 run, 4.1% and 1.2% of the gas had v_y greater than 3σ_y and 4σ_y (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 σ_y  ∼  100 km s⁻¹ 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 principle, this could be resolved if outflows were driven primarily from localized regions that significantly exceed the local Eddington ratio. However, in an RTI-dominated picture, the bulk of the gas can be in a quasi-equilibrium turbulent state with 〈v_y〉  ∼  0, but the presence of low-density channels with larger than average F_ry 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 (see Socrates & Sironi 2013), although more realistic calculations are required to assess the effectiveness of this mechanism.

Another question is whether the inhomogeneities driven by the RTI have a notable impact on the spectrum emitted by the dusty gas. The simulation domain remains very optically thick to UV radiation, so unless the young stars are very near the photosphere, the inhomogeneities will not significantly alter the spectrum at these wavelengths. However, the inhomogeneities could, in principle, allow emission from hotter regions to escape, increasing the mid-infrared flux. We consider this possibility by using our short characteristics solver to generate a spectrum from a single snapshot of our T3_F0.5 run at t = 75t_*. For this exercise, we assume radiative equilibrium and only consider the dust absorption opacity. We model the frequency-dependent dust opacity using a model with R_V = 5.5 and assuming a gas to dust mass ratio of 105 (Weingartner & Draine 2001; Draine 2003). Although we do not have detailed knowledge of the dust opacity for ULIRGs and high-redshift galaxies, a higher value of R_V is typical for giant molecular clouds in our galaxy and these models have previously been used provide satisfactory models of star forming galaxy SEDs (e.g., Chakrabarti & McKee 2005).

Figure 15 shows the resulting spectrum compared to a calculation with the same mass surface density, but, assuming uniform temperature model with T = 97 K, which approximately corresponds to the average photospheric temperature in the simulation. Some enhancement over the isothermal model is apparent at all wavelengths, because radiation coming from deeper, hotter regions can escape, but this enhancement is larger at shorter wavelengths. The spectrum peaks at shorter wavelengths than are typically observed in ULIRG spectral energy densities, but this follows from the simulation's assumption of a fairly large photospheric temperature. The assumed temperature might be characteristic of hotter regions within a system like Arp 220, and the emission at ≲ 10 μm would not be out of line with observed systems. Therefore, it does not seem to be the case that the inhomegeneities will allow too much short wavelength flux to emerge.

Zoom In Zoom Out Reset image size

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 $t_{\rm g}^{-1}\,{\sim}\,\sqrt{2\pi g /\lambda _{ x}}$, if λ_x 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 the order of unity unless radiation pressure is much greater that gas pressure. For our simulations, this timescale evaluates to $t_{\rm g} = t_* \sqrt{\lambda _{ x}/2\pi h_*}$ and we have t_g  ∼  6t_* for λ_x  ∼  256h_*. 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 t = 25t_* (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 T3_F0.5 calculations except that f_* = 1.25. In the standard runs, with κ_R∝T², 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 κ_R∝T² opacity law are essential to the dynamics, providing a clear indication that the RTI is the dominant instability mechanism.

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. The primary question of interest is what is the correct 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)

$$\begin{equation} \frac{d p}{d t} = (1 + \eta \tau _{\rm IR}) \frac{L}{c}, \end{equation} \tag{ 23 }$$

where dp/dt is an estimate of the momentum transferred to the interstellar gas, τ_IR is an estimate of the integrated optical depth through the dusty gas at infrared wavelengths, L is luminosity of a star (or star cluster), and η is an ad hoc reduction factor (η ≲ 1) included to account for inhomogeneities in the surrounding gas.⁸ Since τ_V corresponds to a volume-weighted average total optical depth, one can estimate τ_IR ≃ τ_V, while τ_F represents the effective optical depth for momentum exchange. Therefore, we estimate the η in Equation (23) as the ratio τ_F/τ_V = 0.87, 0.68, and 0.4 for runs T1_0.5, T3_0.5, and T10_F0.5, respectively. Here we have time-averaged the runs t = 30t_* to the end of the simulation. Values of η ≃ 0.4–0.87 are in the range considered by Hopkins et al. (2011) in their calculations.

Our values of τ_F ≃ 6 and 19 for T3_F0.5 and T10_F0.5 (respectively) roughly agree with KT12's f_trap = 5 and 22 for the same runs (since τ_F = f_trap + 1). However, in their discussion, KT12 argue that Hopkins et al. (2011) and others likely overestimate the enhancement due to photon trapping. This may be, in part, because values of η > 1 were considered. It also lies partly in how one estimates τ_IR and partly in the dependence of τ_V and τ_F on τ_* and f_{E, *}. First, KT12 associate the previous estimates of τ_IR with κ(T_mp)Σ_*, where T_mp is the temperature at the base of the simulation domain. For T3_F0.5, this quantity is a factor of ∼1.5 larger than τ_V, which we associate with τ_IR.

We can gain insight into the scalings of τ_V, τ_F, and η as functions of τ_* and f_{E, *} by assuming that all simulations will approach a characteristic horizontally average vertical profile. At late times our simulations can be fit approximately by⁹

$$\begin{equation} \overline{T}^4=T_*^4 3 \left(\frac{2}{3} +\overline{\tau }_{\rm F} \right), \end{equation} \tag{ 24 }$$

where quantities with an overbar correspond to y-dependent (2D) or z-dependent (3D), horizontally averaged quantities. For example, $\overline{\rho } = L_x^{-1} \int \rho \; dx$ in 2D and $\overline{\rho } = (L_x L_y)^{-1} \int \int \rho \; dx dy$ in 3D. The exceptions are

$$\begin{eqnarray} \overline{\tau }_{\rm F}(y) & \equiv & \int ^{L_y}_y \frac{\overline{\kappa _{\rm R} \rho F_{ry}}}{\overline{F_{ry}}} \; dy^{\prime }, \nonumber \\ \overline{\tau }_{\rm V}(y) & \equiv &\int ^{L_y}_y \overline{\kappa _{\rm R} \rho } \; dy^{\prime }, \end{eqnarray} \tag{ 25 }$$

with y replaced by z in the 3D case.

To a reasonable approximation, we find that

$$\begin{equation} \overline{\tau }_{\rm F}\simeq \left\lbrace \begin{array}{@{}lr}\overline{\tau }_{\rm V}, & \; \overline{\tau }_{\rm V} < 1, \\ \eta \overline{\tau }_{\rm V}, & \; \overline{\tau }_{\rm V} \gtrsim 1, \end{array}\right. \end{equation} \tag{ 26 }$$

with a constant η ≃ τ_F/τ_V. Using the fact that

$$\begin{equation} \frac{\partial \overline{\tau }}{\partial \Sigma } = \kappa _{\rm R,*} \left(\frac{\overline{T}}{T_*}\right)^2, \end{equation} \tag{ 27 }$$

with $d \Sigma = - \overline{\rho } dy$, we find

$$\begin{equation} \overline{\tau }_{\rm V} = \frac{3}{4}\eta (\kappa _{\rm R,*} \Sigma )^2 +\sqrt{2} \kappa _{\rm R,*} \Sigma. \end{equation} \tag{ 28 }$$

By definition, $\overline{\tau }_{\rm V}(0)=\tau _{\rm V}$ and Σ(0) = Σ_* = τ_*/κ_{R, *} so

$$\begin{equation} \tau _{\rm V} = \frac{3}{4}\eta \tau _*^2 +\sqrt{2} \tau _*. \end{equation} \tag{ 29 }$$

If we assume that all simulations yield f_{E, V} ≈ 1, we can infer that

$$\begin{equation} \tau _{\rm F} \approx \frac{\tau _*}{f_{\rm E,*}}=\frac{c g}{F_*} \Sigma _*, \end{equation} \tag{ 30 }$$

which is equivalent to KT12's Equation (43) with 〈f_E〉 = 1. Thus, for a given input flux F_* and constant gravitational acceleration g, we can define a characteristic opacity

$$\begin{equation} \kappa _{\rm E}=c g/F_*, \end{equation} \tag{ 31 }$$

which corresponds to the opacity where the Eddington ratio is unity. Then the optical depth τ_F = κ_EΣ_* can be used in place of ητ_IR in Equation (23). In our calculations, κ_E is only slightly larger than κ_R near the photosphere for most of our simulations, which is in approximate agreement with conclusion of Krumholz & Thompson (2013) that the photospheric opacity is closer to the relevant opacity than the midplane value for evaluating momentum feedback. However, κ_E could, in principle, significantly exceed the photospheric opacity for systems with lower f_{E, *} and higher τ_*.

We can estimate for η = τ_F/τ_V by combining Equations (29) and (30) to obtain

$$\begin{equation} \eta = \frac{2}{3\tau _*} \left(\sqrt{2+3\tau _*/f_{\rm E,*}}-\sqrt{2} \right), \end{equation} \tag{ 32 }$$

Finally, we can use Equation (32) in Equation (29) to obtain

$$\begin{equation} \tau _{\rm V} = \frac{\tau _*}{\sqrt{2}} \left[1+\left(1 + \frac{3 \tau _*}{2f_{\rm E,*}} \right)^{1/2}\right]. \end{equation} \tag{ 33 }$$

Table 2 summarizes the time-averaged values for 〈τ_F〉, 〈τ_V〉, and 〈η〉 for runs T1_F0.5, T3_F0.5, and T10_F0.5. The quantities in parentheses are the estimates provided by Equations (30), (32), and (33). The estimates prove to be fairly accurate, particularly for the τ_* = 3, 10 runs. For T1_F0.5, which fails to maintain f_{E, V} ≳ 1, we see that both η and τ_F are slightly low, as anticipated. Therefore, we expect the estimates provided here to be useful predictions of how momentum feedback will scale with τ_* and f_{E, *} (or, alternatively, Σ_* and F_*, given g), as long as the system reaches an equilibrium with an Eddington ratio near unity.

The upshot of this analysis is that κ_E would be the most suitable opacity to use in determining the trapping factor if the system being modeled has values of T_* and Σ_* that allow κ_R to consistently reach or exceed κ_E at the base of the domain. A seemingly necessary (but not sufficient) condition for this is that τ_V > τ_F. Combining Equations (30) and (33), this becomes

$$\begin{equation} \tau _* > \frac{4}{3}\left(\frac{1}{f_{\rm E,*}}-\sqrt{2} \right). \end{equation} \tag{ 34 }$$

Note that all of the simulations except T10_F0.02 obey this constraint, but T1_F0.5 has τ_V barely exceeding τ_F and τ_F ≲ κ_EΣ_* so it would seem this relation is only an approximate guide.

For our simulations with f_{E, *} = 0.5, κ_E corresponds to ∼4 cm² g⁻¹. More generally, one can estimate characteristic values for κ_E using Equation (31) with g = 2πGΣ_*/f_gas and $F_*=\epsilon c^2 \dot{\Sigma }_*$ to obtain

$$\begin{equation} \kappa _{\rm E} \simeq \frac{2 \pi G \Sigma _*}{c f_{\rm gas} \epsilon \dot{\Sigma }_*}, \end{equation} \tag{ 35 }$$

where is the fraction of rest mass energy converted into radiation by stars, $\dot{\Sigma }_*$ is the star-formation rate per unit area, and f_gas is the gas fraction of the total (stellar plus gas) mass. Scaling this relation in terms of characteristic values: ₋₃ = /10⁻³, Σ₀ = Σ_*/(1 g cm⁻²), $\dot{\Sigma }_{*,3}=\dot{\Sigma }_*/ (10^3 \, {\rm M_\odot \, yr}^{-1} \, {\rm kpc}^{-2})$, we find

$$\begin{equation} \kappa _{\rm E} \simeq 2.1 \; {\rm cm^2 \; g^{-1}} \epsilon ^{-1}_{-3} \dot{\Sigma }_{*,3}^{-1} f_{\rm gas}^{-1} \Sigma _0. \end{equation} \tag{ 36 }$$

Given the uncertainties and approximations involved, this estimate is in rough agreement with the values of opacity assumed by previous authors (e.g., 5 cm² g⁻¹; Hopkins et al. 2011).

Finally, we note that all of this analysis and discussion assumes that only radiation pressure and gravitational forces play a role. In real systems, mechanical feedback from stellar winds, supernovae, cosmic rays, magnetic fields, photoionization and photoelectric heating, etc. 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 f_{E, V} toward a value ∼1 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.