Free Access
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/0004-6361/201117808
Published online 19 January 2012

© 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.01cm2   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” lc 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 Rcavity. 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   Myr-1, and the temperature at the location Rcavity 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 length-scale. Hence, we can be sure that the length scale of cooling is resolved.

thumbnail Fig. 1

The cooling length scale lc (dashed line) of the thermal dust emission in the gray approximation as a function of the actual location of the cavity shell Rcavity. 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 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.

thumbnail 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 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 Ftot is divided into the flux F(ν,r) of the frequency-dependent stellar irradiation and the flux of the frequency-averaged thermal radiation energy density Fwhere Eq. (1)denotes the evolution of the thermal radiation energy density ER. The factor depends only on the ratio of internal to radiation energy and contains the specific heat capacity, cV, 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 =  −DER). 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 (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 (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 rmin = 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 rmax = 0.1 pc, and the inner core radius is fixed to rmin = 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 (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 (T0 = 20K) pre-stellar core of gas and dust. The initially constant dust-to-gas mass ratio is chosen to be Mdust/Mgas = 0.01. Dynamically, the model describes a so-called quiescent collapse scenario without initial turbulent motion (ur = uθ = 0) and the core is initially in slow solid-body rotation . The total mass Mcore 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).

Table 1

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 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.

thumbnail Fig. 3

Simulation snapshots of the gas density for one of the FLD (left panels) and one of the ray-tracing + 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 Mcore = 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 radiation-pressure-dominated cavities in the simulations are available as online material, too.

Open with DEXTER

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 Rcavity 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.

thumbnail Fig. 4

Size of the cavity Rcavity 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 Rcavity 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 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 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 Mcore = 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 Mcore = 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 Mcore = 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 Mcore = 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.

thumbnail Fig. 5

Gas density ρ (upper left panel), temperature T (upper right panel), radial velocity v (lower left panel), and radial mass-loss 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 Mcore = 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 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 Mdust/Mgas(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.

thumbnail 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 density-dependent 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 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 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.

thumbnail 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 proto-star, 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 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   Myr-1.

7.2. Radiative acceleration

To compute the radiative acceleration arad 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 (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 Rcavity 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 (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.

thumbnail Fig. 8

Radiative and gravitational acceleration at a position Rcavity = 2000 AU (left panel) and Rcavity = 10   000 AU (right panel) above the massive proto-star 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 “Ray-Tracing” 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 Rcavity above the star (8)This radiative acceleration has to be compared with the gravitational attraction agrav of the dust grains at the cavity shell position Rcavity 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 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 Rcavity ≫ R that the distance Rcavity to the shell is much larger than the stellar radius R(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 . 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(Rcavity) 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 Rcavity of the cavity shell material. In the RT case, the opacity in regions cooler than the dust evaporation temperature is independent of the location Rcavity, i.e. the radiative acceleration scales with in the same way as the gravitational acceleration. The ratio 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 (14)

7.3. Results

In the following, we compare the accelerations depending on the remaining variables M and Rcavity. 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 (Mcore = 50   M) case, which are not detected in the Mcore = 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 Rcavity = 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 × 104   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 (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 Rcavity = 1400 AU, one obtains an expansion timescale of tcavity ≈ 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 Rcavity = 0.1 pc at the outer boundary of our computational domain, the expansion timescale is tcavity ≈ 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 Rcavity ≈ 104 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   Myr-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   Myr-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.

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

  1. Arce, H. G., Shepherd, D., Gueth, F., et al. 2007, Protostars and Planets V, 245 [Google Scholar]
  2. Banerjee, R., & Pudritz, R. E. 2006, ApJ, 641, 949 [NASA ADS] [CrossRef] [Google Scholar]
  3. Banerjee, R., & Pudritz, R. E. 2007, ApJ, 660, 479 [NASA ADS] [CrossRef] [Google Scholar]
  4. Bate, M. R. 2010, MNRAS, 404, L79 [NASA ADS] [Google Scholar]
  5. Beltran, M. T., Girart, J. M., & Estalella, R. 2006, A&A, 457, 865 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  6. Beuther, H., Linz, H., Bik, A., Goto, M., & Henning, T. 2010, A&A, 512, A29 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  7. Bjorkman, J. E., & Wood, K. 2001, ApJ, 554, 615 [NASA ADS] [CrossRef] [Google Scholar]
  8. Bodenheimer, P., Yorke, H. W., Rozyczka, M., & Tohline, J. E. 1990, ApJ, 355, 651 [NASA ADS] [CrossRef] [Google Scholar]
  9. Boley, A. C., Durisen, R. H., Nordlund, Å., & Lord, J. 2007, ApJ, 665, 1254 [NASA ADS] [CrossRef] [Google Scholar]
  10. Bonnell, I. A., & Bate, M. R. 2006, MNRAS, 370, 488 [NASA ADS] [CrossRef] [Google Scholar]
  11. Bruenn, S. 1985, ApJS, 58, 771 [NASA ADS] [CrossRef] [Google Scholar]
  12. Commercon, B., Hennebelle, P., Audit, E., Chabrier, G., & Teyssier, R. 2010, A&A, 510, L3 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  13. Cunningham, A. J., Klein, R. I., Krumholz, M. R., & McKee, C. F. 2011, ApJ, 740, 107 [NASA ADS] [CrossRef] [Google Scholar]
  14. Dale, J. E., & Bonnell, I. A. 2011, MNRAS, 414, 321 [NASA ADS] [CrossRef] [Google Scholar]
  15. Dullemond, C. P., & Dominik, C. 2004, A&A, 417, 159 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  16. Dullemond, C. P., & Turolla, R. 2000, A&A, 360, 1187 [NASA ADS] [Google Scholar]
  17. Edgar, R., & Clarke, C. 2003, MNRAS, 338, 962 [NASA ADS] [CrossRef] [Google Scholar]
  18. Efstathiou, A., & Rowan-Robinson, M. 1990, MNRAS, 245, 275 [NASA ADS] [Google Scholar]
  19. Fallscheer, C., Beuther, H., Zhang, Q., Keto, E. R., & Sridharan, T. K. 2009, A&A, 504, 127 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  20. Fendt, C., Vaidya, B., Porth, O., & Nezami, S. S. 2011, Jets at all Scales, 275, 383 [NASA ADS] [Google Scholar]
  21. Hennebelle, P., Commercon, B., Joos, M., et al. 2011, A&A, 528, 72 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  22. Hosokawa, T., & Omukai, K. 2009, ApJ, 691, 823 [NASA ADS] [CrossRef] [Google Scholar]
  23. Isella, A., & Natta, A. 2005, A&A, 438, 899 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  24. Jacquet, E., & Krumholz, M. R. 2011, ApJ, 730, 116 [NASA ADS] [CrossRef] [Google Scholar]
  25. Klahr, H., & Kley, W. 2006, A&A, 445, 747 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  26. Klahr, H., Henning, T., & Kley, W. 1999, ApJ, 514, 325 [NASA ADS] [CrossRef] [Google Scholar]
  27. Kley, W. 1989, A&A, 208, 98 [NASA ADS] [Google Scholar]
  28. Krumholz, M. R., Klein, R. I., & McKee, C. F. 2007, ApJ, 656, 959 [NASA ADS] [CrossRef] [Google Scholar]
  29. 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]
  30. Kuiper, R., Klahr, H., Beuther, H., & Henning, T. 2010a, ApJ, 722, 1556 [NASA ADS] [CrossRef] [Google Scholar]
  31. Kuiper, R., Klahr, H., Dullemond, C. P., Kley, W., & Henning, T. 2010b, A&A, 511, A81 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  32. Kuiper, R., Klahr, H., Beuther, H., & Henning, T. 2011, ApJ, 732, 20 [NASA ADS] [CrossRef] [Google Scholar]
  33. Laor, A., & Draine, B. T. 1993, ApJ, 402, 441 [NASA ADS] [CrossRef] [Google Scholar]
  34. Levermore, C. D., & Pomraning, G. C. 1981, ApJ, 248, 321 [NASA ADS] [CrossRef] [Google Scholar]
  35. McKee, C. F., & Ostriker, E. C. 2007, ARA&A, 45, 565 [NASA ADS] [CrossRef] [Google Scholar]
  36. Mignone, A., Bodo, G., Massaglia, S., et al. 2007, ApJS, 170, 228 [NASA ADS] [CrossRef] [Google Scholar]
  37. Murray, S. D., Castor, J. I., Klein, R. I., & McKee, C. F. 1994, ApJ, 435, 631 [NASA ADS] [CrossRef] [Google Scholar]
  38. Nakano, T. 1989, ApJ, 345, 464 [NASA ADS] [CrossRef] [Google Scholar]
  39. Pascucci, I., Wolf, S., Steinacker, J., et al. 2004, A&A, 417, 793 [NASA ADS] [CrossRef] [EDP Sciences] [MathSciNet] [Google Scholar]
  40. Peters, T., Banerjee, R., Klessen, R. S., et al. 2010, ApJ, 711, 1017 [NASA ADS] [CrossRef] [Google Scholar]
  41. Peters, T., Banerjee, R., Klessen, R. S., & Mac Low, M.-M. 2011, ApJ, 729, 72 [NASA ADS] [CrossRef] [Google Scholar]
  42. Pinte, C., Harries, T. J., Min, M., et al. 2009, A&A, 498, 967 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  43. Pollack, J. B., Hollenbach, D., Beckwith, S. V. W., et al. 1994, ApJ, 421, 615 [NASA ADS] [CrossRef] [Google Scholar]
  44. Price, D. J., & Bate, M. R. 2009, MNRAS, 398, 33 [NASA ADS] [CrossRef] [Google Scholar]
  45. 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]
  46. Quanz, S. P., Beuther, H., Steinacker, J., et al. 2010, ApJ, 717, 693 [NASA ADS] [CrossRef] [Google Scholar]
  47. Schartmann, M., Krause, M., & Burkert, A. 2011, MNRAS, 415, 741 [NASA ADS] [CrossRef] [Google Scholar]
  48. Shepherd, D., & Churchwell, E. 1996, ApJ, 472, 225 [NASA ADS] [CrossRef] [Google Scholar]
  49. Spitzer, L. 1978 (New York: Wiley-Interscience), 333 [Google Scholar]
  50. Steinacker, J., Henning, T., Bacmann, A., & Semenov, D. 2003, A&A, 401, 405 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  51. Tanaka, K. E. I., & Nakamoto, T. 2011, ApJ, 739, L50 [NASA ADS] [CrossRef] [Google Scholar]
  52. Turner, N. J., Quataert, E., & Yorke, H. W. 2007, ApJ, 662, 1052 [NASA ADS] [CrossRef] [Google Scholar]
  53. van Leer, B. 1979, JCP, 32, 101 [NASA ADS] [CrossRef] [Google Scholar]
  54. Wolfire, M. G., & Cassinelli, J. P. 1986, ApJ, 310, 207 [NASA ADS] [CrossRef] [Google Scholar]
  55. Wolfire, M. G., & Cassinelli, J. P. 1987, ApJ, 319, 850 [NASA ADS] [CrossRef] [Google Scholar]
  56. Yorke, H. W., & Sonnhalter, C. 2002, ApJ, 569, 846 [NASA ADS] [CrossRef] [Google Scholar]
  57. Zhang, Q., Hunter, T. R., Brand, J., et al. 2001, ApJ, 552, L167 [NASA ADS] [CrossRef] [Google Scholar]
  58. 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

Table 1

Overview of simulations performed.

All Figures

thumbnail Fig. 1

The cooling length scale lc (dashed line) of the thermal dust emission in the gray approximation as a function of the actual location of the cavity shell Rcavity. The solid line denotes the resolution of our numerical grid.

Open with DEXTER
In the text
thumbnail 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
thumbnail Fig. 3

Simulation snapshots of the gas density for one of the FLD (left panels) and one of the ray-tracing + 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 Mcore = 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 radiation-pressure-dominated cavities in the simulations are available as online material, too.

Open with DEXTER
In the text
thumbnail Fig. 4

Size of the cavity Rcavity 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 Rcavity 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
thumbnail Fig. 5

Gas density ρ (upper left panel), temperature T (upper right panel), radial velocity v (lower left panel), and radial mass-loss 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 Mcore = 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
thumbnail 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 density-dependent evaporation of dust grains using the low density of the polar cavities.

Open with DEXTER
In the text
thumbnail 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 proto-star, 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
thumbnail Fig. 8

Radiative and gravitational acceleration at a position Rcavity = 2000 AU (left panel) and Rcavity = 10   000 AU (right panel) above the massive proto-star 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 “Ray-Tracing” 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 (full-text 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 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.