The effect of the virial state of molecular clouds on the influence of feedback from massive stars

TL;DR

This study uses SPH simulations to examine how the virial state of molecular clouds affects feedback from massive stars, revealing modest reductions in star formation efficiency and significant unbinding of material, primarily driven by momentum transfer rather than thermal expansion.

Contribution

It demonstrates that feedback effects can be effectively modeled as steady momentum injection, highlighting the role of virial ratios and escape velocities in feedback efficiency.

Findings

Star formation efficiency reduced by 30-50% due to feedback.

Unbound material ranges from 20-60%, especially in low escape velocity clouds.

Feedback effects are mainly due to momentum transfer, not thermal expansion.

Abstract

A set of Smoothed Particle Hydrodynamics simulations of the influence of photoionising radiation and stellar winds on a series of 10$^{4}$M$_{\odot}$ turbulent molecular clouds with initial virial ratios of 0.7, 1.1, 1.5, 1.9 and 2.3 and initial mean densities of 136, 1135 and 9096\,cm$^{-3}$ are presented. Reductions in star formation efficiency rates are found to be modest, in the range $30\%-50\%$ and to not vary greatly across the parameter space. In no case was star formation entirely terminated over the $\approx3$\,Myr duration of the simulations. The fractions of material unbound by feedback are in the range $20-60\%$, clouds with the lowest escape velocities being the most strongly affected.\\ Leakage of ionised gas leads to the HII regions rapidly becoming underpressured. The destructive effects of ionisation are thus largely not due to thermally--driven expansion of the HII regions, but to momentum transfer by photoevaporation of fresh material. Our simulations have similar global ionisation rates and we show that the effects of feedback upon them can be adequately modelled as a steady injection of momentum, resembling a momentum--conserving wind.

Tables

Table 1. Table 1: Initial properties of clouds. The meaning of the columns is as follows: 1: run identifier; 2: initial radius (pc); 3: initial virial parameter; 4: initial mean molecular hydrogen number density (cm -3 ); 5: initial root–mean–square 3D turbulent velocity (km s -1 ); 6: initial mean temperature (K); 7: initial mean Mach number; 8: initial freefall time at the mean density (Myr); 9: initial escape velocity (km s -1 ).

Run	R0(pc)	${\alpha }_{\mathrm{vir},0}$	$⟨n{\left(H_2\right)}_{,0}⟩$ (cm-3)	vRMS,0(km s-1)	$⟨T_0⟩\left(K\right)$	$⟨{ℳ}_0⟩$	tff,0 (Myr)	$v_{\mathrm{ESC},0}$(km s-1)
r07i	10.0	0.7	136	2.1	52	4.6	2.6	2.9
r07j	5.0	0.7	1135	3.0	30	8.5	0.9	4.1
r11o	5.0	1.1	1135	3.8	30	10.8	0.9	4.1
r15l	5.0	1.5	1135	4.4	30	12.5	0.9	4.1
r15m	2.5	1.5	9096	6.1	18	22.4	0.4	5.8
r19s	5.0	1.9	1135	4.9	30	13.9	0.9	4.1
r19t	2.5	1.9	9096	6.9	18	25.3	0.4	5.8
r23p	2.5	2.3	9096	7.6	18	27.9	0.4	5.8
r23q	5.0	2.3	1135	5.4	30	14.9	0.9	4.1

Equations

$$\displaystyle\begin{array}[]{rlrl}\gamma&=0.75;&\hfill\hfil&\rho\leq\rho_{1}\\ \gamma&=1.0;&\rho_{1}\leq&\rho\leq\rho_{2}\\ \gamma&=1.4;&\hfill\rho_{2}\leq&\rho\leq\rho_{3}\\ \gamma&=1.0;&\hfill\hfil&\rho\geq\rho_{3},\\ \end{array}$$

$${\dot{M}}_{\mathrm{leak}}\approx 4\pi R_{\mathrm{b}}^2{\rho }_{\mathrm{i}}{f}_{\mathrm{leak}}{c}_{\mathrm{II}}$$

$${\rho }_{\mathrm{i}}{k}_{\mathrm{B}}{T}_{\mathrm{i}}/{\mu }_{\mathrm{i}}=\frac{(1-f_{\mathrm{leak}}){\dot{M}}_{\mathrm{w}}{v}_{\infty }}{(1-f_{\mathrm{leak}})4\pi R_{\mathrm{b}}^2}$$

$${\dot{M}}_{\mathrm{i}}\approx f_{\mathrm{leak}}\frac{{\dot{M}}_{\mathrm{w}}{v}_{\infty }{c}_{\mathrm{II}}{\mu }_{\mathrm{i}}}{k_{\mathrm{B}}{T}_{\mathrm{i}}}$$

$$p_{\mathrm{tot}}=\sum _{\mathrm{i}}{m}_{\mathrm{i}}|{\mathbf{v}}_{\mathrm{i}}|,$$

$$p_{\mathrm{rad}}=\sum _{\mathrm{i}}{m}_{\mathrm{i}}\frac{{\mathbf{v}}_{\mathrm{i}}.{\mathbf{r}}_{\mathrm{i}}}{|{\mathbf{r}}_{\mathrm{i}}|},$$

$$p_{\mathrm{azi}}=\sum _{\mathrm{i}}{m}_{\mathrm{i}}{\left({\mathbf{v}}_{\mathrm{i}}.{\mathbf{v}}_{\mathrm{i}}-{\left[\frac{{\mathbf{v}}_{\mathrm{i}}.{\mathbf{r}}_{\mathrm{i}}}{|{\mathbf{r}}_{\mathrm{i}}|}\right]}^2\right)}^{\frac{1}{2}}$$

$$p_{\mathrm{phot}}(t)=2M_{\mathrm{i}}(t)c_{\mathrm{II}},$$

$$p_{\mathrm{grav}}(t)={\int }_0^t\frac{G{M}_{\mathrm{hm}}(t^{\prime }{)}^2}{R_{\mathrm{hm}}(t^{\prime }{)}^2}\mathrm{d}{t}^{\prime },$$

$$\frac{\mathrm{d}}{\mathrm{d}t}(M_{\mathrm{s}}{v}_{\mathrm{s}})=2⟨{\dot{M}}_{\mathrm{i}}⟩c_{\mathrm{II}},$$

$$M_{\mathrm{s}}={\int }_0^{R}4\pi r^2\rho (r)\mathrm{d}r=\frac{4\pi {\rho }_0}{(3+\alpha )r_0^{\alpha }}{R}^{3+\alpha }$$

