© ESO, 2012

1. Introduction

Outflows are a unique characteristic of the star formation process in general. The strong radiation pressure provided by massive luminous stars plays an important role in the dynamics of these outflows. Observations show that massive star-forming regions involve several large-scale outflows and their interaction with the interstellar material depicts an important part of the complex morphology of the cluster center (e.g. Shepherd & Churchwell 1996; Zhang et al. 2001; Beltran et al. 2006; Zhang et al. 2007; Fallscheer et al. 2009; Beuther et al. 2010; Quanz et al. 2010). Non-radiative MHD simulations (Banerjee & Pudritz 2006, 2007; Hennebelle et al. 2011) have shown that cavities of lower density potentially form before the onset of star formation during the collapse of magnetized clouds. The treatment of the radiation field is the next step in investigating – especially the thermodynamics of – the cavities.

From the perspective of numerical development, the implementation and improvement of radiation transport schemes in (magneto-)hydrodynamics simulations is both a challenging and timely problem in modern astrophysics, especially in the field of (massive) star formation, in which the radiative feedback onto the environment plays a crucial role (Yorke & Sonnhalter 2002; Krumholz et al. 2007, 2009; Price & Bate 2009; Commercon et al. 2010; Bate 2010; Kuiper et al. 2010a, 2011).

The hybrid radiation transport scheme applied in this paper was previously used in our studies of the radiation pressure barrier in the formation of massive stars (Kuiper et al. 2010a) and the angular momentum transport in disks around massive stars (Kuiper et al. 2011). We introduced the derivation of the solving algorithm as well as its numerical implementation into a three-dimensional (3D) (magneto-)hydrodynamics code in Kuiper et al. (2010b). The solver algorithm superposes two radiation fields: a frequency-dependent ray-tracing (RT) component representing the stellar irradiation field and a frequency-averaged (gray) flux-limited diffusion (FLD) component representing the thermal dust emission. These splitting schemes were previously used in one-dimensional (1D) simulations (Wolfire & Cassinelli 1986, 1987; Edgar & Clarke 2003) and in two-dimensional (2D) hydrostatic disk atmosphere models (Murray et al. 1994) .

The RT methods are commonly used in radiative transfer models without hydrodynamical motion (e.g. Efstathiou & Rowan-Robinson 1990; Steinacker et al. 2003). Alternative methods for these kinds of radiative transfer models include the short characteristics method (Dullemond & Turolla 2000) and Monte Carlo radiative transfer (Bjorkman & Wood 2001; Dullemond & Dominik 2004). Comparison benchmark studies of these methods were presented in Pascucci et al. (2004) and Pinte et al. (2009). These methods are able to solve the radiation transport problem to high accuracy for a fixed static configuration, but they require much CPU time, which yields low efficiency in time-dependent hydrodynamics simulations.

In this sense, the FLD approximation (Kley 1989; Bodenheimer et al. 1990; Klahr et al. 1999; Klahr & Kley 2006) provides a fast method to determine the temperature evolution accurately in the optically thick (diffusion) and optically thin (free-floating) limit. However, the approximation becomes inaccurate in multi-dimensional anisotropic problems as well as in transition regions from optically thin to thick (shown e.g. in Boley et al. 2007). One example of this transition region is the cavity shell investigated in this paper. The FLD approximation in the frequency-averaged (gray) limit remains the default technique in modern multi-dimensional radiation hydrodynamics. In an exceptional case, a frequency-dependent FLD solver was used in Yorke & Sonnhalter (2002).

From the theoretical point of view, the stability of radiation-pressure-dominated cavities has a direct impact on the star formation efficiency and accretion onto a luminous massive star. Krumholz et al. (2009) proposed that massive stars can grow beyond their Eddington limit because these radiation-pressure-dominated cavity shells are subject to the so-called “radiative Rayleigh-Taylor instability” allowing additional mass to be fed onto the central accretion disk. Following the approach of Nakano (1989) and Yorke & Sonnhalter (2002), we proposed the feeding of a massive star beyond its Eddington limit by disk accretion only (Kuiper et al. 2010a, 2011). While the anisotropy of the thermal radiation field sufficiently diminishes the radiation pressure onto the disk accretion flow near the midplane, the radiation pressure in the polar direction is able to launch an outflow, forming a stable cavity.

The removal of a substantial fraction of the initial core mass by the radiative launching of outflows influences the star formation efficiency in the central core region. These feedback effects are proposed to play an important role in explaining the low star-formation efficiencies observed (McKee & Ostriker 2007). Simulations of HII regions in massive star forming regions by Dale & Bonnell (2011) and close to massive stars forming in the cluster center by Peters et al. (2010) suggest that ionized gas fills preexisting voids and bubbles. Long-lived cleared polar cavities therefore would make the expansion of HII regions and their interaction with the stellar cluster environment easier. The feedback of the energy injection from jets and outflows onto the stellar cluster formation is briefly discussed in Bonnell & Bate (2006).

In this study, we investigate different implementation methods for direct stellar irradiation feedback (see Sect. 3) and their influence on the stability of radiation pressure dominated cavities. We present in Sect. 4 the qualitative outcome and Sects. 5 and 6 the quantitative analyses of the simulations. The observed difference in the radiative acceleration of the cavity shell depending on the applied radiation transport method is analytically derived in Sect. 7. Finally, a comprehensive comparison section (Sect. 8) to previous simulations, analytic work, and observations as well as a brief summary (Sect. 9) are provided.

2. Physical and numerical model

2.1. Physics included

We compare the effect of two different radiation transport methods on stellar radiation feedback. The goal is to investigate the stability of the shell surrounding the radiation-pressure-dominated cavity. We emphasize that in these simulations we do not aim to provide a complete description of an outflow region around a massive star.

A list of potential limitations and caveats includes:

• 1.

  The collimation effects by magnetic fields – which are certainly important for the morphology of the cavity – are neglected. Recent results of MHD simulations of jet formation are summarized by Fendt et al. (2011). Potentially, the relative importance of radiative or magnetic forces in jet formation can be distinguished observationally based on morphology: a magnetically launched fast jet would be more collimated than a radiation-pressure-driven wide-angle outflow (see e.g. Arce et al. 2007; Pudritz et al. 2007, and citations therein).

• 2.

  Another example of the influence of magnetic fields, the so-called photon bubble instability (Turner et al. 2007) could potentially diminish the radiation pressure and could therefore alter the morphology of a cavity filled with magnetized gas.

• 3.

  Evaporation of dust is included, but otherwise the dust is assumed to be perfectly coupled to the gas.

• 4.

  The dust opacity clearly dominates the absorption in any dusty region, but in the dust-free zone around the massive star up to the dust sublimation front the (remaining) gas opacity is set to be constant κ_gas = 0.01cm²   g^-1. Such a constant gas opacity of course denotes only a first order approximation to the complex molecular line absorption features. An optically thick gas disk could e.g. amplify the flashlight effect (Tanaka & Nakamoto 2011) and therefore increase the radiative flux into the polar regions.

• 5.

  Since the cavity growth is driven by radiation pressure, the quantitative details of the launching process and the growth phase depend on the stellar evolution model as well as the dust model. Future studies will attempt to establish the details of this dependence.

• 6.

  The low-density cavity is potentially enriched, but certainly influenced, by a stellar wind from the central massive star, which is neglected herein. Early proto-stellar winds could, e.g., lead to an early formation of rather collimated polar cavities (Cunningham et al. 2011). The radiation pressure at later epochs in such a configuration could be channeled into pre-existing cavities, leading to smaller opening angles for the radiation pressure dominated cavities.

• 7.

  In Peters et al. (2010, 2011), ionizing radiation was proposed to dominate the cavity dynamics in regions of massive star formation.

2.2. Numerical resolution

