The early impact of ionizing radiation on forming molecular clouds
As part of the SILCC-ZOOM project we present our first sub-parsec resolution radiation-hydrodynamic simulations of two molecular clouds self-consistently forming from a turbulent, multi-phase ISM. The clouds have similar initial masses of few , escape velocities of 5 km s, and a similar initial energy budget. We follow the formation of star clusters with a sink based model and the impact of radiation from individual massive stars with the tree-based radiation transfer module TreeRay. Photo-ionizing radiation is coupled to a chemical network to follow gas heating, cooling and molecule formation and dissociation. For the first 3 Myr of cloud evolution we find that the overall star formation efficiency is considerably reduced by a factor of 4 to global cloud values of 10 % as the mass accretion of sinks that host massive stars is terminated after 1 Myr. Despite the low efficiency, star formation is triggered across the clouds. Therefore, a much larger region of the cloud is affected by radiation and the clouds begin to disperse. The time scale on which the clouds are dispersed sensitively depends on the cloud substructure and in particular on the amount of gas at high visual extinction. The damage of radiation done to the highly shielded cloud (MC) is delayed. We also show that the radiation input can sustain the thermal and kinetic energy of the clouds at a constant level. Our results strongly support the importance of ionizing radiation from massive stars for explaining the low observed star formation efficiency of molecular clouds.
keywords:hydrodynamics; methods: numerical; ISM: kinematics and dynamics; ISM: clouds; stars: formation
Molecular clouds (MC) condense out of the diffuse, interstellar medium (ISM). These dense regions host filamentary substructures of molecular gas with an atomic envelope (André et al., 2014; Dobbs et al., 2014; Klessen & Glover, 2016). Massive stars form in infrared dark clouds, which are the densest parts of MCs (Goldreich & Kwan, 1974; Lada & Lada, 2003; Rathborne et al., 2006; Klessen, 2011; Ragan et al., 2012). During their lifetime, massive stars emit ionizing radiation and eject high-velocity winds, which result in the deposition of momentum, kinetic, and thermal energy in the ISM and change the chemical composition. The underlying physical processes are collectively termed stellar feedback, i.e. stellar winds (Castor et al., 1975; Weaver et al., 1977; Wünsch et al., 2011), ionizing radiation (Spitzer, 1978; Hosokawa & Inutsuka, 2006; Dale et al., 2012; Walch et al., 2012), radiation pressure (Krumholz & Matzner, 2009; Fall et al., 2010; Murray et al., 2010), and supernovae (SNe, Sedov 1958; Ostriker & McKee 1988; Walch & Naab 2015; Körtgen et al. 2016). Feedback modifies density structures, counteracts the gravitational collapse, interrupts mass accretion, and directly influences the cycle of star formation. However, the detailed impact on molecular cloud evolution is still a matter of discussion (Whitworth, 1979; Krumholz, 2006; Krumholz et al., 2009; Walch et al., 2012; Dale, 2015). It seems clear that stellar feedback can change the local and global multi-phase structure of the ISM with dramatic consequences for star formation (Naab & Ostriker, 2017).
MCs are complex. Observations indicate that they are embedded in their galactic environment (Mac Low et al., 2004) and coupled to large scale (some 100 pc) motions (Hughes et al., 2013; Colombo et al., 2014). Galactic turbulent velocity fields seem to be inherited (Brunt et al., 2009) with consequences for the star formation rate (Rey-Raposo et al., 2015). Hence, it is likely that the cloud properties are already imprinted during early formation and continuously reshaped by physical processes on galactic scales (Dobbs et al., 2012; Walch et al., 2015; Girichidis et al., 2016; Seifried et al., 2017; Rey-Raposo et al., 2017). This also suggests that the possible support of MCs by internal stellar feedback is highly variable and depends on the cloud structure. Analytical models which are usefully guiding our theoretical understanding, may not fully reflect the complexity of self-consistently evolving MCs (Matzner, 2002).
Early studies treat MCs in isolated environments to investigate gravitational collapse and implications for the star formation rate (Shu, 1977; Foster & Chevalier, 1993; Hetem & Lepine, 1993; Klessen et al., 2000; Dale et al., 2005; Gavagnin et al., 2017). Follow-up studies started to investigate connections to the surrounding ISM with idealised gas replenishing scenarios such as colliding flows (Heitsch et al., 2005; Vázquez-Semadeni et al., 2007, 2010) or cloud-cloud collisions (Whitworth et al., 1994; Inoue & Fukui, 2013; Dobbs et al., 2015; Balfour et al., 2015). In galactic-scale simulations (de Avillez & Breitschwerdt, 2005; Slyz et al., 2005; Joung & Mac Low, 2006; Hill et al., 2012; Kim et al., 2013; Hennebelle & Iffrig, 2014; Smith et al., 2014a; Walch et al., 2015; Dobbs, 2015; Girichidis et al., 2016) the statistical properties of MCs are analysed and the global importance of individual feedback processes estimated (Girichidis et al., 2016; Padoan et al., 2016; Gatto et al., 2017; Peters et al., 2017; Padoan et al., 2017; Kim & Ostriker, 2018). Recent progress in computational performance enables us to simulate a galactic-scale environment and simultaneously increase the spatial and time resolution in forming molecular clouds. This technique is referred to as a zoom-in simulation (Clark et al., 2012; Bonnell et al., 2013; Smith et al., 2014b; Dobbs, 2015; Butler et al., 2017; Ibáñez-Mejía et al., 2017; Kuffmeier et al., 2017; Pettitt et al., 2017; Seifried et al., 2017; Nordlund et al., 2017). The advantage is that large scale influences (e.g. SN shocks) are propagated down to MC-scales and cloud formation can be studied in a self-consistently evolved environment.
The impact of stellar feedback strongly depends on the mass of the star, hence the UV-luminosity (Geen et al., 2018), and its environment. A massive star with 23 M emits a factor of 100 less energy in a wind than it releases in radiative energy (Matzner, 2002) but higher/lower mass stars have stronger/weaker winds relative to radiation. Furthermore, stellar winds are inefficiently coupled to dense environment (Haid et al., 2018). Therefore, in massive MCs ( 10 M), the impact of stellar winds seems negligible (Dale et al., 2014; Geen et al., 2015; Ngoumou et al., 2015; Howard et al., 2017). However, in low-mass MCs ( 10 M), winds are able to reshape the clouds, ablate dense material, and even drive gas out of the clouds through low density channels (Rogers & Pittard, 2013). Winds are also more important than radiation if the environment of the massive star is already warm or hot, because in this case the radiation does not couple to the surrounding ISM and the radiative energy cannot be deposited in the gas (i.e. low coupling efficiency; Haid et al., 2018). Ionizing radiation also struggles to impact bound, massive MCs (Dale et al., 2012, 2013), while clouds with the sound speed of the photo-ionized gas being similar to the escape velocity can be dispersed completely within a few Myr (Walch et al., 2012). In any case, both processes shape the environment for the final SN explosions to leak out, thereby dispersing the clouds effectively (Harper-Clark & Murray, 2009; Pittard, 2013; Rosen et al., 2014; Gatto et al., 2017; Naab & Ostriker, 2017; Peters et al., 2017; Wareing et al., 2017).
The observed star formation in MCs is low with only a few percent of gas that is converted into stars during one free-fall time (Zuckerman & Evans, 1974; Evans et al., 2009; Murray, 2011). This inefficiency suggests that processes inside a cloud oppose the gravitational collapse. Stellar feedback is discussed to be an internal driver of supersonic turbulence with a velocity dispersion of a few km s (Mac Low & Klessen, 2004; Mac Low et al., 2004; Mellema et al., 2006; Walch et al., 2012). However, numerical simulations fail to reproduce this low level of star formation (Klessen et al., 2000; Vázquez-Semadeni et al., 2003; Dale et al., 2014).
Therefore, two aspects of the interaction of stellar feedback with the MC environment remain a matter of discussion. Is star formation limited to the low observed values of a few percent as a consequence of internal feedback processes? What is the role of MC substructure and filling factor on the coupling efficiencies of stellar winds and ionizing radiation (Haid et al., 2018)?
To address this questions, we present three-dimensional, radiation-hydrodynamic adaptive mesh refinement (AMR) simulations of MCs as part of a SN-driven, multi-phase ISM in a piece of a galactic disc (within the SILCC project; Walch et al., 2015; Girichidis et al., 2016). We apply a zoom-in technique to follow the formation and evolution of two MCs with total gas masses of a few 10 M with an effective resolution of 0.122 pc (Seifried et al., 2017). Sink particles are integrated by our novel predictor-corrector scheme (Dinnbier et al., in prep.). With a model of star cluster formation within sink particles, we couple ionizing radiation to the ambient medium (Haid et al., 2018). The radiation is treated by the novel, tree-based radiative transfer scheme TreeRay (Wünsch et al. 2018, in prep.) based on the tree solver for gravity and diffuse radiation implemented in Flash (Wünsch et al., 2018). For now, we neglect stellar winds, as their contribution in the early, dense phase of MCs is likely subordinate to ionizing radiation (Dale et al., 2014; Haid et al., 2018). We focus on the interplay of ionizing radiation and the particular MC morphology and star formation efficiency.
The paper is organized as follows. In Section 2, we present the numerical method. In Section 3, we give an overview of the simulation setup. We discuss the morphological impact of ionizing radiation in Section 4. We depict the effect of radiative feedback on the environment around the stellar component (Section 5). The evolution of global cloud properties is shown in Section 6 and the differences between the two clouds are discussed in Section 7. Finally, we conclude in Section 8.
2 Numerical method
We use the three-dimensional AMR magneto-hydrodynamics code (MHD) FLASH 4 (Fryxell et al., 2000; Dubey et al., 2008) with the directionally split, Bouchut HLL5R solver (Bouchut et al., 2007, 2010; Waagan, 2009; Waagan et al., 2011) including self-gravity, a chemical network to follow molecule formation and dissociation, the novel radiative transfer module TreeRay, sink particles, and the stellar evolution of massive stars.
2.1 Sink particles
Sink particles represent the unresolved formation of stars or clusters by gravitational collapse. In the simulations, we use a new particle module (Dinnbier et al., in prep.) which uses a Hermite predictor-corrector integrator and is coupled to the Barnes-Hut tree (Wünsch et al., 2018). The sink formation and accretion criteria are the same as in Federrath et al. (2010). In this work, sink particles represent star clusters (hereafter also simply called sinks) within which multiple massive stars (hereafter also stars) can form. For further information on the cluster sink implementation we refer to Gatto et al. (2017).
A sink particle can only be formed in a computational cell and followed through the computational domain if the harbouring cell lives on the highest refinement level (smallest spatial resolution) in the AMR hierarchy. The accretion radius is set to = = 0.31 pc. We further demand that the gas within is Jeans unstable, is in a converging flow and represents a local gravitational potential minimum (Federrath et al., 2010). Under the assumption of an isothermal gas with a temperature = 100 K, we derive the density threshold above which sinks can form, = 1.110 g cm following the Jeans criterion:
Sink particles accrete gas. A fraction of the accreted gas is turned into massive stars by means of the star cluster subgrid model. Assuming a Kroupa stellar initial mass function (IMF), one new massive star (9 M 120 M) is randomly sampled for every 120 M accreted on a sink (Kroupa, 2001). We assume the Salpeter slope of -2.35 in the high mass regime of the IMF (Salpeter, 1955). Each sink with a mass can contain stars with individual initial masses and individual stellar evolutions tracks (see Ekström et al. 2012, Gatto et al. 2017, Peters et al. 2017 and references therein). We refer to the number of massive stars, , in a sink as the active stellar component, , with . The residual gas is converted into low-mass stars, which are not recorded individually as they currently provide no feedback to the surrounding medium.
Each sink particle is subject to the gravitational attraction of the gas and the other sink particles. Their trajectories are integrated by a predictor-corrector scheme, which is inspired by the two nested fourth-order Hermite predictor-corrector integrators used in the Nbody6 code (Makino, 1991; Makino & Aarseth, 1992; Aarseth, 1999, 2003). Here, the outer (regular) integrator takes into account the slowly varying force due to the gas, while the inner (irregular) integrator takes into account the fast varying force due to the other sink particles. It is an anlogue to the Ahmad-Cohen scheme (Ahmad & Cohen, 1973) where the division to regular and irregular forces is based on the kind of interaction (gas or sink particles) instead of physical proximity. The regular time-step corresponds to the hydrodynamical time-step, while the irregular time-step is calculated according to the standard formula (Aarseth, 2003)
and then quantised to bins differing by factor of 2 in time. We set the constant for integration to be . The quantities a, , and are the acceleration and the higher time derivatives acting on the particle due to the other sink particles. The scheme uses the softening kernel described in Monaghan & Lattanzio (1985) with softening length corresponding to = 0.31 pc at the highest refinement level. Likewise, gas is attracted by sink particles, which are placed to the tree to facilitate the force evaluation. We present the detailed description of the sink particle integrator as well as numerical tests in Dinnbier et al. (in prep.).
2.2 Ionizing radiation and radiative heating
The transfer of ionizing radiation is calculated by a new module for the Flash code called TreeRay. It is an extension of the Flash tree solver described in Wünsch et al. (2018). TreeRay uses the octal-tree data structure constructed and updated at each time step by the tree solver and shares it with the Gravity (calculates gas self-gravity, see Wünsch et al. 2018), Optical-Depth (calculates the optical depth and parameters for the total, H and CO shielding, see Walch et al. 2015; Wünsch et al. 2018), and EUV modules. The latter is the new module which calculates the local flux of ionizing radiation. Here, we only give basic information about TreeRay; a detailed description alongside with accuracy and performance tests will be presented in Wünsch et al. (2018, in prep). TreeRay has already been benchmarked in Bisbas et al. (2015) and applied in homogeneous media (Haid et al., 2018).
Each node of the octal-tree represents a cuboidal collection of grid cells and stores the total gas mass contained in it, masses of individual chemical species and the position of the mass centre. In addition to that, TreeRay stores for each node the total amount of the radiation luminosity generated inside the node, the radiation energy flux passing through the node, and the node volume. Before the tree is traversed for each grid cell (called target cell), a system of rays pointing from the target cell to different directions is constructed. The directions are determined by the HEALPIX algorithm (Górski et al., 2005), which tessellates the unit sphere into elements of equal spatial angle. We use = 48. Each ray is then divided into segments with lengths increasing linearly with the distance from the target cell. In this way, the segment lengths correspond approximately to the sizes of the nodes interacting with the target cell during the tree walk if the Barnes-Hut (BH) criterion for node acceptance is used. Here, we use the BH criterion with an opening angle of . When the tree is traversed, node densities, radiation luminosities, and energy fluxes are mapped onto the ray according to the node and the volume belonging to the ray segment.
Finally, after the tree walk, the one-dimensional radiative transport equation is solved using the On-the-Spot approximation along each ray using the case B recombination coefficient with the temperature dependence in the range of = [5000, 20000] K given by (Draine, 2011)
The radiative transfer equation along the ray towards the target cell is given by
where is the received flux in the target cell, is the number of segments along a ray, is the emission coefficient in segment , the total flux coming into segment , the source in a segment if a source exists, the distance from the segment to the target cell, and the volume of the segment. As the radiation flux passing through a given segment from different directions has to be taken into account, the solution has to be searched for iteratively, repeating the whole process of tree construction, tree walk and solving the radiation transport equation until the maximum relative error drops below 0.01. To speed up convergence, we use the result of the previous hydrodynamic time-step as the radiation field typically changes only slightly between times-steps, in most cases only one or two iterations are needed in each time step.
We use the prescription given in Gatto et al. (2017) and Peters et al. (2017) to simulate the evolution, and in particular the radiative energy output, of massive stars using the Geneva stellar tracks from the zero-age main sequence to the Wolf-Rayet phase (Kudritzki & Puls, 2000; Markova et al., 2004; Markova & Puls, 2008; Puls et al., 2008; Ekström et al., 2012). An initial proto-stellar phase is not included. The corresponding time evolution of the radiative luminosity is shown in Fig. 1 for three stars with = 12 (red), 23 (black), and 60 M (blue). For the later discussion, we include dashed horizontal lines that correspond to the luminosities used in Dale et al. (2012), Dale et al. (2014) based on stellar models of Diaz-Miller et al. (1998).
The aforementioned stellar tracks provide the time-dependent number of Lyman continuum photons, , and the effective stellar temperature (Peters et al., 2017). In TreeRay, this information is processed to get the average excess photon energy, between = 13.6 eV and the average photon frequency , by assuming a black-body spectrum for each star and integrating it in the Lyman continuum (Rybicki & Lightman, 2004). Note that, since we only consider the radiative transfer in a single energy band (all photons in the Lyman continuum), we do not distinguish between the direct ionization of H and H, as necessary for detailed models of photon-dominated regions (Röllig et al., 2007; Baczynski et al., 2015).
We calculate the heating rate, , in the ionization-recombination equilibrium with (Tielens, 2005)
where is the photon flux, the hydrogen photo-ionization cross-section, the hydrogen number density, and is the Planck constant. The ionization heating rate and number of ionizing photons are provided to the Chemistry module (see Section 2.3), where the temperature is self-consistently increased by balancing heating and cooling processes, the mean hydrogen ionization state is updated using the given photo-ionization rate (Haid et al., 2018) and CO is dissociated. In ionization-recombination equilibrium, an Hii region develops around the sink particle with interior temperatures between 7000 – 9000 K. In homogeneous media, this is well explained by the ionization of the Strömgren sphere followed by the Spitzer expansion (Strömgren, 1939; Spitzer, 1978; Hosokawa & Inutsuka, 2006). However, the equilibrium temperature strongly depends on the density of the ionized gas within the Hii region and can be significantly lower in young, embedded Hii regions, which are still quite dense (see Section 5.2) .
2.3 Gas cooling, heating and chemistry
We include a simple chemical network, which is explained in detail in Walch et al. (2015). It is based on Glover & Mac Low (2007a, b); Glover et al. (2010) and Nelson & Langer (1997) to follow the abundances of seven chemical species: molecular, atomic and ionized hydrogen as well as carbon monoxide, ionized carbon, atomic oxygen and free electrons (H, H, H, CO, C, O, e). The gas has solar metallicity (Sembach et al., 2000) with fixed elemental abundances of carbon, oxygen and silicon ( = 1.14 10, = 3.16 10, = 1.5 10) and the dust-to-gas mass ratio is set to 0.01. We include a background interstellar radiation field (ISRF) of homogeneous strength G = 1.7 (Habing, 1968; Draine, 1978). So far, TreeRay does not treat the FUV regime. The effect of radiation in the FUV energy band will be discussed in a follow-up paper. For the cloud dynamics, we still expect photoionization to be the dominant process (Peters et al., 2010a; Walch et al., 2012; Baczynski et al., 2015) with typical temperatures around 8000 K almost independently of gas densities. FUV radiation is considered to be important in photo-dissociation regions which are forming ahead of the ionization shock fronts. As we show in Section 5.2, a hypothetical FUV field of 1000 G increases the gas temperature in such dense (10 g cm) photo-dissociation regions to a few 100 K at most. Therefore, the predicted dynamical effect resulting from the FUV heating is considered to be negligible in respect to the EUV heating.
The ISRF is attenuated in shielded regions depending on the column densities of total gas, H, and CO. Thus, we consider dust shielding and molecular (self-) shielding for H and CO (Glover et al., 2010) by calculating the shielding coefficients with the TreeRay Optical-Depth module (Wünsch et al., 2018). From the effective column density in each cell the visual extinction, , is calculated by
where the total gas column density is given by = where is the surface density, the mean molecular weight and the proton mass.
For gas with temperatures above 10 K we model the cooling rates according to Gnat & Ferland (2012) in collisional ionization equilibrium. Non-equilibrium cooling (also for Lyman ) is followed at lower temperatures through the chemical network. Within the Hii region, we neglect both C and O cooling because these species are predominantly in a higher ionization state. Heating rates include the photoelectric effect, cosmic ray ionization with a rate of = 310 s, X-ray ionization by Wolfire et al. (1995), and photo-ionization heating (see Section 2.2).
3 Simulation setup
3.1 The SILCC simulation
The SILCC simulation (Walch et al., 2015; Girichidis et al., 2016) is the basic setup and is used to self-consistently study the evolution of the SN-driven multi-phase ISM. The computational domain with an extent of 500 pc x 500 pc 5 kpc has a disc midplane with galactic properties at low red-shift similar to the solar neighbourhood. The boundary conditions for the gas are periodic in - and - direction and outflow in -direction. For gravity the boundary conditions are periodic in - and - direction and isolated in -direction (see Wünsch et al. 2018 for mixed gravity boundary conditions. The base grid resolution, denoted as l = 5 in the following, is = 3.9 pc.
Initially, the disc has a gas surface density of = 10 M pc and the density profile follows a Gaussian distribution in the vertical direction
where the scale height of the gas = 30 pc and the mid-plane density = 9 10 g cm. The initial temperature of the gas near the mid-plane is 4500 K and the disc is made up of H and C. The density is floored to 10 g cm with a temperature of 4 10 K in the gas at high altitudes above and below the galactic plane.
The simulation includes a static background potential to model the old and inactive stellar component in the disc, which is modelled as an isothermal sheet with a stellar surface density = 30 M pc and a scale height of 100 pc (Spitzer, 1942). This static potential is added to the gravitational potential of the self-gravitating gas that is calculated in every time-step.
For the first = 11.9 Myr from the start of the simulations, the development of a multi-phase ISM is driven by SN explosions. Therefore, we inject SNe at a fixed rate of 15 Myr. The SN rate bases on the Kennicutt-Schmidt relation (Schmidt, 1959; Kennicutt, 1998) and a standard initial mass function for . We use mixed SN driving, where SNe explode in density peaks and at random positions by an equal share (1:1 ratio). In -direction, the positioning of the random SNe is weighted with a Gaussian distribution with a scale height of 50 pc for (see Walch et al. 2015 and Girichidis et al. 2016 for more details).
A single SN injects 10 erg of energy. Whether the energy is injected in the form of internal energy or momentum depends on the ability to resolve the Sedov-Taylor radius. The spherical injection region has a minimum radius of four grid cells. In case the density is low, the Sedov-Taylor radius is resolved and the energy thermally injected. If it is unresolved, the temperature in the region is raised to = 10 K and the SN bubble is momentum-driven (Blondin et al., 1998; Gatto et al., 2015; Walch et al., 2015; Haid et al., 2016).
3.2 Initial conditions for this work: the zoom-in simulations
We refer to the re-simulation of selected clouds with a higher spatial resolution as zoom-in. Two different molecular clouds, MC and MC, are selected from the SILCC-ZOOM simulation (Seifried et al., 2017). The selected domains (see Table 1) are traced back in time to properly model the formation process from the beginning. The zoom-in starts at = 11.9 Myr where SN driving is suspended and the typical number densities in the selected zoom-in regions do not exceed some 10 cm.
During the zoom-in simulation, the resolution is gradually increased in both MCs. Starting from the SILCC base grid resolution of = 3.9 pc (l = 5), we allow adaptive refinement down to = 0.122 pc (l = 10). Two refinement criteria are used. The refinement on the second derivative of gas densities, which picks up density fluctuations, is limited to a maximum refinement of 0.5 pc. Further refinement to the smallest depends on the local Jeans length, , which is computed for each cell. We require that is resolved with at least 16 cells in each spatial dimension, otherwise we refine.
The zoom-in is not carried out in a single time-step. Starting from the base grid resolution, we increase the refinement step-by-step and require about 200 time steps in between two steps. On the one hand, this choice allows the relaxation of the gas to prevent filamentary grid artefacts, which appear in case of an instantaneous zoom-in (Seifried et al., 2017). On the other hand, it avoids the formation of large-scale, rotating, disc-like structures in case of a slower zoom-in in which case compressive motions are dissipated too efficiently. The zoom-in simulation reaches the highest refinement level at = 13.2 Myr.
3.3 Simulation overview
In this work, we continue from the zoom-in simulation at = 13.2 Myr and allow for the formation of cluster sink particles. From this time we start two simulations. The reference run, ZI_NOFB, does not include any stellar feedback. In the second simulation, ZI_RAD, the forming, active stellar component provides ionizing radiation.
Note that each simulation contains both molecular clouds, MC and MC with similar volumes , in which the zoom-in is enabled, and total gas masses within 10 percent. Table 1 summarizes the initial properties of the clouds. Each cloud develops its own sink evolution, i.e. star formation history. The first sink in MC forms at = 13.51 Myr and in MC at = 13.40 Myr. As radiation feedback sets in, the clouds start to evolve differently. Therefore, we define six times relative to , which we use for the analysis of the clouds as:
The subscript indicates the time in Myr after , i.e. = 0 Myr refers to and = 3.0 Myr after . A second subscript is used to indicate the respective MC.
Fig. 2 shows the total gas column density, , in the --plane in MC (top) and MC (bottom) at . To obtain , we integrate the density along the -direction within the volume (see Table 1). The black frames indicate the central subregions (later referred to as center of volume, CoV) to be shown in more detail in Fig. 3 and Fig. 4 with properties summarized in Table 1.
4 The morphology of the molecular clouds
The masses and volumes of MC and MC are comparable (see Table 1). However, their formation out of the turbulent, multi-phase ISM leads to different morphologies (see Fig. 2). MC contains a highly collimated, dense ( 500 M pc), T-shaped filament, where the bar is one horizontal structure with extended ends. The vertical trunk is divided into two, roughly parallel substructures. Each of the dense filaments is surrounded by an ’envelope’ with intermediate column densities ( 5 – 50 M pc). To the bottom, right a low column density ( 0.05 – 0.5 M pc) cavity is situated, which originates from a previous SN explosion outside of the cloud. In MC the main filamentary structure is vertically elongated and less condensed with a central, hub-like condensation. Qualitatively, the surface density maps of MC1 and MC2 span the same dynamic range.
In Fig. 3 and Fig. 4 we show the time evolution (from top to bottom) of the column densities of MC and MC without (leftmost column) and with radiative feedback (2 column). Sink particles without and with active stellar components are indicated with circles and stars, respectively with their age indicated by a second color scheme ranging from 0 to 3 Myr. The third and 4 column show the column densities of molecular hydrogen and ionized hydrogen for the runs with radiation. Here, we only show the central (40 pc x 40 pc) subregion of the clouds (see Table 1, bottom, subscript ’CoV’) and therefore the column densities are obtained from the integration along the -direction over the corresponding 40 pc, henceforth . By comparing the maximum values of in Fig. 2 with Fig. 3 and Fig. 4, one can see that only lower density gas from the fore- and background has been cut.
Without feedback (runs ZI_NOFB), gravity is further condensing the initial structures while the lower column density gas surrounding the main filaments is accreted. In MC the gas is gravitationally collapsing but the global structure of the cloud does not change significantly and is still recognizable at . In both clouds, sink formation occurs in the densest filament(s) and its debris.
Radiative feedback (runs ZI_RAD), does not significantly alter the global dynamics of MC during the first 2 Myr but more filamentary substructures appear while the existing substructures seem to be locally supported against gravitational collapse. The dense regions are puffed up by the expanding radiative shocks. Multiple radiation driven, partly or fully embedded bubbles develop (see in the right columns of Fig. 3 and Fig. 4). Some active sinks do not form a noticeable bubble of ionized hydrogen, in particular if the contained massive stars that are less massive than 20 M. During the last Myr, the clouds decompose and a variety of filamentary substructures evolve into all directions. The envelope is widened and heated gas is expelled into the cavity. Star formation takes place in the bar and its remnants but primarily in the central dense clump (compare to in the third column of Fig. 3). In MC, the bottom half of the cloud forms massive stars quickly, while the upper half forms only low-mass and hence inactive sink particles. Feedback from the bottom half disrupts the cloud into an upper, crescent-shaped filament (compare to in Fig. 4, third column) and some left-over, dispersed gas below (see at in Fig. 4, forth column) . At later time, the emerging feedback triggers a second generation (y = 195 pc, z = 5 pc) of massive stars in the upper part. The low-density envelope is replenished with expelled gas.
MC and MC with radiative feedback show a significant difference in morphology. The first cloud evolves into one massive structure with multiple embedded Hii regions and is surrounded by a low density envelope, which is only slowly evolving. The central structure hosts almost all stars and star formation continues. The second cloud is partly destroyed by a rapidly forming first generation of stars in the lower cloud filament. A new generation of stars is triggered in the upper part of the central subregion. That demonstrates that not only the mass (which is roughly similar for both clouds) but also the morphology prior to stellar feedback influences its impact. We investigate this further in Section 6.
5 Cloud environments with massive stars
In runs ZI_NOFB, the total numbers of sinks in MC and MC at time is 39 and 19 with masses of 1.810 and 1.510 M, respectively. In ZI_RAD, the two clouds host 31 sinks with 5900 M and 23 sinks with 3300 M, respectively (see top panel of Fig. 18). Hence in MC more sinks with a smaller average mass per sink particle are formed than in MC. The different fragmentation properties of the two clouds is caused by the different cloud substructure. In MC, the total number of sinks is slightly increased by radiative feedback, although the mass in sinks is dramatically reduced. This shows that radiative feedback may regulate star formation and, at the same time, trigger star formation. A more detailed investigation of triggered star formation is postponed to a follow-up paper.
The IMF for massive stars in MC (red) and MC (black) is shown in Fig. 5 for the simulations ZI_NOFB (thick, transparent lines) and ZI_RAD (thin, opaque lines) at time . The blue line indicates the Salpeter slope of the IMF in the high mass range proportional to (Salpeter, 1955). In the runs ZI_RAD, 31 and 23 massive stars form with 830 and 480 M in the total volume of MC and MC, respectively. Within the period of 3 Myr, a small number of massive stars is forming, which leaves the high-mass end of the IMF under-sampled. Hypothetical sampling of low-mass stars ( 9 M) from the residual sink mass results in a well-represented low-mass end of the Kroupa IMF. The corresponding formation history of the massive stars is shown in Fig. 20.
Sink particles accrete gas from their environment as long as the gas is e.g. gravitationally bound. When accretion stops the sink particle has reached its maximum mass, . It is useful to investigate the accretion time, , which is the time elapsed from the formation of the sink until the maximum sink mass is reached, which quantitatively demonstrates the impact of radiative feedback on the local star formation rate. In Fig. 6, we show as a function of for the simulations of MC (red) and MC (black) with (ZI_RAD; full markers) and without radiative feedback (ZI_NOFB; open markers) within . Transparent markers are sinks with masses below the massive star formation limit of 120 M. Blue crosses indicate that the accretion onto the corresponding sink has stopped. In ZI_NOFB, the accretion times stretch over a wider temporal range and sinks grow to a few 1000 M because accretion cannot be halted. In ZI_RAD, the accretion time is less than 1 Myr. Sink particles with masses above the star formation threshold not only stop their own accretion but effect or even interrupt the mass accretion of any nearby companion. This results in a large fraction of sinks that remain below the star formation mass threshold. We expect that short accretion times are accompanied by a drastic change in the environmental density of the sink particles as a function of time. This is investigated in the following Section.
5.1 Environmental densities
Over the life time of massive stars, their environmental densities are continuously changing. The ambient density determines the impact of radiative feedback (Haid et al., 2018). Right after the star a star is born the surrounding gas is typically so dense, that the young Hii region is confined (Wood & Churchwell, 1989; Peters et al., 2010b). The bubble then expands hydro-dynamically, while more gas is ionized. However, once the environmental density has significantly decreased (after the ionized gas has leaked out of the bubble or the star has moved out of the dense star-forming filament) the ionization front spreads out and the impact of radiative feedback is not locally confined.
In Fig. 7, we show the cumulative distribution of the number of sinks (top, ) and the cumulative sink mass (bottom, ) as a function of the environmental density obtained by averaging the ambient density of each sink in a sphere with a radius of 1 pc, . We show MC (left) and MC (right) with (ZI_RAD; thin lines) and without feedback (ZI_NOFB; thick lines) at two times: the formation time of a sink particle is denoted with (red), and the end time (black). The green vertical line indicates the sink formation threshold density, .
First of all, when comparing ZI_RAD and ZI_NOFB in MC, we can see that with radiative feedback more sinks are formed, which hints towards triggered star formation. However, the higher number of sink contains a significantly lower total mass (see bottom panels), which indicates that feedback limits the accretion onto star forming dense regions (see also Fig. 6). Furthermore, we can see that the sinks’ environmental densities are severely changed by radiative feedback. All sinks are born in very dense gas ( g cm) and without feedback most of them also stay there (modulo some wandering off a bit). This implies that there is a large enough gas reservoir to feed the sink particles for the simulated time, even if their mass has grown significantly. With radiative feedback, however, the distribution is significantly shifted towards lower densities at compared to . Even though there is still a number of sinks at g cm, these are mostly the young sink particles which did not have time to disperse their environment. In MC about 60 percent of all sinks are surrounded by gas with g cm. The dispersal of MC has progressed farther and 80 percent of all sinks are found at g cm.
Fig. 7 shows that the environmental densities for more than 90 and 60 percent of the sink particles in MC, respectively MC lie between 10 g cm and 10 g cm at . These are conditions, where stellar winds were shown not to be important (Geen et al., 2015; Haid et al., 2018). Therefore, we do not include this additional feedback process in this work. Nevertheless, the environmental densities are continuously reduced by ionizing radiation at later stages () stellar wind might become important.
5.2 The multi-phase evolution
We show the mass-weighted (color) density-temperature and density-pressure (pressure over the Boltzmann constant) distributions for both clouds in Fig. 8 and Fig. 9, respectively. Note that, according to the Jeans criterion, the depicted gas density is fully resolved, even at the high density end. For each cloud (MC, top panels; MC, bottom panels) we show different times , , and from top to bottom. To guide the eye, the black lines show the thermal equilibrium curves calculated using a stand-alone version of the chemistry module with increasing G of 1.7 (solid), 17 (dashed), 170 (dash-dotted), and 1700 (dotted) in units of Habing fields. The thermal equilibrium curve can be assumed to be the transition to the CO dominated gas (Röllig et al., 2007). In each row, the left column shows the runs without radiative feedback (ZI_NOFB) and the three panels to the right base on run ZI_RAD. The second panel shows the total gas within , while the third and forth panel show gas above and below a visual extinction, 1 mag, which is computed self-consistently for every cell in the computational domain using the TreeRay OpticalDepth module (see Section 2.3). The markers indicate the average environmental density, temperature, and pressure of sinks without (circles) and with active stellar components (stars) within a sphere of radius 1 pc around each sink particle, , , and . The numbers in the lower left corners indicate the mass within at given time or the fraction of mass at high and low relative to , respectively. It can be seen that MC has a significantly higher fraction of shielded gas than MC (for further analysis see Section 6.1). The sinks are shown in the respective high/low panels depending on their average within the surrounding 1 pc radius.
The gas distributions of runs ZI_NOFB (Fig. 8 and 9, left column) are more or less constant in time and follow the computed equilibrium curves. For g cm, most of the gas is more deeply embedded and cools down to 10 K. For runs with radiative feedback, the phase diagrams change significantly as a lot of gas is lifted above the equilibrium curve towards high temperatures and pressures. Several new horizontal branches become apparent in the temperature-density diagram. Young and deeply embedded Hii regions first appear in the 1 mag distribution (see Fig. 17 for the ionization state of this gas). With time, these embedded bubbles grow, burst out of the dense filament and leak into the low density environment. Fully developed Hii regions are heated up to temperatures between 8000 K and 10000 K occurring at 1 mag.
Inactive stars are usually deeply embedded inside the cloud (at high density and low temperature), unless they reside near an active sink which influences their environment. The Hii regions seem to expand until the pressure gradient between ambient medium and Hii region across their outer boundary is diminished, which can be seen from Fig. 9 where the gas inside the Hii region joins the rising equilibrium pressure branch of the warm ambient medium.
6 The interaction between molecular clouds and radiative feedback
MC and MC were chosen as two clouds with similar initial parameters (see Seifried et al., 2017, and Table 1). Nevertheless, the clouds evolve differently in the presence of radiative feedback (see Fig. 3 and Fig. 4) where MC seems less affected than MC. In this Section, we discuss the physical property of the cloud, i.e. the local extinction, which we ultimately (after a careful and extensive analysis) identify to be responsible for the apparent differences. Next, we discuss the energy content and the star formation properties of both clouds.
6.1 Extinction matters!
Fig. 10 shows the fraction of cloud mass constrained by different extinction thresholds in MC (red) and MC (black) in the total cloud ( /, top) and the central subregion ( /, bottom) as a function of time in simulations ZI_RAD (thin) and ZI_NOFB (thick). Note that the evolutions of and are shown in Fig. 16. We evaluate the mass with the extinction below (dashed) and above (solid) a visual extinction of = 1 mag. For both simulations, ZI_NOFB and ZI_RAD, the evolutions are similar for the first 1.5 Myr. In the total domain of MC and MC, most mass resides at 1 mag with 0.7 of the total mass, M. Hence, only a small fraction of the gas is well-shielded with 1 mag and the well-shielded mass fraction is higher in MC than in MC by 30 percent at up to 80 percent at and finally become similar at . In the central subregions, the evolutions of the well-shielded gas follow the larger volume with higher initial fractions of 60 percent and to 50 percent of M. During the evolution, gas is dispersed by feedback and the fraction of A 1 mag dominates (compare to Fig. 8 and Fig. 9).
In Fig. 11, we show the mass-weighted (top) and volume-weighted (bottom) density (left) and column density PDF (right) in simulation ZI_RAD at (top panels) and (bottom panel). The PDF includes gas within for MC (red) and MC (black). Dotted lines indicate the density PDF of gas with mag. This is calculated for every cell in the computational domain via the TreeRay OpticalDepth module. We find that the mass-weighted and volume-weighted total gas density distributions of MC and MC are similar over a wide range of densities and at and . The column density PDFs (right column) of the clouds are slightly different: While MC dominates in the high regime (), MC hosts more gas at lower column densities.
The differences between the two clouds become apparent when inspecting the high- gas (blue, green). MC has more mass in gas at high than MC and this gas occupies a larger fraction of the cloud volume. Also, there is basically no difference between the two time steps and for MC, apart from the very high density tail in the mass-weighted density PDF, which forms at . This indicates that the gas in MC is still relatively confined at . On the other hand, MC has less mass at high and this mass occupies a smaller fraction of the cloud volume. Also, the mass and volume fractions of gas at high are clearly decreasing as a function of time, which is a clear sign of cloud dispersal.
Overall we find that the impact of radiative feedback is very sensitive to the detailed cloud substructure. In this regard, even if the volume density distributions are similar, the distribution of the three-dimensional, visual extinction may be different and these differences are enhanced when the cloud is exposed to radiative feedback.
6.2 Energy evolution
In Fig. 12, we show the evolution of the internal (red, ) and kinetic (black, ) energy of the gas for the simulations ZI_RAD (thin) and ZI_NOFB (thick) for MC (top) and MC (bottom). Note that we only consider the gas while the contribution to from sink particles is neglected. The initial kinetic and internal energies are similar for both MCs with 310 erg and 210 erg. Both clouds are initially bound with virial parameters of 0.72 and 0.89 as calculated for in . In ZI_NOFB, the energies in both clouds decrease as no feedback energy is injected. With radiative feedback the internal energy increases following the formation of massive stars. Radiative feedback also clearly enhances the kinetic energy content of the clouds, i.e. it drives turbulence (see e.g. Gritschneder et al., 2009; Walch et al., 2012).
In both clouds, we see jumps in by up to 50 percent. This behaviour is linked to the confinement of radiative bubbles. Initially, they are embedded in dense structures, only a small volume is affected and the radiative impact is delayed. But as soon as the Hii regions open, radiation and ionized material leak out and increase the internal energy in a larger domain. As the average rate of ionizing photons in the central volume is comparable for both clouds with a factor of 2 difference (see Fig. 20), the final energtic states are similar.
6.3 Star formation
A common way to separate the diffuse ISM from the dense gas is to choose a density threshold. With the subscript ’100’ we refer to gas with number density 100 cm ( 3.8410 g cm). To trace predominantly molecular gas, we use the subscript ’H2’, which means that the mass fraction of H in every cell is equal or greater than 50 percent (see Seifried et al., 2017). A general way to indicate either of the two constraints is a subscript ’x’.
The instantaneous star formation rate surface density assumes that all gas that is accreted onto a sink particle is immediately forming an ensemble of low and high mass stars. It is defined as
where = 0.1 Myr, is the area of the cloud and the mass accretion rate of the sinks over a time period (Matzner & McKee, 2000; Gatto et al., 2017). is calculated from the mass-weighted radius, which we calculate from the distance of all cells above a given threshold relative to the center of mass in the volume .
In Fig. 13, we depict for MC (top) and MC (bottom) for the simulations ZI_RAD (thin) and ZI_NOFB (thick). We compute the mass-weighted cloud area from all cells with 100 cm. The time averaged values for the simulation ZI_NOFB are 1.9 and 1.3 M Myr pc as well as 0.6 and 0.3 M Myr pc in run ZI_RAD for MC and MC, respectively. This shows that radiative feedback reduces the star formation rate surface density by a factor of 4. Similar values are obtained for a radius constrained with molecular hydrogen dominated gas.
The star formation efficiency per free-fall time (SFE) is the dimensionless ratio of the mass in stars that forms within a free-fall time divided by the mass of the cloud with x = [100, H2]. In this paper, the mass in stars is equivalent to the mass in the sink particles, = . The free-fall time is given with , where . This gives
In Fig. 14, we show for simulation ZI_RAD in MC (top) and MC (bottom) for (red) and (black). The time-averaged values are shown as horizontal, solid lines. The horizontal, dashed lines are the time-averaged values from the corresponding simulations without radiative feedback, ZI_NOFB. The time-averaged free-fall times are 2.7, 2.4 Myr in MC and 2.8, 1.8 Myr in MC for the thresholds x = [100, H2], respectively. For the number density threshold and the H-based threshold , the time averaged SFE are 9 and 13 percent in MC and 5 and 6 percent in MC with radiative feedback, and 29 and 40 percent in MC, respectively 31 and 37 percent in MC without feedback. Thus, in the simulations ZI_NOFB the average values are 4 times higher and in agreement with the findings for . In general, we obtain somewhat higher average values due to short episodes of high star formation, although the SFE regularly drops down to the 5-percent regime (green, shaded area). Note that SFE is larger due to a smaller mass that is available for star formation (see Eq. 10).
Values for SFE are observed for low-mass clouds to be around a few percent (Krumholz & Tan, 2007; Evans et al., 2009). In clouds which are more massive and/or have longer free-fall times, the efficiency can increase up to 30 percent (Murray, 2011). Isolated, bound MCs in numerical simulations by Dale et al. (2012, 2014) show SFE without and with feedback of 16 and 11 percent, which indicates that radiative feedback is inefficient in regulating star formation. This contradicts our findings, where ionizing radiation reduces the SFE on average by a factor of 4. Similar values are found by Howard et al. (2016). One reason for the inefficiency of radiative feedback in Dale et al. (2012) and Dale et al. (2014) can be found in the underlying model of Diaz-Miller et al. (1998), which systematically underestimates the ionizing luminosities. These differ from the presented model by up to a factor of 10 lower values with increasing stellar mass (see Fig. 1).
Fig. 15 compares the SFE obtained from the simulations with resolved observations of Milky Way molecular clouds, kpc-scale observations of Local Group galaxies, and from unresolved observations of both disk and starburst galaxies in the local universe and at high redshift published in Lada et al. (2010) and Heiderman et al. (2010). We relate the star formation rate surface density with the surface density over the free-fall time derived for MC (red) and MC (black) for the simulations ZI_RAD (full markers) and ZI_NOFB (empty markers) at . For each cloud, we only consider gas within its V above the number density threshold 100 cm. The black line and grey shaded area shows the fitted behaviour found in Krumholz et al. (2012, 2013) and the associated uncertainty.
The SFE obtained for the ZI_NOFB runs are too high and clearly offset from the observed relation. However, both clouds with radiative feedback are right on top of the relation, with MC at slightly lower than MC. The two observed points which sit directly on top of our results correspond to the Taurus and Ophiuchus MCs (points near MC) and Lupus 3 (point on top of MC result). This result is reassuring because these observed clouds lie in the solar neighbourhood - the environment simulated here - and have total masses and other physical properties that are comparable to our simulated clouds.
7 Discussion: Differences between MC and MC
MC and MC condense out of the same multi-phase ISM. They were selected to have similar initial properties like masses around 10 M, similar volumes, comparable kinetic and internal energies, and similar virial parameters of 0.72 and 0.89. However, during the evolution under the influence of ionizing radiation both clouds seem to diverge in respect to their morphologies, meaning that most of the gas in MC remains (see Fig. 3) while MC is almost fully dispersed at (see Fig. 4).
In Fig. 7, we show that the environmental densities of sink particles gradually decrease with age and particularly the sinks in MC are embedded in low-density media. Together with Fig. 11, where we show that this cloud has much less well-shielded gas with 1 mag, we interpret that MC has more deeply embedded dense structures and a thicker envelope. The density-temperature distribution of the central region of the clouds (see Fig. 8) with 1 mag indicates that some sources are deeply embedded in these well-shielded regions. Therefore, radiative feedback is confined to small bubbles in MC. The radiative impact is delayed until the radiative bubbles open into the ambient medium, ionized gas and radiation leak out and induce kinetic motions.
It is important to mention, that the emitted radiative energy is similar in both central regions (see Fig. 20 for the rate of UV-photons). There are also massive stars forming in the rest of the cloud. The most relevant stars have high masses. In our simulations, almost all of those are situated far away from the central subregion. With distance and with decreasing mass the UV-photon rates per volume decrease and easily drop below the rates form a 9 M, hence are considered as minor. The effect of feedback from massive stars outside the subregions on the dense structures is minor. Otherwise it should also be visible in the mass evolution of well shielded gas (see Fig. 10 and compare ZI_NOFB with ZI_RAD) and in the density-temperature distribution, but both remain almost unchanged (see Fig 8 and Fig. 9).
Concerning the star formation in both clouds, ionizing radiation is able to lower the star formation rate surface density by a factor of 4. The star formation efficiency constrained by gas above 100 cm is found to be on average 5 – 9 percent in both clouds (see Fig. 14 and Fig. 13). In ZI_NOFB, a few, massive sink particles evolve, whereas in ZI_RAD in cloud MC the sink masses are significantly reduced but their number increased. Star formation is triggered by radiative feedback.
The comparison of MC and MC shows that, not only the cloud masses (Dale et al., 2012, 2013), the corresponding luminosities (Geen et al., 2018) and escape velocities influence the impact from radiative feedback, but that the initial cloud sub-structure significantly determines the cloud evolution. The initial conditions are imprinted during the formation process of the cloud (Brunt et al., 2009; Rey-Raposo et al., 2017). We find that the fully three-dimensional shielding properties determine the time scales of molecule formation (Seifried et al., 2017) as well as the time scales for cloud dispersal (this paper).
In this paper, we investigate the impact of ionizing radiation feedback from massive stars in the early evolution of MCs up to 3 Myr. We perform hydrodynamic simulations with the AMR code FLASH 4 and include the novel radiative transfer scheme TreeRay, which is coupled to a chemical network to treat the effect of ionizing radiation. We self-consistently follow the formation of two, initially bound MCs from a SN-driven, multi-phase ISM down to a resolution of 0.122 pc within the SILCC-ZOOM project. We allow for sink particle formation on the highest refinement level. In the simulation ZI_RAD, ionizing radiation is coupled to massive stars. Simulation ZI_NOFB is the reference run without feedback process. In the following we list the main conclusions.
Despite the similar initial masses of the two MCs, the morphological evolutions under the influence of ionizing radiations is different. In MC, a central blob of gas remains, whereas a part of MC is fully dispersed. We show that this difference is linked to the mass of internal substructures of dense and well shielded gas, which embeds, and delays radiative feedback. The substructures are imprinted during the formation of the clouds.
We show that the total gas density PDFs are nearly identical for the different MCs. However, the well-shielded gas ( 1 mag) reveals cloud-specific properties. In this work, MC shows more volume-filling gas at intermediate densities, i.e. a thicker envelope surrounding the densest, star-forming filaments. These are responsible for sustaining the cloud structure despite the impact of radiative feedback. MC contains less of this gas, hence becomes more easily dispersed.
In the simulation, we find for some massive stars that the environmental densities are high and the radiative bubbles are embedded in substructures. When the radiative bubble opens into the ambient medium, internal energy and hot, ionized gas is released and the embedded phase is terminated. This behaviour is reflected by small jumps in the internal energy evolution. In this phase the ionized gas inside the Hii region is heated to the prototypical 8000 K.
Star formation can be regulated by radiative feedback. In simulations ZI_RAD the star formation efficiency drops by a factor of 4 compared to the ZI_NOFB in both clouds. The star formation efficiency in gas with densities above 100 cm lies at 5 – 9 percent. This indicates that internal morphologies regulate the impact of photo-ionizing radiation, hence the star formation.
Without feedback, a few sink particles accrete a significant fraction of the cloud mass. Radiative feedback significantly reduces the sink mass and instead may increase its number. This seems to be triggered star formation, even though the overall star formation efficiency is so severely reduced.
When comparing with observational data, we find that our two clouds, which were simulated using typical solar neighbourhood conditions, lie on top of the results derived for Taurus, Ophiuchus, and Lupus 3, three low-mass, star forming, nearby molecular clouds with similar total masses.
SH, SW, DS and FD acknowledge the support by the Bonn-Cologne Graduate School for physics and astronomy which is funded through the German Excellence Initiative. SH, SW, and DS also acknowledge funding by the Deutsche Forschungsgemeinschaft (DFG) via the Sonderforschungsbereich SFB 956 ”Conditions and Impact of Star Formation” (subproject C5). SH, SW, DS, FD and TN acknowledge the support by the DFG Priority Program 1573 ”The physics of the interstellar medium”. SH and SW acknowledge funding by the European Research Council through ERC Starting Grant No. 679852 ”RADFEEDBACK”. TN acknowledges support from the DFG cluster of excellence ”Origin and Structure of the Universe”. R.W. acknowledges support by the Albert Einstein Centre for Gravitation and Astrophysics via the Czech Science Foundation grant 14-37086G and by the institutional project RVO:67985815 of the Academy of Sciences of the Czech Republic. The software used in this work was developed in part by the DOE NNSA ASC- and DOE Office Science ASCR-supported FLASH Center for Computational Science at University of Chicago. The authors gratefully acknowledge the Gauss Centre for Supercomputing e.V. (www.gauss-centre.eu) for funding this project (pr62su) by providing computing time on the GCS Supercomputer SuperMUC at Leibniz Supercomputing Centre (www.lrz.de). We thank the YT-PROJECT community (Turk et al., 2011) for the YT analysis package, which we used to analyse and plot most of the data. We thank the anonymous referee for the constructive input.
- Aarseth (1999) Aarseth, S. J. 1999, PASP, 111, 1333
- Aarseth (2003) Aarseth, S. J. 2003, Gravitational N-Body Simulations, 430
- Ahmad & Cohen (1973) Ahmad, A. & Cohen, L. 1973, Journal of Computational Physics, 12, 389
- André et al. (2014) André, P., Di Francesco, J., Ward-Thompson, D., et al. 2014, Protostars and Planets VI, 27
- Baczynski et al. (2015) Baczynski, C., Glover, S. C. O., & Klessen, R. S. 2015, MNRAS, 454, 380
- Balfour et al. (2015) Balfour, S. K., Whitworth, A. P., Hubber, D. A., & Jaffa, S. E. 2015, MNRAS, 453, 2471
- Bisbas et al. (2015) Bisbas, T. G., Haworth, T. J., Williams, R. J. R., et al. 2015, MNRAS, 453, 1324
- Blondin et al. (1998) Blondin, J. M., Wright, E. B., Borkowski, K. J., & Reynolds, S. P. 1998, ApJ, 500, 342
- Bonnell et al. (2013) Bonnell, I. A., Dobbs, C. L., & Smith, R. J. 2013, MNRAS, 430, 1790
- Bouchut et al. (2007) Bouchut, F., Klingenberg, C., & Waagan, K. 2007, Numerische Mathematik, 108, 7, 10.1007/s00211-007-0108-8
- Bouchut et al. (2010) Bouchut, F., Klingenberg, C., & Waagan, K. 2010, Numerische Mathematik, 115, 647, 10.1007/s00211-010-0289-4
- Brunt et al. (2009) Brunt, C. M., Heyer, M. H., & Mac Low, M.-M. 2009, A&A, 504, 883
- Butler et al. (2017) Butler, M. J., Tan, J. C., Teyssier, R., et al. 2017, ApJ, 841, 82
- Castor et al. (1975) Castor, J., McCray, R., & Weaver, R. 1975, ApJ, 200, L107
- Clark et al. (2012) Clark, P. C., Glover, S. C. O., Klessen, R. S., & Bonnell, I. A. 2012, MNRAS, 424, 2599
- Colombo et al. (2014) Colombo, D., Hughes, A., Schinnerer, E., et al. 2014, ApJ, 784, 3
- Dale (2015) Dale, J. E. 2015, New A Rev., 68, 1
- Dale et al. (2005) Dale, J. E., Bonnell, I. A., Clarke, C. J., & Bate, M. R. 2005, MNRAS, 358, 291
- Dale et al. (2012) Dale, J. E., Ercolano, B., & Bonnell, I. A. 2012, MNRAS, 424, 377
- Dale et al. (2013) Dale, J. E., Ngoumou, J., Ercolano, B., & Bonnell, I. A. 2013, MNRAS, 436, 3430
- Dale et al. (2014) Dale, J. E., Ngoumou, J., Ercolano, B., & Bonnell, I. A. 2014, MNRAS, 442, 694
- de Avillez & Breitschwerdt (2005) de Avillez, M. A. & Breitschwerdt, D. 2005, A&A, 436, 585
- Diaz-Miller et al. (1998) Diaz-Miller, R. I., Franco, J., & Shore, S. N. 1998, ApJ, 501, 192
- Dobbs (2015) Dobbs, C. L. 2015, MNRAS, 447, 3390
- Dobbs et al. (2014) Dobbs, C. L., Krumholz, M. R., Ballesteros-Paredes, J., et al. 2014, Protostars and Planets VI, 3
- Dobbs et al. (2012) Dobbs, C. L., Pringle, J. E., & Burkert, A. 2012, MNRAS, 425, 2157
- Dobbs et al. (2015) Dobbs, C. L., Pringle, J. E., & Duarte-Cabral, A. 2015, MNRAS, 446, 3608
- Draine (1978) Draine, B. T. 1978, ApJS, 36, 595
- Draine (2011) Draine, B. T. 2011, Physics of the Interstellar and Intergalactic Medium (Princeton Series in Astrophysics)
- Dubey et al. (2008) Dubey, A., Fisher, R., Graziani, C., et al. 2008, in Astronomical Society of the Pacific Conference Series, Vol. 385, Numerical Modeling of Space Plasma Flows, ed. N. V. Pogorelov, E. Audit, & G. P. Zank, 145
- Ekström et al. (2012) Ekström, S., Georgy, C., Eggenberger, P., et al. 2012, A&A, 537, A146
- Evans et al. (2009) Evans, II, N. J., Dunham, M. M., Jørgensen, J. K., et al. 2009, ApJS, 181, 321
- Fall et al. (2010) Fall, S. M., Krumholz, M. R., & Matzner, C. D. 2010, ApJ, 710, L142
- Federrath et al. (2010) Federrath, C., Banerjee, R., Clark, P. C., & Klessen, R. S. 2010, ApJ, 713, 269
- Foster & Chevalier (1993) Foster, P. N. & Chevalier, R. A. 1993, ApJ, 416, 303
- Fryxell et al. (2000) Fryxell, B., Olson, K., Ricker, P., et al. 2000, ApJS, 131, 273
- Gatto et al. (2015) Gatto, A., Walch, S., Low, M.-M. M., et al. 2015, MNRAS, 449, 1057
- Gatto et al. (2017) Gatto, A., Walch, S., Naab, T., et al. 2017, MNRAS, 466, 1903
- Gavagnin et al. (2017) Gavagnin, E., Bleuler, A., Rosdahl, J., & Teyssier, R. 2017, MNRAS, 472, 4155
- Geen et al. (2015) Geen, S., Rosdahl, J., Blaizot, J., Devriendt, J., & Slyz, A. 2015, MNRAS, 448, 3248
- Geen et al. (2018) Geen, S., Watson, S. K., Rosdahl, J., et al. 2018, MNRAS, 481, 2548
- Girichidis et al. (2016) Girichidis, P., Walch, S., Naab, T., et al. 2016, MNRAS, 456, 3432
- Glover et al. (2010) Glover, S. C. O., Federrath, C., Mac Low, M.-M., & Klessen, R. S. 2010, MNRAS, 404, 2
- Glover & Mac Low (2007a) Glover, S. C. O. & Mac Low, M.-M. 2007a, ApJS, 169, 239
- Glover & Mac Low (2007b) Glover, S. C. O. & Mac Low, M.-M. 2007b, ApJ, 659, 1317
- Gnat & Ferland (2012) Gnat, O. & Ferland, G. J. 2012, ApJS, 199, 20
- Goldreich & Kwan (1974) Goldreich, P. & Kwan, J. 1974, ApJ, 189, 441
- Górski et al. (2005) Górski, K. M., Hivon, E., Banday, A. J., et al. 2005, ApJ, 622, 759
- Gritschneder et al. (2009) Gritschneder, M., Naab, T., Walch, S., Burkert, A., & Heitsch, F. 2009, ApJ, 694, L26
- Habing (1968) Habing, H. J. 1968, Bull. Astron. Inst. Netherlands, 19, 421
- Haid et al. (2016) Haid, S., Walch, S., Naab, T., et al. 2016, MNRAS, 460, 2962
- Haid et al. (2018) Haid, S., Walch, S., Seifried, D., et al. 2018, MNRAS, 478, 4799
- Harper-Clark & Murray (2009) Harper-Clark, E. & Murray, N. 2009, ApJ, 693, 1696
- Heiderman et al. (2010) Heiderman, A., Evans, II, N. J., Allen, L. E., Huard, T., & Heyer, M. 2010, ApJ, 723, 1019
- Heitsch et al. (2005) Heitsch, F., Burkert, A., Hartmann, L. W., Slyz, A. D., & Devriendt, J. E. G. 2005, ApJ, 633, L113
- Hennebelle & Iffrig (2014) Hennebelle, P. & Iffrig, O. 2014, A&A, 570, A81
- Hetem & Lepine (1993) Hetem, Jr., A. & Lepine, J. R. D. 1993, A&A, 270, 451
- Hill et al. (2012) Hill, A. S., Joung, M. R., Mac Low, M.-M., et al. 2012, ApJ, 750, 104
- Hosokawa & Inutsuka (2006) Hosokawa, T. & Inutsuka, S.-i. 2006, ApJ, 646, 240
- Howard et al. (2017) Howard, C., Pudritz, R., & Klessen, R. 2017, ApJ, 834, 40
- Howard et al. (2016) Howard, C. S., Pudritz, R. E., & Harris, W. E. 2016, MNRAS, 461, 2953
- Hughes et al. (2013) Hughes, A., Meidt, S. E., Colombo, D., et al. 2013, ApJ, 779, 46
- Ibáñez-Mejía et al. (2017) Ibáñez-Mejía, J. C., Mac Low, M.-M., Klessen, R. S., & Baczynski, C. 2017, ApJ, 850, 62
- Inoue & Fukui (2013) Inoue, T. & Fukui, Y. 2013, ApJ, 774, L31
- Joung & Mac Low (2006) Joung, M. K. R. & Mac Low, M.-M. 2006, ApJ, 653, 1266
- Kennicutt (1998) Kennicutt, Jr., R. C. 1998, ApJ, 498, 541
- Kim & Ostriker (2018) Kim, C.-G. & Ostriker, E. C. 2018, ApJ, 853, 173
- Kim et al. (2013) Kim, C.-G., Ostriker, E. C., & Kim, W.-T. 2013, ApJ, 776, 1
- Klessen (2011) Klessen, R. S. 2011, in EAS Publications Series, Vol. 51, EAS Publications Series, ed. C. Charbonnel & T. Montmerle, 133–167
- Klessen & Glover (2016) Klessen, R. S. & Glover, S. C. O. 2016, Star Formation in Galaxy Evolution: Connecting Numerical Models to Reality, Saas-Fee Advanced Course, Volume 43. ISBN 978-3-662-47889-9. Springer-Verlag Berlin Heidelberg, 2016, p. 85, 43, 85
- Klessen et al. (2000) Klessen, R. S., Heitsch, F., & Mac Low, M.-M. 2000, ApJ, 535, 887
- Körtgen et al. (2016) Körtgen, B., Seifried, D., Banerjee, R., Vázquez-Semadeni, E., & Zamora-Avilés, M. 2016, MNRAS, 459, 3460
- Kroupa (2001) Kroupa, P. 2001, MNRAS, 322, 231
- Krumholz (2006) Krumholz, M. R. 2006, ApJ, 641, L45
- Krumholz et al. (2012) Krumholz, M. R., Dekel, A., & McKee, C. F. 2012, ApJ, 745, 69
- Krumholz et al. (2013) Krumholz, M. R., Dekel, A., & McKee, C. F. 2013, ApJ, 779, 89
- Krumholz et al. (2009) Krumholz, M. R., Klein, R. I., McKee, C. F., Offner, S. S. R., & Cunningham, A. J. 2009, Science, 323, 754
- Krumholz & Matzner (2009) Krumholz, M. R. & Matzner, C. D. 2009, ApJ, 703, 1352
- Krumholz & Tan (2007) Krumholz, M. R. & Tan, J. C. 2007, ApJ, 654, 304
- Kudritzki & Puls (2000) Kudritzki, R.-P. & Puls, J. 2000, ARA&A, 38, 613
- Kuffmeier et al. (2017) Kuffmeier, M., Haugbølle, T., & Nordlund, Å. 2017, ApJ, 846, 7
- Lada & Lada (2003) Lada, C. J. & Lada, E. A. 2003, ARA&A, 41, 57
- Lada et al. (2010) Lada, C. J., Lombardi, M., & Alves, J. F. 2010, ApJ, 724, 687
- Mac Low et al. (2004) Mac Low, M.-M., de Avillez, M. A., & Korpi, M. J. 2004, in Astrophysics and Space Science Library, Vol. 315, How Does the Galaxy Work?, ed. E. J. Alfaro, E. Pérez, & J. Franco, 339
- Mac Low & Klessen (2004) Mac Low, M.-M. & Klessen, R. S. 2004, Reviews of Modern Physics, 76, 125
- Makino (1991) Makino, J. 1991, ApJ, 369, 200
- Makino & Aarseth (1992) Makino, J. & Aarseth, S. J. 1992, PASJ, 44, 141
- Markova & Puls (2008) Markova, N. & Puls, J. 2008, A&A, 478, 823
- Markova et al. (2004) Markova, N., Puls, J., Repolust, T., & Markov, H. 2004, A&A, 413, 693
- Matzner (2002) Matzner, C. D. 2002, ApJ, 566, 302
- Matzner & McKee (2000) Matzner, C. D. & McKee, C. F. 2000, ApJ, 545, 364
- Mellema et al. (2006) Mellema, G., Arthur, S. J., Henney, W. J., Iliev, I. T., & Shapiro, P. R. 2006, ApJ, 647, 397
- Monaghan & Lattanzio (1985) Monaghan, J. J. & Lattanzio, J. C. 1985, A&A, 149, 135
- Murray (2011) Murray, N. 2011, ApJ, 729, 133
- Murray et al. (2010) Murray, N., Quataert, E., & Thompson, T. A. 2010, ApJ, 709, 191
- Naab & Ostriker (2017) Naab, T. & Ostriker, J. P. 2017, ARA&A, 55, 59
- Nelson & Langer (1997) Nelson, R. P. & Langer, W. D. 1997, ApJ, 482, 796
- Ngoumou et al. (2015) Ngoumou, J., Hubber, D., Dale, J. E., & Burkert, A. 2015, ApJ, 798, 32
- Nordlund et al. (2017) Nordlund, Å., Ramsey, J. P., Popovas, A., & Kuffmeier, M. 2017, ArXiv e-prints
- Ostriker & McKee (1988) Ostriker, J. P. & McKee, C. F. 1988, Reviews of Modern Physics, 60, 1
- Padoan et al. (2017) Padoan, P., Haugbølle, T., Nordlund, Å., & Frimann, S. 2017, ApJ, 840, 48
- Padoan et al. (2016) Padoan, P., Pan, L., Haugbølle, T., & Nordlund, Å. 2016, ApJ, 822, 11
- Peters et al. (2010a) Peters, T., Klessen, R. S., Mac Low, M.-M., & Banerjee, R. 2010a, ApJ, 725, 134
- Peters et al. (2010b) Peters, T., Mac Low, M.-M., Banerjee, R., Klessen, R. S., & Dullemond, C. P. 2010b, ApJ, 719, 831
- Peters et al. (2017) Peters, T., Naab, T., Walch, S., et al. 2017, MNRAS, 466, 3293
- Pettitt et al. (2017) Pettitt, A. R., Tasker, E. J., Wadsley, J. W., Keller, B. W., & Benincasa, S. M. 2017, MNRAS, 468, 4189
- Pittard (2013) Pittard, J. M. 2013, MNRAS, 435, 3600
- Puls et al. (2008) Puls, J., Vink, J. S., & Najarro, F. 2008, A&A Rev., 16, 209
- Ragan et al. (2012) Ragan, S., Henning, T., Krause, O., et al. 2012, A&A, 547, A49
- Rathborne et al. (2006) Rathborne, J. M., Jackson, J. M., & Simon, R. 2006, ApJ, 641, 389
- Rey-Raposo et al. (2017) Rey-Raposo, R., Dobbs, C., Agertz, O., & Alig, C. 2017, MNRAS, 464, 3536
- Rey-Raposo et al. (2015) Rey-Raposo, R., Dobbs, C., & Duarte-Cabral, A. 2015, MNRAS, 446, L46
- Rogers & Pittard (2013) Rogers, H. & Pittard, J. M. 2013, MNRAS, 431, 1337
- Röllig et al. (2007) Röllig, M., Abel, N. P., Bell, T., et al. 2007, A&A, 467, 187
- Rosen et al. (2014) Rosen, A. L., Lopez, L. A., Krumholz, M. R., & Ramirez-Ruiz, E. 2014, MNRAS, 442, 2701
- Rybicki & Lightman (2004) Rybicki, G. B. & Lightman, A. P. 2004, Radiative Processes in Astrophysics (Wiley-Vch Verlag)
- Salpeter (1955) Salpeter, E. E. 1955, ApJ, 121, 161
- Schmidt (1959) Schmidt, M. 1959, ApJ, 129, 243
- Sedov (1958) Sedov, L. I. 1958, Reviews of Modern Physics, 30, 1077
- Seifried et al. (2017) Seifried, D., Walch, S., Girichidis, P., et al. 2017, MNRAS, 472, 4797
- Sembach et al. (2000) Sembach, K. R., Howk, J. C., Ryans, R. S. I., & Keenan, F. P. 2000, ApJ, 528, 310
- Shu (1977) Shu, F. H. 1977, ApJ, 214, 488
- Slyz et al. (2005) Slyz, A. D., Devriendt, J. E. G., Bryan, G., & Silk, J. 2005, MNRAS, 356, 737
- Smith et al. (2014a) Smith, R. J., Glover, S. C. O., Clark, P. C., Klessen, R. S., & Springel, V. 2014a, MNRAS, 441, 1628
- Smith et al. (2014b) Smith, R. J., Glover, S. C. O., & Klessen, R. S. 2014b, MNRAS, 445, 2900
- Spitzer (1978) Spitzer, L. 1978, Physical processes in the interstellar medium (A Wiley-Interscience Publication)
- Spitzer (1942) Spitzer, Jr., L. 1942, ApJ, 95, 329
- Strömgren (1939) Strömgren, B. 1939, ApJ, 89, 526
- Tielens (2005) Tielens, A. G. G. M. 2005, The Physics and Chemistry of the Interstellar Medium (Cambridge, UK: Cambridge University Press)
- Truelove et al. (1997) Truelove, J. K., Klein, R. I., McKee, C. F., et al. 1997, ApJ, 489, L179
- Turk et al. (2011) Turk, M. J., Smith, B. D., Oishi, J. S., et al. 2011, ApJS, 192, 9
- Vázquez-Semadeni et al. (2003) Vázquez-Semadeni, E., Ballesteros-Paredes, J., & Klessen, R. S. 2003, ApJ, 585, L131
- Vázquez-Semadeni et al. (2010) Vázquez-Semadeni, E., Colín, P., Gómez, G. C., Ballesteros-Paredes, J., & Watson, A. W. 2010, ApJ, 715, 1302
- Vázquez-Semadeni et al. (2007) Vázquez-Semadeni, E., Gómez, G. C., Jappsen, A. K., et al. 2007, ApJ, 657, 870
- Waagan (2009) Waagan, K. 2009, Journal of Computational Physics, 228, 8609
- Waagan et al. (2011) Waagan, K., Federrath, C., & Klingenberg, C. 2011, Journal of Computational Physics, 230, 3331
- Walch et al. (2015) Walch, S., Girichidis, P., Naab, T., et al. 2015, MNRAS, 454, 238
- Walch & Naab (2015) Walch, S. & Naab, T. 2015, MNRAS, 451, 2757
- Walch et al. (2012) Walch, S. K., Whitworth, A. P., Bisbas, T., Wünsch, R., & Hubber, D. 2012, MNRAS, 427, 625
- Wareing et al. (2017) Wareing, C. J., Pittard, J. M., & Falle, S. A. E. G. 2017, MNRAS, 470, 2283
- Weaver et al. (1977) Weaver, R., McCray, R., Castor, J., Shapiro, P., & Moore, R. 1977, ApJ, 218, 377
- Whitworth (1979) Whitworth, A. 1979, MNRAS, 186, 59
- Whitworth et al. (1994) Whitworth, A. P., Bhattal, A. S., Chapman, S. J., Disney, M. J., & Turner, J. A. 1994, MNRAS, 268, 291
- Wolfire et al. (1995) Wolfire, M. G., Hollenbach, D., McKee, C. F., Tielens, A. G. G. M., & Bakes, E. L. O. 1995, ApJ, 443, 152
- Wood & Churchwell (1989) Wood, D. O. S. & Churchwell, E. 1989, ApJS, 69, 831
- Wünsch et al. (2011) Wünsch, R., Silich, S., Palouš, J., Tenorio-Tagle, G., & Muñoz-Tuñón, C. 2011, ApJ, 740, 75
- Wünsch et al. (2018) Wünsch, R., Walch, S., Dinnbier, F., & Whitworth, A. 2018, MNRAS
- Zuckerman & Evans (1974) Zuckerman, B. & Evans, II, N. J. 1974, ApJ, 192, L149
Appendix A Mass evolution
Fig. 16 shows the time evolution of the total mass in MC (red) and MC (black) for simulation ZI_RAD (thin) and ZI_NOFB (thick) in the total volume, (solid), and the central subregion, (dashed).
Fig. 17 shows the density-temperature distribution of MC (top subpanel) and MC (bottom subpanel) at (top) and (bottom). The molecular hydrogen fraction (left) and ionized hydrogen fraction (right) is indicated by color.
Appendix B Stellar mass evolution
Fig. 18 shows the time evolution of the sink mass, , (top) and the total mass of massive stars, , (bottom) in the total volume of MC (red) and MC (black) for simulation ZI_NOFB (thick, only top) and ZI_RAD (thin). The initial phase of shows oscillations, which is due to two particles, which move out of the domain.
Fig. 20 shows the time evolution of the initial stellar mass (top), the distance (center) to the center of the small, central subregion of cloud MC (red) and MC (black) and the luminosity of Lyman continuum photons (bottom). Full markers and solid lines and empty markers and dashed lines indicate that the star is located in the central subregion or the rest of the cloud, respectively.
Fig. 19 shows the ratio between the stellar and sink mass at the time of the stellar formation, / , in MC (red) and MC (black) for simulation ZI_RAD in the total cloud (solid) and the central subregions (dashed). The mass fractions of the high mass range in respect to the underlying Kroupa IMF lies at 18 percent (green horizontal line). Hence, a stellar population (in a cloud) which satisfies this ratio represents the IMF well. Ratios above and below indicate that massive stars over-, respectively under-represented the IMF. After an initial massive star deficit, both central subregions and MC2 are well sampling the IMF. The high values in MC1 are caused by an initial 100 M star.
Appendix C Other Projections
Fig. 21 and 22 as well as Fig. 23, 22 show the time evolution of the gas column density in the --plane and in the - plane for the central subregion (, see Fig. 2 and Table 1) of MC and MC,respectively. The two columns show simulations ZI_NOFB (left) and ZI_RAD (right), respectively, at times to (from top to bottom). Circles indicate sink particles without an active stellar component, i.e. without massive stars. Star-shaped markers are cluster sink particles with active stellar feedback. The color of the markers indicate the sink age.