$$n_{\mathrm{PE}}=2/(4+\alpha ),$$

$$k_{\mathrm{PE}}={\left[\frac{(4+\alpha )(3+\alpha )r_0^{\alpha }⟨{\dot{M}}_{\mathrm{i}}⟩c_{\mathrm{II}}}{4\pi {\rho }_0}\right]}^{\frac{1}{4+\alpha }}.$$

Full text

Leakage of ionised gas leads to the HII regions rapidly becoming underpressured. The destructive effects of ionisation are thus largely not due to thermally–driven expansion of the HII regions, but to momentum transfer by photoevaporation of fresh material. Our simulations have similar global ionisation rates and we show that the effects of feedback upon them can be adequately modelled as a steady injection of momentum, resembling a momentum–conserving wind.

††pagerange: LABEL:firstpage–References††pubyear: 2006

1 Introduction

Star formation is one of the most important subprocesses in the cycling of baryonic matter in the interstellar medium (ISM). Star formation acts as a sink of molecular material, and a source of ionised gas, metal–enriched gas, photons, cosmic rays, momentum and energy. It is one of the primary drivers of the evolution and dynamics of the ISM, and indeed the IGM, via stellar–feedback driven galactic outflows (e.g. Slyz et al., 2005; Dubois & Teyssier, 2008; Ostriker & Shetty, 2011; Girichidis et al., 2016; Gatto et al., 2016).

The various stellar feedback mechanisms influence both the process of star formation itself on the scale of molecular cores, clumps, and clouds, and the larger–scale dynamics of the ISM on the scales of spiral arms or galactic disks. While supernova are almost certainly the most important stellar contribution at these larger scales, many studies have shown that the effects of supernovae are modulated by the environments in which they explode (e.g. Slyz et al., 2005; Hennebelle & Iffrig, 2014; Übler et al., 2014; Gatto et al., 2015). The influence of supernovae thus depends on the prior action of other feedback mechanisms, in particular photoionisation (e.g. Whalen et al., 2008; Walch & Naab, 2015), and/or winds (e.g. Rogers & Pittard, 2013; Fierlinger et al., 2016) and, especially in very dense systems, radiation pressure (e.g. Fall et al., 2010; Skinner & Ostriker, 2015; Kim et al., 2016).

Ionisation also has more subtle large–scale influence owing to the fact that substantial quantities of ionising photons evidently escape from young clusters and their host clouds. This enabled massive stars to contribute to the reionisation of the Universe at high redshifts (e.g. Madau et al., 1999; Petkova & Springel, 2011; Pawlik et al., 2015), and supplies photons which help to maintain the thick layers of ionised gas observed to clothe galactic disks, such as the Milky Way’s Reynolds Layer (e.g. Reynolds et al., 1995; Wood & Mathis, 2004; Barnes et al., 2014).

Both supernovae and photoionisation are themselves indirectly influenced by feedback acting on still smaller scales, such as jets and outflows and thermal accretion feedback, since these regulate the rate of formation of the stars which are the sources of ionising photons (e.g. Krumholz et al., 2012; Federrath, 2015). While different types of feedback surely act at different scales, they are thus all connected.

In this paper, we examine the effects of ionisation and wind feedback from massive stars acting on molecular cloud scales. The clouds themselves are modelled simply as initially–smooth spheres with Gaussian density profiles and imposed turbulent velocity fields which are allowed to form a few massive stars. At this point each simulation is forked. A control simulation is allowed to continue as before, while in a dual–feedback counterpart, winds and photoionisation from the massive stars are enabled.

It is well known that molecular clouds have a range of virial ratios (e.g. Dobbs et al., 2011) and our previous work (e.g. Dale et al., 2014) focussed on clouds with initial virial ratios of either 0.7 or 2.3. These values exclude approximately the 20% most bound and 20% least bound of the clouds in the Heyer et al. (2009) sample analysed by Dobbs et al. (2011), and thus roughly the bracket the ‘typical’ 60% of clouds around the median virial ratio of $\approx 1.5$. In this work, we fill in the gap between these two bracketing values of the virial parameter by complementing our original calculations with clouds whose initial value of $\alpha$ is 1.1, 1.5 and 1.9, in order to examine the influence of the clouds’ initial energetic states on their response to feedback in a more graduated fashion.

The initial density and virial state of a cloud are likely the two most important quantities deciding how fast it will form stars, and we here model clouds with three different initial densities of 136, 1135 or 9096 cm*-3*. Raising the density decreases the freefall time, which straightforwardly increases the absolute rate at which the cloud converts gas to stars, as measured by the star–formation efficiency rate, $SFER$, the rate at which the star formation efficiency increases. However, the star formation efficiency rate per freefall time, $SFER_{\rm ff}$ should not depend on the density. Observations (e.g. Krumholz & Tan, 2007; Feldmann & Gnedin, 2011) suggest that the star formation rate per freefall may be constant over a wide range of densities. Results presented in Dale et al. (2014) (Figure 16 in that paper) imply that $SFER_{\rm ff}$ indeed correlates weakly with cloud density, but is a function of cloud virial ratio, with the more bound simulations exhibiting $SFER_{\rm ff}$ more than a factor of three higher than the less bound clouds. This reflects the role of turbulence in resisting, preventing or even reversing cloud collapse.

We focus only on lower–mass 104M*⊙* objects to make the simulation set as homogeneous as possible, with the object of modelling a representative sample of clouds of this mass covering a range of virial parameters and densities. The outline for the remainder of this paper is as follows. Section 2 briefly recapitulates our numerical methods, which are identical to previous papers. The results of the simulations are presented in Section 3 and discussed in Section 4. Our conclusions are summarised in Section 5.

2 Numerical methods

We use a hybrid N–body/Smoothed Particle Hydrodynamics (SPH) code based on that described by Benz (1990) and implementing point mass sink particles to treat star formation as in Bate et al. (1995). We present models of 104M*⊙* clouds using $10^{6}$ SPH particles, giving a mass resolution of $\approx 0.5-1.0$M*⊙. Sink particle accretion radii are fixed at 0.005 pc and the minimum creation density is $7\times 10^{7}$cm-3*. Sink particle mergers are not permitted and sink–sink gravitational interactions are smoothed at the accretion radius. We model the clouds as initially smooth spheres of gas with a mild Gaussian density profile seeded with a three–dimensional turbulent velocity field. The turbulence is initially divergence–free, satisfying $P(k)\propto k^{-4}$, and decays freely during the simulations.