To obtain unambiguous results, the relevant physical length scales have to be resolved by the numerical simulations. Therefore, we discuss these scales with respect to the resolution of our radiation hydrodynamical simulations. We demonstrate that the simulations resolve all relevant length scales sufficiently.

A study of a radiative Rayleigh-Taylor instability in the shell surrounding the outflow cavity needs to resolve the wavelengths of this instability along the shell discontinuity. At small wavelengths, the radiative Rayleigh-Taylor instability is likely to be suppressed by diffusion, but at large wavelengths – comparable to the physical size of the cavity – the instability occurs. Quantitatively, the shortest unstable wavelength for instance is determined to be 1000 AU for a 10 M_⊙ and 10 000 AU for a 100 M_⊙ star (Jacquet & Krumholz 2011). The polar resolution rΔθ of the spherical grid in our simulations grows linearly with the radius and is fixed to rΔθ ≲ 0.1r. The grid size along the shell discontinuity is typically a factor of ten smaller than the shortest unstable wavelength. Hence, the length scale of the radiative Rayleigh-Taylor instability is fully resolved.

Two other important length scales are related to the radiative properties of the shell. The “cooling length scale” l_c is given by an optical depth of τ = 1 in the long wavelength regime, and the “irradiation length scale” l_∗ is given by an optical depth of τ = 1 of the broad stellar irradiation spectrum. Figure 1 shows the cooling length scale in the gray approximation as a function of the location of the cavity shell R_cavity. The shell density is chosen to be as depicted in Fig. 5. The opacities are computed as the Rosseland mean opacities of Laor & Draine (1993). The stellar temperature is taken from the Hosokawa & Omukai (2009) tracks for a 20 M_⊙ star and an accretion rate of Ṁ = 10^-3   M_⊙yr^-1, and the temperature at the location R_cavity is computed with a slope $\mathit{T}\mathrm{\propto }{\mathit{R}}_{\mathrm{cavity}}^{-0.5}$ in the gray approximation. For simplification, the evaporation of dust grains is neglected, but would result in an even larger cooling length-scale. Hence, we can be sure that the length scale of cooling is resolved.

The absorption length scale l_∗ of the stellar irradiation spectrum is much smaller than the cooling length scale for thermal dust emission. Using again the opacities of Laor & Draine (1993), Fig. 2 shows the absorption length scale of the whole stellar irradiation spectrum for different gas densities. For the aforementioned shell density of and a typical stellar temperature at the interesting point in time at the potential onset of instability of about 10   000 K  < T_∗ < 30   000 K, the absorption length scale of the photons with frequencies higher than the peak of the stellar black-body spectrum is about 10 to 100 AU; these photons are most likely to be absorbed in the first grid cell on top of the cleared cavity. For lower frequency photons, the absorption length scale smoothly rises up to the 1000 AU scale and is most likely to be resolved. In the long wavelength regime of the spectrum, the results of the previous paragraph for the thermal dust emission of course apply.

The morphology in the radial direction during the cavity growth phase can be distinguished into three important regimes: the inner low-density cleared cavity, the shell of swept-up material on top of this cavity, and the remnant of the in-falling envelope on larger scales. The resolved transition of the gas density, temperature, and velocity in these regimes is depicted in the context of the following analyses in Figs. 3 (two-dimensional snapshots) and 5 (along the polar axis). Even the steep gradient of the density at the inner rim of the cavity shell (Fig. 5, upper left panel) is resolved by four grid cells, and the smooth decrease to larger radii by dozens of grid cells. Therefore, the shell morphology is clearly resolved in the simulations.

2.3. Axial symmetry

The simulations in this paper are performed assuming axial symmetry. We are convinced that this is not a critical limitation of the result we present. In 2D simulations with the same initial conditions, but different methods for stellar radiation feedback, we observe strong differences in the resulting radiation feedback onto the shell on top of the cleared cavity. These differences do not depend on the detailed morphology of this shell and therefore the conclusion is also applicable to non-axially symmetric, 3D configurations.

In our 2D simulation results with the FLD approximation, the cavity shell undergoes an instability; this result matches the 3D simulation result of Krumholz et al. (2009). In our 2D simulations as well as in our 3D simulation (Kuiper et al. 2011), including the RT step for direct irradiation feedback, the radiation-pressure-dominated cavity remains stable.

3. Method

3.1. Equations

To follow the motion of the gas, we solve the equations of compressible self-gravity hydrodynamics, including shear-viscosity as described in Kuiper et al. (2010a). For this task, we use the open source MHD code Pluto (Mignone et al. 2007). This set of equations is coupled to the radiation transport in the medium. In contrast to our previous studies, we investigate the influence of two different methods to determine the radiative flux, namely our hybrid radiation transport scheme (ray-tracing + flux-limited diffusion) and flux-limited diffusion only.

3.1.1. Ray-tracing + flux-limited diffusion

For simulations labeled “RT+FLD”, we use our hybrid radiation transport module (see Kuiper et al. 2010b). The total radiation energy density F_tot is divided into the flux F_∗(ν,r) of the frequency-dependent stellar irradiation and the flux of the frequency-averaged thermal radiation energy density F$\begin{array}{ccc}{\mathit{\partial }}_{\mathit{t}}{\mathit{E}}_{\mathrm{R}}\mathrm{+}{\mathit{f}}_{\mathrm{c}}\mathbf{\nabla }\mathrm{·}{F}& \mathrm{=}& \mathrm{-}{\mathit{f}}_{\mathrm{c}}\left(\mathbf{\nabla }\mathrm{·}{F}\mathrm{\ast }\mathrm{-}{\mathit{Q}}^{\mathrm{+}}\right)\mathit{,}\\ {F}\mathrm{\ast }\mathrm{(}\mathit{\nu ,}{\mathit{r}}^{\mathrm{)}}\mathit{/}{F}\mathrm{\ast }\mathrm{(}\mathit{\nu ,}{\mathit{R}}_{\mathrm{\ast }}\mathrm{)}& \mathrm{=}& {\mathrm{(}{\mathit{R}}_{\mathrm{\ast }}\mathit{/}{\mathit{r}}^{\mathrm{)}}}^{\mathrm{2}}\exp \mathrm{(}\mathrm{-}\mathit{\tau }{\mathrm{(}\mathit{\nu ,}{\mathit{r}}^{\mathrm{)}}}^{\mathrm{)}}\mathit{,}\end{array}$where Eq. (1)denotes the evolution of the thermal radiation energy density E_R. The factor ${\mathit{f}}_{\mathrm{c}}\mathrm{=}{{}^{\mathrm{(}}\mathit{c}_{\mathrm{V}}\mathit{\rho }\mathit{/}\mathrm{4}\mathit{a}{\mathit{T}}^{\mathrm{3}}\mathrm{+}{\mathrm{1}}^{\mathrm{)}}}^{-1}$ depends only on the ratio of internal to radiation energy and contains the specific heat capacity, c_V, and the radiation constant, a. The source term Q ⁺  depends on the physics included and contains any additional energy source such as hydrodynamical compression  − P∇·u and viscous heating. We solve Eq. (1)using the so-called flux-limited diffusion approximation, in which the flux is set proportional to the gradient of the radiation energy density (F =  −D∇E_R). The diffusion constant D = λc/ρκ_R depends on the flux limiter λ and the Rosseland mean opacity κ_R. The quantity c denotes the speed of light in a vacuum. We use the flux limiter proposed by Levermore & Pomraning (1981) and neglect scattering.

