Massive star feedback in clusters: variation of the FUV interstellar radiation field in time and space

TL;DR

This study uses advanced radiation-hydrodynamics simulations to analyze how a massive star's feedback affects a turbulent cloud's evolution and the variability of the FUV radiation field impacting nearby objects.

Contribution

It introduces a detailed model of radiative feedback from a massive star in a turbulent cloud, including the first measurement of FUV field variation in a dynamic cluster environment.

Findings

75% of cloud mass disperses within 4.3 Myr

Ionized mass is lower in higher mass clouds despite larger ionized volume

FUV radiation field varies significantly over time, affecting protoplanetary disc evolution

Abstract

We investigate radiative feedback from a 34 M$_\odot$ star in a $10^4$ M$_\odot$ turbulent cloud using three-dimensional radiation-hydrodynamics (RHD) models. We use Monte Carlo radiative transfer to accurately compute photoionization equilibrium and radiation pressure, with multiple atomic species and silicate dust grains. We include the diffuse radiation field, dust absorption/re-emission, and scattering. The cloud is efficiently dispersed, with 75 per cent of the mass leaving the (32.3 pc)$^3$ grid within 4.3 Myr (1.1 $t_{ff}$). This compares to all mass exiting within 1.6 Myr (0.74 $t_{ff}$) in our previously published $10^3$ M$_\odot$ cloud. At most 20 per cent of the mass is ionized, compared to 40 per cent in the lower mass model, despite the ionized volume fraction being 80 per cent in both, implying the higher mass cloud is more resilient to feedback. The total Jeans-unstable mass increases linearly up to 1500 M$_\odot$ before plateauing after 2 Myr, corresponding to a core formation efficiency of 15 per cent. We also measure the time-variation of the far-ultraviolet (FUV) radiation field, $G_0$, impinging on other cluster members, taking into account for the first time how this changes in a dynamic cluster environment with intervening opacity sources and stellar motions. Many objects remain shielded in the first 0.5 Myr whilst the massive star is embedded, after which $G_0$ increases by orders of magnitude. Gas motions later on cause comparable drops which happen instantaneously and last for $\sim$ 1 Myr before being restored. This highly variable UV field will influence the photoevaporation of protoplanetary discs near massive stars.

Figures

Equations

$$n(a)\propto a^{-q}$$

$$\begin{split}G_{0}&=\frac{1}{H}\int_{$912\text{\,}\mathrm{\SIUnitSymbolAngstrom}$}^{$2400\text{\,}\mathrm{\SIUnitSymbolAngstrom}$}4\pi J_{\lambda}\operatorname{d\!}\lambda\\ &=\frac{1}{H}\frac{\epsilon}{\Delta tV}\sum_{\lambda=$912\text{\,}\mathrm{\SIUnitSymbolAngstrom}$}^{$2400\text{\,}\mathrm{\SIUnitSymbolAngstrom}$}\ell\end{split}$$

$${\tau }_{\text{FUV}}=\int {\kappa }_{\text{FUV}}\rho ds$$

$$\kappa_{\textrm{FUV}}=\frac{\int_{$912\text{\,}\mathrm{\SIUnitSymbolAngstrom}$}^{$2400\text{\,}\mathrm{\SIUnitSymbolAngstrom}$}(\kappa_{\textrm{abs}}+\kappa_{\textrm{sca}})\operatorname{d\!}\lambda}{\int_{$912\text{\,}\mathrm{\SIUnitSymbolAngstrom}$}^{$2400\text{\,}\mathrm{\SIUnitSymbolAngstrom}$}\operatorname{d\!}\lambda}$$

$${\lambda }_J=\sqrt{\frac{\pi c_s^2}{G\rho }}$$

Full text

††pubyear: 2018††pagerange: Massive star feedback in clusters: variation of the FUV interstellar radiation field in time and space–Massive star feedback in clusters: variation of the FUV interstellar radiation field in time and space

1 Introduction

Star formation occurs in giant molecular clouds (GMCs) ranging in mass from $\sim 10^{4}$ to ${10}^{6}\text{\,}\mathrm{M_{\astrosun}{}}$, with radii of a few to $100\text{\,}\mathrm{p}\mathrm{c}$ (Solomon et al., 1987; Heyer et al., 2009). The dynamics in these GMCs are likely dominated by massive O stars via radiative feedback (Matzner, 2002) occurring in two wavelength regimes: extreme-ultraviolet (EUV) photons with energy greater than $13.6\text{\,}\mathrm{eV}$ can ionize atomic hydrogen to form an H ii region. Far-ultraviolet (FUV) radiation between 5 and $13.6\text{\,}\mathrm{eV}$, to which ionized gas is optically thin, is absorbed by molecules, leading to them being photodissociated into their constituent atoms (Stecher & Williams, 1967). Furthermore, FUV absorption by dust grains can heat the grains themselves or heat the gas via photoelectric ejection (Draine, 1978; Wolfire et al., 1995; Hollenbach & Tielens, 1999). Massive stars are therefore surrounded by an H ii region of ionized gas at ${10}^{4}\text{\,}\mathrm{K}$, bounded by a photodissociation region (PDR) of atomic gas between 100 to $1000\text{\,}\mathrm{K}$. This is finally followed by molecular gas at $10\text{\,}\mathrm{K}$.