As in our earlier work, the thermodynamics of the cold gas is treated in a simplified fashion using a piecewise barotropic equation of state from Larson (2005), defined so that $P=k\rho^{\gamma}$, where

[TABLE]

and $\rho_{1}=5.5\times 10^{-19}{\rm g\ cm}^{-3},\rho_{2}=5.5\times 10^{-15}{\rm gcm}^{-3},\rho_{3}=2\times 10^{-13}{\rm g\ cm}^{-3}$. At low densities, $\gamma$ is less than unity, implicitly ensuring that the gas in this regime has temperatures in excess of the canonical cloud temperature of $10$K. Dust cooling of the gas is modelled implicitly by the isothermal $\gamma=1.0$ segment, and the $\gamma=1.4$ segment represents the regime where dense collapsing cores become adiabatic. The final isothermal phase of the equation of state is simply in order to allow sink-particle formation to occur. The minimum gas temperature is set to 7.5K, fixing the relation between $\rho$ and $T$. All our simulated clouds have initial average densities $<\rho_{1}$, so that they lie in the line–cooling regime with mean gas temperatures in excess of 10K.

The individual initial properties of the simulated clouds are given in Table 1. Cloud radii are set to 2.5, 5.0 or 10.0 pc. The magnitudes of the turbulent velocities are normalised to control the virial ratios of the clouds. The gravitational potentials of each one are computed exactly using the SPH code’s gravity tree and the velocity fields set to give virial ratios (defined as the simple ratio of total kinetic energy in the centre–of–mass frame to total gravitational binding energy $E_{\rm kin}/|E_{\rm grav}|=\alpha_{\rm vir}=0.7,1.1,1.5,1.9,2.3$). Clouds with $\alpha_{\rm vir}=0.7$ and 2.3 have already been presented in previous papers (Dale et al., 2012, 2013a).

Calculations are evolved until three sinks with masses in excess of $20$M*⊙* have formed at which point they are forked into a control simulation which proceeds as before, and a dual–feedback run with ionisation and wind feedback enabled. Sinks exceeding the mass threshold are assigned ionising photon and wind momentum fluxes as described in Dale et al. (2012) and Dale et al. (2013c) respectively. The simulation pairs are then evolved for as close as is practicable to 3 Myr to determine the effects of pre–supernova O–star feedback.

The two feedback mechanisms considered – photoionising radiation and main–sequence stellar winds – are modelled in the same fashion as in previous papers. A ‘Strömgren volume’ technique is used to locate the ionisation fronts by performing reverse–ray–tracing on all gas particles and computing a discretised approximation to the directionally–dependent Strömgren integral given in Dale et al. (2007) to determine if sufficient ionising photons reach each particle during the current timestep to ionise them (if they are neutral), or to keep them ionised otherwise. Ionised particles are maintained at the canonical 104 K. If an ionised particle is deprived of photons, it is allowed to recombine on the appropriate timescale, after which it descends a cooling curve. The original algorithm of Dale et al. (2007) was updated to allow for the presence of multiple ionising sources, as described in Dale et al. (2012). The contribution of the recombination rate of a given intermediate particle to the Strömgren integral from a given source to a given target particle is multiplied by the ratio of the ionising flux reaching the intermediate particle from that source to the total flux arriving at that particle from all sources. Several iterations are then required to solve for the complete radiation field.

Winds are modelled approximately as described in Dale et al. (2013c). Unlike, for example, Rogers & Pittard (2013), we do not inject hot gas into the simulations, since this is technically difficult in SPH. We instead inject only the momentum, or ram–pressure, carried by the stellar winds, by locating the gas particles nearest each massive star which are struck first by large numbers of ‘momentum packets’ emitted by the sources in a Monte Carlo fashion. The quantities of momentum received by each of these particles are then converted to forces which are added to the other forces acting upon them. The thermal energies of the gas particles impacted by the winds are not modified. The assumption inherent in this method is that, when the fast–moving stellar winds collide with the ISM, the thermal energy generated by the resulting shocks is radiated away instantaneously. This clearly means that we are modelling the lower limit of the effects of the winds. While this is clearly an oversimplification, more sophisticated simulations (e.g Rogers & Pittard, 2013) have shown that much of the hot gas injected by winds in reality in fact leaks out of the host clouds, advecting the thermal energy along with it. Observational studies of several young star–forming regions by Rosen et al. (2014) also suggest that only a small fraction of the total injected wind energy is converted into work done on the surrounding clouds, and that much is lost by hydrodynamic leakage.

To make it easier to analyse the new simulations and the old synoptically, we adopt a uniform naming convention, where each run is given a four character name. The first (redundant) character is ‘r’, the middle two are digits representing the initial virial ratio of the cloud (07, 11, 15, 19 or 23), and the fourth is an identifying letter preserving the names of the older simulations, and adding new letters for the new runs.

3 Results

3.1 General morphology of clouds

Figure 1 shows column–density images viewed along the $z$–axis of all nine clouds at the same 25$\times$25 pc scale, with sink particles shown as white dots, just before feedback is initiated, showing the environments that the combined HII regions and wind bubbles first encounter. The most bound clouds are shown at the top left and least bound at the bottom right. The clouds all show to some degree a centrally–condensed filamentary structure, with the main concentrations of stars being around filament junctions or hubs near the centres of clouds and therefore at the deepest parts of their potentials. This configuration resembles that commonly observed in real molecular clouds (e.g. Myers, 2009; Schneider et al., 2012), and is therefore a good starting point for the study of the influence of feedback.

It should be noted, however, that a more realistic starting point would be a cloud formed in a larger (galactic) scale simulation, suitably re–resolved, as has been done recently by (e.g Rey-Raposo et al., 2015; Smilgys & Bonnell, 2016). While this is undeniable, it is not obvious that arbitrarily rescaling a cloud drawn from a larger simulation in order to study the influence of a particular parameter is an appropriate thing to do. The flexibility offered by the initial conditions used here therefore compensates for their simplicity.

Figure 2 shows the end results of the control–fork simulations, where feedback is neglected. In general, most of the clouds have expended to some degree, owing to the fact that most have them have virial ratios in excess of unity. It is also apparent that the clouds’ central regions have collapsed, resulting in more sharply–defined filaments, and greater concentrations of both stars and gas towards the cloud centres. This is largely a result of the more rapid draining of turbulent energy from the denser central regions of the clouds, as discussed in Dale et al. (2013a).