Equation (2)calculates the flux of the frequency-dependent stellar irradiation in a ray-tracing step. The first factor on the right hand side describes the geometrical decrease in the flux proportional to r^-2. The second factor describes the absorption of the stellar light as a function of the optical depth τ(ν,r) = κ(ν)ρ(r)r depending on the frequency-dependent mass absorption coefficients κ(ν). For this purpose, we use tabulated dust opacities by Laor & Draine (1993), including 79 frequency bins, and calculate the local evaporation temperature of the dust grains by using the formula of Isella & Natta (2005). The flux at the inner radial boundary is given by the luminosity L_∗, temperature T_∗, and radius R_∗ of the forming star. For this purpose, we use tabulated stellar evolutionary tracks for accreting high-mass stars, calculated by Hosokawa & Omukai (2009). The gas and dust temperature T is finally calculated in equilibrium with the combined stellar irradiation and thermal radiation field $\mathit{a}{\mathit{T}}^{\mathrm{4}}\mathrm{=}{\mathit{E}}_{\mathrm{R}}\mathrm{+}\frac{\mathit{\kappa }\mathrm{(}{\mathit{\nu }}^{\mathrm{)}}}{{\mathit{\kappa }}_{\mathrm{P}}\mathrm{\left(}\mathit{T}\mathrm{\right)}}\frac{\mathrm{|}{F}\mathrm{\ast }\mathrm{|}}{\mathit{c}}$(3)with the Planck mean opacities κ_P.

Numerical details, test cases, including a comparison of gray and frequency-dependent irradiation, as well as performance studies of the hybrid radiation transport scheme are summarized in Kuiper et al. (2010b). The viscosity prescription as well as the tabulated dust and stellar evolution model are presented in Kuiper et al. (2010a).

3.1.2. Flux-limited diffusion only

For simulations labeled “FLD”, we do not ray-trace the stellar irradiation, but add the luminosity of the forming star as a source term to the right hand side of the FLD equation $\begin{array}{ccc}{\mathit{\partial }}_{\mathit{t}}{\mathit{E}}_{\mathrm{R}}\mathrm{+}{\mathit{f}}_{\mathrm{c}}\mathbf{\nabla }\mathrm{·}{F}& \mathrm{=}& {\mathit{f}}_{\mathrm{c}}\left({\mathit{L}}_{\mathrm{\ast }}\mathit{\delta }\mathrm{\left(}{r}\mathrm{\right)}\mathrm{+}{\mathit{Q}}^{\mathrm{+}}\right)\mathit{.}\end{array}$(4)Since the forming star is enclosed in the inner sink cell in our grid in spherical coordinates, the luminosity of the central star is implemented as a Dirichlet boundary condition at the inner rim of the computational domain. The radiation energy at this inner boundary is therefore given by the stellar surface temperature and the assumption that the medium between the stellar surface R_∗ and the inner rim at r_min = 10 AU is optically thin.

The hybrid radiation transport scheme (Sect. 3.1.1) as well as the FLD approach were tested in detail in comparison with a Monte Carlo solution using the RADMC code (Dullemond & Dominik 2004) or a short characteristics method using the RADICAL code (Dullemond & Turolla 2000) in a setup including a central star, an accretion disk, and an envelope as described in Pascucci et al. (2004) and Pinte et al. (2009).

3.2. Numerical configuration

The simulations are performed on a time-independent grid in spherical coordinates with a logarithmically stretched radial coordinate axis. The outer core radius is fixed to r_max = 0.1 pc, and the inner core radius is fixed to r_min = 10 AU. The accurate size of this inner sink cell was determined in a parameter scan presented in Kuiper et al. (2010a), Sect. 5.1. The polar angle extends from 0° to 90° assuming midplane symmetry. The grid consists of 64    ×    16 grid cells, i.e. the highest resolution of the non-uniform grid is chosen to be $\mathrm{\Delta }\mathit{r}\mathrm{\times }\mathit{r}\mathrm{\Delta }\mathit{\theta }\mathrm{=}\mathrm{1.27}\mathrm{AU}\mathrm{\times }\mathrm{1.04}\mathrm{AU}$(5)around the forming massive star. The resolution decreases logarithmically in the radial outward direction proportional to the radius. The radially inner and outer boundary of the computational domain are semi-permeable walls, i.e. the gas is allowed to leave the computational domain (by accretion onto the central star or due to radiative and centrifugal forces over the outer boundary), but cannot enter this domain. This outer boundary condition allows the mass reservoir for stellar accretion to be controlled by the initial choice of the mass of the pre-stellar core.

The Pluto code uses high-order Godunov solver methods to compute the hydrodynamics, i.e. it uses a shock capturing Riemann solver within a conservative finite volume scheme. Our default configuration consists of a Harten-Lax-Van Leer Riemann solver and a so-called “minmod” flux limiter using piecewise linear interpolation and a Runge-Kutta 2 time-integration, also known as the predictor-corrector-method; for comparison we also refer to van Leer (1979). Therefore, the total difference scheme is accurate to second order in time and space.

The internal iterations of the implicit solver of the FLD in Eqs. (1)or (4)is stopped at an accuracy of the resulting temperature distribution of either ΔT/T ≤ 10^-3 or ΔT ≤ 0.1 K. The internal iterations of the implicit solver for Poisson’s equation is stopped at an accuracy of the resulting gravitational potential of ΔΦ/Φ ≤ 10^-5.

3.3. Initial conditions

We start from a cold (T₀ = 20K) pre-stellar core of gas and dust. The initially constant dust-to-gas mass ratio is chosen to be M_dust/M_gas = 0.01. Dynamically, the model describes a so-called quiescent collapse scenario without initial turbulent motion (u_r = u_θ = 0) and the core is initially in slow solid-body rotation ${\mathrm{(}}^{}\mathrm{\left|}{u}{\mathit{\phi }}_{}\mathrm{\right|}\mathit{/}\mathit{R}\mathrm{=}{\mathrm{\Omega }}_{\mathrm{0}}\mathrm{=}\mathrm{5}\mathrm{\ast }{\mathrm{10}}^{-13}\mathrm{Hz}{\mathrm{)}}^{}$. The total mass M_core in the computational domain is chosen to be 50 or 100 M_⊙. In addition to the influence of the different radiation transport methods described in the last sections, we check the dependence of the stability of the radiation pressure dominated cavity on the initial density slope of the pre-stellar core. We chose the initial density slope to be either ρ ∝ r^-1.5 or ρ ∝ r^-2. An overview of the runs is given in Table 1. For a more comprehensive parameter scan of the initial conditions of this pre-stellar core collapse model, we refer to Kuiper et al. (2010a, 2011).

4. Qualitative simulation results

We summarize the qualitative outcome, the stability, and the overall morphology of the radiation pressure dominated cavities, depending on the different initial conditions and the radiation transfer method. In the six simulations, the radiation pressure is generally high enough to launch an outflow and form a cleared polar cavity.

In the case of the RT+FLD radiation transport method, these outflow cavities, once launched, rapidly grow in their extent up to the outer computational domain at 0.1 pc away from the central star. In contrast to this monotonic growth, the outflow cavities in the FLD radiation transport scheme stop their increase along the polar axis at an extent of the order of several 100 to 1400 AU depending on the initial conditions. This stopping of their growth is followed by a penetration of the gas mass on top of the outflow cavity formed so far, i.e. the outflow cavity collapses. During a follow-up second and third epoch of outflow launching, the resulting cavities grow in their extent, but on much longer timescales than in the RT+FLD case. The epochs of frozen cavity growth in the FLD-only simulations suggest that (parts of) the top layer of the outflow cavity are in marginal Eddington equilibrium, i.e. the radiation pressure force is in equilibrium with the stellar gravity. The qualitative outcome of the simulations (depending on whether an extended epoch of marginal Eddington equilibrium occurs) is summarized in Table 1, Col. 5.