The thermal pressure gradient between the warm ionized gas and the cold neutral surroundings causes the H ii region to expand into the neutral medium. The dispersal of gas via expanding H ii regions has been well studied using analytical and numerical models, including two-dimensional studies of champagne flows (Whitworth, 1979; Tenorio-Tagle, 1979), filamentary configurations (Bodenheimer et al., 1979), and power-law density fields (Franco et al., 1990). These studies showed how ionizing feedback from massive stars can efficiently disperse their host molecular clouds, even with a small star formation efficiency (SFE) of a few per cent, as is observed in the Galaxy (Lada & Lada, 2003). However, the problem becomes more complex in three-dimensional, inhomogeneous gas distributions, as shown by Dale et al. (2005); unlike previous works, they concluded that clouds may not necessarily be destroyed by ionizing radiation, as a small fraction of the gas may effectively shield the rest of the cloud, carrying significant amounts of energy out of the system instead of distributing it evenly. Parameter studies by Dale & Bonnell (2011) and Dale et al. (2012, 2013) found that the degree of dispersal was closely coupled to the initial conditions of their simulations, such as cloud mass (see also Howard et al., 2017). Dispersal has also been shown to depend on initial gas surface density (Kim et al., 2018), morphology (Geen et al., 2018), and cluster luminosity (Geen et al., 2016, 2018). One problem is that even the most recent models simplify the radiative transfer, for example by using the on-the-spot approximation for recombination (thus neglecting the ionizing photons re-emitted by the gas), by using simple two-step temperature schemes for neutral and ionized hydrogen, or by neglecting dust microphysics such as absorption and scattering. These have a non-negligible effect on dynamics and morphology (see e.g. Ercolano & Gritschneder, 2011; Haworth & Harries, 2012; Haworth et al., 2015) and therefore must be accounted for to accurately model gas dispersal by radiative feedback.

In addition to dispersal on a grand scale, there is also the thermodynamical evolution of individual structures inside the H ii region. Proplyds – neutral clumps of density $\sim${10}^{5}\text{,}{\mathrm{cm}}^{-3}$$, size $\sim$ 100 to ${10}^{5}\text{\,}\mathrm{A}\mathrm{U}$, with ionized boundaries, and containing a circumstellar or protoplanetary disc – are seen in the vicinity of O stars, most notably near $\theta^{1}$ Ori C in the Orion Nebula Cluster (O’dell et al., 1993; Bally et al., 2000; O’dell, 2001). They have also been seen in other H ii regions such as NGC 3603 (Brandner et al., 2000), the Carina Nebula (Smith et al., 2003), Cygnus OB2 (Wright et al., 2012), and even around B stars with weaker UV fields (Kim et al., 2016). EUV photons ionize the outermost layer of a proplyd, resulting in a photoevaporative flow of ionized gas into the diffuse environment. FUV photons are able to penetrate further in, photodissociating and heating the molecular disc surface, which drives a neutral atomic wind from the PDR (as modelled by Johnstone et al. 1998 and measured spectroscopically by Henney & O’Dell 1999).

Many observational studies of different star-forming regions have found a positive correlation between disc frequency/size/mass and projected distance from ionizing sources (O’dell, 1998; Guarcello et al., 2007, 2010; Guarcello et al., 2016; Ansdell et al., 2017; Eisner et al., 2018). Hydrodynamic models indeed show a correlation between extent of disc photoevaporation and $G_{0}$ (Adams et al., 2004; Clarke, 2007; Anderson et al., 2013; Haworth et al., 2017; Winter et al., 2018; Haworth et al., 2018), where $G_{0}$ is the FUV flux in units of the Habing (1968) field ($1.63\text{\times}{10}^{-3}\text{\,}\mathrm{erg}\text{\,}{\mathrm{s}}^{-1}\text{\,}{\mathrm{cm}}^{-2}$). Combined with simulations which show that disc destruction may occur primarily through external photoevaporation rather than stellar encounters (Scally & Clarke, 2001; Winter et al., 2018), the interstellar radiation field (ISRF) appears to be of crucial importance in determining the evolution of disc properties. The above studies explore the parameter space of $G_{0}$ from a solar neighbourhood value of 1 to an inner ONC value of $10^{6}$, depending on an assumed distance to an FUV source. Real proplyds and discs are located within gaseous, dusty, and turbulent clouds. Massive stars disperse material, as described above, and this alters the column density between them and any proplyds. Furthermore, the cluster members have their own motions, meaning the separations also change. Stellar velocity dispersions are of the order $1\text{\,}\mathrm{km}\text{\,}{\mathrm{s}}^{-1}$ ($=$1\text{,}\mathrm{pc}\text{,}{\mathrm{Myr}}^{-1}$$), proplyds can be propelled to $5\text{\,}\mathrm{km}\text{\,}{\mathrm{s}}^{-1}$, and the ionized sound speed is $10\text{\,}\mathrm{km}\text{\,}{\mathrm{s}}^{-1}$ – therefore, over the course of a proplyd lifetime, there will be non-negligible changes in the incident flux due to the kinematics of both the stars and the gas; both effects combined will vary the geometric dilution ($\propto r^{-2}$ for separation $r$) and the FUV attenuation ($\propto\textrm{e}^{-\tau}$ for optical depth $\tau$) as a function of time and distance. This is the behaviour we aim to describe in this paper.

In Ali et al. (2018, hereafter Paper I), we modelled a ${10}^{3}\text{\,}\mathrm{M_{\astrosun}{}}$ cloud with a $34\text{\,}\mathrm{M_{\astrosun}{}}$ massive star, using a detailed radiative transfer scheme to compute photoionization and radiation pressure – this included the diffuse radiation field, scattering and absorption by dust, and calculations of thermal equilibrium with cooling from recombination lines, forbidden lines, and free-free continuum. In this paper, we apply the same scheme to a 10-times more massive cloud, using the same initial mean surface density and stellar luminosity. We outline the numerical methods in section 2 and the initial setup in section 3. We present and discuss the results in section 4 (for the dynamics) and section 5 (for the FUV ISRF). Finally, we draw our conclusions in section 6.

2 Numerical methods