Figure 3 shows the end states of the dual–feedback simulations. Note that Figures 1, 2 and 3 show all the clouds at the same 25$\times$25 pc scale, with the same column–density colour scheme. Comparison of Figures 2 and 3 shows that feedback has clearly had a profound effect, both on the structure of the gas and on the distribution of the stars, in all simulations. Many of the clouds feature well–cleared but irregularly–shaped bubbles. The degree to which the clouds are cleared out appears to vary with the cloud size, in that the larger clouds appear to have had larger fractions of their volumes evacuated by a few coherent bubbles, whereas the smaller clouds, particularly Runs r07j and r11o and r15m, exhibit a jumble of smaller bubbles separated by surviving lanes of cold gas. Since these clouds are all the same mass, smaller size implies higher density, and it therefore appears that the denser clouds are more difficult for feedback to clear effectively.

Comparing Figures 2 and 3 also gives the general impression that the feedback calculations generally have stars spread out over a larger fraction of the cloud volumes, and less condensed in clusters. This is broadly consistent with the conclusion drawn in Dale et al. (2013b) that feedback redistributes star formation. We will return to this issue in detail in a companion paper.

3.2 Basic diagnostic properties

In Figure 4 we plot the time–evolution of a series of basic global diagnostic quantities for all the models, using a rainbow colour scheme where the more strongly bound clouds are towards the red end of the spectrum and the less bound towards the violet end, as an aid to interpretation. Solid lines represent control simulations (where appropriate) and dashed lines dual–feedback simulations.

The first two panels are concerned with the mass unbound from the clouds. The first panel shows for both control and dual–feedback simulations the total quantity of matter with positive energy in the cloud centre–of–mass frame as a function of time. Note that in all control simulations, some non–zero quantity of mass is always unbound. The amount of such material increases with virial ratio, as expected, but in a given simulation, does not increase very much with time over the timescales considered here, even in the most strongly unbound and rapidly expanding clouds.

Turning to the dual feedback simulations, clearly substantial amounts of additional mass are unbound by feedback, but no clouds are destroyed outright. In the second panel of Figure 4, we plot the unbound mass fractions in the dual–feedback calculations ‘corrected’ by subtracting at each timestep the fraction of material unbound in the corresponding control simulation. This gives an idea of how much extra material is expelled from the clouds by feedback, on top of what is already unbound due to the remnant turbulent velocity field and the cloud expansion driven by it. We see that 30–60$\%$ of cloud–mass is expelled by this definition, increasing to 40–80$\%$ if the lines were extrapolated to 3Myr. Since the lines of different colours are well mixed, there is no apparent correlation between the virial states of the clouds and how much harm is done to them by feedback.

The third panel of Figure 4 depicts the increase in star formation efficiency with time for all runs. The star formation efficiencies of the simulations span values between 10 and 50$\%$ and some of the control simulations would likely achieve efficiencies of $\approx 80\%$ if they were allowed to continue to 3Myr. In all simulations, feedback reduces the star formation rates and efficiencies, but not by large factors, in no case by more than a factor of two. There is once again no apparent correlation with the virial parameter.

The fourth panel of Figure 4 shows the increase in the numbers of stars with time for all calculations. Not only is there no connection between the number of stars each cloud produces and its virial state, but the effect of feedback on the number of stars formed can be in either direction, reflecting the multiple roles of feedback discussed in previous papers (e.g Dale et al., 2013b) of triggering, aborting and redistributing star formation. Recall that the effect of feedback on the total stellar mass is always negative.

The fifth and sixth panels of Figure 4 show the evolution of the 50$\%$ Lagrange radii and the three–dimensional velocity dispersion within this volume. All the clouds save r07j expand to some degree, with expansion in the control simulations being driven purely by the fading turbulent velocity fields, so that the more unbound clouds expand by larger factors. The expansion is, not surprisingly, substantially faster in the dual–feedback simulations. The velocity dispersions of the more bound control simulations grow with time for the more bound clouds, reflecting the gravitational collapse in the cloud cores, which also causes the star formation rates in these clouds to hold steady or accelerate somewhat. In the less bound clouds, the velocity dispersion declines as the clouds’ turbulence dies away and decreasing density leads to a steady diminution in star formation activity. In the dual–feedback simulations, the velocity dispersion increases very sharply as feedback is enabled, and then plateaus or even declines as the HII region expansion is slowed by mass–loading.

3.3 Variations across parameter space

It is instructive to plot the values at common epochs of some of the quantities discussed above in the [escape velocity:turbulent velocity dispersion] parameter space in which these simulations were conceived, and against the initial virial ratio of the clouds, which we do in Figures 5, 6 and 7. The escape velocity on the $x$–axis and the initial turbulent velocity dispersion on the $y$–axis are both normalised by the ionised sound speed of 11.3 km s*-1*. The different plot symbols represent different cloud initial densities as a reminder of the variation of this parameter – circles denote $\langle n(H_{2})_{,0}\rangle$=136 cm*-3*, squares $\langle n(H_{2})_{,0}\rangle$=1135 cm*-3* and squares $\langle n(H_{2})_{,0}\rangle$=9096 cm*-3*.

In Figure 5, we plot by symbol size and for the control simulations, the mean star formation rate per Myr (left panel) and per freefall time (centre panel) from the point where the simulations are forked until the ends of the control runs over the parameter space. The right panel instead plots $SFER_{\rm ff}$ directly against the cloud initial virial ratios.

The star formation efficiency rates per Myr and per freefall time vary by roughly an order of magnitude amongst the simulations. In general, the clouds with the highest densities and the lowest virial ratios form stars fastest. The relation between $SFER_{\rm ff}$ and $\alpha$ appears very close to a linear decline, albeit with some scatter.

Figure 6 depicts the same quantities shown using the same scales, but for the dual–feedback simulations in the period between the initiation of feedback and the ends of the simulations. In general, the same trends in these quantities across the parameter space and with the virial ratio are repeated here, except that the star formation efficiency rates, measured either in absolute times or per freefall times, are scaled down by factors of $\approx 30-50\%$. In the absence of feedback, most simulations achieve $SFER_{\rm ff}$ of $>5\%$ and two simulations push this over 10%. None of the dual–feedback simulations reach such high $SFER_{\rm ff}$ and all but two (Runs r07j and r11o) exhibit $SFER_{\rm ff}<5\%$. This demonstrates that, in these smaller clouds, it is possible for the combined effects of winds and ionisation alone to reduce star formation efficiency rates per freefall time to values factors of only a few larger than those observed. Similar values of $SFER_{\rm ff}$ were recently reported by Federrath (2015), who modelled the combined effects of driven turbulence, jets and magnetic fields.

