H ii regions: Witnesses to massive star formation
We describe the first three-dimensional simulation of the gravitational collapse of a massive, rotating molecular cloud that includes heating by both non-ionizing and ionizing radiation. These models were performed with the FLASH code, incorporating a hybrid, long characteristic, ray tracing technique. We find that as the first protostars gain sufficient mass to ionize the accretion flow, their H ii regions are initially gravitationally trapped, but soon begin to rapidly fluctuate between trapped and extended states, in agreement with observations. Over time, the same ultracompact H ii region can expand anisotropically, contract again, and take on any of the observed morphological classes. In their extended phases, expanding H ii regions drive bipolar neutral outflows characteristic of high-mass star formation. The total lifetime of H ii regions is given by the global accretion timescale, rather than their short internal sound-crossing time. This explains the observed number statistics. The pressure of the hot, ionized gas does not terminate accretion. Instead the final stellar mass is set by fragmentation-induced starvation. Local gravitational instabilities in the accretion flow lead to the build-up of a small cluster of stars, all with relatively high masses due to heating from accretion radiation. These companions subsequently compete with the initial high-mass star for the same common gas reservoir and limit its mass growth. This is contrary to the classical competitive accretion model, where the massive stars are never hindered in growth by the low-mass stars in the cluster. Our findings show that the most significant differences between the formation of low-mass and high-mass stars are all explained as the result of rapid accretion within a dense, gravitationally unstable, ionized flow.
Massive stars influence the surrounding universe far out of proportion to their numbers through ionizing radiation, supernova explosions, and heavy element production. Their formation requires the collapse of massive interstellar gas clouds with accretion rates exceeding M yr (Beuther et al., 2002; Beltrán et al., 2006) to reach their final masses before exhausting their nuclear fuel (Keto & Wood, 2006). Massive stars can ionize the gas around them, forming H ii regions, that traditionally have been modelled as spherical bubbles expanding in a uniform medium (Spitzer, 1978). This approach, however, fails to explain their observed numbers and morphologies (Wood & Churchwell, 1989; Kurtz et al., 1994). A more recent alternative picture models H ii regions as the ionized section of the accretion flow that feeds the young massive star (Keto, 2002, 2007). The rich diversity of observed H ii regions thus bears witness to the complexity of high-mass star formation.
H ii regions form around accreting protostars once they exceed M, equivalent to a spectral type of early B. Thus, accretion and ionization must occur together in the formation of massive stars. The pressure of the 10 K ionized gas far exceeds that in the 10 K accreting molecular gas, producing feedback in the form of ionized outflows (Keto, 2002, 2003, 2007).
High-mass stars form in denser and more massive cloud cores (Motte et al., 2008) than their low-mass counterparts (Myers et al., 1986). High densities also result in local gravitational instabilities in the accretion flow, resulting in the formation of multiple additional stars (Klessen & Burkert, 2000; Kratter & Matzner, 2006). Young massive stars are almost always observed to have companions (Ho & Haschick, 1981), and the number of their companions significantly exceeds those of low-mass stars (Zinnecker & Yorke, 2007). Such companions influence subsequent accretion onto the initial star (Krumholz et al., 2009a). Observations show an upper mass limit of about M. It remains unclear whether limits on internal stability or termination of accretion by stellar feedback determines the value of the upper mass limit (Zinnecker & Yorke, 2007).
Around the most luminous stars, the outward radiation pressure can counterbalance the inward gravitational attraction. A spherically symmetric calculation of radiation pressure on dust yields equality at just under 10 M (Wolfire & Cassinelli, 1987). However, the dust opacity is wavelength dependent, the accretion is non-spherical, the mass-luminosity ratio is different for multiple companions than for a single star, and the momentum of the accretion flow must be reversed, rather than just achieving static force balance (Larson & Starrfield, 1971; Kahn, 1974; Yorke & Krügel, 1977; Nakano et al., 1995; Sigalotti et al., 2009). Observations (Chini et al., 2004; Patel et al., 2005; Keto & Wood, 2006; Beltrán et al., 2006) provide evidence for the presence of all these mitigating factors, and numerical experiments combining some of these effects (Yorke & Sonnhalter, 2002; Krumholz et al., 2007b) confirm their effectiveness, showing that radiation pressure is not dynamically significant below the Eddington limit.
The most significant differences between massive star formation and low-mass star formation seem to be the clustered nature of star formation in dense accretion flows and the ionization of these flows. We present the first three-dimensional simulations of the collapse of a molecular cloud to form a cluster of massive stars that include ionization feedback, allowing us to study these effects simultaneously. During recent years the problem of massive star formation was tackled in a number of three-dimensional numerical studies (see e.g. Klessen et al., 2009, for a recent review). However, none of these simulations included the self-consistent ionization feedback from massive protostars that dynamically condense out of gravitationally unstable regions. Dale et al. (2005) performed the first hydrodynamic calculation of a star cluster forming region to include a simplified treatment of ionization feeback, but only from a single source that was inserted into a pre-existing simulation and whose ionizing flux did not depend on the mass of the growing (proto-)star. Also, heating by the accretion luminosity and from the stellar photosphere is neglected in these calculations, but is included in ours. Other calculations that incorporate ionizing radiation focus on triggered star formation and turbulence by external sources (Dale et al., 2007a, b; Gritschneder et al., 2009). Using a flux-limited diffusion approximation, Krumholz et al. (2007b) incorporated radiation feedback in their calculations of massive star formation. These calculations showed that radiation pressure from non-ionizing radiation cannot halt accretion onto the massive protostar and therefore is not able to set an upper mass limit of such stars (see also Krumholz et al., 2009a). Nevertheless, these simulations lack the important feedback from ionizing radiation distinctive of massive star formation.
In the following, we present our numerical method (§ 2) and give a detailed analysis of the simulation results. We discuss the accretion history of the stellar cluster (§ 3), the fragmentation of the rotationally flattened structure, and the generation of ionization-driven bipolar outflows (§ 4), and compare synthetic radio continuum and line observations with observed data (§ 5). Finally, we discuss the relation of our numerical model to previous work (§ 6), and the relevance of our findings for massive star formation (§ 7).
2 Numerical Method
We use a modified version of the FLASH code (Fryxell et al., 2000), an adaptive-mesh code that integrates the compressible gas dynamic equations with self-gravity and radiation feedback. Our modifications include the propagation of both ionizing and non-ionizing radiation to follow heating, though not radiation pressure, a refinement criterion to resolve the Jeans length (Banerjee et al., 2004), various cooling processes of relevance for protostellar collapse (Banerjee et al., 2004; Banerjee & Pudritz, 2006; Banerjee et al., 2006); and Lagrangian sink particles (Banerjee et al., 2009; Federrath et al., 2009). Here we report on our new modifications to the code.
To propagate radiation, we use the hybrid characteristics raytracing module originally developed by Rijkhorst et al. (2006). Parallel raytracing on a block-structured adaptive mesh requires the solution of an inverse problem: Given the position in the domain, the corresponding block identifier is needed. This problem was originally approached by storing the block identifier in one single large array corresponding to a fully refined domain. Each entry in this array could be mapped to a (potential) cell in the domain and vice versa, so that the block identifier could be obtained by looking up the stored value. This solution is straightforward and fast, but it violates the principles of adaptive mesh refinement and becomes extremely memory consuming at high resolution. Collapse simulations are impossible using this method.
Instead of creating one large array, we use the fact that FLASH stores hierarchical information about its block structure. Every block in the adaptive-mesh hierarchy has information about its parent and child blocks as well as on its neighbors at the same level of refinement. Valid data is stored for blocks of highest refinement only, which are called leaf blocks. Neighboring leaf blocks can differ only by one level of refinement. This hierarchical information makes it possible to perform a tree walk in the adaptive-mesh hierarchy.
Starting from a block whose identifier is known, we determine the direction in which the ray leaves the block. If there is a neighboring block on the same level of refinement, it can either still be a leaf block or have one level of children. In the former case the new block is already found, in the latter case the child blocks must be checked. If there is no neighboring block at the same refinement level, the new block must be the neighbor of the parent of the current block in this direction.
Although the method is more complicated than the original one, it is equally fast since the hierarchical data can be effectively stored and communicated among processors.
The radiation module has been tested for its capability of accurately casting shadows (Rijkhorst et al., 2006) and tracing R-type ionization fronts in a cosmological setting (Iliev et al., 2006). For simulations of massive star formation, the expansion of D-type ionization fronts should be adequately modeled. We have tested the ionization physics in our model against an approximate solution found by Spitzer (1978) for the expansion of a D-type ionization front into a homogeneous medium. Given the Strömgren radius and the sound speed of the ionized gas , the radius of the ionization front at time is given by
A comparison of our numerical simulation with this analytic solution is plotted in Fig. 1. We show the data from two different simulations, both of which had initial density g cm, ionizing luminosity s, and were run with a maximum resolution of pc. In the first run, we used an isothermal equation of state () and included no cooling processes. In the second, we set and used the cooling curve from Dalgarno & McCray (1972). In both simulations, the ionization fraction and temperature of the ionized gas was determined self-consistently by solving a rate equation for the ionization fraction as discussed in Rijkhorst et al. (2006).
We compare the numerical result in each case with the result from equation (1) using the appropriate sound speed in the ionized gas , where is Boltzmann’s constant and is the mean mass per particle. The agreement of each run with the analytical solution and with each other is acceptable. The decent agreement even in the presence of cooling demonstrates that numerical diffusion from the cold, dense, neutral shell into the hot, ionized interior does not lead to significant unphysical cooling of the hot interior.
We show a slice of the expanding H ii region in figure 2. The whole box shown has an effective numerical resolution of grid points. Small deviations from a perfectly spherical shape are visible. This is because the sphere is initially not very well resolved, so that dynamical processes inside the H ii region can lead to an amplification of grid effects. Resolution studies in two dimensions that we have performed show that this effect diminishes with increasing resolution.
To trace the collapse of individual fragments, we use a sink particle technique that we developed for the FLASH code (Federrath et al., 2009) (see also Bate et al. (1995)). Here, sink particles are created if the local density exceeds a critical value and the region within the sink particle radius is gravitationally bound and collapsing. For the high resolution simulations (Run A and Run B), we use g cm and AU, for the lower resolution simulations (Run Ca and Run Cb), we take g cm and AU. The sink particles gain mass through accretion of overdense gas within the accretion radius that is gravitationally bound to it. The particles interact gravitationally with the gas and with themselves and are free to move within the simulation box, independent of the underlying grid structure, i.e. they are treated fully Lagrangian.
We use the properties of these Lagrangian sink particles as sources for the radiation module. A prestellar model is used to determine luminosity and effective temperature of the source as function of protostellar mass and accretion rate. We set the stellar luminosity by a zero-age main sequence (Paxton, 2004) and the accretion luminosity by interpolation from a more detailed model (Hosokawa & Omukai, 2009).
The photoionization rate and the photoionization heating rate are set by the specific mean intensity along a ray
Here, is the distance from the source, the radius of the star, the effective temperature of the star, the frequency, the speed of light, Planck’s constant, Boltzmann’s constant and the optical depth for the ionizing radiation, which is calculated using the absorption coefficient of atomic hydrogen (Rijkhorst et al., 2006). Hence, the strength of the ionizing radiation is totally determined by and . We use a zero-age main sequence (ZAMS) model derived from the freely available stellar evolution code EZ (Paxton, 2004) to calculate these quantities.
To also include heating by non-ionizing radiation, we have extended the raytracing module to calculate Planck mean opacities following the Krumholz et al. (2007b) interpolation of data from Pollack et al. (1994). We use the non-relativistic approximation for the dust heating term (see e.g., Krumholz et al., 2007a)
with the Planck mean opacity , the gas density and the total radiation energy . In the raytracing approximation, the heating due to stellar radiation is thus given by
with the Stefan-Boltzmann constant and the Planck mean opacity .
In addition to the stellar heating, we also include a model for accretion heating. To this end, we assume that the potential energy of the gas which is released during accretion is fully converted into radiation at the accretion radius , so that the accretion luminosity is given by
with Newton’s constant , the protostellar mass and the accretion rate . Since we do not account for diffuse radiation, the accretion heating also originates from the sink particles and leads to another heating term of the form
To determine the accretion radius and the effective temperature of the accretion luminosity, we interpolate the results of the prestellar evolution model of Hosokawa & Omukai (2009) from the current accretion rate and protostellar mass.
2.2 Initial Conditions
We start our calculations with a molecular cloud with a mass of M and an initial temperature of K. The cloud has a flat inner region with pc radius, surrounded by a region in which the density drops as . The core density is g cm. The whole simulation box has an edge length of pc. The core is in solid body rotation with a ratio of rotational to gravitational energy . We then run the simulation allowing sink particles to continue forming for kyr after the formation of the first one. Note that the total simulation time is Myr.
These initial conditions represent an extrapolation from lower mass prestellar cores, which are well described by density profiles similar to Bonnor-Ebert spheres (see e.g., Pirogov, 2009), towards higher masses and larger scales. After a simulation runtime of kyr, the flat inner part has totally disappeared, and a centrally condensed structure with a peak density of g cm and a radius of pc has emerged. This is the stage at which infrared dark clouds are typically observed (see e.g., Beuther & Henning, 2009).
In some runs, we choose to suppress secondary sink formation and instead use a density-dependent temperature floor to prevent runaway collapse of dense blobs of gas. We need to resolve the local Jeans length with cells (Truelove et al., 1997). To do so, we introduce a dynamical temperature floor
where is the mean molecular weight, the proton mass, and the cell size.
We consider two full resolution simulations with a cell size at the highest refinement level of 98 AU (11 refinement levels). In one, Run A, we artificially suppress disk fragmentation using the dynamical temperature floor, while in the other, Run B, we permit it and dynamically form secondary sink particles. In addition, we ran two half resolution simulations with maximum resolution of 196 AU (10 refinement levels) and suppressed secondary sink formation. One of these runs, Run Ca, starts with identical initial conditions to Runs A and B, while Run Cb has an additional -perturbation of of the gas density (Boss & Bodenheimer, 1979). The high-resolution stellar cluster simulation Run B reaches grid cells by the end of the simulation, and required CPU hours to complete. An overview of the runs with their different model parameters is given in Tab. 1.
|Run A||98 AU||no|
|Run B||98 AU||yes|
|Run Ca||196 AU||no|
|Run Cb||196 AU||no|
We first examine the accretion histories of our different models. Figure 3 shows that in Run A, with only one sink particle allowed to form, nothing halts accretion onto the central protostar. It continues to grow at an average rate of Myr until we stop the calculation when the star has reached M. The increasingly massive star ionizes the surrounding gas, raising it to high pressure. This gas breaks out above and below the disk plane, but it cannot halt mass growth through the disk mid-plane .
We ran our lower resolution simulations with a single sink particle, Runs Ca and Cb, until the sink particles reached masses of M and M without any reduction in the mass accretion rate. Figure 4 shows the accretion history for these simulations. One can clearly see that neither the additional perturbation nor the change in resolution has any significant influence on the overall accretion behavior. Runs Ca and Cb accrete at mean rates of Myr and Myr, almost identical to the value for Run A, showing that this quantity is well converged. These simulations with a single sink particle strongly suggest that ionization feedback alone does not limit massive star accretion.
In contrast, in Run B, where multiple sink particles are allowed to form, two subsequent sink particles form and begin accreting soon after the first one, and many more follow within the next yr (see Fig. 3). By that time the first sink has accreted M. Within another yr seven further fragments have formed, with masses ranging from M to M while the first three sink particles have masses between M and M, all within a radius of pc from the most massive object. Accretion by the secondary sinks terminates the mass growth of the central massive sinks. Material that moves inwards through the disk driven by gravitational torques accretes preferentially onto the sinks at larger radii (Bate, 2002). Eventually, hardly any gas makes it all the way to the center to fall onto the most massive objects. This fragmentation-induced starvation prevents any star from reaching a mass M in this case.
The accretion behavior in Run B contrasts sharply with competitive accretion models (Bonnell et al., 2001, 2004), which have no mechanism to turn off accretion onto the most massive stars. In these models, material falls all the way down to the massive stars sitting in the center of the gravitational potential which thereby take away the gas from the surrounding low-mass stars. In our fragmentation-induced starvation scenario, exactly the opposite happens.
The different behaviour of the multiple sink simulation (Run B) and the single sink simulation (Run A) can be understood by looking at slices of density in the disk plane, as shown in Figures 5 and 6. In Run B, secondary sink particles form in the disk filaments, which remove gas from the common reservoir. Consequently, the most massive stars cannot maintain their high accretion rates. It is this halting of the accretion flow that allows the ionizing radiation to create a bubble around the massive protostars. In contrast, the sink particle in Run A is embedded in an accretion flow at all times. The ionizing radiation does blow away a significant fraction of gas from the sink particle, but it is not able to do so in all directions. There is always a region in the disk midplane around the sink particle where the gas is dense enough that evaporation by the ionizing radiation remains insufficient to stop the accretion flow.
4 Bipolar Outflows
In all runs, the H ii regions are trapped in the disk plane but drive a bipolar outflow perpendicular to the disk. The highly variable rate of accretion onto protostars as they pass through dense filaments causes fast ionization and recombination of large parts of the interior of the perpendicular outflow. The H ii regions around the massive protostars do not uniformly expand, but instead fluctuate rapidly in size, shape and flux.
The differences between the single and multiple fragment Runs A and Run B are also evident from the characteristics of the bipolar outflows. Figures 7 and 8 show density slices perpendicular to the disk plane for both simulations. Figure 7 shows the outflows driven by the most massive sink particle in Run B. The efficient shielding by the filaments and the motion of the sink particle hinder the thermal pressure of the ionized gas from driving a symmetric bipolar outflow at early times. The shielding of the ionizing radiation leads to a drop in the thermal pressure which drives and supports the outflow, and thus the outflow can be quenched again. At later times, the strong accretion flow onto the sink particle stops, and this allows the sink particle to drive a much larger outflow, which in the end has the form of a bubble.
The situation for the sink particle in Run A is presented in Figure 8. The general properties of the outflow are very similar, but the big difference is that here the accretion flow never stops. On the contrary, Figure 3 shows that the accretion rate increases towards the end of the simulation. This strong accretion flow fully absorbs the ionizing radiation, so the sink particle is unable to build up a bubble, even though it has more than three times the mass of the most massive star in Run B.
Observed outflows (Arce et al., 2007) from O-type protostars with ages yr tend to be much wider than those from low- or even intermediate-mass protostars. When resolved, these outflows show a disordered appearance that may be the outcome of multiple powering sources (including H ii regions) and non-steady driving. The pressure-driven outflows produced by H ii regions in our simulations reproduce these broad, disordered outflows well. Collimated outflows from less massive protostars without H ii regions are likely magnetically driven, though (see e.g., Banerjee & Pudritz, 2006). A more detailed model linking these two types of outflows is out of the scope of this paper.
5 Comparison to Observations
We can directly compare our models with radio observations of free-free continuum, hydrogen recombination lines, and NH inversion lines by generating synthetic maps.
We generate radio continuum maps from our numerical data by integrating the radiative transfer equation in the Rayleigh-Jeans limit while neglecting scattering (Kraus, 1966; Gordon & Sorochenko, 2002). To this end, we calculate for every cell in the domain the free-free absorption coefficient for atomic hydrogen
where is the number density of free electrons, the gas temperature and the frequency. The optical depth at distance is then
Without scattering, the radiative transfer equation for the brightness temperature can be directly integrated and yields
If is the solid angle subtended by the beam, the resulting flux density is
Following the algorithm described in Mac Low et al. (1991), we convolve the resulting image with the beam width and add some noise according to the telescope parameters.
For a time-series comparison at high temporal resolution, we select a few time intervals in the multi-sink run (Run B) during which the continuum emission from a particular, well-isolated, H ii region varies strongly. We then produced 2 cm continuum maps every yr from the simulation output and measured the properties of the region, including the observational scale length , the flux , and the accretion rate of the protostar that powers the H ii region. is defined as the equivalent diameter of a circle with the same area as the emission. is obtained by integrating the intensity of the maps over solid angle. is the accretion rate of the sink particle.
We generate synthetic observations of line emission using our radiative transfer code MOLLIE (Keto, 1990). This code generally uses the Accelerated Lambda iteration algorithm (Rybicki & Hummer, 1991) although in this investigation we compute the level populations of NH assuming local thermodynamic equilibrium (LTE) and assume optically thin conditions for the H53 recombination line. The computation is done on 4 nested grids with spatial scales ranging from 80 to 640 AU and covering 0.2 pc. The maps of the NH and H53 velocities show the first moment, in the case of NH integrating over the main hyperfine line of the (3,3) transition.
Our simulated observations of radio continuum emission reproduce the morphologies reported in surveys of ultracompact H ii regions (Wood & Churchwell, 1989; Kurtz et al., 1994). Figure 9 shows typical images from Run B, and the Online Material contains a movie of the entire run. The H ii regions in our model fluctuate rapidly between different shapes while accretion onto the protostar continues. When the gas reservoir around the two most massive stars is exhausted, their H ii regions merge into a compact H ii region similar to the extended envelopes often observed around ultracompact H ii regions (Kim & Koo, 2001; Kurtz et al., 1999). Our results suggest that the lifetime problem of ultracompact H ii regions (Wood & Churchwell, 1989) is only apparent. Since H ii regions embedded in accretion flows are continuously fed, and since they flicker with variations in the flow rate, their size does not depend on their age until late in their lifetimes.
The same mechanisms that cause the different accretion behavior (§ 3) and bipolar outflows (§ 4) for Run A and Run B also lead to differences in the H ii region morphologies. One the one hand, the permanent accretion flow in Run A quenches the H ii region in the disk plane. On the other hand, the higher stellar mass in Run A leads to larger H ii regions perpendicular to the disk plane. Since the second effect dominates, the bottom line is that the H ii regions in Run A are generally larger than in Run B. Hence, multiple sink formation is not only necessary to model the gas fragmentation correctly, it is also crucial to reproduce the large relative fraction of spherical H ii regions observed in surveys (Wood & Churchwell, 1989; Kurtz et al., 1994). A detailed investigation of the H ii region morpholgies, including a comprehensive statistical analysis for both simulations, will be presented in a follow-up paper.
To compare more directly to observations of the time variability of H ii regions, we analyzed a few time intervals of interest at a time resolution of yr. We find that when the accretion rate to the star powering the H ii region has a large, sudden increase, the ionized region shrinks, and then slowly re-expands. This agrees with the contraction, changes in shape, or anisotropic expansion observed in radio continuum observations of ultracompact H ii regions over intervals of 10yr (Franco-Hernández & Rodríguez, 2004; Rodríguez et al., 2007; Galván-Madrid et al., 2008). Figure 10 shows that the sudden accretion of large amounts of material is accompanied by a fast decrease in the observed size and flux of the H ii region. In the left panel, the H ii region is initially relatively large, and accretion is almost shut off. A large (from 0 to yr), sudden accretion event causes the H ii region to shrink and decrease in flux. The star at this moment has a mass of M. In the right panel, the star has a larger mass (M), the H ii region is initially smaller, and the star is constantly accreting gas. The ionizing-photon flux appears to be able to ionize the infalling gas stably, until a peak in the accretion rate by a factor of three and the subsequent continuous accretion of gas makes the H ii region shrink and decrease in flux. The H ii region does not shrink immediately after the accretion peak because the increase is relatively mild and the geometry of the infalling gas permitted ionizing photons to escape in one direction. Our results show that observations of large, fast changes in ultracompact H ii regions (Franco-Hernández & Rodríguez, 2004; Rodríguez et al., 2007; Galván-Madrid et al., 2008) are controlled by the accretion process. Both the scale length and flux decrease at rates of 5–7 % per year in our model, agreeing well with observed fluctuations of 2–9 % per year (Franco-Hernández & Rodríguez, 2004; Galván-Madrid et al., 2008). Shortly after the minimum values are reached, the H ii regions re-expand, on timescales yr.
We can compare our models to a well-studied example of an ultracompact H ii region, W51e2 (Zhang et al., 1998; Keto & Klaassen, 2008). Although we did not choose parameters specifically to model this region, the comparison is nevertheless revealing. In Figure 11 we show simulated and observed maps of NH emission, 1.3 cm thermal continuum emission, and the H53 radio recombination line. The simulated maps were made at a time when a 20 M star has formed in Run B.
In the simulated observations, the origin of the coordinates is referenced to the location of the 20 M protostar at the center of the accretion flow as marked. The accretion flow and disk are viewed edge-on. The spatial scale in the observations (lower panels) is 5100 AU per arc sec assuming W51e2 is at a distance of 5.1 kpc (Xu et al., 2009). The color bar on the right of each figure shows the molecular velocities in km s. The velocities of the observations include the LSR velocity of W51e2, approximately 57 kms. The white contour levels are 50% through 90% of the peak brightness temperature of 71 K (upper left) and 0.1, 0.2 and 0.3 Jy beam kms (lower left). The red contour levels are 30%, 70%, and 95% of the peak 1.3 cm continuum brightness temperature of 8902 K (upper left) and 0.07, 0.14, and 0.22 Jy/beam (lower left). The white contours on the H53 observations (lower-right) show the 7 mm continuum emission at 2, 4, 10, 30, 50, 70, and 90% of the peak emission of 0.15 Jy beam. The molecular line observations are from Zhang et al. (1998). The H53 observations are from Keto & Klaassen (2008). Both observations have a beamsize of about 1. Coordinates are in the B1950 epoch.
The brightest NH emission indicates the dense accretion disk surrounding the most massive star in the model, one of several within the larger-scale rotationally flattened flow. The disk shows the signature of rotation, a gradient from redshifted to blueshifted velocities across the star. A rotating accretion flow is identifiable in the observations, oriented from the SE (red velocities) to the NW (blue velocities) at a projection angle of 135 east of north (counterclockwise).
The 1.3 cm radio continuum traces the ionized gas, which in the model expands downward perpendicularly to the accretion disk, down the steepest density gradient. As a result, the simulated map shows the brightest radio continuum emission just off the mid-plane of the accretion disk, separated from the central star, rather than surrounding it spherically. In the observations of W51e2, continuum emission is indeed offset from the accretion disk traced in ammonia. The NH(3,3) in front of the H ii region is seen in absorption and red-shifted by its inward flow toward the protostar. The density gradient in the ionized flow determines the apparent size of the H ii region. Therefore the accretion time scale determines the age of the H ii region rather than the much shorter sound-crossing time of the H ii region.
Photoevaporation of the rotationally flattened molecular accretion flow supplies the ionized outflow. Therefore, the ionized gas rotates as it flows outward, tracing a spiral. An observation that only partially resolves the spatial structure of the ionized flow sees a velocity gradient oriented in a direction between that of rotation and of outflow, as shown in the simulated observation. The observed H53 recombination line (Keto & Klaassen, 2008) in Figure 11 indeed shows a velocity gradient oriented between the directions of rotation and the outflow.
In this work, we focus on the impact of feedback from radiative heating by both ionizing and non-ionizing radiation during massive star formation. Our work is complementary to recent simulations that neglect ionizing radiation but include other effects.
For example, we use a ray-tracing algorithm to describe ionizing and non-ionizing radiation while Yorke & Sonnhalter (2002); Krumholz et al. (2007b, 2009b); Price & Bate (2009) rely on flux-limited diffusion to follow the radiation field. Like us, these models include heating from accretion luminosity. They neglect heating by ionizing radiation, but do include radiation pressure and scattering effects. Radiation pressure in particular becomes important at the 10 AU scale. These models therefore simulate rather smaller scales than we do, starting from a core of only pc radius, as opposed to our more massive clump with and a pc radius. As a result, our calculations are better geared towards describing the cluster-forming scale as opposed to modeling in detail the mass growth of individual protostars through their inner accretion disks. Yorke & Sonnhalter (2002); Krumholz et al. (2007b, 2009b); Sigalotti et al. (2009) show that radiation pressure cannot stop accretion onto the massive protostar and is dynamically unimportant except for massive star clusters near the Galactic center or in starburst galaxies (Krumholz & Matzner, 2009). Krumholz et al. (2007b) initialized a microturbulent flow within their dense core, but Krumholz et al. (2009b) found that it made little difference by the time material had collapsed to the size of the disk. In addition, Price & Bate (2009) include magnetic fields and find that they provide some support for low-density gas, thus reducing the amount of matter available for accretion. Banerjee & Pudritz (2006), Machida et al. (2005), Hennebelle & Teyssier (2008), Hennebelle & Fromang (2008) and Hennebelle & Ciardi (2009) furthermore show that the presence of magnetic fields can dramatically alter the shape and size of accretion disks and strongly suppress the formation of close binary systems. These last results, may, however, be an artifact of the neglect of turbulent and ambipolar diffusion within the disks in these models. Dale & Bonnell (2008) and Wang et al. (2009), on the other hand, focused on mechanical feedback through winds and bipolar outflows, without and with magnetic fields respectively, but neglect radiative heating entirely. They find that outflows further limit accretion, but do not fundamentally change the picture of cluster formation in a gravitationally unstable accretion flow. Ultimately, the quantitative characterization of the star formation efficiency will require treating the effects of radiation, magnetic fields, and outflows simultaneously.
As noted above, we have neglected a number of physical processes, including radiation pressure, scattering, magnetic fields, winds and bipolar outflows, and turbulent initial conditions. In addition, we here discuss other approximations that we have made.
Our prestellar model is designed such that the number of UV photons emitted by the protostar is overestimated for three reasons. First, we use a ZAMS model to set the ionizing luminosity, while the models of Hosokawa & Omukai (2009) suggest that the protostellar radii of massive stars can be larger than the ZAMS value, leading to a smaller ionizing luminosity compared to the ZAMS value. Second, we neglect the aborption of UV photons by dust inside the H ii region. Thus, the sizes of our H ii regions represent upper limits to the real values. Third, the dense inner disk around each protostar is subsumed in the sink particle, allowing radiation to propagate along the disk midplane to the edge of the sink particle, where densities will be lower.
However, as discussed in § 3, we find that, although we overestimate the strength of the ionizing radiation, it is unable to stop accretion onto a single massive star. Similarly, in § 5 we showed that the strength of the accretion flow restricts the expansion of the H ii region, which solves the lifetime problem of ultracompact H ii regions (Wood & Churchwell, 1989). Both results are robust in the sense that the impact of the ionizing radiation on the accretion flow and the sizes of the H ii regions would even be smaller if we used a more detailed modeling of the protostellar evolution and environment, and included dust absorption.
We also note that, although we nominally use a simplified treatment of radiative losses in the optically thick regime, this will not have an impact on the dynamics of the simulated system, as we expect that the gas will become optically thick at densities of g cm (e.g. Larson, 1969), which is two orders of magnitude above the threshold density for sink creation.
Our results suggest that the accretion flow onto massive stars cannot be stopped by ionizing radiation (§ 3). Instead, the growth of massive stars is halted by competition from lower-mass companions that form by gravitational fragmentation in the dense, rotationally flattened, accretion flow. This process, which we call fragmentation-induced starvation, limits the maximum stellar mass of the highest-mass star in our simulations.
Our models lead to a self-consistent picture of massive star formation where massive stars form through (unsteady) mass growth within a gravitationally unstable accretion flow. Radiation feedback from the accretion luminosity and the protostellar surface heats up the environment in which the first massive protostar forms (Krumholz et al., 2007b). This increases the local Jeans mass, so that subsequent fragments in the neigborhood of this star are relatively massive, naturally resulting in a cluster of massive stars.
Unlike in the competitive accretion model (Bonnell & Bate, 2006) where the most massive protostar continues to accrete from the common gas reservoir, in our model accretion is shut off by the subsequent fragmentation of the gas because the lower-mass companions intercept material that would otherwise accrete onto the central star, starving it of material. Fragmentation-induced starvation will only occur in regions with high enough density for many massive fragments to form despite accretion heating. In our models, the most massive star starts to form early and completes accretion while the whole cluster is still growing, contrary to standard competitive accretion models in which the most massive star grows over the whole cluster evolution timescale. On the other hand, competition for a common reservoir continues to define the top end of the initial mass function, unlike in models that rely on the clump mass spectrum to define the stellar initial mass function.
Our models show (§ 5) that the H ii regions formed around massive young stars flicker but do not grow systematically while heavy accretion from the protostellar core continues. This behavior seems able to explain the many puzzles raised by observations of ultracompact H ii regions. The chaotic interaction between gravitationally unstable gas filaments in the accretion disk and ionizing radiation from the young stars creates a great variety of different ultracompact H ii region morphologies. We find shell-like, core-halo, cometary, spherical and irregular morphologies in a single simulation, depending on the current flow field and viewing angle. Our models show behavior consistent with observations of shrinking H ii regions with changes in size or flux of a few percent per year, as a consequence of denser material being accreted the star and shielding the ionizing radiation. Shielding by the disk also produces uncollimated bipolar outflows in our models (see § 4). These results lead to an evolution scenario for H ii regions that is substantially different from the classical picture of monotonically expanding hot bubbles. The surprising dynamics of H ii region flickering resolves the lifetime problem for ultracompact H ii regions, since the size of the H ii region remains unrelated to the age of the protostar, as long as fast accretion continues. A detailed comparison of simulated observations of our model with the ultracompact H ii region W51e2 demonstrates that our model reproduces distinctive features of massive star forming regions such as H ii regions expanding asymmetrically away from the protostar down the steepest density gradient, and outward-twisting, spiraling, ionized flows. All these observations can be understood as consequences of ionization within a gravitationally unstable accretion flow.
- affiliation: Centro de Radioastronomía y Astrofísica, UNAM, A.P. 3-72 Xangari, Morelia 58089, Mexico
- affiliation: Academia Sinica Institute of Astronomy and Astrophysics, P.O. Box 23-141, Taipei 106, Taiwan
- Arce, H. G., Shepherd, D., Gueth, F., Lee, C.-F., Bachiller, R., Rosen, A., & Beuther, H. 2007, in Protostars and Planets V, ed. B. Reipurth, D. Jewitt, & K. Keil (The University of Arizona Press), 245–260
- Banerjee, R. & Pudritz, R. E. 2006, Astrophys. J., 641, 949
- Banerjee, R., Pudritz, R. E., & Anderson, D. W. 2006, Mon. Not. R. Astron. Soc., 373, 1091
- Banerjee, R., Pudritz, R. E., & Holmes, L. 2004, Mon. Not. R. Astron. Soc., 355, 248
- Banerjee, R., VVázquez-Semadeni, E., Hennebelle, P., & Klessen, R. S. 2009, Mon. Not. R. Astron. Soc., 398, 1082
- Bate, M. R. 2002, Mon. Not. R. Astron. Soc., 314, 33
- Bate, M. R., Bonnell, I. A., & Price, N. M. 1995, Mon. Not. R. Astron. Soc., 277, 363
- Beltrán, M. T., Cesaroni, R., Codella, C., Testi, L., Furuya, R. S., & Olmi, L. 2006, Nature, 443, 427
- Beuther, H. & Henning, T. 2009, Astron. Astrophys., 503, 859
- Beuther, H., Schilke, P., Sridharan, T. K., Menten, K. M., Walmsley, C. M., & Wyrowski, F. 2002, Astron. Astrophys., 383, 892
- Bonnell, I. A. & Bate, M. R. 2006, Mon. Not. R. Astron. Soc., 370, 488
- Bonnell, I. A., Bate, M. R., Clarke, C. J., & Pringle, J. E. 2001, Mon. Not. R. Astron. Soc., 323, 785
- Bonnell, I. A., Vine, S. G., & Bate, M. R. 2004, Mon. Not. R. Astron. Soc., 349, 735
- Boss, A. P. & Bodenheimer, P. 1979, Astrophys. J., 234, 289
- Chini, R., Hoffmeister, V., Kimeswenger, S., Nielbock, M., Nürnberger, D., Schmidtobreick, L., & Sterzik, M. 2004, Nature, 429, 155
- Dale, J. E. & Bonnell, I. A. 2008, Mon. Not. R. Astron. Soc., 391, 2
- Dale, J. E., Bonnell, I. A., Clarke, C. J., & Bate, M. R. 2005, Mon. Not. R. Astron. Soc., 358, 291
- Dale, J. E., Bonnell, I. A., & Whitworth, A. P. 2007a, Mon. Not. R. Astron. Soc., 375, 1291
- Dale, J. E., Clark, P. C., & Bonnell, I. A. 2007b, Mon. Not. R. Astron. Soc., 377, 534
- Dalgarno, A. & McCray, R. A. 1972, Ann. Rev. Astron. Astrophys., 10, 375
- Federrath, C., Clark, P. C., Banerjee, R., & Klessen, R. S. 2009, in preparation
- Franco-Hernández, R. & Rodríguez, L. F. 2004, Astrophys. J., 604, L105
- Fryxell, B., Olson, K., Ricker, P., Timmes, F. X., gale, M., Lamb, D. Q., MacNeice, P., Rosner, R., Truran, J. W., & Tufo, H. 2000, Astrophys. J. Suppl. Ser., 131, 273
- Galván-Madrid, R., Rodríguez, L. F., Ho, P. T. P., & Keto, E. 2008, Astrophys. J., 674, L33
- Gordon, M. A. & Sorochenko, R. L. 2002, Radio Recombination Lines (Kluwer Academic Publishers Group)
- Gritschneder, M., Naab, T., Walch, S., Burkert, A., & Heitsch, F. 2009, Astrophys. J., 694, L26
- Hennebelle, P. & Ciardi, A. 2009, Astron. Astrophys., 506, L29
- Hennebelle, P. & Fromang, S. 2008, Astron. Astrophys., 477, 9
- Hennebelle, P. & Teyssier, R. 2008, Astron. Astrophys., 477, 25
- Ho, P. T. P. & Haschick, A. D. 1981, Astrophys. J., 248, 622
- Hosokawa, T. & Omukai, K. 2009, Astrophys. J., 691, 823
- Iliev, I. T., Ciardi, B., Alvarez, M. A., Maselli, A., Ferrara, A., Gnedin, N. Y., Mellema, G., Nakamoto, T., Norman, M. L., Razoumov, A. O., Rijkhorst, E.-J., Ritzerveld, J., Shapiro, P. R., Susa, H., Umemura, M., & Whalen, D. J. 2006, Mon. Not. R. Astron. Soc., 371, 1057
- Kahn, F. D. 1974, Astron. Astrophys., 37, 149
- Keto, E. 2002, Astrophys. J., 580, 980
- —. 2003, Astrophys. J., 599, 1196
- —. 2007, Astrophys. J., 666, 976
- Keto, E. & Klaassen, P. 2008, Astrophys. J., 678, L109
- Keto, E. & Wood, K. 2006, Astrophys. J., 637, 850
- Keto, E. R. 1990, Astrophys. J., 355, 190
- Kim, K.-T. & Koo, B.-C. 2001, Astrophys. J., 549, 979
- Klessen, R. S. & Burkert, A. 2000, Astrophys. J. Suppl. Ser., 128, 287
- Klessen, R. S., Krumholz, M. R., & Heitsch, F. 2009, arxiv:0906.4452v1
- Kratter, K. M. & Matzner, C. D. 2006, Mon. Not. R. Astron. Soc., 373, 1563
- Kraus, J. D. 1966, Radio Astronomy (McGraw-Hill, Inc.)
- Krumholz, M. R., Klein, R. I., & McKee, C. F. 2007a, Astrophys. J., 667, 626
- —. 2007b, Astrophys. J., 656, 959
- Krumholz, M. R., Klein, R. I., McKee, C. F., Offner, S. S. R., & Cunningham, A. J. 2009a, Science, 323, 754
- —. 2009b, Science, 323, 754
- Krumholz, M. R. & Matzner, C. D. 2009, Astrophys. J., 703, 1352
- Kurtz, S., Churchwell, E., & Wood, D. O. S. 1994, Astrophys. J. Suppl. Ser., 91, 659
- Kurtz, S. E., Watson, A. M., Hofner, P., & Otte, B. 1999, Astrophys. J., 514, 232
- Larson, R. B. 1969, Mon. Not. R. Astron. Soc., 145, 271
- Larson, R. B. & Starrfield, S. 1971, Astron. Astrophys., 13, 190
- Mac Low, M.-M., Van Buren, D., Wood, D. O. S., & Churchwell, E. 1991, Astrophys. J., 369, 395
- Machida, M. N., Matsumoto, T., Hanawa, T., & Tomisaka, K. 2005, Mon. Not. R. Astron. Soc., 362, 382
- Motte, F., Bontemps, S., Schneider, N., Schilke, P., & Menten, K. M. 2008, in Astronomical Society of the Pacific Conference Series, Vol. 387, Massive Star Formation: Observations Confront Theory, ed. H. Beuther, H. Linz, & T. Henning, 22–29
- Myers, P. C., Dame, T. M., Thaddeus, P., Cohen, R. S., Silverberg, R. F., Dwek, E., & Hauser, M. G. 1986, Astrophys. J., 301, 398
- Nakano, T., Hasegawa, T., & Norman, C. 1995, Astrophys. J., 450, 183
- Patel, N. A., Curiel, S., Sridharan, T. K., Zhang, Q., Hunter, T. R., Ho, P. T. P., Torrelles, J. M., Moran, J. M., Gómez, J. F., & Anglada, G. 2005, Nature, 437, 109
- Paxton, B. 2004, Publ. Astron. Soc. Pac., 116, 699
- Pirogov, L. E. 2009, Astron. Rep., 53, 1127
- Pollack, J. B., Hollenbach, D., Beckwith, S., Simonelli, D. P., Roush, T., & Fong, W. 1994, Astrophys. J., 421, 615
- Price, D. J. & Bate, M. R. 2009, Mon. Not. R. Astron. Soc., 398, 33
- Rijkhorst, E.-J., Plewa, T., Dubey, A., & Mellema, G. 2006, Astron. Astrophys., 452, 907
- Rodríguez, L. F., Gómez, Y., & Tafoya, D. 2007, Astrophys. J., 663, 1083
- Rybicki, G. B. & Hummer, D. G. 1991, Astron. Astrophys., 245, 171
- Sigalotti, L. D. G., de Felice, F., & Daza-Montero, J. 2009, Astrophys. J., 707, 1438
- Spitzer, Jr., L. 1978, Physical Processes in the Interstellar Medium (John Wiley & Sons, Inc.)
- Truelove, J. K., Klein, R. I., McKee, C. F., Hollman II, J. H., Howell, L. H., & Greenough, J. A. 1997, Astrophys. J., 489, L179
- Wang, P., Li, Z.-Y., Abel, T., & Nakamura, F. 2009, arxiv:0908.4129v1
- Wolfire, M. G. & Cassinelli, J. P. 1987, Astrophys. J., 319, 850
- Wood, D. O. S. & Churchwell, E. 1989, Astrophys. J. Suppl. Ser., 69, 831
- Xu, Y., Reid, M. J., Menten, K. M., Brunthaler, A., Zheng, X. W., & Moscadelli, L. 2009, Astrophys. J., 693, 413
- Yorke, H. W. & Krügel, E. 1977, Astron. Astrophys., 54, 183
- Yorke, H. W. & Sonnhalter, C. 2002, Astrophys. J., 569, 846
- Zhang, Q., Ho, P. T. P., & Ohashi, N. 1998, Astrophys. J., 494, 636
- Zinnecker, H. & Yorke, H. W. 2007, Ann. Rev. Astron. Astrophys., 45, 481