Issue 
A&A
Volume 537, January 2012



Article Number  A122  
Number of page(s)  13  
Section  Interstellar and circumstellar matter  
DOI  https://doi.org/10.1051/00046361/201117808  
Published online  19 January 2012 
On the stability of radiationpressuredominated cavities^{⋆}
^{1} Argelander Institute for Astronomy, Bonn University, Auf dem Hügel 71, 53121 Bonn, Germany
email: kuiper@astro.unibonn.de
^{2} Max Planck Institute for Astronomy, Königstuhl 17, 69117 Heidelberg, Germany
Received: 2 August 2011
Accepted: 8 November 2011
Context. When massive stars exert a radiation pressure onto their environment that is higher than their gravitational attraction (superEddington condition), they launch a radiationpressuredriven outflow, which creates cleared cavities. These cavities should prevent any further accretion onto the star from the direction of the bubble, although it has been claimed that a radiative RayleighTaylor instability should lead to the collapse of the outflow cavity and foster the growth of massive stars.
Aims. We investigate the stability of idealized radiationpressuredominated cavities, focusing on its dependence on the radiation transport approach used in numerical simulations for the stellar radiation feedback.
Methods. We compare two different methods for stellar radiation feedback: gray fluxlimited diffusion (FLD) and raytracing (RT). Both methods are implemented in our selfgravity radiation hydrodynamics simulations for various initial density structures of the collapsing clouds, eventually forming massive stars. We also derive simple analytical models to support our findings.
Results. Both methods lead to the launch of a radiationpressuredominated outflow cavity. However, only the FLD cases lead to prominent instability in the cavity shell. The RT cases do not show such instability; once the outflow has started, it precedes continuously. The FLD cases display extended epochs of marginal Eddington equilibrium in the cavity shell, making them prone to the radiative RayleighTaylor instability. In the RT cases, the radiation pressure exceeds gravity by 1–2 orders of magnitude. The radiative RayleighTaylor instability is then consequently suppressed. It is a fundamental property of the gray FLD method to neglect the stellar radiation temperature at the location of absorption and thus to underestimate the opacity at the location of the cavity shell.
Conclusions. Treating the stellar irradiation in the gray FLD approximation underestimates the radiative forces acting on the cavity shell. This can lead artificially to situations that are affected by the radiative RayleighTaylor instability. The proper treatment of direct stellar irradiation by massive stars is crucial for the stability of radiationpressuredominated cavities.
Key words: stars: winds, outflows / stars: massive / stars: formation / radiative transfer / hydrodynamics / methods: numerical
Movies are available in electronic form at http://www.aanda.org
© 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 starforming regions involve several largescale 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). Nonradiative 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 threedimensional (3D) (magneto)hydrodynamics code in Kuiper et al. (2010b). The solver algorithm superposes two radiation fields: a frequencydependent raytracing (RT) component representing the stellar irradiation field and a frequencyaveraged (gray) fluxlimited diffusion (FLD) component representing the thermal dust emission. These splitting schemes were previously used in onedimensional (1D) simulations (Wolfire & Cassinelli 1986, 1987; Edgar & Clarke 2003) and in twodimensional (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 & RowanRobinson 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 timedependent 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 (freefloating) limit. However, the approximation becomes inaccurate in multidimensional 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 frequencyaveraged (gray) limit remains the default technique in modern multidimensional radiation hydrodynamics. In an exceptional case, a frequencydependent FLD solver was used in Yorke & Sonnhalter (2002).
From the theoretical point of view, the stability of radiationpressuredominated 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 radiationpressuredominated cavity shells are subject to the socalled “radiative RayleighTaylor 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 starformation 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. Longlived 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 radiationpressuredominated 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 radiationpressuredriven wideangle 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 socalled 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 dustfree zone around the massive star up to the dust sublimation front the (remaining) gas opacity is set to be constant κ_{gas} = 0.01cm^{2} 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 lowdensity cavity is potentially enriched, but certainly influenced, by a stellar wind from the central massive star, which is neglected herein. Early protostellar 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 preexisting 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 RayleighTaylor 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 RayleighTaylor 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 RayleighTaylor 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 in the gray approximation. For simplification, the evaporation of dust grains is neglected, but would result in an even larger cooling lengthscale. Hence, we can be sure that the length scale of cooling is resolved.
Fig. 1 The cooling length scale l_{c} (dashed line) of the thermal dust emission in the gray approximation as a function of the actual location of the cavity shell R_{cavity}. The solid line denotes the resolution of our numerical grid. 

Open with DEXTER 
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 blackbody 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 lowdensity cleared cavity, the shell of sweptup material on top of this cavity, and the remnant of the infalling 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 (twodimensional 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.
Fig. 2 The absorption length scale l_{∗} regarding the stellar irradiation spectrum as a function of frequency ν or wavelength λ, respectively. From top to bottom, the different solid lines denote a density of , respectively. The two horizontal lines mark the highest and lowest resolution of the spherical grid. The right vertical axis shows the scale for the stellar spectra: from top to bottom, the dashed lines denote black body spectra of T_{∗} = 100 000,30 000,10 000,and3000 K. The effect of potential dust evaporation is not included here, but would lead to a higher resolution of the numerical grid in terms of the irradiation spectrum. 

Open with DEXTER 
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 nonaxially 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 radiationpressuredominated cavity remains stable.
3. Method
3.1. Equations
To follow the motion of the gas, we solve the equations of compressible selfgravity hydrodynamics, including shearviscosity 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 (raytracing + fluxlimited diffusion) and fluxlimited diffusion only.
3.1.1. Raytracing + fluxlimited 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 frequencydependent stellar irradiation and the flux of the frequencyaveraged thermal radiation energy density Fwhere Eq. (1)denotes the evolution of the thermal radiation energy density E_{R}. The factor 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 socalled fluxlimited 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 frequencydependent stellar irradiation in a raytracing 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 frequencydependent 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 highmass 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 (3)with the Planck mean opacities κ_{P}.
Numerical details, test cases, including a comparison of gray and frequencydependent 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. Fluxlimited diffusion only
For simulations labeled “FLD”, we do not raytrace the stellar irradiation, but add the luminosity of the forming star as a source term to the right hand side of the FLD equation (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 timeindependent 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 nonuniform grid is chosen to be (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 semipermeable 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 prestellar core.
The Pluto code uses highorder 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 HartenLaxVan Leer Riemann solver and a socalled “minmod” flux limiter using piecewise linear interpolation and a RungeKutta 2 timeintegration, also known as the predictorcorrectormethod; 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_{0} = 20K) prestellar core of gas and dust. The initially constant dusttogas mass ratio is chosen to be M_{dust}/M_{gas} = 0.01. Dynamically, the model describes a socalled quiescent collapse scenario without initial turbulent motion (u_{r} = u_{θ} = 0) and the core is initially in slow solidbody rotation . 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 prestellar 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 prestellar core collapse model, we refer to Kuiper et al. (2010a, 2011).
Overview of simulations performed.
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 followup 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 FLDonly 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.
Fig. 3 Simulation snapshots of the gas density for one of the FLD (left panels) and one of the raytracing + FLD (right panels) runs for three different points in time during the launch of the radiation pressure dominated cavities. The data is taken from the runs with initial core mass M_{core} = 50 M_{⊙} and a radial density slope of β = −2. The spatial section of the RT+FLD snapshots increases with time to follow the rapid expansion of the outflow cavity. Animations of the launch of the radiationpressuredominated cavities in the simulations are available as online material, too. 

Open with DEXTER 
In addition to the different growth rates, the radiationpressuredriven 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 sweptup 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 raytracing method, the isotropic stellar irradiation flux directly impinges onto the sweptup mass on top of the polar cavity and pushes the mass to larger radii. The resulting largescale 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 FLDonly 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 RayleighTaylor instability.
5. Quantitative analysis of the cavity growth
We analyze quantitatively the timedependent 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.
Fig. 4 Size of the cavity R_{cavity} as a function of time t for the three different initial conditions. All runs were performed using the RT+FLD hybrid radiation transport scheme (solid lines) and the FLD scheme (dashed lines). The lowest and highest extent of R_{cavity} of 10 AU and 0.1 pc, respectively, are given by the inner and outer rim of the computational domain. 

Open with DEXTER 
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 raytracing scheme (labeled “RT+FLD”) than simply included in the fluxlimited 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 with FLD, the cavity shell region has been proposed to be subject to the radiative RayleighTaylor 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 sweptup mass in the cavity shell. The peaks in the density distributions correspond to and in the RT+FLD and FLD cases, respectively.
Fig. 5 Gas density ρ (upper left panel), temperature T (upper right panel), radial velocity v (lower left panel), and radial massloss rate (lower right panel) as a function of height z above the star along the polar axis. Data is for the run with initial core mass M_{core} = 100 M_{⊙} and radial density slope β = −2. The snapshot in time (t = 16.7 kyr) is chosen in such a way that the cavity shell has arrived at the same location in both simulations (RT+FLD and FLD only). The lowest density of in the upper left panel is given by the numerical floor value of the hydrodynamics solver. 

Open with DEXTER 
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 sweptup 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 sweptup 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 sweptup mass on top of the cavity are embedded in a local radiation field, whereas the RT+FLD method takes into account the (frequencydependent) 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 multidimensional environment with a nonisotropic optical depth, this assumption of the FLD approximation breaks down. Including higherorder 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 taytracing equation, which takes into account the correct photon propagation direction, i.e. the stellar irradiation flux directly impinges on the sweptup mass shell at the top of the optically thin cavity.
7. Analytical determination of the acceleration caused by stellar irradiation in the raytracing approach and the fluxlimited 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 sweptup 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 sweptup mass at the top of the optically thin outflow cavity. The dust grains in the sweptup material therefore absorb the stellar flux following the frequencydependent 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 frequencydependent opacity table (Laor & Draine 1993). To obtain the absorption coefficient in a given grid cell, these opacities are multiplied by the local dusttogas 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 powerlaw 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 sweptup 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.
Fig. 6 Rosseland mean opacities κ_{R} (left panel) and Planck mean opacities κ_{P} (right panel) per gram gas as a function of the local radiation temperature T. The label “FLD 1” denotes the opacity description by Krumholz et al. (2007), Eq. (11) therein. The label “FLD 2” denotes the opacities used in Kuiper et al. (2010a, 2011) as well as in this paper; these opacities are computed based on the densitydependent evaporation of dust grains using the low density of the polar cavities. 

Open with DEXTER 
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 frequencydependent 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 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 are multiplied by the evaporation probability, and, therefore, also depend on the local gas density and temperature.
Fig. 7 Planck mean opacities per gram gas as a function of the radiation temperature T_{∗} of the stellar photosphere as used in the RT step of the hybrid radiation transport scheme in Kuiper et al. (2010a, 2011) and in this paper. The upper horizontal axis shows the mass of the protostar, which corresponds to the stellar surface temperature T_{∗}, taken from the Hosokawa & Omukai (2009) evolutionary track for a constant accretion rate of Ṁ = 10^{3} M_{⊙} yr^{1}. 

Open with DEXTER 
The mass of the protostar 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 higherorder treatment of stellar irradiation via RT results in a stable cavity shell with radial velocities up to v ≳ 100 km s^{1} and massloss rates of Ṁ ≈ 2 × 10^{4} M_{⊙}yr^{1}.
7.2. Radiative acceleration
To compute the radiative acceleration a_{rad} of the sweptup 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 wellknown Eddington limit. We start with the formula for radiative acceleration (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 sweptup 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 (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 sweptup 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 reemit 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.
Fig. 8 Radiative and gravitational acceleration at a position R_{cavity} = 2000 AU (left panel) and R_{cavity} = 10 000 AU (right panel) above the massive protostar as a function of the stellar mass M_{∗} for three different radiation transport approaches: The label “FLD 1” denotes the FLD approximation with the opacity description of Krumholz et al. (2007). The label “FLD 2” refers to the FLD approximation with the opacities used in Kuiper et al. (2010a, 2011) and in this paper. The label “RayTracing” refers to the RT step as in the hybrid radiation transport scheme. 

Open with DEXTER 
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 (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_{∗}(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 largescale 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 protostar 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_{∗}(10)where β_{opacity} represents the exponent of the slope of the frequencydependent opacities at longer wavelengths. The assumption of gray absorption would, e.g., imply β_{opacity} = 0 and, therefore, results in a temperature slope of . 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 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 in the same way as the gravitational acceleration. The ratio of these accelerations, therefore, is scalefree and depends only on the protostellar properties leading to the wellknown formula for the generalized Eddington limit (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 protostar 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 superEddington (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 superEddington 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 lowmass (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 superEddington, 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 protostellar 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 RayleighTaylor 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_{⊙} protostar is superEddington by a factor of 20 and increases at large radii up to 100 for a protostar 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 reemission 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 selfgravity radiation hydrodynamics simulation of the collapse of a 100 M_{⊙} prestellar 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 prestellar core is initially in solidbody rotation without any turbulent motion. The model describes a socalled monolithic collapse scenario. The applied radiation transport method is the gray fluxlimited fiffusion approximation. These properties of this configuration are the same as in our simulation run “CFLD”. 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 “radiationfilled bubbles” are blown into the environment of the forming massive star. To clarify our terminology, these “radiationfilled bubbles” correspond to the “radiationpressuredominated 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 RayleighTaylor 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 lowdensity 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 cavingin of material – the continuing evolution of their simulation differs in terms of morphology from our simulations because of the evolving nonaxisymmetric structures. We comment further on this difference after the following paragraph on the RT + FLD simulations.
In our simulations using the frequencydependent 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 cavingin 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 frequencydependent absorption coefficients as demonstrated in the previous sections.
The cavity is affected by RayleighTaylorlike 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 raytracing 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 (raytracing + FLD), the cavity remains stable;
in our 3D simulation (Kuiper et al. 2011), using a frequencydependent RT step to account for stellar irradiation, the outflow cavity contains strong nonaxisymmetric features, but no instability is detected.
8.2. Comparison with frequencydependent FLD simulations
In Yorke & Sonnhalter (2002), the authors presented six selfgravity radiation hydrodynamics simulations of the collapse of a prestellar 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 prestellar core is initially in solidbody 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 frequencydependent 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^{4} L_{⊙}. No polar cavity forms before the end of the simulation at 110 kyr.
In the three simulations with frequencydependent 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 radiationpressuredominated cocoonlike 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 raytracing 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 highfrequency 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 longwavelength part of the stellar spectrum explicitly can accelerate the layers on top of the sweptup 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 radiationpressuredominated cavities for the adiabatic approximation. They provided growth times of the radiative RayleighTaylor 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 (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 RayleighTaylor 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 RayleighTaylor 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 RayleighTaylor 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 RayleighTaylor 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^{4} 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 nonadiabatic regime even out to an extent of 0.1 pc for a stellar mass of M_{∗} ≤ 40M_{⊙}. In the nonadiabatic regime, the radiative RayleighTaylor instability is damped by diffusion. The quantitative change caused by the damping remains unclear, and the growth rates of the radiative RayleighTaylor instability in the nonadiabatic regime have not yet been derived.
8.4. Comparison with observations
As mentioned in the introduction, radiationpressuredominated cleared cavities or optically thin bubbles of the discussed sizes around massive protostars have not yet been observed. Similar cavities are also proposed to form during a nonradiative but magnetized prestellar 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 protostar. However, largescale polar cavities could potentially diminish the extinction of the envelope at poleon 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 radiationpressuredominated outflows investigated in this paper include a shell velocity on the order of v ≈ 100kms^{1}, as well as a massloss rate of roughly Ṁ ≈ 10^{4} M_{⊙}yr^{1} (see Fig. 5).
9. Summary
We have performed six selfgravity radiation hydrodynamics simulations of a collapsing prestellar core including the launch of a radiationpressuredominated 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 frequencydependent 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 prestellar 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 longterm 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 frequencydependent 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 RayleighTaylor 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 dusttogas 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 highresolution observations of the inner core regions, the velocity and the massloss rate of these radiationpressuredominated outflows are given: in the simulations, we find a shell velocity on the order of v ≈ 100kms^{1} and a massloss 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 RayleighTaylor instability in the radiationpressuredominated cavities around forming highmass stars.
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 proofreading the manuscript.
References
 Arce, H. G., Shepherd, D., Gueth, F., et al. 2007, Protostars and Planets V, 245 [Google Scholar]
 Banerjee, R., & Pudritz, R. E. 2006, ApJ, 641, 949 [NASA ADS] [CrossRef] [Google Scholar]
 Banerjee, R., & Pudritz, R. E. 2007, ApJ, 660, 479 [NASA ADS] [CrossRef] [Google Scholar]
 Bate, M. R. 2010, MNRAS, 404, L79 [NASA ADS] [Google Scholar]
 Beltran, M. T., Girart, J. M., & Estalella, R. 2006, A&A, 457, 865 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Beuther, H., Linz, H., Bik, A., Goto, M., & Henning, T. 2010, A&A, 512, A29 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bjorkman, J. E., & Wood, K. 2001, ApJ, 554, 615 [NASA ADS] [CrossRef] [Google Scholar]
 Bodenheimer, P., Yorke, H. W., Rozyczka, M., & Tohline, J. E. 1990, ApJ, 355, 651 [NASA ADS] [CrossRef] [Google Scholar]
 Boley, A. C., Durisen, R. H., Nordlund, Å., & Lord, J. 2007, ApJ, 665, 1254 [NASA ADS] [CrossRef] [Google Scholar]
 Bonnell, I. A., & Bate, M. R. 2006, MNRAS, 370, 488 [NASA ADS] [CrossRef] [Google Scholar]
 Bruenn, S. 1985, ApJS, 58, 771 [NASA ADS] [CrossRef] [Google Scholar]
 Commercon, B., Hennebelle, P., Audit, E., Chabrier, G., & Teyssier, R. 2010, A&A, 510, L3 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Cunningham, A. J., Klein, R. I., Krumholz, M. R., & McKee, C. F. 2011, ApJ, 740, 107 [NASA ADS] [CrossRef] [Google Scholar]
 Dale, J. E., & Bonnell, I. A. 2011, MNRAS, 414, 321 [NASA ADS] [CrossRef] [Google Scholar]
 Dullemond, C. P., & Dominik, C. 2004, A&A, 417, 159 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Dullemond, C. P., & Turolla, R. 2000, A&A, 360, 1187 [NASA ADS] [Google Scholar]
 Edgar, R., & Clarke, C. 2003, MNRAS, 338, 962 [NASA ADS] [CrossRef] [Google Scholar]
 Efstathiou, A., & RowanRobinson, M. 1990, MNRAS, 245, 275 [NASA ADS] [Google Scholar]
 Fallscheer, C., Beuther, H., Zhang, Q., Keto, E. R., & Sridharan, T. K. 2009, A&A, 504, 127 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Fendt, C., Vaidya, B., Porth, O., & Nezami, S. S. 2011, Jets at all Scales, 275, 383 [NASA ADS] [Google Scholar]
 Hennebelle, P., Commercon, B., Joos, M., et al. 2011, A&A, 528, 72 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hosokawa, T., & Omukai, K. 2009, ApJ, 691, 823 [NASA ADS] [CrossRef] [Google Scholar]
 Isella, A., & Natta, A. 2005, A&A, 438, 899 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Jacquet, E., & Krumholz, M. R. 2011, ApJ, 730, 116 [NASA ADS] [CrossRef] [Google Scholar]
 Klahr, H., & Kley, W. 2006, A&A, 445, 747 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Klahr, H., Henning, T., & Kley, W. 1999, ApJ, 514, 325 [NASA ADS] [CrossRef] [Google Scholar]
 Kley, W. 1989, A&A, 208, 98 [NASA ADS] [Google Scholar]
 Krumholz, M. R., Klein, R. I., & McKee, C. F. 2007, ApJ, 656, 959 [NASA ADS] [CrossRef] [Google Scholar]
 Krumholz, M. R., Klein, R. I., McKee, C. F., Offner, S. S. R., & Cunningham, A. J. 2009, Science, 323, 754 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Kuiper, R., Klahr, H., Beuther, H., & Henning, T. 2010a, ApJ, 722, 1556 [NASA ADS] [CrossRef] [Google Scholar]
 Kuiper, R., Klahr, H., Dullemond, C. P., Kley, W., & Henning, T. 2010b, A&A, 511, A81 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kuiper, R., Klahr, H., Beuther, H., & Henning, T. 2011, ApJ, 732, 20 [NASA ADS] [CrossRef] [Google Scholar]
 Laor, A., & Draine, B. T. 1993, ApJ, 402, 441 [NASA ADS] [CrossRef] [Google Scholar]
 Levermore, C. D., & Pomraning, G. C. 1981, ApJ, 248, 321 [NASA ADS] [CrossRef] [Google Scholar]
 McKee, C. F., & Ostriker, E. C. 2007, ARA&A, 45, 565 [NASA ADS] [CrossRef] [Google Scholar]
 Mignone, A., Bodo, G., Massaglia, S., et al. 2007, ApJS, 170, 228 [NASA ADS] [CrossRef] [Google Scholar]
 Murray, S. D., Castor, J. I., Klein, R. I., & McKee, C. F. 1994, ApJ, 435, 631 [NASA ADS] [CrossRef] [Google Scholar]
 Nakano, T. 1989, ApJ, 345, 464 [NASA ADS] [CrossRef] [Google Scholar]
 Pascucci, I., Wolf, S., Steinacker, J., et al. 2004, A&A, 417, 793 [NASA ADS] [CrossRef] [EDP Sciences] [MathSciNet] [Google Scholar]
 Peters, T., Banerjee, R., Klessen, R. S., et al. 2010, ApJ, 711, 1017 [NASA ADS] [CrossRef] [Google Scholar]
 Peters, T., Banerjee, R., Klessen, R. S., & Mac Low, M.M. 2011, ApJ, 729, 72 [NASA ADS] [CrossRef] [Google Scholar]
 Pinte, C., Harries, T. J., Min, M., et al. 2009, A&A, 498, 967 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Pollack, J. B., Hollenbach, D., Beckwith, S. V. W., et al. 1994, ApJ, 421, 615 [NASA ADS] [CrossRef] [Google Scholar]
 Price, D. J., & Bate, M. R. 2009, MNRAS, 398, 33 [NASA ADS] [CrossRef] [Google Scholar]
 Pudritz, R. E., Ouyed, R., Fendt, C., & Brandenburg, A. 2007, in Protostars and Planets V, ed. B. Reipurth, D. Jewitt, & K. Keil, 277 [Google Scholar]
 Quanz, S. P., Beuther, H., Steinacker, J., et al. 2010, ApJ, 717, 693 [NASA ADS] [CrossRef] [Google Scholar]
 Schartmann, M., Krause, M., & Burkert, A. 2011, MNRAS, 415, 741 [NASA ADS] [CrossRef] [Google Scholar]
 Shepherd, D., & Churchwell, E. 1996, ApJ, 472, 225 [NASA ADS] [CrossRef] [Google Scholar]
 Spitzer, L. 1978 (New York: WileyInterscience), 333 [Google Scholar]
 Steinacker, J., Henning, T., Bacmann, A., & Semenov, D. 2003, A&A, 401, 405 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Tanaka, K. E. I., & Nakamoto, T. 2011, ApJ, 739, L50 [NASA ADS] [CrossRef] [Google Scholar]
 Turner, N. J., Quataert, E., & Yorke, H. W. 2007, ApJ, 662, 1052 [NASA ADS] [CrossRef] [Google Scholar]
 van Leer, B. 1979, JCP, 32, 101 [NASA ADS] [CrossRef] [Google Scholar]
 Wolfire, M. G., & Cassinelli, J. P. 1986, ApJ, 310, 207 [NASA ADS] [CrossRef] [Google Scholar]
 Wolfire, M. G., & Cassinelli, J. P. 1987, ApJ, 319, 850 [NASA ADS] [CrossRef] [Google Scholar]
 Yorke, H. W., & Sonnhalter, C. 2002, ApJ, 569, 846 [NASA ADS] [CrossRef] [Google Scholar]
 Zhang, Q., Hunter, T. R., Brand, J., et al. 2001, ApJ, 552, L167 [NASA ADS] [CrossRef] [Google Scholar]
 Zhang, Q., Hunter, T. R., Beuther, H., et al. 2007, ApJ, 658, 1152 [NASA ADS] [CrossRef] [Google Scholar]
Online material
Movie 1 (Access here)
Movie 2 (Access here)
Movie 3 (Access here)
All Tables
All Figures
Fig. 1 The cooling length scale l_{c} (dashed line) of the thermal dust emission in the gray approximation as a function of the actual location of the cavity shell R_{cavity}. The solid line denotes the resolution of our numerical grid. 

Open with DEXTER  
In the text 
Fig. 2 The absorption length scale l_{∗} regarding the stellar irradiation spectrum as a function of frequency ν or wavelength λ, respectively. From top to bottom, the different solid lines denote a density of , respectively. The two horizontal lines mark the highest and lowest resolution of the spherical grid. The right vertical axis shows the scale for the stellar spectra: from top to bottom, the dashed lines denote black body spectra of T_{∗} = 100 000,30 000,10 000,and3000 K. The effect of potential dust evaporation is not included here, but would lead to a higher resolution of the numerical grid in terms of the irradiation spectrum. 

Open with DEXTER  
In the text 
Fig. 3 Simulation snapshots of the gas density for one of the FLD (left panels) and one of the raytracing + FLD (right panels) runs for three different points in time during the launch of the radiation pressure dominated cavities. The data is taken from the runs with initial core mass M_{core} = 50 M_{⊙} and a radial density slope of β = −2. The spatial section of the RT+FLD snapshots increases with time to follow the rapid expansion of the outflow cavity. Animations of the launch of the radiationpressuredominated cavities in the simulations are available as online material, too. 

Open with DEXTER  
In the text 
Fig. 4 Size of the cavity R_{cavity} as a function of time t for the three different initial conditions. All runs were performed using the RT+FLD hybrid radiation transport scheme (solid lines) and the FLD scheme (dashed lines). The lowest and highest extent of R_{cavity} of 10 AU and 0.1 pc, respectively, are given by the inner and outer rim of the computational domain. 

Open with DEXTER  
In the text 
Fig. 5 Gas density ρ (upper left panel), temperature T (upper right panel), radial velocity v (lower left panel), and radial massloss rate (lower right panel) as a function of height z above the star along the polar axis. Data is for the run with initial core mass M_{core} = 100 M_{⊙} and radial density slope β = −2. The snapshot in time (t = 16.7 kyr) is chosen in such a way that the cavity shell has arrived at the same location in both simulations (RT+FLD and FLD only). The lowest density of in the upper left panel is given by the numerical floor value of the hydrodynamics solver. 

Open with DEXTER  
In the text 
Fig. 6 Rosseland mean opacities κ_{R} (left panel) and Planck mean opacities κ_{P} (right panel) per gram gas as a function of the local radiation temperature T. The label “FLD 1” denotes the opacity description by Krumholz et al. (2007), Eq. (11) therein. The label “FLD 2” denotes the opacities used in Kuiper et al. (2010a, 2011) as well as in this paper; these opacities are computed based on the densitydependent evaporation of dust grains using the low density of the polar cavities. 

Open with DEXTER  
In the text 
Fig. 7 Planck mean opacities per gram gas as a function of the radiation temperature T_{∗} of the stellar photosphere as used in the RT step of the hybrid radiation transport scheme in Kuiper et al. (2010a, 2011) and in this paper. The upper horizontal axis shows the mass of the protostar, which corresponds to the stellar surface temperature T_{∗}, taken from the Hosokawa & Omukai (2009) evolutionary track for a constant accretion rate of Ṁ = 10^{3} M_{⊙} yr^{1}. 

Open with DEXTER  
In the text 
Fig. 8 Radiative and gravitational acceleration at a position R_{cavity} = 2000 AU (left panel) and R_{cavity} = 10 000 AU (right panel) above the massive protostar as a function of the stellar mass M_{∗} for three different radiation transport approaches: The label “FLD 1” denotes the FLD approximation with the opacity description of Krumholz et al. (2007). The label “FLD 2” refers to the FLD approximation with the opacities used in Kuiper et al. (2010a, 2011) and in this paper. The label “RayTracing” refers to the RT step as in the hybrid radiation transport scheme. 

Open with DEXTER  
In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.