As shown by Figure 4, in no case is star formation in whole clouds terminated, so that these clouds would continue towards ever larger star formation efficiencies absent some other effect which turns star formation off. However, feedback redistributes star formation to the peripheries of the systems and many of the simulations exhibit stellar clusters/subclusters which are devoid of gas and which are not accreting from the remains of their host clouds, so star formation within them has been locally terminated. This is an encouraging sign from the point of view of reproducing observations, since observations of very young massive clusters empty of gas and without ongoing star formation abound (e.g. Bastian et al., 2014; Hollyhead et al., 2015), and such systems appear to become gas–free on timescales of a few Myr. It should be noted, however, that most of the clusters studied in the cited works are much more massive than those produced here, which are at most a few$\times 10^{3}$M*⊙*.

In Figure 7, we plot in the same parameter space the ionisation fractions, unbound mass fractions, and the corrected unbound mass fractions. All quantities are extrapolated to a common time of 3 Myr after the initiation of feedback. We see that there is very little variation in the first of these, with $\approx 10\%$ of the mass of all clouds being ionised. The variation in the total unbound mass is also modest, with all simulations finishing with unbound mass fractions of 40–80%, and hints that, for a given escape velocity, clouds with higher turbulent velocities finish with larger unbound mass fractions. The corrected unbound mass fractions show a similar degree of variation and there is a trend that clouds with higher escape velocities suffer less damage, as we found in Dale et al. (2014).

4 Discussion

A striking result from Figure 7 is that the rate at which fresh neutral gas is ionised is very similar in all simulations. We suggest here a possible explanation for this observation. For most of the duration of the simulations, the configuration of the gas can be thought of as approximating a leaky bubble whose inner surface is being continuously ionised by OB stars located near the centre of the bubble, while ionised gas continually escapes through the holes in the walls. If the bubble subtends a fraction $(1-f_{\rm leak})$ of the sky as seen from the ionising sources, the rate at which ionised gas escapes from the bubble is approximately

[TABLE]

The density of the ionised gas can be inferred by assuming that is in pressure balance with the momentum flux of the stellar winds also emanating from the massive stars, so that

[TABLE]

If we now assert that the HII region is in a steady state, so that the rate at which material is being ionised is equal to the rate at which it is leaking from the bubble, we can combine the above two expressions to obtain

[TABLE]

If we insert typical numbers of $f_{\rm leak}\approx 0.5$, $\dot{M}_{\rm w}\sim 10^{-6}$M*⊙* yr*-1*, $v_{\infty}\approx 2\times 10^{3}$km s*-1*, $\mu_{\rm i}=2.3m_{\rm H}$, $T_{\rm i}=10^{4}$K, we obtain an ionisation rate of $\sim 3\times 10^{-4}$ $M_{\odot}$ yr*-1*, corresponding to a mass of $\sim 10^{3}$ M*⊙* over 3 Myr, which is close to the values we obtain from the simulations.

Walch et al. (2012) model ionising feedback in $10^{4}$ M*⊙* non–centrally–condensed fractal clouds of initial radius 6.4 pc and an initial mean density of $\approx 190$ cm*-3*. Their ionising sources have a common photon luminosity of $Q_{\rm H}=10^{49}$ s*-1*, comparable to the total ionising luminosities of our simulations. They find global ionisation rates rising to and plateauing at $\sim 1\times 10^{-3}$ $M_{\odot}$ yr*-1*, with relatively little difference between clouds with fractal dimensions in the range 2.0–2.8. Their mean cloud density is lower than all but that in Run r07i, but even that simulation does not achieve the ionisation rates seen by Walch et al. (2012). Taken in conjunction with our findings, their results imply that a turbulent density and velocity field is more difficult to ionise than a fractal density field of comparable mean density, which may well be due to the resistance to HII region expansion offered by the ram pressures in the turbulent flows in the neutral gas.

4.1 Momentum feedback from photoionisation

In Dale et al. (2013c) and Dale et al. (2014), we showed that the holes though which the ionised gas leaks cause the thermal pressure in the HII region to drop to much lower values than would obtain in a sealed bubble (see figures 18 and 5 in those papers respectively, and pertinent discussion). This being the case, the principal transfer of momentum to the cold gas due to ionisation at late stages of the simulations is not due to thermal pressure of the HII region, but to the rocket effect of new gas being ionised. The transfer of momentum resulting from heating $M_{\rm i}=1\times 10^{3}$ M*⊙* to $c_{\rm II}\approx 10$ km s*-1* is $2M_{\rm i}c_{\rm II}\approx 2\times 10^{4}$ M*⊙* km s*-1* (Matzner, 2002). We therefore examine in detail the uptake of momentum by the clouds.

4.2 Distribution of momentum

We compute for each particle the modulus of the total momentum, the signed value of the the total radial momentum, and modulus of the azimuthal (i.e. non–radial) momentum in the clouds’ centre of mass frames, and sum these quantities, distinguishing between cold neutral particles and hot ionised particles, as

[TABLE]

Note that, as defined, $p_{\rm tot}\neq p_{\rm rad}+p_{\rm azi}$.

In Figure 8 we plot the evolution of these quantities in time from three simulations covering the range of virial parameters under study here. Note that the control runs do not contain any hot gas. Unsurprisingly, as the virial ratio increases from run r07i through run r15l to run r23q, the momentum contained in the cold gas in the control simulation makes a larger and larger contribution to the total gas momentum in the corresponding feedback simulations, to the point where the momentum stored in the cold gas and hot gas in the dual–feedback run r23q are nearly the same at the end of the simulation.

4.3 Momentum loss and wastage

Substantial quantities of momentum are stored in azimuthal motions. While the amount of azimuthal momentum initially also increases strongly in the dual–feedback simulations, this quantity tends to decline as the simulations progress, in contrast to the radial momentum which always increases. Given that stellar winds and HII region expansion, be it thermally or photevaporatively driven, should principally drive radial motions, the driving of azimuthal flows demands some explanation.