In addition to the different growth rates, the radiation-pressure-driven outflow cavities also develop pronounced differences in their morphology depending on the applied radiation transport method. As an example, the different morphologies produced by the different radiation transport methods are visualized in three snapshots in time during the onset of the outflow launching in Fig. 3. In the case of FLD, the outflow is easily stopped by the infalling matter along the polar axis, while in the RT+FLD case the outflow cavity grows rapidly with time. In the FLD approximation, the radiative flux tends to follow a path that minimizes the optical depth and hence avoids the swept-up mass on top of the cavity. This avoidance is alleviated by the centrifugal forces, which diminish the gravitational attraction in regions slightly above the disk. In the ray-tracing method, the isotropic stellar irradiation flux directly impinges onto the swept-up mass on top of the polar cavity and pushes the mass to larger radii. The resulting large-scale morphology of the cavity is far more isotropic. The opening angle of the cavity is determined by the inner disk structure.

In the following Sects. 5 and 6, we analyze quantitatively the outcome of the simulations presented here. We check our hypothesis of epochs of marginal Eddington equilibrium in the FLD-only runs and determine via analytical estimates (Sect. 7), why the cavity shells in the RT+FLD case do not undergo these epochs and therefore remain stable with respect to the radiative Rayleigh-Taylor instability.

5. Quantitative analysis of the cavity growth

We analyze quantitatively the time-dependent extent of the outflow cavity. The radial extent R_cavity of the outflow cavities above the central star is determined as the extent of the cleared cavity along the polar axis, as visualized in Fig. 3. The resulting cavity radii as a function of time are shown for the three different initial conditions as well as the different radiation transport methods in Fig. 4.

For all initial conditions, the extent of the radiation pressure dominated cavities increases far more rapidly, if the stellar source of radiation is treated via a ray-tracing scheme (labeled “RT+FLD”) than simply included in the flux-limited diffusion solver (labeled “FLD”). If using the FLD approximation, the outflow is launched a little bit earlier in time and after its initial growth phase undergoes an epoch of marginal stability, in which the radiation pressure force along the polar axis seems to be balanced by gravity; the cavity growth along the polar axis is first stopped and then reversed. Under these conditions of marginal Eddington equilibrium ${}^{\mathrm{(}}\mathit{f}_{\mathrm{rad}}\lesssim {\mathit{f}}_{\mathrm{grav}}{\mathrm{)}}^{}$ with FLD, the cavity shell region has been proposed to be subject to the radiative Rayleigh-Taylor instability (Jacquet & Krumholz 2011). As a consequence, subsequent growth phases of the outflows in the FLD runs depend on the details of the subgrid models; the dynamics of the marginally balanced cavity shells are strongly influenced by the stellar evolution model and the dust model, particularly the treatment of sublimation and evaporation.

In contrast to this epoch of marginal Eddington equilibrium, the extent of the outflow cavity in the RT + FLD simulations increases rapidly in time. With the exception of the M_core = 50   M_⊙ case, in which a correspondingly lower mass star (yielding much lower luminosity) is formed, the outflow cavity increases in size monotonically in the RT + FLD simulations.

The difference in the growth rate for the various radiation transport schemes is most prominent in the M_core = 100   M_⊙ and ρ ∝ r^-2 case (Fig. 4 left panel), i.e. in the case of the initially highest mass in the central core region, allowing for a rapid formation and growth of the central star. In the M_core = 50   M_⊙ case (Fig. 4 middle panel), the much lower stellar luminosity results in an initially unstable cavity region even in the RT+FLD case (but with a smaller extent than 100 AU). In the ρ ∝ r^-1.5 case (Fig. 4 right panel), the mass is initially more concentrated at larger radii and hence the cavity shell has to sweep up far more mass during its increase. Therefore, the growth phase of the cavity takes much longer than in the cases of initially steeper density profiles.

6. Quantitative analysis of the cavity shell morphology and dynamics

To unveil the physical background of the qualitative and quantitative difference of the outflow cavity structure in the simulations using either the RT+FLD or the FLD method, we now investigate the morphology and dynamics of the cavity shell depending on the radiation transport method. In this section, we analyze quantitatively the difference in the morphology and the dynamics of the cavity shell in the RT+FLD and FLD cases. We focus on the data of the simulation with an initial core mass M_core = 100   M_⊙ and a radial density slope of β =  −2.

The cavity in the FLD case increases in its first growth phase up to roughly 1400 AU before the expansion along the polar axis stops and the mass flow is reversed by gravity. We analyze the gas dynamics at the point in time when in both simulations (with RT+FLD and with FLD only) the expansion of the outflow cavity has arrived at the same location, namely at t = 16.7 kyr. The gas density, temperature, radial velocity, and the radial mass loss rate along the polar axis at this time are shown in Fig. 5. The upper left panel of this figure shows that there is almost no difference in the radial extent, compactness, and peak density of the swept-up mass in the cavity shell. The peaks in the density distributions correspond to and in the RT+FLD and FLD cases, respectively.

The upper right panel of Fig. 5 shows that the RT+FLD run results in a continuously slightly higher temperature in the cleared cavity, the cavity shell, and beyond. The temperature distribution in the pure FLD run stays 24% to 45% below the temperature distribution of the RT+FLD run.

In contrast to these similarities, the two lower panels of Fig. 5 highlight the striking difference of both runs in their radial velocities and mass flux along the polar axis. In the RT+FLD run, the swept-up shell material moves into the interstellar medium with a speed slightly higher than 100 km s^-1, which is roughly a factor of three higher than in the corresponding FLD run. At larger radii (r > 3000 AU), the velocity slope is still dominated by gravitational infall and hence is independent of the radiation transport method.

Up to the onset of the Eddington equilibrium epoch in the FLD run, the mass of the central star in both simulations (RT+FLD and FLD only) closely match, leading to the same strength of gravitational attraction for the matter in the cavity and shell towards the direction of the star. Furthermore, the density structure is roughly the same (Fig. 5, upper left panel) and the deviation in temperature within the cavity is smaller than 24–45% (Fig. 5, upper right panel), hence the slightly higher thermal pressure cannot be the main driver of the enormous difference in the radial velocities and mass flux (Fig. 5, lower panels).

As a consequence, the radiation pressure acting on the swept-up mass on top of the outflow cavity has to differ in the two radiation transport methods. First, the FLD approximation assumes that the dust grains of the swept-up mass on top of the cavity are embedded in a local radiation field, whereas the RT+FLD method takes into account the (frequency-dependent) absorption probability caused by direct stellar irradiation. The resulting difference in the radiative acceleration at the top of the outflow cavity is analytically estimated in the following section.

Secondly, since the FLD approximation is mathematically a moment method, the derivation of the FLD approximation includes the integral over the angular distribution of the radiative flux; hence the emitted photons of the central star do not move along straight rays until they are absorbed. The radiative flux in the FLD approximation instead follows a path that minimizes the optical depth. The FLD method was originally introduced to describe spherically symmetric flows, hence the integral over the propagation direction did not result in unphysical behavior (see e.g. Bruenn 1985). In a multi-dimensional environment with a non-isotropic optical depth, this assumption of the FLD approximation breaks down. Including higher-order terms in the derivation would minimize this inaccuracy and potentially yield a sufficient tracing of the correct photon path. However, in the RT+FLD scheme the irradiation by the central star is computed via a tay-tracing equation, which takes into account the correct photon propagation direction, i.e. the stellar irradiation flux directly impinges on the swept-up mass shell at the top of the optically thin cavity.

7. Analytical determination of the acceleration caused by stellar irradiation in the ray-tracing approach and the flux-limited diffusion approximation