We use the Monte Carlo (MC) radiative transfer (RT) and hydrodynamics (HD) code, torus (described in detail by Harries et al., 2019). We use the same treatment as in Paper I, but provide a summary here. The MCRT technique follows the Lucy (1999) method, whereby energy packets representing photons are propagated through a medium and undergo scattering and absorption/re-emission events. This includes the radiation field from stars as well as the diffuse field from gas and dust. Photon wavelengths are interpolated from 1000 logarithmically spaced bins between ${10}^{2}\text{\,}$ and ${10}^{7}\text{\,}\mathrm{\SIUnitSymbolAngstrom}$. Silicate dust grains (Draine & Lee, 1984) are distributed according to a standard ISM power-law density relationship (Mathis et al., 1977)

[TABLE]

using grain sizes $a$ between 0.1 to $1\text{\,}\mathrm{\SIUnitSymbolMicro m}$ and a power-law index $q$ of 3.5, giving a median grain size of $0.12\text{\,}\mathrm{\SIUnitSymbolMicro m}$. The dust-to-gas mass ratio is kept constant at 0.01. The dust temperature is calculated from the balance between photon absorption and thermal Planck emission. This is calculated separately from the gas temperature, but an additional term accounts for collisional heat exchange between the two (Hollenbach & McKee, 1979). For the gas, we calculate photoionization balance for H i–ii, He i–iii, C i–iv, N i–iii, O i–iii, Ne i–iii, and S i–iv. Thermal equilibrium for the gas temperature takes into account heating due to photoionization of H and He, and cooling from H and He recombination lines, collisionally excited metal forbidden lines, and free–free continuum (Haworth & Harries, 2012).

Radiative transfer is calculated before every hydrodynamics step, which evolves isothermally. Self-gravity is included, using a V-cycling multigrid method to solve Poisson’s equation with Dirichlet boundary conditions based on a multipole expansion. Radiation pressure is calculated using the momentum-transfer algorithm presented in Harries (2015). Sink particles (see Harries, 2015) are included but do not accrete material or form on-the-fly (see section 3 for the formation method); however, they do evolve in mass, effective temperature, and luminosity using Schaller et al. (1992) evolutionary tracks. O stars emit radiation according to ostar2002 (Lanz & Hubeny, 2003) spectra, while later-type stars use Kurucz (1991) spectra. They move due to N-body interactions and the mutual gravitational force from gas. These are used as tracer particles to follow the behaviour near the possible locations of actual stars or proplyds.

The FUV radiation field is calculated in each cell as part of the MCRT algorithm, using

[TABLE]

given in units of the Habing (1968) field, $H=$1.63\text{\times}{10}^{-3}\text{,}\mathrm{erg}\text{,}{\mathrm{s}}^{-1}\text{,}{\mathrm{cm}}^{-2}$$. $J_{\lambda}$ is the mean intensity, $V$ is the cell volume, and $\ell$ is the path length travelled between events (scattering/absorption/cell-boundary crossing) by a MC photon packet with energy $\epsilon$ during time step $\Delta t$. This is calculated in every cell at all time steps. To track how it can vary at the location of a particular sink particle (representing a star or proplyd), we follow $G_{0}$ in the cell containing that sink. When a sink moves from one cell to another, the intensity may jump up or down due to the discrete nature of the grid structure (and hence the optical depth). Therefore, in the results presented in section 5, we average around a radius of 2.5 cells; if accretion were enabled, this would represent the computational accretion radius, but is simply a convenience factor here. For section 5, we integrate the optical depth along the ray joining two sinks as a measure of the dust attenuation; that is,

[TABLE]

where $\rho$ is density and $\kappa_{\textrm{FUV}}$ is a representative average of the dust opacity in each cell given by

[TABLE]

for absorption and scattering opacities $\kappa_{\textrm{abs}}$ and $\kappa_{\textrm{sca}}$, respectively. It should be noted that the radiative transfer uses the actual opacities, not this average value.

3 Initial conditions

We model a ${10}^{4}\text{\,}\mathrm{M_{\astrosun}{}}$ cloud which is initially spherical with a mean surface density of $0.01\text{\,}\mathrm{g}\text{\,}{\mathrm{cm}}^{-2}$ and radius $R_{s}=$8.41\text{,}\mathrm{pc}$$. The mass volume density is uniform in the inner core ($4.9\text{\times}{10}^{-22}\text{\,}\mathrm{g}\text{\,}{\mathrm{cm}}^{-3}$ between $0>r>R_{s}/2$), then falls off as $r^{-1.5}$ (reaching $1.7\text{\times}{10}^{-22}\text{\,}\mathrm{g}\text{\,}{\mathrm{cm}}^{-3}$ at $r=R_{s}$). The density outside the cloud is 1 per cent of the density at the sphere boundary. These conditions are similar to Paper I, where the sphere mass and radius were ${10}^{3}\text{\,}\mathrm{M_{\astrosun}{}}$ and $2.66\text{\,}\mathrm{pc}$, respectively. The mean surface densities of both spheres are the same. We refer to the ${10}^{4}\text{\,}\mathrm{M_{\astrosun}{}}$ cloud presented in this paper as the ‘high-mass cloud’, and the ${10}^{3}\text{\,}\mathrm{M_{\astrosun}{}}$ cloud in Paper I as the ‘low-mass cloud’.

The free-fall time associated with a uniform sphere with density $3M_{s}/4\pi R_{s}^{3}=$3\text{\times}{10}^{-22}\text{,}\mathrm{g}\text{,}{\mathrm{cm}}^{-3} is $\langle t_{\textrm{ff}}\rangle=$3.86\text{\,}\mathrm{Myr} (where $M_{\textrm{s}}$ is the total sphere mass). The temperature is $10\text{\,}\mathrm{K}$ everywhere for both gas and dust until radiation sources are switched on. The grid size from end to end is $32.3\text{\,}\mathrm{pc}$, approximately 4 times the sphere radius, giving a linear resolution of $0.13\text{\,}\mathrm{pc}$ per cell with $256^{3}$ cells. The grid structure is Cartesian, uniform, and fixed.