We suggest two mechanisms for driving azimuthal motions, both of which lead to momentum that could drive cloud destruction to effectively be wasted:

(i) the feedback sources are not all at the same location, so that some component of the flows they drive much be non–radial with respect to the cloud, and some of the momentum they inject is consumed when oppositely–directed flows collide (illustrated in the left panel of Figure 9);

(ii) the complex structure of the gas means that momentum imparted to the cold gas by photoevaporation need not be parallel to the direction of photon propagation, which also results in the collisions of flows (illustrated by comparing the centre and right panels of Figure 9, which show the momentum transfer from a smooth and a structured shell, respectively).

This latter point is similar to that adduced in the driving of the ionisation–front instability (e.g Giuliani, 1979; Whalen & Norman, 2008). Here, however, the complex structures of the bubbles is more likely to be due to the complex density and velocity fields in which the HII regions are expanding.

In Figures 10 and 11, we show qualitatively that these two momentum–loss mechanisms do indeed operate in our simulations. Figure 10 shows an early stage of the r23q simulation, where two large bubbles driven by different subclusters are interacting with each other (the evolved remains of these bubbles are visible in the bottom–centre panel of Figure 3) The blue arrows show the velocity field of the cold gas driven by the bubbles, and it is clear that the velocities are oppositely–directly in the region squeezed between the bubbles.

Figure 11 instead shows a region on the edge of the left–hand bubble in the same simulation at a later time, with blue arrows again depicting the velocity field of the cold gas. It is clear from the crossing of the arrows at many locations that the generally radial flows have azimuthal components which often result in their convergence, leading to partial loss of the momentum injected by photoevaporation on the inside of the bubble walls.

4.4 Isolating momentum injected by feedback

As with several other quantities discussed here, it is useful to compute the differences in the values of the momenta derived above between the dual–feedback and control simulations, to isolate the changes in cloud evolution due to feedback alone. In Figure 12, we plot these relative quantities by subtracting the relevant momentum components of the control simulation from those of the dual–feedback simulation at each timestep, resulting in the total (blue solid lines), outward radial (green solid lines) and azimuthal (teal solid lines) momentum injected by feedback. For comparison, we also compute estimates for the cumulative momentum input from photoevaporation $p_{\rm phot}(t)$ (black dashed lines) as

[TABLE]

and from radial gravitational forces $p_{\rm grav}(t)$ (grey dashed lines) as

[TABLE]

where $M_{\rm hm}$ is half the total mass of stars and cold gas, $R_{\rm hm}$ is the radius from the system centre of mass containing $M_{\rm hm}$, and $t=0$ is the instant when feedback begins.

Figure 12 confirms that a large fraction of the momentum content of the dual–feedback simulations is stored in the expanding hot gas, and shows that the momentum so stored exceeds that taken up by the cold gas by factors of a few. The momentum contained in the expanding ionised gas exceeds that injected by photoevaporation, implying that the hot gas gains momentum from additional sources. The obvious candidates are thermal pressure forces, and momentum input from the stellar winds, which impact the hot gas lining the bubble interiors. We confirm that this is likely the case in Figure 13, which depicts the evolution of the momentum–differences in the ionisation–only r07i calculation. The momentum of the hot gas is clearly much more strongly influenced by the winds than that of the cold gas, so that the main effect of the winds is to help flush hot gas out of the cloud interiors.

Much of the momentum which has found its way into the simulations is thus contained in hot gas which has expanded to large radii and is no longer able to influence the evolution of the dense cold gas. The momentum of the cold gas at later times is instead better traced by the momentum input from photoevaporation, implying that this is indeed the main mechanism by which feedback influences the cold gas. However, at earlier times, the cold gas momentum does exceed $p_{\rm phot}(t)$, indicating that thermal pressure does contribute to the cold gas momentum budget at early times when the bubbles are still partially sealed. We find that the momentum uptake in the cold gas in these simulations (and the others, which are not shown in the interests of saving space) does indeed reach a common value of $\approx 2\times 10^{4}$ M*⊙* km s*-1*, as suggested above by multiplying the typical mass of ionised gas by twice the ionised sound speed.

We also see that $p_{\rm grav}$, while close to $p_{\rm phot}(t)$ in Runs r07i and r23q, is substantially larger in Run r15l (and in other simulations, particularly r07j and r15m, where the discrepancy is even greater). At first sight, this implies that momentum injected by gravitational forces in some of the simulations should completely overwhelm the effects of feedback. However, while the clouds for which this is case have lower unbound mass fractions, these fractions are still substantial ($>20$% in all cases), so that the fact that $p_{\rm grav}>p_{\rm phot}(t)$ does not imply that feedback is unable to unbind any material from a cloud. The reason that this criterion is unreliable is that both gravity and photoevaporation inject momentum in a highly non–uniform and in homogenous manner and comparing the gross quantities of momentum input can be misleading.

It is instructive to compare the total estimated momentum input from photoionisation with that from other feedback mechanisms. Winds are included in these simulations in a simple way as an injection of momentum. The most massive stars in our simulations have wind mass loss rates $\dot{M}_{\rm w}\approx 1\times 10^{-6}$ M*⊙* yr*-1*, and terminal velocity $v_{\infty}\approx 2\times 10^{3}$ km s*-1*. Over 3 Myr, the momentum injected by such a star is $\approx 6\times 10^{3}$ M*⊙* km s*-1*, a factor of a few less than that introduced by photoevaporation.

We have not included radiation pressure in our simulations, but these same massive stars have bolometric luminosities of $L_{*}\sim 10^{5}$ L*⊙, which would give rise to a corresponding radiative momentum flux $L_{*}/c\approx 6\times 10^{3}$ M⊙* km s*-1*, which is comparable to that injected by the winds, as might be expected. The momentum actually absorbed by the cloud depends critically upon the radiative trapping factor, which is still a matter of some controversy. Krumholz & Matzner (2009) estimate it to be 2–3, while Skinner & Ostriker (2015) measure from their radiation transport calculations a value of barely over unity at most, and substantially less than that throughout most of their simulated clouds. The upper end of the range suggested by Krumholz & Matzner (2009) would make radiation pressure of roughy equal importance as photoevaration, in terms of momentum deposition.

For the sake of completeness, we also note that a supernova explosion in which 10 M*⊙* of ejecta carry 1051 erg also carries $\approx 3\times 10^{4}$ M*⊙* km s*-1* of momentum. We conclude from the foregoing that momentum injected by photoevaporation is likely to be dominant at least until the detonation of supernovae and that, even after a few of these events, it will still constitute a major contribution to the momentum budget of GMCs which are in the process of disruption.