Our aforementioned analysis indicates that the difference between the cavity shell stability in the two radiation transport methods is due to the higher radiative acceleration of the swept-up material when the RT approach is used for the stellar irradiation instead of the FLD approximation. The use of RT guarantees that the stellar radiation directly interacts with the swept-up mass at the top of the optically thin outflow cavity. The dust grains in the swept-up material therefore absorb the stellar flux following the frequency-dependent absorption coefficient κ(ν). If we integrate over frequencies, the dust grains absorb the stellar spectrum with respect to the Planck mean opacity κ_P(T_∗) at a temperature T_∗ of the stellar photosphere. In contrast, the FLD approximation assumes that dust grains are embedded in a locally isotropic radiation field and that their absorption behavior is based on the Rosseland mean opacity κ_R(T) at the local radiation temperature T at the top of the cavity. Computing the Planck or Rosseland mean opacity for a given temperature T mostly results in a difference of only a factor of two. However, computing the Planck mean opacity with the stellar temperature T_∗ instead of the local radiation temperature T leads to enormous differences. In the next subsection, we determine and illustrate the difference between these opacities.

7.1. Opacities

For the FLD part of our hybrid radiation transport scheme, we compute the Rosseland mean opacities κ_R(T) directly from the frequency-dependent opacity table (Laor & Draine 1993). To obtain the absorption coefficient in a given grid cell, these opacities are multiplied by the local dust-to-gas mass ratio M_dust/M_gas(x). In our models, the evaporation temperature of dust grains is calculated depending on the local gas density by applying the formula of Isella & Natta (2005), which represents a power-law fit to the Pollack et al. (1994) data. For the opacities in the simulation of Krumholz et al. (2009), the authors use linear regressions to the Pollack et al. (1994) data in a small couple of intervals; the density dependence of the evaporation of dust grains is not taken into account.

The Rosseland κ_R(T) and Planck κ_P(T) mean opacities used in our simulations and those applied in Krumholz et al. (2007), Eq. (11), are visualized in Fig. 6. Owing to the similarities between our and the Krumholz et al. opacities, the resulting radiation pressure in the cavity and shell should be comparable in the runs using the FLD radiation transport method. If there were a difference at all, the higher mean opacities in our simulations for T > 600 K would lead to a slightly higher radiation pressure onto the swept-up mass shell than in the simulations of Krumholz et al. If our opacities in the FLD run lead to an unstable cavity shell, the opacities used by Krumholz et al. have to yield the same outcome.

In our RT+FLD runs, the FLD approximation is only used to determine the thermal radiation by dust grains. In addition, the irradiation of the stellar luminosity up to the first absorption event (i.e. up to an optical depth of τ(ν) ≈ 1) is computed by a RT step that appropriately resembles the propagation and absorption of the stellar photons. The dust grains absorb the star light according to the frequency-dependent absorption coefficient κ(ν). Although this frequency dependence is fully covered in our RT solver (we use 79 frequency bins in the case of the Laor & Draine 1993, opacities), in this derivation – for simplicity – we use the gray approximation by computing the corresponding Planck mean opacities κ_P(T_∗) that depend on the stellar surface temperature T_∗. Fig. 7 shows the resulting Planck mean opacities ${\mathit{\kappa }}_{\mathrm{P}}^{\mathrm{0}}\mathrm{\left(}{\mathit{T}}_{\mathrm{\ast }}\mathrm{\right)}$ as a function of the stellar surface temperature T_∗. To obtain the absorption coefficient κ_P(T_∗) in a given grid cell of the computational domain for these stellar photons, the opacities ${\mathit{\kappa }}_{\mathrm{P}}^{\mathrm{0}}\mathrm{\left(}{\mathit{T}}_{\mathrm{\ast }}\mathrm{\right)}$ are multiplied by the evaporation probability, and, therefore, also depend on the local gas density and temperature.

The mass of the proto-star during the onset of the cavity shell instability in the FLD run is about 20M_⊙ and higher. Comparing the absorption coefficient of the dust grains in the shell on top of the optically thin cavity in the FLD approximation (Fig. 6) with the RT method (Fig. 7), immediately illustrates that the FLD approximation underestimates the absorption coefficient tremendously (at least by a factor of 20–50).

In the following subsections, we show that this underestimated absorption coefficient for direct stellar irradiation in the FLD method leads to epochs of marginal Eddington equilibrium in the cavity shell, allowing the shell to become unstable, whereas the higher-order treatment of stellar irradiation via RT results in a stable cavity shell with radial velocities up to v ≳ 100   km   s^-1 and mass-loss rates of Ṁ ≈ 2 × 10^-4   M_⊙yr^-1.

7.2. Radiative acceleration

To compute the radiative acceleration a_rad of the swept-up material at the top of the optically thin cavity, we derive analytic expressions for the different radiation transport approaches based on formulae similar to the well-known Eddington limit. We start with the formula for radiative acceleration ${a}\mathrm{rad}\mathrm{=}\mathit{\kappa }\frac{{F}}{\mathit{c}}$(6)for the radiative flux F and the speed of light c. Since we only wish to estimate the radiation pressure by the first absorption of the (isotropic) stellar light and the cavity is treated as optically thin up to the location R_cavity of the swept-up material on top of it, we compute the radiative acceleration onto dust grains only along the direction of the polar axis (1D) and express the radiative flux by $\mathit{F}\mathrm{\left(}{\mathit{R}}_{\mathrm{cavity}}\mathrm{\right)}\mathrm{=}\frac{{\mathit{L}}_{\mathrm{\ast }}}{\mathrm{4}\mathit{\pi } {\mathit{R}}_{\mathrm{cavity}}^{\mathrm{2}}}\mathrm{·}$(7)Equation (7) is valid either for the assumption of an optically thin cavity or – owing to the conservation of momentum and energy – as long as the swept-up mass in the cavity shell can be treated as spherically symmetric. By means of the absorption of stellar photons, the dust grains heat up and then re-emit photons at longer wavelengths. Integrating over all photons penetrating a shell at the radius r and the momentum gained by the dust grains ensures that there is an excellent conservation of momentum.

Inserting Eq. (7) into Eq. (6) results in a radiative acceleration that depends on the absorption coefficient κ, the stellar luminosity L_∗, and the cavity height R_cavity above the star ${\mathit{a}}_{\mathrm{rad}}\mathrm{=}\mathit{\kappa }\frac{{\mathit{L}}_{\mathrm{\ast }}}{\mathrm{4}\mathit{\pi } \mathit{c} {\mathit{R}}_{\mathrm{cavity}}^{\mathrm{2}}}\mathrm{·}$(8)This radiative acceleration has to be compared with the gravitational attraction a_grav of the dust grains at the cavity shell position R_cavity of the central stellar mass M_∗${\mathit{a}}_{\mathrm{grav}}\mathrm{=}\frac{\mathit{G} {\mathit{M}}_{\mathrm{\ast }}}{{\mathit{R}}_{\mathrm{cavity}}^{\mathrm{2}}}\mathit{,}$(9)where G is the gravitational constant. The gravity of the less massive circumstellar accretion disk is neglected here. In observations and simulations, even the mass of the accretion disk plus the rotating large-scale torus further away from the cavity shell is less than or comparable to the mass of the central star; hence neglecting this mass reduces the stellar gravity by only a factor of two at most. To minimize the numbers of independent variables, we derive the stellar luminosity L_∗ from the stellar mass M_∗ by assuming that the proto-star follows the evolutionary track of Hosokawa & Omukai (2009) with a constant accretion rate of Ṁ = 10^-3   M_⊙   yr^-1. The assumption of this particular evolutionary track does of course influence the exact quantitative comparison of the radiative versus the gravitational acceleration, but the qualitative outcome does not change. The critical factor in this comparison is the strength of the absorption coefficient κ.

In the FLD approximation, the absorption coefficient κ is given by the Rosseland mean opacity κ_R(T) at the local radiation temperature T of the cavity shell.