As in Paper I, the sphere evolves under self-gravity and a seeded turbulent velocity field for 0.75 $\langle t_{\textrm{ff}}\rangle$ without stars. We use the same random Gaussian turbulent velocity field as Bate et al. (2002), with a power spectrum $P(k)\propto k^{-4}$ for wavenumber $k$, such that the kinetic energy equals the gravitational potential energy, i.e. the virial parameter $\alpha_{\textrm{vir}}\equiv 2E_{\textrm{kin}}/E_{\textrm{grav}}=2$. At 0.75 $\langle t_{\textrm{ff}}\rangle$, a $33.7\text{\,}\mathrm{M_{\astrosun}{}}$ star is placed in the cloud’s most massive core. 28 other stars are inserted according to a probability density function (PDF) proportional to star formation rate $\dot{M}(\textbf{{r}})$ at some position r; that is, $p(\textbf{{r}})\propto\dot{M}(\textbf{{r}})\propto\rho(\textbf{{r}})/t_{\textrm{ff}}\propto\rho(\textbf{{r}})^{1.5}$. We keep the same number of stars and stellar masses as Paper I in order to compare how the same luminosity affects a more massive cloud. This means the imposed star formation efficiency in the high-mass model is 1 per cent as opposed to 10 per cent in the low-mass model. The initial velocity of each star is the velocity of the gas in the cell containing that star. The initial positions (plotted in projection in Fig. 1) are not necessarily the same as in Paper I, as they are placed using the aforementioned PDF without any further constraints; however, since the turbulence in the initial hydrodynamics-only phase results in a similar morphology at 0.75 $\langle t_{\textrm{ff}}\rangle$, when stars are created, the most massive star ends up being placed in roughly the same area of the cloud. This keeps the initial conditions broadly consistent, but just scaled up in cloud mass. The parameter space of cloud mass, radius, and velocity dispersion is shown in Fig. 2, alongside observations of GMCs in the Galaxy by Heyer et al. (2009) and the suite of simulations by Dale et al. (2012).

4 Dynamics

In this section we present and discuss the results of the RHD simulations for the ${10}^{4}\text{\,}\mathrm{M_{\astrosun}{}}$ cloud. One includes photoionization and radiation pressure, and one has photoionization only. The models are evolved for $\sim$ $4\text{\,}\mathrm{Myr}$ (1 $\langle t_{\textrm{ff}}\rangle$), with $t=0$ corresponding to the initiation of feedback. Results for the ${10}^{3}\text{\,}\mathrm{M_{\astrosun}{}}$ model (‘low-mass cloud’) are presented in Paper I; this was evolved under feedback for $1.6\text{\,}\mathrm{Myr}$ (0.74 $\langle t_{\textrm{ff}}\rangle$).

4.1 Bulk grid properties

Fig. 3 shows the grid properties as a function of time. The first panel plots the total mass, which starts just under $1.2\text{\times}{10}^{4}\text{\,}\mathrm{M_{\astrosun}{}}$ (including the low-density material outside the ${10}^{4}\text{\,}\mathrm{M_{\astrosun}{}}$ turbulent sphere); between 1 and $4\text{\,}\mathrm{Myr}$, there is a steady decrease as mass leaves the grid. The total mass flux at the grid boundaries, shown in the second panel of Fig. 3, peaks between 1.5 and $2.5\text{\,}\mathrm{Myr}$ (0.39 and 0.65 $\langle t_{\textrm{ff}}\rangle$ respectively) with a value of $4.7\text{\times}{10}^{-3}\text{\,}\mathrm{M_{\astrosun}{}}\text{\,}{\mathrm{yr}}^{-1}$. This is more than double the mass flux in the $1000\text{\,}\mathrm{M_{\astrosun}{}}$ model presented in Paper I. The shape of the curve is essentially the same, with the first phase being dominated by the removal of low-density material, peaking at the crossing time for ionized gas travelling from the centre of the grid to the edge, with a sound speed $\approx$10\text{,}\mathrm{km}\text{,}{\mathrm{s}}^{-1}$$; the decreasing second phase is overlaid with sharp, short-lived spikes corresponding to the removal of dense clumps.

At $1.75\text{\,}\mathrm{Myr}$ (0.45 $\langle t_{\textrm{ff}}\rangle$), the highest value of ionized mass is reached, with $1600\text{\,}\mathrm{M_{\astrosun}{}}$, or about 17 per cent, of the mass being ionized; this proportion is half that of the low-mass model, which peaked at 40 per cent, again at the time of peak mass flux. This is despite about 80 per cent of the volume being ionized, which is comparable to the low-mass model. The conclusion reached in that model is even more so the case here – most of the mass remains in small, dense, neutral areas.

The final frame of Fig. 3 shows the total mass in Jeans-unstable cells.It should be noted that dust-heating by (F)UV absorption raises the temperature from $10\text{\,}\mathrm{K}$ to a high of $\sim$ $100\text{\,}\mathrm{K}$, which results in a higher Jeans mass and therefore greater resistance against fragmentation. The figure represents an upper limit to the amount of mass that could fragment to form stars. Only a third of a core mass may actually be converted into stars (Alves et al., 2007). Furthermore, to prevent artificial numerical fragmentation, the Jeans length

[TABLE]

should be resolved by at least four grid cells (Truelove et al., 1997). This limit corresponds to a density $\sim${10}^{-20}\text{,}\mathrm{g}\text{,}{\mathrm{cm}}^{-3}$$.

The total Jeans-unstable mass increases linearly to $1500\text{\,}\mathrm{M_{\astrosun}{}}$ before plateauing at around $2\text{\,}\mathrm{Myr}$ (0.5 $\langle t_{\textrm{ff}}\rangle$). Dense clumps start leaving the grid after this stage, which results in the unstable mass decreasing. In terms of the original cloud mass, this represents a core formation efficiency of 15 per cent, or a star formation efficiency of 5 per cent. This order of magnitude is in line with observations in the Galaxy (Lada & Lada, 2003). It also agrees with the simulation by Geen et al. (2017), labelled ‘L’, which has a similar initial mass and surface density, and includes ionization feedback – Geen et al. (2017) found that the SFE for that cloud was closer to observations of local star-forming regions compared to denser models, where SFEs approached 100 per cent.