4.5 A simple model of photovaporation–driven evolution

The HII regions in these clouds are under–pressured due to leakage of hot gas from the bubbles, and the principal transfer of momentum to the clouds is via photoevaporation. Over the course of the simulations, the quantities of mass ionised are remarkably similar (as shown in Figure7, left panel) and the total momentum injected is also. One could in principle then regard the dynamical effect of photoionisation on the cold gas as similar to that of a steady momentum–driven wind.

If one treats the clouds as swept–up shells of mass $M_{\rm s}$ expanding at $v_{\rm s}$ and driven by an injection of momentum resulting from a photevaporation flow at their inner surfaces, one may follow a similar procedure to that used in the analysis of momentum–driven winds in Dale et al. (2013c) and write

[TABLE]

where $\langle\dot{M}_{\rm i}\rangle$ denotes the mean rate at which gas is ionised, and $c_{\rm II}$ is the sound speed in this material.

If the clouds have density profiles described by $\rho(r)=\rho_{0}(r/r_{0})^{\alpha}$, the shell mass is given by

[TABLE]

for $\alpha\geq-2$. Inserting this expression into Equation 14, setting $v_{\rm s}=\dot{R}$ and making the customary assumption (e.g. Lamers & Cassinelli, 1999) that $R(t)=k_{\rm PE}t^{n_{\rm PE}}$, $k$ being a constant and the subscript PE denoting photoevaporation, we obtain

[TABLE]

and

[TABLE]

While leakage of the ionised gas from the cloud is essential to deriving this expression since it allows the assumption of momentum conservation to be made, the hydrodynamic leakage factor as defined in Dale et al. (2014) does not appear in any of the expressions. Any dependence on the leakage factor would enter into the determination of the quantity $\langle\dot{M}_{\rm i}\rangle$, which is here treated as a constant.

In order to test the fidelity of the relations derived above, we evaluate $R(t)=k_{\rm PE}t^{n_{\rm PE}}$ making the following assumptions:

(i) $\alpha=-2$;

(ii) $r_{0}=R_{0}/100$;

(iii) $\rho_{0}$ is set by integrating $\rho(r)$ from $r=0$ to $r=R_{0}$ and setting the total cloud mass to $M_{\rm c}=10^{4}$M*⊙*;

(iv) $\langle\dot{M}_{\rm i}\rangle=\bar{f_{\rm i}}M_{\rm c}/t_{\rm SN}$ with $\bar{f_{\rm i}}$ the typical final global ionisation fraction, taken to be 0.1, and $t_{\rm SN}$ the timescale over which feedback acts, taken to be 3 Myr, giving $\langle\dot{M}_{\rm i}\rangle\approx 3.3\times 10^{-4}$M*⊙* yr*-1*.

Under these assumptions, note that $R(t)\propto t$.

To compare with the behaviour of the simulations, we identify the $N$ cold gas particles nearest the most massive star and compute their mean radius from the star and the interquartile range of this quantity. Following the time evolution of these quantities traces the expansion of the main bubble in each simulation. We do not use the ionisation fronts or the ionised gas to trace the bubble expansion because of the very irregular shape of the ionisation fronts and the leakage of ionised gas from the clouds to large radii. The choice of $N$ is not obvious, so we repeat the analysis using 0.1% and 1% of the $10^{6}$ particles in each simulation.

We plot the results for all simulations in Figure 14. Red coloured areas show the mean radius and its interquartile range for $N=10000$, and blue areas the same quantities for $N=1000$. The black solid line denotes the model. We note that the choice of $N$ does not affect the inferred bubble radius greatly. We also repeated the analysis using the median radius of the $N$ particles rather than the mean, but the resulting differences were negligible.

Overall, the ability of the model to explain the evolution of the simulations is quite good. For most of the evolution of most of the simulations, the assertion that $R(t)\propto t$ is well reproduced, and the gradient (i.e. the value of $k_{\rm PE}$) is also well modelled in most simulations. The denser clouds – r07j, r15m and r23p – are those in which the agreement is poorest. There are many reasons which could be advocated for this, perhaps the most obvious being that the model is spherically symmetric and assumes that the feedback–driving stars are located at $r=0$, which is not strictly true in any simulation and is least accurate in the denser clouds which are very irregularly cleared by feedback. However, we consider that the agreement shown in Figure 14 establishes that the assertion that the clouds are largely impacted by momentum imparted by photoevaporation flows is a reasonable one. We do not propose to try to make the model more complex in the hope of achieving closer congruence with the SPH simulations.

4.6 Observations of photoevaporation

These simulations suggest that photoevaporation is a process of prime importance in the disruption of molecular clouds and it is of paramount importance to test this assertion observationally by measuring photoevaporation rates directly. This can be done by measuring the electron density close to the ionisation front, in order to trace the most recently ionised gas, assuming that the gas streams away at the typical ionised sound speed of $\approx 10$ km s*-1*, and multiplying by a suitably–measured area. The advent of advanced integral–field units has recently made it possible to measure photoevaporation rates over significant areas of star–forming regions. McLeod et al. (2015) and McLeod et al. 2016 (submitted) recently employed MUSE at the VLT to measure the electron density via the $S^{23}$ sulphur line ratio for pillar structures respectively in M16, the Carina Nebula and NGC 3603. Their results clearly demonstrate the importance of photoevaporation in sculpting the remains of the natal molecular clouds in these systems, and they highlight that photoevaporation rates even for the isolated portions of the clouds represented by the pillars range over four orders of magnitude.

4.7 Physics missing from this work

There are several potentially important physical processes absent from the simulations presented here. Feedback from low– and intermediate–mass stars is neglected. As pointed out by Offner et al. (2009), this is the only form of feedback operative between the onset of star formation and the birth of the first OB–type stars. Roughly speaking, such feedback falls into two classes – thermal feedback from protstellar contraction, deuterium burning and the thermalisation of gravitational potential energy, and momentum–feedback from jets. The effects of the former are largely to suppress fragmentation and control the initial mass function (e.g Bate, 2009; Offner et al., 2009; Urban et al., 2010). The latter form is more effective in limiting the overall growth in stellar mass (e.g Hansen et al., 2012; Federrath, 2015). This is of particular interest in attacking the persistent problem of simulations forming stars too quickly. It is possible that including low–mass stellar feedback in addition to the forms of feedback modelled here would produce star formation rates more compatible with those observed.