In the following, we compute the local gas temperature T at the top of the cavity by the formula of Spitzer (1978) for an optically thin stellar environment with the prerequisite R_cavity ≫ R_∗ that the distance R_cavity to the shell is much larger than the stellar radius R_∗$\mathit{T}\mathrm{\left(}{\mathit{R}}_{\mathrm{cavity}}\mathrm{\right)}\mathrm{=}{\left(\mathrm{0.5} \frac{{\mathit{R}}_{\mathrm{\ast }}}{{\mathit{R}}_{\mathrm{cavity}}}\right)}^{\frac{\mathrm{2}}{\mathrm{4}\mathrm{+}{\mathit{\beta }}_{\mathrm{opacity}}}}\mathrm{\times }{\mathit{T}}_{\mathrm{\ast }}\mathit{,}$(10)where β_opacity represents the exponent of the slope of the frequency-dependent opacities at longer wavelengths. The assumption of gray absorption would, e.g., imply β_opacity = 0 and, therefore, results in a temperature slope of $\mathit{T}\mathrm{\propto }{\mathit{R}}_{\mathrm{cavity}}^{-0.5}$. The slope for the opacities of Laor & Draine (1993) is β_opacity ≈ 2.05.

We can express the radiative acceleration in the FLD case using Eq. 8 with the Rosseland mean opacities of Fig. 6 at the local radiation temperature T(R_cavity) computed via Eq. (10). The radiative acceleration in the RT case is computed by using Eq. (8) with the corresponding Planck mean opacities (including the local dust evaporation probability) of Fig. 7 at the stellar temperature T_∗, which is for a specific stellar mass M_∗ also taken from the evolutionary track of Hosokawa & Omukai (2009). The final formulae for the accelerations are given by $\begin{array}{ccc}{\mathit{a}}_{\mathrm{rad}}^{\mathrm{FLD}}& \mathrm{=}& {\mathit{\kappa }}_{\mathrm{R}}\mathrm{\left(}\mathit{T}\mathrm{\right(}{\mathit{R}}_{\mathrm{cavity}}\mathrm{\left)}\mathrm{\right)}\mathrm{\times }\frac{{\mathit{L}}_{\mathrm{\ast }}\mathrm{\left(}{\mathit{M}}_{\mathrm{\ast }}\mathrm{\right)}}{\mathrm{4}\mathit{\pi } \mathit{c} {\mathit{R}}_{\mathrm{cavity}}^{\mathrm{2}}}\mathit{,}\\ {\mathit{a}}_{\mathrm{rad}}^{\mathrm{RT}}& \mathrm{=}& {\mathit{\kappa }}_{\mathrm{P}}\mathrm{\left(}{\mathit{T}}_{\mathrm{\ast }}\mathrm{\right(}{\mathit{M}}_{\mathrm{\ast }}\mathrm{\left)}\mathrm{\right)}\mathrm{\times }\frac{{\mathit{L}}_{\mathrm{\ast }}\mathrm{\left(}{\mathit{M}}_{\mathrm{\ast }}\mathrm{\right)}}{\mathrm{4}\mathit{\pi } \mathit{c} {\mathit{R}}_{\mathrm{cavity}}^{\mathrm{2}}}\mathit{,}\\ {\mathit{a}}_{\mathrm{grav}}& \mathrm{=}& \frac{\mathit{G} {\mathit{M}}_{\mathrm{\ast }}}{{\mathit{R}}_{\mathrm{cavity}}^{\mathrm{2}}}\mathrm{·}\end{array}$The two remaining independent parameters are the actual stellar mass M_∗ and the location R_cavity of the cavity shell material. In the RT case, the opacity in regions cooler than the dust evaporation temperature is independent of the location R_cavity, i.e. the radiative acceleration scales with $\mathit{a}\mathrm{\propto }{\mathit{R}}_{\mathrm{cavity}}^{-2}$ in the same way as the gravitational acceleration. The ratio ${\mathit{a}}_{\mathrm{rad}}^{\mathrm{RT}}\mathit{/}{\mathit{a}}_{\mathrm{grav}}$ of these accelerations, therefore, is scale-free and depends only on the proto-stellar properties leading to the well-known formula for the generalized Eddington limit $\frac{{\mathit{L}}_{\mathrm{\ast }}}{{\mathit{M}}_{\mathrm{\ast }}}\mathrm{=}\frac{\mathrm{4}\mathit{\pi } \mathit{G} \mathit{c}}{{\mathit{\kappa }}_{\mathrm{P}}\mathrm{\left(}{\mathit{T}}_{\mathrm{\ast }}\mathrm{\right)}}\mathrm{·}$(14)

7.3. Results

In the following, we compare the accelerations depending on the remaining variables M_∗ and R_cavity. Therefore, Fig. 8 compares the radiative accelerations in both the FLD approximation and the RT approach to gravity at two different locations from the proto-star as a function of the stellar mass M_∗. The steep increase in the radiative forces at roughly 6–8 M_⊙ marks the point in time at which the stellar luminosity in the Hosokawa & Omukai (2009) tracks starts to dominate over the accretion luminosity.

The results of the analytical estimate derived in the previous section and visualized in Fig. 8 fully support the outcome of the radiation hydrodynamical simulations: In the RT case, the cavity shell becomes super-Eddington (ratio of radiative to gravitational acceleration of unity) for a 10 M_⊙ star and beyond. Afterwards, the radiative acceleration increases rapidly to be one or two orders of magnitude higher than gravity. The cavity shell is highly super-Eddington in all RT+FLD simulations. Furthermore, the continuous enhancement in radiative acceleration from a 10 M_⊙ up to a 26 M_⊙ star explains the difference between the first and second outflow launch in the low-mass (M_core = 50   M_⊙) case, which are not detected in the M_core = 100   M_⊙ cases, yielding higher accretion rates and more rapid stellar evolution.

At large radii, also in the FLD runs, the cavity shell is finally super-Eddington, but even for a 30 M_⊙ star only by a factor of a few. Most of the time, the cavity shell in the FLD case is in marginal Eddington equilibrium. In the runs for the FLD approximation, the cavity shell can be estimated to be in marginal Eddington equilibrium for a proto-stellar mass of roughly 20 M_⊙ and an extent of the cavity of roughly 2000 AU (left panel of Fig. 8). As a consequence, the launched cavity shell stops its expansion along the polar axis at a height R_cavity = 2000 AU, where it undergoes an epoch of marginal Eddington equilibrium, i.e. radiative forces are balanced by gravity. In this situation, these cavity shells are supposed to become radiative Rayleigh-Taylor unstable (Jacquet & Krumholz 2011). This is in full agreement with the FLD simulation results herein.

This behavior changes drastically when the direct absorption of stellar light is taken into account. In the simulations, which treat the stellar irradiation with a RT step, the dust grains of the cavity shell absorb the stellar light according to its frequency. In this way, the acceleration of the shell turns out to be one to two orders of magnitude higher than in the FLD approximation. In the analytical estimate of the RT case, the environment of the 20 M_⊙ proto-star is super-Eddington by a factor of 20 and increases at large radii up to 100 for a proto-star of 30 M_⊙.

In the simulations, the difference in the radiation transport methods exists only up to an optical depth of τ(ν) ≈ 1 in the cavity shell, i.e. up to the location where all stellar photons are absorbed. In our hybrid radiation transport scheme, the re-emission of the photons by dust grains (mostly at longer wavelengths in the infrared regime) is computed in the FLD solver step. The corresponding length scale is given by l_∗ = (κ_∗(ν)   ρ_cavity)^-1 and therefore depends on the frequency ν of the photons as well as on the actual shell density ρ_cavity. This absorption length scale l_∗ is shown as a function of frequency ν for different densities ρ of the cavity shell in Fig. 2.

8. Comparisons

8.1. Comparison with a 3D gray FLD simulation