The bulk grid properties for the model without radiation pressure are in general similar to the combined feedback model, although there are minor differences. For example, there is a difference of $300\text{\,}\mathrm{M_{\astrosun}{}}$ for the total Jeans-unstable mass. Since the inferred core- or star-formation efficiency is a rough estimate, this difference may not be significant, but indicates the possibility of radiation pressure producing denser structures.

Fig. 5 shows column density histograms for the first $4\text{\,}\mathrm{Myr}$. Compared to the low-mass cloud in Paper I, they have similar shapes but are more homogeneous with time, only significantly shifting to lower surface densities beyond $3\text{\,}\mathrm{Myr}$. This model starts off with more material in the low-density end, at $\Sigma=$e-4.5\text{,}\mathrm{g}\text{,}{\mathrm{cm}}^{-2}, whereas the low-mass cloud stops at $\Sigma=${10}^{-4}\text{\,}\mathrm{g}\text{\,}{\mathrm{cm}}^{-2}; therefore this model requires less of a left-ward shift to reach the $\Sigma=${10}^{-5}\text{,}\mathrm{g}\text{,}{\mathrm{cm}}^{-2} peak by the end of the simulation. There is very little variation during the first $1.5\text{\,}\mathrm{Myr}$ (the total run time for the low-mass model), other than the spreading out of the low-density spike ($\sim${10}^{-4}\text{\,}\mathrm{g}\text{\,}{\mathrm{cm}}^{-2}) corresponding to the diffuse gas outside the cloud.

4.2 Dispersal

As in Paper I, the dense filaments on the left side of the cloud are more resistant to being dispersed than the right side. In particular, a cavity is easily blown out in the bottom right quadrant of the $x$-$y$ plane in the first $\mathrm{Myr}$, allowing ionized material to stream out. This is also the case at $(x=0,y>0)$ and $(x\approx$3\text{,}\mathrm{pc}$,y>0)$. It is not surprising that these locations are where the H ii breaks out, as this also happens in the low-mass cloud – both have a similar gravo-turbulent structure and a massive star in similar locations.

However, in this massive cloud, the dense filaments on the left side ($x<0$ in the $x$-$y$ plane) have more time to accumulate mass before the H ii region can reach them. Indeed, the filaments are still collapsing towards the centre while feedback progresses. This makes them less prone to dispersal, whether by ionization heating or the rocket effect – by $t=$4\text{,}\mathrm{Myr}$$, many of the globule–pillar objects are still present inside a $8\text{\,}\mathrm{pc}$ radius around the massive star, which corresponds to the radius of the entire grid in Paper I.

The initial conditions are similar to run I of Dale et al. (2012) (the closest point to this model in Fig. 2). Although their model starts off with a smoother density distribution than our model, a similar butterfly-shaped H ii region is produced – this is the only model in their suite of conditions to create this morphology. Run I showed the greatest effect in terms of structure and dynamics, with 60 per cent of the gas being unbound. Their total ionization fraction did not exceed 10 per cent, whereas the model here approaches 20 per cent, even with similar ionizing fluxes ($\sim${10}^{49}\text{,}{\mathrm{s}}^{-1}$$).

Dale et al. (2013) summarised their parameter study of GMC masses between $10^{4}$ to ${10}^{6}\text{\,}\mathrm{M_{\astrosun}{}}$, stating that the ionized mass fraction is around 5 to 10 percent regardless of the cloud properties. This disagrees with the models presented here and in Paper I, which approach fractions of 20 per cent and 40 per cent, respectively; that is, they are both higher and they vary with cloud mass. This is despite their models having ionizing luminosities which increase with cloud mass. One contributor to the difference between models may be the diffuse radiation field, which the Dale et al. models not include. As shown by Ercolano & Gritschneder (2011) and Haworth & Harries (2012), this can penetrate into what would otherwise be shadowed regions behind dense gas, creating changes in morphology such as detaching pillar heads. If this makes the environment more permeable to ionizing photons, the total ionized mass could be increased. This bears more investigation as the method for forming stars differs between the two studies, as does the gas distribution when stars begin to radiate, and these differences can also change the effectiveness of feedback (Geen et al., 2018). Dale et al. (2013) also find a strong dependence on dispersal efficiency with cloud escape velocity. $v_{\textrm{esc}}$ for the cloud presented here is still very much below the ionized sound speed ($\sim$3\text{,}\mathrm{km}\text{,}{\mathrm{s}}^{-1}$$ compared to $10\text{\,}\mathrm{km}\text{\,}{\mathrm{s}}^{-1}$, respectively), but the full parameter space is yet to be explored.

The ineffectiveness of radiation pressure as a dispersive feedback mechanism is in agreement with the analytical model of Fall et al. (2010), who showed that a ${10}^{4}\text{\,}\mathrm{M_{\astrosun}{}}$ cloud requires at least ten times greater surface density than the cloud we model here for radiation pressure to become dominant. In their numerical models, Kim et al. (2018) found that surface density was the important parameter in determining the feedback efficiency; photoionization dominated over radiation pressure for similar initial conditions to ours, and in general the latter only became relevant once the surface density exceeded a few times that of our cloud. Skinner & Ostriker (2015) showed that the radiation pressure from dust-processed infrared photons was ineffective even in high-mass, high-density conditions ($M>${10}^{5}\text{,}\mathrm{M_{\astrosun}{}}$,\Sigma>$0.3\text{,}\mathrm{g}\text{,}{\mathrm{cm}}^{-2}$$), except for dust opacities corresponding to metallicities higher than found in Galactic GMCs (e.g. a few times Solar). This latter regime was studied by Tsang & Milosavljević (2018) in a ${10}^{7}\text{\,}\mathrm{M_{\astrosun}{}}$ turbulent box, where radiation pressure only modestly reduced the star formation efficiency and rate, and did not disperse enough gas to stop star formation completely (helped by the porous gas distribution, in contrast with the Fall et al. (2010) model). In short, the relatively low surface density and metallicity of our cloud (as well as the low luminosity) results in radiation pressure being unimportant as a source of feedback as it pertains to gas dispersal. For future models, we aim to include sink particle formation and accretion and therefore will be able to comment on how these feedback mechanisms affect star formation directly.