Magnetic fields are likewise not treated here. These also have several possible modes of influence. At early stages in GMC evolution, they effectively provide an additional (non–isotropic) pressure, which slows cloud collapse somewhat (e.g Price & Bate, 2008). Magnetic fields may then provide an additional brake on the global star formation rate which, in combination with that due to low–mass stellar feedback, may buy more time for the massive stars to actually disrupt clouds and terminate star formation.

However, magnetic fields will also likely influence the severity of high–mass stellar feedback. Krumholz et al. (2007) and Peters et al. (2011) showed that smooth magnetic fields suppresses the expansion of HII regions, although Arthur et al. (2011) found that a disordered field had a much more modest influence, at least during the earlier stages of HII region expansion. Perhaps of more relevance here, Gendelev & Krumholz (2012) examined the effects of uniform fields on blister–type HII regions or champagne flows. They found that the magnetic field reduced the hydrodynamic kinetic energy of the flow, but that the field was capable of storing a quantity of energy which far outweighed this reduction. However, it is not clear how this energy could be later tapped. The effects of magnetic fields on very leaky bubbles growing in turbulent clouds, such as those presented here, is rather unclear.

There has been considerable recent interest in modelling radiation pressure, in particular in dense GMCs. Skinner & Ostriker (2015) find that radiation pressure from star clusters formed in clouds of 106 M*⊙* and 10 pc in radius, corresponding to a mean density of $\approx 5\times 10^{3}$ cm*-3*. They find that radiation pressure is able to terminate star formation and disrupt the clouds. Their models are less dense than our Runs r15l, r19t and r23p, implying that these simulations might be sensitive to radiation pressure, although the simple analytic calculation given above implied that the radiative trapping factor would need to be $>3$ for radiation pressure to be more importance than photoevaporation, whereas Skinner & Ostriker (2015) find that the trapping factor is rarely in excess of unity. The importance of radiation pressure on low–mass clouds such as those examined here is thus not obvious.

The expansion of HII regions and wind bubbles, particularly in declining density profiles, is likely to be subject to several hydrodynamic instabilities, such as the Rayleigh–Taylor, Vishniac (Vishniac, 1983), Franco (Franco et al., 1990), and thin–shell (Elmegreen, 1994) instabilities. These may be poorly resolved in SPH, since small–scale perturbations are smeared out by the smoothing process, retarding the growth of the instabilities.

Dale et al. (2009) demonstrated, using detailed comparisons between the SPH code used for the calculations presented here, the FLASH (Fryxell et al., 2000) grid code, and the analysis of (Elmegreen, 1994), that SPH codes in fact are able to model the thin–shell instability. However, in Dale et al. (2013c) (section 5.3) we showed that the steep radial density profiles of our clouds is likely to suppress this instability.

Of more relevance here, Franco et al. (1990) showed that HII region expansion in a density profile steeper than $\rho(r)\propto r^{-3/2}$ is guaranteed to be unstable because the expansion is unable to sweep up enough material to contain the ionisation front. The azimuthally–averaged density profiles of our model clouds are in the range $\rho(r)\propto r^{-3/2--2}$. Even if the clouds were initially smooth, one would thus expect them to be susceptible to this instability. However, the clouds are strongly inhomogeneous before the onset of feedback. This is only likely to make the Franco instability more violent, since it offers directions in which the ionisation front is able to accelerate more or less than average. The front therefore tends to burst out in directions in which the local radial density profile is steepest, or in which there is the least material to sweep up.

Our simulations do clearly exhibit this phenomenon, as shown by the early HII region morphologies depicted in Figure 1 of Dale et al. (2014). The ionisation fronts rapidly escape to large radii and ionised gas penetrates large fractions of the cloud volumes. The rapidity of this behaviour is likely to lessen the importance of other types of instabilities by washing out their effects. The existence of ionised gas throughout much of the cloud volume allows radiative feedback to influence the dynamics of a much large fraction of the cold gas in the clouds than would be allowed by the slower expansion of confined swept–up shells. Boneberg et al. (2015) showed that the velocity fields in the photoionised clouds do indeed possess some characteristics of turbulence, consistent with driving on scales of several pc.

5 Conclusions

We have presented simulations of the effects of ionisation and wind feedback on a set of nine 104 M*⊙* turbulent molecular clouds covering a range of initial virial ratios between 0.7 and 2.3 in order to evaluate the influence of the cloud virial ratio on the evolution of the clouds. Our main conclusions are as follows:

(i) Star formation efficiency rates measured absolutely and per freefall time vary considerably amongst the model clouds, being highest for clouds with the shortest freefall times and lowest virial ratios. The star formation rate per freefall time in particular declines linearly with increasing virial parameter. Feedback, in the main, does not change these relations substantially, but lowers the star formation efficiency rates by similar modest factors of 30–50% across the parameter space. This results in seven of the nine clouds exhibiting $SFER_{\rm ff}<0.05$.

(ii) Feedback fails to terminate star formation in any cloud, although several clouds sport (sub)clusters which are free of gas and not accreting, and where star formation has been locally stopped.

(iii) The fractions of the cloud masses unbound by feedback vary somewhat more across the parameter space, being confined to the range $20-60\%$, clouds with the lowest escape velocities being those most strongly affected, in agreement with our previous studies.

(iv) Depressurisation of the HII regions in these low–mass clouds by hydrodynamic leakage (Dale et al., 2013c, 2014) allows the influence of photoionisation to be adequately explained as being primarily doe to momentum deposition via the rocket effect, despite the fact that HII regions are generally regarding as a form of thermal feedback. Momentum deposition is non–uniform and its efficiency is somewhat limited by interactions of separate bubbles, and by crossing of flows at irregular ionisation fronts, leading to the deposition of momentum in azimuthal directions, increasing the probability of momentum cancellation in opposing flows.

(v) Since all the clouds exhibit similar time–averaged ionisation rates, the momentum injection rates are approximately constant and cloud evolution may be modelled as being due to momentum–driven bubbles expanding in power–law density profiles. The density profiles of our model clouds are well approximated by $\rho\propto r^{-2}$, resulting the expansion of the bubbles being linear in time.

6 Acknowledgements

This research was supported by the DFG cluster of excellence ‘Origin and Structure of the Universe’ (JED).