Krumholz et al. (2009) presented a self-gravity radiation hydrodynamics simulation of the collapse of a 100 M_⊙ pre-stellar core. The outer core radius was chosen to be 0.1 pc and the density profile declines in proportion to r^-1.5. The initial isothermal temperature of the core is 20 K. The pre-stellar core is initially in solid-body rotation without any turbulent motion. The model describes a so-called monolithic collapse scenario. The applied radiation transport method is the gray flux-limited fiffusion approximation. These properties of this configuration are the same as in our simulation run “C-FLD”. The equations are solved on a 3D adaptive mesh refinement grid in cartesian coordinates. Densities above the Jeans density on the finest grid level are represented by sink particles.

During the simulation, bipolar “radiation-filled bubbles” are blown into the environment of the forming massive star. To clarify our terminology, these “radiation-filled bubbles” correspond to the “radiation-pressure-dominated cavities” and the “bubble wall” is called a “cavity shell” in this paper. At an extent of roughly from 1200 to 2000 AU (from Krumholz et al. 2009, Fig. 1 therein) the cavity shell is subject to a “radiative Rayleigh-Taylor instability”, meaning that the radiation pressure expands the optically thin region only at specific solid angles, while at other solid angles the concentrated mass load on top of the cavity shell is able to penetrate into the low-density region.

This mechanism of shell instability, the focus of radiation pressure onto solid angles with lower optical depth, and the collapse of condensed material from the infalling envelope at other solid angles, precisely matches the outcome of our simulations, if the FLD approximation is used for the stellar radiation feedback. After the epoch of instability – the caving-in of material – the continuing evolution of their simulation differs in terms of morphology from our simulations because of the evolving non-axisymmetric structures. We comment further on this difference after the following paragraph on the RT + FLD simulations.

In our simulations using the frequency-dependent RT step for the stellar irradiation, the mass load on top of the cavity shell causes a dependence of the optical depth on the solid angle (owing to the decreasing centrifugal forces towards the pole). However, in these simulations, the caving-in of material at solid angles of high optical depth is prohibited; this is most likely caused by the treatment of the direct stellar irradiation by means of a RT approach preserving the natural isotropy of the irradiation up to the first absorption layer. Furthermore, the RT approach accounts for the respective frequency-dependent absorption coefficients as demonstrated in the previous sections.

The cavity is affected by Rayleigh-Taylor-like instabilities, when the FLD approximation is used for the radiation transport. This implication is verified by the axial symmetric simulations performed herein as well as by the 3D simulation in Krumholz et al. (2009). Improving the radiation transport scheme by incorporating a ray-tracing of the direct stellar irradiation corrects for this behavior and leads to a stable cavity growth. This implication is verified by the axial symmetric simulations performed herein as well as by the 3D simulation in Kuiper et al. (2011).The largest differences between the numerical treatment in the simulation of Krumholz et al. (2009) and the simulations with similar initial conditions presented herein are twofold: the radiation transport approaches and our assumed axial symmetry. As mentioned above, the simplification to axial symmetry increases the difficulty of a comparison of epochs after the occurrence of the instability, but the overall simulation results indicate that the outcome, regardless of whether the instability occurs, does not depend on the axial symmetric approach:

• the 2D simulations in the FLD approximation match the outcome of the 3D simulation in the FLD approximation, namely the instability of the cavity;

• in the 2D simulations with the improved radiation transport scheme (ray-tracing + FLD), the cavity remains stable;

• in our 3D simulation (Kuiper et al. 2011), using a frequency-dependent RT step to account for stellar irradiation, the outflow cavity contains strong non-axisymmetric features, but no instability is detected.

This implies that the stability of the outflow cavity is strongly effected by the direct stellar irradiation of the massive star, independent of the symmetry.

8.2. Comparison with frequency-dependent FLD simulations

In Yorke & Sonnhalter (2002), the authors presented six self-gravity radiation hydrodynamics simulations of the collapse of a pre-stellar core. The initial core mass in the simulations was chosen to be 30, 60, and 120 M_⊙ with an outer core radius of 0.05, 0.1, and 0.2 pc. The density profile drops proportional to r^-2. Since in these simulations an additional inflow of mass at the outer boundary into the computational domain is allowed, the total mass reservoir of the forming star is not limited to these initial core masses. The initial isothermal temperature of the core is 20 K. The pre-stellar core is initially in solid-body rotation without any turbulent motion. The equations are solved on a static nested grid of three levels in cylindrical coordinates assuming axial symmetry and midplane symmetry. The forming massive star is represented by a sink cell at the origin of the domain. The radiation transport method used is the gray as well as the frequency-dependent FLD approximation.

The outflow properties in the simulations with gray FLD and initial core masses of 30 and 120 M_⊙ are neither visualized nor discussed in Yorke & Sonnhalter (2002). In the simulation with gray FLD and an initial core mass of 60 M_⊙, the massive star stops its mass growth at M_∗ = 20.7   M_⊙ and maintains a luminosity of L_∗ = 5.2 × 10⁴   L_⊙. No polar cavity forms before the end of the simulation at 110 kyr.

In the three simulations with frequency-dependent FLD, polar cavities are formed. In the 30 M_⊙ case, the cavity interacts with the infalling envelope and collapses within a few thousand years. In the two higher mass cases, the infall even at larger radii is reversed shortly after the launch of the outflow cavity; therefore the cavity does not interact with a gravitationally driven infall anymore and simply expands in time. The 60 M_⊙ simulation run was stopped after an evolutionary age of 45 kyr. Nevertheless, in the 120 M_⊙ case, the radiation pressure eventually evacuates the stellar environment in all directions, including the disk midplane. After its first expansion phase, this radiation-pressure-dominated cocoon-like structure decreases again, especially in the direction of the highest optical depth (the midplane), similar to the behavior of the expansion of the cavity shell in the FLD simulation by Krumholz et al. (2009) as well as in our FLD simulations. The expansion velocity of the cavities of v ≳ 10kms^-1 corresponds more to the velocities observed in our FLD simulation (v = 6 − 30kms^-1) than those using the ray-tracing method (v ~ 100kms^-1).

The frequency dependence of the FLD solver (as in our RT method) might also be important. Using the FLD approximation in the gray limit leads to a downgrading of the radiation field throughout the cavity region, i.e. the dust grains in the cavity shell absorb the stellar flux with the Rosseland mean opacity at the local radiation temperature of the cavity shell, which does not resemble the broad stellar irradiation spectrum. By dividing the stellar flux into several frequency bins and treating the flux of each bin in the FLD approximation, one also can ensure that the high-frequency part of the stellar spectrum is absorbed with the appropriate absorption coefficient. In addition to a more correct description of the absorption behavior in the cavity shell, the RT approach provides a most realistic representation of the photon paths emitted at the stellar photosphere. Owing to the integral over the solid angle in the FLD approximation, the resulting flux does not include the correct propagation direction.

Moreover, including the long-wavelength part of the stellar spectrum explicitly can accelerate the layers on top of the swept-up cavity shell, which is optically thick for the short wavelengths only. Furthermore, Schartmann et al. (2011) demonstrated that in numerical simulations using a RT approach for the radiation feedback, the acceleration of a gas clump by radiation pressure against a gravitational potential is inversely proportional to the optical depth of the clump, as expected by analytic estimates similar to the computation of the Eddington limit.

8.3. Comparison with an analytic stability analysis

Jacquet & Krumholz (2011) derived a criterion for the stability of radiation-pressure-dominated cavities for the adiabatic approximation. They provided growth times of the radiative Rayleigh-Taylor instability in cavity shells around massive stars.