4.3 Temperature and electron density

Fig. 6 and Fig. 7 show the average temperature and electron density, respectively, of the ionized gas, using two weighted averages: (a) $w=n_{e}(n_{\text{H\,{ii}}}+n_{\text{He\,{ii}}})\approx n_{e}^{2}$; (b) cell mass if hydrogen in the cell is more than 90 per cent ionized or $w=0$ if less. The results follow a similar pattern as the low-mass cloud, with both weighted-averages of the temperature decreasing down to $8000\text{\,}\mathrm{K}$ by the end of the $4\text{\,}\mathrm{Myr}$. The electron density goes down to $0.3\text{\,}{\mathrm{cm}}^{-3}$, again closely matching the low-mass cloud. This average density corresponds to a proton mass density of $5\text{\times}{10}^{-25}\text{\,}\mathrm{g}\text{\,}{\mathrm{cm}}^{-3}$, which is five times more than the simulation floor density of ${10}^{-25}\text{\,}\mathrm{g}\text{\,}{\mathrm{cm}}^{-3}$. The temperature and density (and thus thermal pressure) by the end of the simulation are as expected for the warm ionized medium (Wolfire et al., 1995).

5 FUV interstellar radiation field

Fig. 8 shows $G_{0}$ as a function of distance from the most massive star, at eight snapshots in time. Stars are labelled by a number going from 1 to 29 (where 1 is the most massive star, and 29 is the least massive). The corresponding (initial) stellar masses are labelled in Fig. 9. The straight line in Fig. 8 shows the purely geometric dilution of $G_{0}$ from the cell containing the massive star. There is less scatter after the first Myr, as the massive star disperses gas from the cluster, removing sources of FUV opacity, resulting in most of the stars falling along the straight line. There are some exceptions – stars 2 to 8 (B stars from $11\text{\,}\mathrm{M_{\astrosun}{}}$ to $3.6\text{\,}\mathrm{M_{\astrosun}{}}$) produce their own modest FUV luminosity which raises them above the massive star dilution line. In particular, star 2 has its own radius of influence ($\sim$0.5\text{,}\mathrm{pc}), in which some of the lower-mass stars (e.g. stars 19, 20, and 26) receive higher fluxes. Beyond this, star 1 dominates. At $t=$3\text{\,}\mathrm{Myr}, star 10 is a factor of a few below the dilution line, due to a higher optical depth between it and star 1; this effect is more clearly seen in the plots of $G_{0}$ vs time, as discussed below.

The time-variation of the flux is plotted for every star in Fig. 9, alongside the integrated optical depth to the $34\text{\,}\mathrm{M_{\astrosun}{}}$ star, which dominates the emission. The left axis shows $G_{0}$, where the solid line plots the spatial average within a radius of 2.5 cells as it varies with time; the dashed line shows the time-average value (with filled boundaries for the standard deviation). The right axis shows $\exp(-\tau_{\textrm{FUV}})$ integrated along a ray between star 1 and the star of interest, where $\tau_{\textrm{FUV}}$ is given in equation 2.3. This is the path between the central cell of each sphere, so does not take into account any spatial averaging, but is a measure of the density moving into the ray directly.

The results can be roughly split into three categories: stars which start off shielded but then become exposed with $G_{0}$ rapidly jumping up; stars which experience shielding later on such that $G_{0}$ either dips temporarily, or remains constant despite getting closer to the massive star; stars where the separation is the dominant factor in determining the variation in $G_{0}$. These categories are not exclusive. In the following sections, we describe the former two categories,

5.1 Sudden illumination

Star 7 begins $4.25\text{\,}\mathrm{p}\mathrm{c}$ away from the massive star (star 1), then moves to $6\text{\,}\mathrm{p}\mathrm{c}$ over $4\text{\,}\mathrm{Myr}$. However, its value of $G_{0}$ is lowest at the beginning, despite this being its closest approach; it takes $0.75\text{\,}\mathrm{Myr}$ to rise by a factor of $\sim 3$ to $G_{0}=200$. This occurs as the massive star disperses material around itself, allowing its radiation field to reach more stars in the cluster with less intervening opacity. The column density between the two stars reduces around the $4\text{\,}\mathrm{Myr}$ mark, which results in a corresponding increase in $G_{0}$, despite the stars getting further away.

A similar evolution occurs for other stars, such as star 4 and star 15, where the initial exposure takes $0.5\text{\,}\mathrm{Myr}$. The latter’s position is in one of the cavities blown out by the expanding H ii region (the bottom right quadrant of the plots shown in Fig. 4), which decreases in column density. This means it does not take as long for $G_{0}$ to rise compared to star 7. It remains in the cavity with little intervening material, so there are no rapid rises or dips after the maximum value of $G_{0}=90$ is reached. Stars 22 and 23 follow a similar pattern, with the latter positioned in the upper right quadrant cavity.

$G_{0}$ for star 16 rises by an order of magnitude to 150 in the first $0.45\text{\,}\mathrm{Myr}$, before slowly decreasing down to 110 at $t=$0.8\text{,}\mathrm{Myr}$$, then jumping up to 280 in just $0.12\text{\,}\mathrm{Myr}$. This is despite the star starting off at a distance of $3\text{\,}\mathrm{pc}$ and only getting further away. The delayed exposure is due to its location near the dense core in which star 1 is embedded. As that core expands and gets eroded, star 16 receives more flux; but due to the inhomogeneous distribution of gas, combined with the motion of the stars themselves, the column density increases, causing the slightly lower flux. Once this is itself dispersed, the star is more fully exposed.

This is even more the case for star 20, which is initially only $1.2\text{\,}\mathrm{pc}$ away from star 1. $G_{0}$ increases by a factor of 6 in the first $\mathrm{Myr}$, reaching a maximum value of 5700.

The evolution of star 21 is somewhat more complex, as it is located further away – starting at $7\text{\,}\mathrm{pc}$ and increasing to $13\text{\,}\mathrm{pc}$. As the H ii region expands, it pushes material in the way, so the illumination of star 21 occurs in a series of three steps up to $1.45\text{\,}\mathrm{Myr}$.

Star 29 remains embedded in a dense filament (and later globule), which shields it up until $2.5\text{\,}\mathrm{Myr}$, when the gas is photoevaporated away and the star becomes exposed. $G_{0}$ rises from 2 to 920 from $t=2.45$ to $3.15\text{\,}\mathrm{Myr}$ – an increase of almost three orders of magnitude in the space of $0.7\text{\,}\mathrm{Myr}$.

5.2 Temporary occultation

The wide separation between stars 21 and 1 (7 to $13\text{\,}\mathrm{pc}$) means there is a greater chance for material to intervene: at $1.68\text{\,}\mathrm{Myr}$, this causes a dip by a factor of 2 from $G_{0}=40$ to 20. This is short-lived, as it rises after another $0.1\text{\,}\mathrm{Myr}$.

Star 5 is surrounded by a moving clump of gas which becomes aligned such that it blocks light from star 1. A column density snapshot is shown in Fig. 11. At $1.5\text{\,}\mathrm{Myr}$, the flux decreases from 420 to 10, before being restored $0.85\text{\,}\mathrm{Myr}$ later. This happens again around $3.5\text{\,}\mathrm{Myr}$, as $G_{0}$ drops from just over 400 to 7, a factor of 60 decrease. After $0.2\text{\,}\mathrm{Myr}$, it then rapidly jumps up to 750 as the gas is blown away. This, combined with the narrower separation between the two stars, results in a final $G_{0}$ of 1000. Star 10 follows a similar pattern to star 5 as it located near the same clump of gas. The occultation at $t=$3\text{,}\mathrm{Myr}$$ can be observed in the corresponding snapshot of Fig. 8, which shows $G_{0}$ as a function of distance.

5.3 Implications