This analytic analysis fully supports the outcome of our investigations. From our estimate of the accelerations in Sect. 7, we compute the expansion timescale of the cavity shell to be ${\mathit{t}}_{\mathrm{cavity}}\mathrm{=}\sqrt{\frac{\mathrm{2} {\mathit{R}}_{\mathrm{cavity}}}{{\mathit{a}}_{\mathrm{rad}}\mathrm{-}{\mathit{a}}_{\mathrm{grav}}}}\mathrm{·}$(15)In the case of FLD simulations, the radiation pressure onto the cavity shell is in marginal equilibrium with gravity (E ≈ 1) for long epochs in time. The corresponding shell expansion timescale goes to infinity during these epochs, allowing the radiative Rayleigh-Taylor instability to set in. In the case of RT simulations, radiation pressure exceeds gravity by 1–2 orders of magnitude, leading to a much shorter cavity shell expansion timescale. For a stellar mass of 20 M_⊙ and a cavity radius of R_cavity = 1400 AU, one obtains an expansion timescale of t_cavity ≈ 0.68 kyr.

Jacquet & Krumholz (2011) stated that the growth times of the radiative Rayleigh-Taylor instability in cavity shells around massive stars should always be shorter then 1 kyr, but actually the numbers one reads from their plots (Fig. 3 and 4) are instead 6.3–40 kyr depending on the actual stellar mass and the wavelength of the instability. Hence, the expansion timescale of the shell is at least one order of magnitude shorter than the timescale for the growth of the radiative Rayleigh-Taylor instability. Even, for the largest shell radius of R_cavity = 0.1 pc at the outer boundary of our computational domain, the expansion timescale is t_cavity ≈ 10 kyr. During this 10 kyr, the massive star increases its mass approximately from 20 to 30 M_⊙ by disk accretion. Therefore, a cavity shell that has formed around a massive star is expected to be able to expand up to the stellar cluster scale ( ≈ 0.1 pc) before being prone to the radiative Rayleigh-Taylor instability.

One might discuss whether an actual outflow of a massive star might fulfill the conditions for adiabaticity. In Jacquet & Krumholz (2011), the authors assume that the cavity shell is in the “adiabatic regime”. This is estimated in Jacquet & Krumholz (2011) in Eq. (84) under the assumption of a hot cavity shell with T ≲ 1100 K at a very large radius of R_cavity ≈ 10⁴ AU. In the simulation data of Krumholz et al. (2009) and our simulations in this paper, we observe much lower shell temperatures, even at much smaller radii. As a consequence, for the stellar evolutionary tracks of Hosokawa & Omukai (2009) and the shell temperature estimate given in Eq. (10), the cavity shell is in the non-adiabatic regime even out to an extent of 0.1 pc for a stellar mass of M_∗ ≤ 40M_⊙. In the non-adiabatic regime, the radiative Rayleigh-Taylor instability is damped by diffusion. The quantitative change caused by the damping remains unclear, and the growth rates of the radiative Rayleigh-Taylor instability in the non-adiabatic regime have not yet been derived.

8.4. Comparison with observations

As mentioned in the introduction, radiation-pressure-dominated cleared cavities or optically thin bubbles of the discussed sizes around massive proto-stars have not yet been observed. Similar cavities are also proposed to form during a non-radiative but magnetized pre-stellar core collapse as shown in MHD simulations of Banerjee & Pudritz (2006, 2007) and Hennebelle et al. (2011). The observation of cavities, especially in their initial launch phase, are hampered by the extinction (the cavities are still embedded in the infalling envelope) as well as the limited resolution of the inner core structure close to the proto-star. However, large-scale polar cavities could potentially diminish the extinction of the envelope at pole-on viewing angles.

Moreover, a direct comparison of the detailed morphology of the simulated outflow cavities is of course restricted by the limitations and caveats of the numerical simulations discussed in Sect. 2. Nevertheless, the basic properties of the radiation-pressure-dominated outflows investigated in this paper include a shell velocity on the order of v ≈ 100kms^-1, as well as a mass-loss rate of roughly Ṁ ≈ 10^-4   M_⊙yr^-1 (see Fig. 5).

9. Summary

We have performed six self-gravity radiation hydrodynamics simulations of a collapsing pre-stellar core including the launch of a radiation-pressure-dominated outflow cavity by the centrally forming star. For the radiation transport method, three of these simulations were performed using the FLD approximation and three use the hybrid radiation transport scheme introduced in Kuiper et al. (2010b). In the hybrid radiation transport scheme, the stellar irradiation is computed in a frequency-dependent RT step and only the computation of the thermal dust emission makes use of the FLD approximation. The three different initial conditions of the simulations vary in terms of the initial pre-stellar core mass (50 and 100 M_⊙) as well as the initial density slope (ρ ∝ r^-1.5 and ρ ∝ r^-2). The analysis of these simulations indicates that only in the case of the FLD approximation does the cavity shell undergo an epoch of marginal Eddington equilibrium, i.e. radiative forces are balanced by gravity. As a consequence, the expansion of the outflow cavity along the polar axis stops and the infalling envelope material gathered on top of the cavity shell caves in. In contrast to this scenario, simulations using the hybrid radiation transport scheme predict that there are stable and rapidly growing outflow cavities, and that no long-term epochs of marginal Eddington equilibrium occur.

We analyzed the outcome of these simulations in more detail using analytical estimates of the radiative acceleration and the growth timescale of the cavities with respect to the radiation transport method. In contrast to the FLD assumption that dust grains absorb photons according to the local radiation temperature, dust grains in the shell on top of the optically thin cavity are directly irradiated by photons from the luminous massive star. Therefore, including the feedback effect of direct stellar irradiation via a frequency-dependent RT approach leads to the prediction of stronger radiative forces in the cavity shell. Including the direct absorption of stellar photons enhances the acceleration of the shell by between one and two orders of magnitude in comparison with the FLD assumption. This analytic estimate confirms the much higher shell velocity detected in the radiation hydrodynamics simulations including the RT step for stellar radiation feedback. These estimates of the relevant timescales imply that this enhancement in acceleration leads to a cavity expansion time that is much shorter than the timescale needed for the radiative Rayleigh-Taylor instability to occur.

The results of this study does not prohibit in general the existence of (radiative) instabilities in cavity shells. The simulation results for the three different initial conditions show a dependence of the radiative launching and cavity growth on the density distribution of the stellar environment, e.g. initially flatter density profiles lead to larger amounts of mass at larger radii and, therefore, a heavier mass load on top of the cavity shell. The limitations and caveats of the numerical simulations (e.g. the assumptions of no magnetic fields and perfect dust-to-gas coupling) do not allow a complete or even comprehensive description of the complex outflow region of a massive star.

Nevertheless, this study clearly emphasizes the strong influence of the (spectral) stellar irradiation on both the dynamics and the stability of the first absorption layer. It has been demonstrated that a careful numerical treatment of this direct stellar irradiation feedback can be achieved by using a RT approach. For comparison with future high-resolution observations of the inner core regions, the velocity and the mass-loss rate of these radiation-pressure-dominated outflows are given: in the simulations, we find a shell velocity on the order of v ≈ 100kms^-1 and a mass-loss rate of roughly Ṁ ≈ 10^-4   M_⊙yr^-1. The qualitative and quantitative analyses of our numerical simulations, including different initial conditions and radiation transport methods combined with the analytical estimates of the relevant forces and timescales, cast severe doubt on the occurrence of the radiative Rayleigh-Taylor instability in the radiation-pressure-dominated cavities around forming high-mass stars.

Online material

Movie 1 Access here

Movie 2 Access here

Movie 3 Access here

Acknowledgments

We thank Malcolm Walmsley and Harold Yorke for fruitful discussions of the paper. We acknowledge the help of Andrea Mignone, the main developer of the open source MHD code Pluto, as well as Petros Tzeferacos, who implemented the viscosity tensor into Pluto. We thank our colleague Mario Flock for his dedicated participation in the enhancements of our code version and Takashi Hosokawa for contributing the stellar evolutionary tracks. Thanks to Sarah Ragan for proof-reading the manuscript.

References

All Tables

All Figures