The models presented here show how $G_{0}$ evolves in time as both gas and stars move around. During the first $\sim$0.5\text{,}\mathrm{Myr} the most massive star is embedded within a core which is gradually eroded by the expanding H ii region – shells expand outwards and have holes punched through them, allowing radiation to escape into the cluster medium. This means many stars are initially shielded from the massive star’s radiation field, allowing them to evolve in relative isolation, but then become illuminated as material is dispersed. Some parcels of gas may move back in between, causing a drastic decrease in flux. In [Fig. 12](#S5.F12) we show cumulative probability distributions for some representative stars, showing how likely they are to have varying values of $G_{0}$. For example, star 24 (dotted red line) has a near-constant $G_{0}$ over its whole lifetime as it is not obstructed and stays at the same distance from the most massive star; star 26 (dot-dashed black line) is similar with the exception that the distance does vary ($\sim$2\text{\,}\mathrm{pc}) – this causes the 20 per cent variation at high $G_{0}$, Stars 5, 10, and 29 have larger ranges of $G_{0}$, with a preference to be found at the higher end, as their lows are caused by transient shielding.

In real star-forming regions, the interstellar radiation field can truncate and erode protoplanetary discs by photoevaporating them from the outside. Models by Scally & Clarke (2001) and Winter et al. (2018) show that disc destruction is dominated by photoevaporation over stellar encounters, making radiative feedback the primary external factor determining the fate of protoplanetary discs. Discs have lifetimes of the order of a few to $10\text{\,}\mathrm{Myr}$ (Haisch et al., 2001; Mamajek, 2009; Williams & Cieza, 2011), comparable to the main sequence lifetime for O stars. For example, Adams et al. (2004) found that within time-scales of $10\text{\,}\mathrm{Myr}$, discs around solar-mass stars can be severely truncated (down to radii below $15\text{\,}\mathrm{A}\mathrm{U}$) with $G_{0}=$3\text{\times}{10}^{4}\text{,}, while stellar masses below $0.5\text{\,}\mathrm{M_{\astrosun}{}}$ only require $G_{0}=$3\text{\times}{10}^{3}\text{\,}. Such fields are present in the models presented here and in Paper I. Even small values of $G_{0}$ can be significant, as shown by Facchini et al. (2016) who calculated mass-loss rates of ${10}^{-7}\text{\,}\mathrm{M_{\astrosun}{}}\text{\,}{\mathrm{yr}}^{-1}$ with $G_{0}\sim 30$, provided grain growth is taken into account. Disc truncation will by necessity remove the possibility for planets to be formed outside that radius.

For the proplyds, the derived mass-loss rates of $\sim${10}^{-7}\text{,}$$ to ${10}^{-5}\text{\,}\mathrm{M_{\astrosun}{}}\text{\,}{\mathrm{yr}}^{-1}$ imply short lifetimes of ${10}^{4}\text{\,}$ to ${10}^{5}\text{\,}\mathrm{yr}$ (Störzer & Hollenbach, 1999; Henney & O’Dell, 1999; Brandner et al., 2000). This gives rise to the so-called proplyd lifetime problem (see e.g. O’dell, 2001; Scally & Clarke, 2001; Clarke, 2007), as numerous proplyds are still observed in H ii regions. The most commonly proposed solution (assuming the mass measurements are accurate) is that the O stars are very young and have only just begun to shine on nearby objects, and the proplyds will be rapidly eroded away. However, it is unlikely that the observations of several H ii regions where proplyds are found are all capturing a special, short-lived time in their evolution and in the lifetime of the O star. A complementary explanation could lie in the stellar motions and transient shielding found in our model, both early on when the massive star is embedded, and later via occultation by clumps of gas and dust. This may prolong the proplyd lifetimes, as photoevaporation rates would not be constant. A constant $G_{0}$, as used in the aforementioned models, implies a constant separation between disc and external source (as the geometric dilution does not change) as well as constant optical depth between them. These simulations show that the latter effect is particularly important in the first Myr (when the scatter is greatest in Fig. 8, and the extinction factor $\textrm{e}^{-\tau}$ in Fig. 9 is smallest), and afterwards the distance dependence becomes the dominant factor across the ensemble (when the scatter reduces, and $\textrm{e}^{-\tau}\sim 1$). Stars which do not experience any opacity effects can still have large deviations in flux arising from stellar motions (e.g. star 26). In summary, the cluster-gas interactions have a significant influence on the flux impinging on cluster members, and therefore our simulations can provide crucial input parameters for protoplanetary disc models which aim to calculate more realistic photoevaporative mass-loss rates and, when combined with viscous evolutionary models, realistic lifetimes.

6 Summary and conclusions

We have modelled a ${10}^{4}\text{\,}\mathrm{M_{\astrosun}{}}$ cloud with the same surface density and stellar mass function as Ali et al. (2018, Paper I), where a ${10}^{3}\text{\,}\mathrm{M_{\astrosun}{}}$ cloud was investigated. The same numerical methods are present in both models, allowing a comparison to be made between the two mass regimes. In summary:

The higher mass cloud is somewhat more resistant to feedback, with 25 per cent of the initial mass remaining inside $($32.3\text{,}\mathrm{pc}$)^{3}$ after $4.3\text{\,}\mathrm{Myr}$. In the lower-mass model, almost all material is removed from the $($15.5\text{,}\mathrm{pc}$)^{3}$ grid within $1.6\text{\,}\mathrm{Myr}$. 2. 2.

Given that 75 per cent of the mass is removed, photoionization feedback is an effective mechanism for disrupting GMCs. Mass leaves the grid with a peak flux of $4.7\text{\times}{10}^{-3}\text{\,}\mathrm{M_{\astrosun}{}}\text{\,}{\mathrm{yr}}^{-1}$, comparable to the effects of supernovae in the model by Rogers & Pittard (2013). 3. 3.

The total mass in Jeans-unstable cells approaches $1500\text{\,}\mathrm{M_{\astrosun}{}}$, representing a 15 per cent core formation efficiency or potential 5 per cent star formation efficiency. This estimate agrees with observational constraints. 4. 4.

Radiation pressure has a negligible effect on bulk dispersal measures, although there is a minor difference in the total Jeans-unstable mass – there is $300\text{\,}\mathrm{M_{\astrosun}{}}$ more in the combined feedback model than the ionization-only model. This indicates radiation pressure may help to induce core formation. 5. 5.

The time-variation of $G_{0}$ impinging on other cluster members is complex. It depends on the dispersal of the core in which the massive star is embedded (which takes just over $0.5\text{\,}\mathrm{Myr}$, during which objects are shielded), and on occulations by dense ISM material. There can be drastic increases or decreases by orders of magnitude, lasting for durations of $\sim$1\text{,}\mathrm{Myr}$$. This may help to explain the proplyd lifetime problem by temporarily lowering photoevaporation rates. 6. 6.

There is scatter in the relationship between incident $G_{0}$ and distance from the most massive star due to dust extinction. After the central cluster gas is dispersed, the scatter reduces and geometric dilution dominates for most stars. B stars produce a modest FUV field which keeps them above the diluted O-star flux, while occulted stars fall below this level.

The total run time of this model was $4.3\text{\,}\mathrm{Myr}$, about a Myr before the $34\text{\,}\mathrm{M_{\astrosun}{}}$ star will explode as a supernova (SN). A quarter of the gas still remains on the grid, and is unlikely to significantly disperse by the time the SN occurs, as the grid-boundary mass flux falls off over time. However, since this material is located in small, dense filaments and clumps, surrounded by large cavities of diffuse ionized gas, it is unlikely that the SN will transfer a significant proportion of its energy into the gas, whether it occurs now or in a Myr. We aim to investigate this in a future paper.

Acknowledgements

We thank the referee for insightful and constructive comments. We thank Julian Pittard and Tim Naylor for comments which improved an early version of this paper. We also thank Tom Haworth and Matthew Bate for useful discussions. AA is funded by the University of Exeter and a Science and Technology Facilities Council (STFC) studentship. TJH is funded by STFC Consolidated Grant ST/M00127X/1. The calculations for this paper were performed on DiRAC Complexity and DiRAC Data Intensive (DIaL), both at the University of Leicester, and the DiRAC@Durham facility managed by the Institute for Computational Cosmology. These form part of the STFC DiRAC HPC Facility (www.dirac.ac.uk). Complexity was funded by BIS National E-Infrastructure capital grant ST/K000373/1 and STFC DiRAC Operations grant ST/K0003259/1. DIaL was funded by BEIS capital funding via STFC capital grants ST/K000373/1 and ST/R002363/1 and STFC DiRAC Operations grant ST/R001014/1. DiRAC@Durham was funded by BEIS capital funding via STFC capital grants ST/P002293/1, ST/R002371/1 and ST/S002502/1, Durham University and STFC operations grant ST/R000832/1. DiRAC is part of the National e-Infrastructure.

The research data supporting this publication are openly available from the University of Exeter’s institutional repository at https://ore.exeter.ac.uk/repository/.