External Photoevaporation of Protoplanetary Disks: Does Location Matter?

Most stars are observed to form in groups where the stellar density significantly exceeds that in the local solar neighborhood (Lada & Lada 2003; Bressert et al. 2010). In these young star-forming regions, stars are frequently observed to have excess flux in the infrared portion of their spectral energy distributions (Haisch et al. 2001; Richert et al. 2018), which is interpreted as being due to the presence of a protoplanetary disk. Furthermore, such protoplanetary disks have been directly imaged (O'dell & Wen 1994; ALMA Partnership et al. 2015) and can now have their masses (dust, and sometimes gas) measured (e.g., Mann et al. 2014; Ansdell et al. 2017; Eisner et al. 2018).

These protoplanetary disks rapidly evolve (Haisch et al. 2001; Richert et al. 2018), with the fraction of young stars in a given star-forming region hosting a disk dropping from around 80% at 1 Myr to around 10%–20% after 5 Myr (with the caveat that ages of young stars are notoriously difficult to determine; Bell et al. 2013; Soderblom et al. 2014). This either implies that planets are being formed rapidly, and/or a combination of processes is destroying the disks.

Young star-forming regions often host massive stars (>15 M_⊙), whose rapid internal evolution produces intense far-ultraviolet (FUV) and extreme-ultraviolet (EUV) radiation fields. Several authors (Störzer & Hollenbach 1999; Hollenbach et al. 2000; Scally & Clarke 2001; Adams et al. 2004; Fatuzzo & Adams 2008; Haworth et al. 2018b; Concha-Ramírez et al. 2019; Nicholson et al. 2019; Winter et al. 2019; Parker et al. 2021) have shown that these radiation fields can destroy or reduce the masses of protoplanetary disks by photoevaporating the gas (and, to a much lesser extent, dust; Haworth et al. 2018b). Massive stars are inherently rare, but for example in a region containing 1000 low-mass stars one or two stars >20 M_⊙ are expected from randomly sampling the stellar initial mass function (IMF). Several authors (Parker & Goodwin 2007; Maschberger & Clarke 2008) have argued that the only limit on the mass of a star that can form is the molecular cloud mass, so occasionally a low-mass star-forming region containing only tens to hundreds of stars could produce high-mass stars (and observational examples of this appear to exist).

When one or two massive stars are present, they are expected to significantly alter the properties of protoplanetary disks around stars due to photoevaporation within a sphere of influence that usually extends out to ∼0.5 pc (e.g., Scally & Clarke 2001). Many observational studies have measured the masses of protoplanetary disks as a function of projected distance from the photoionizing massive star(s). Many of these studies report a positive correlation between the mass of a protoplanetary disk and the projected distance from the massive star, with the interpretation being that the closer disks are more susceptible to destruction from the radiation fields (Guarcello et al. 2007, 2014; Fang et al. 2012; Mann et al. 2014, 2015; Ansdell et al. 2017; Eisner et al. 2018; van Terwisga et al. 2019).

However, observations are two-dimensional (2D) projections of the inherent three-dimensional (3D) distribution, and so disk-hosting stars that are in the fore- and background of the star-forming region may appear to be close to the massive star(s) in two dimensions. Furthermore, dense, young star-forming regions are often dynamically old. Both theory (e.g., McMillan et al. 2007; Girichidis et al. 2011; Bate 2012; Kuznetsova et al. 2015) and observations (e.g., Gomez et al. 1993; Larson 1995; Cartwright & Whitworth 2004; Gutermuth et al. 2009; Buckner et al. 2019) indicate that star formation results in substructured stellar distributions in star-forming regions. In these regions, the important dynamical timescale is not related to the size scale of the entire region, but rather that of the clumps of (proto)stars in the filaments and substructure (Allison et al. 2010), which are typically of order 0.1 pc (André et al. 2014) and so can be many (local) crossing times old.

In star-forming regions with these initial conditions, the influence of external processes, both from direct interactions (e.g., Parker & Quanz 2012; Zheng et al. 2015; Vincke & Pfalzner 2016) and photoevaporation (e.g., Scally & Clarke 2001; Adams et al. 2004; Parker et al. 2021) have a heightened effect on the disruption of planetary systems. A combination of violent and two-body relaxation constantly processes the population of stars in any fixed field of view, so stars that have their disks destroyed or altered may not stay close to the ionizing sources and, conversely, disk-bearing stars may move from the outskirts to the center.

In this paper, we focus on both of these issues. We first discuss the strength of the correlations reported in the observational data (Section 2) before discussing the effects of 2D projection on a synthetic star cluster with a population of stars that have been affected by photoevaporation (Section 3). We then describe N-body simulations and a post-processing disk photoevaporation analysis (Section 4), before examining the combined effects of projection issues and stellar dynamics in Section 5. We provide a discussion in Section 6 and we conclude in Section 7.

We will compare the protoplanetary disk properties of stars in simulated star-forming regions with several observed samples of disks in the Orion Nebula Cluster (ONC; Mann et al. 2014 and Eisner et al. 2018), the Orion Molecular Cloud (OMC; van Terwisga et al. 2019), NGC 2024 (Mann et al. 2015; van Terwisga et al. 2020) and the σ Orionis star-forming region (Ansdell et al. 2017).

A correlation between the mass of the protoplanetary disk and the projected distance to the most luminous ionizing source was reported in the ONC by Mann et al. (2014), and to a lesser extent by Eisner et al. (2018). A similar correlation was reported for σ Orionis by Ansdell et al. (2017), and van Terwisga et al. (2019) found that the disks closer to the center of the Trapezium cluster were less massive than those further from the center of the Trapezium by a factor of five. Finally, Mann et al. (2015) reported that there is no correlation of disk mass with distance from the ionizing star in NGC 2024, which they attribute being due to the weaker radiation field in NGC 2024 compared to that in the ONC. In these observations, the projected distances range from hundredths of a parsec to several parsecs. Estimates for the ages of these regions vary (e.g., Bell et al. 2013), but generally the ONC is thought to be between 1–2 Myr (Da Rio et al. 2010; Jeffries et al. 2011; Reggiani et al. 2011), NGC 2024 is slightly younger (0.5–1 Myr; Meyer 1996; Levine et al. 2006), and σ Orionis is thought to be slightly older (3–5 Myr; Oliveira et al. 2004).

In Figure 1 we show the the disk masses as a function of projected distance from the brightest ionizing star in each star-forming region. In the data sets of the ONC and the OMC, this corresponds to the distance from the position of θ¹ Ori C, in the σ Orionis star-forming region this corresponds to the distance from the position of the star σ Ori, and in NGC 2024 this corresponds to the distance from the position of IRS 2b. The disk masses are calculated by taking the reported dust masses from sub-millimeter continuum observations and assuming a gas-to-dust ratio of 100.

Zoom In Zoom Out Reset image size

The data from Mann et al. (2014) for the ONC are shown by the orange circles in Figure 1, and those from Eisner et al. (2018) are shown by the green stars. Data for the OMC (van Terwisga et al. 2019) are shown by the gray diamonds. Data for NGC 2024 from Mann et al. (2015) are shown by the raspberry squares, and those from van Terwisga et al. (2020) are shown by the yellow pentagons. Finally, the data for σ Orionis (Ansdell et al. 2017) are shown by the blue triangles, with the darker points indicating an unambiguous detection of gas in the disk.

In order to determine the strength of any correlation in the data sets, we use the Spearman rank coefficient test, which determines whether a correlation exists between two variables. The Spearman test returns a value ρ between −1 and +1 and an associated p-value to assess the significance of the correlation (i.e., the probability of obtaining a correlation with the rank coefficient ρ from random chance, under the null hypothesis that there is no correlation). Positive values for ρ toward unity represent a strong positive correlation, negative values of ρ represent a negative correlation, and ρ values around zero indicate no correlation.

We summarize the Spearman coefficient, ρ, for all of the data sets in Figure 1 in Table 1. Significant (p-values less than 0.01) positive correlations are present in the Mann et al. (2014) ONC sample and the van Terwisga et al. (2020) data for the extended OMC, but not in the Eisner et al. (2018) data for the ONC. The three data sets from Orion probe different physical distance ranges from the ionizing stars at the center of the ONC, and for completeness we calculate the Spearman coefficient for all three data sets. If we combine both the ONC and the OMC data (under the assumption that the photoionizing star(s) at the center of the ONC are responsible for the disk evolution across these different distance scales), the Spearman ρ coefficient is 0.54, with an associated p-value of 7 × 10⁻²⁰. In NGC 2024 the Spearman rank coefficients are negative, although not significantly so. In σ Ori there appears to be a significant positive correlation.

Interestingly, the shapes of the positively correlated distributions are consistent with the idea that the disk masses are proportional to the host star masses. For example, if we draw stellar masses from a standard Maschberger (2013) IMF of the form

$$\begin{eqnarray}&&p(m)\propto {\left(\displaystyle \frac{m}{\mu }\right)}^{-\alpha }{\left(1+{\left(\displaystyle \frac{m}{\mu }\right)}^{1-\alpha }\right)}^{-\beta },\end{eqnarray} \tag{ 1 }$$

and set the disk masses to be 10% of the host star mass, then we produce distributions such as the one shown by the black plus symbols in Figure 1. (Note that in Equation (1), μ = 0.2 M_⊙ is the average stellar mass, α = 2.3 is the Salpeter (1955) power-law exponent for higher-mass stars, and β = 1.4 describes the slope of the IMF for low-mass objects, and we sample masses between 0.1–50 M_⊙.) We draw the same number of stars as in the Mann et al. (2014) sample, i.e., 22. If all distances from the ionizing star are equally probable, the distribution of black points in Figure 1 is simply the stellar mass function placed on its side and reduced by a factor of 10. In this particular realization the Spearman rank coefficient is ρ = 0.46, with a p-value of 3.2 × 10⁻² which, although larger than the p-values calculated for the observed samples, could imply a modest correlation. In 100 repeats of this experiment, we find positive correlations 7% of the time, and negative correlations 5% of the time. We therefore have demonstrated that, before considering any photoevaporation effects, any correlation (if it exists—it is only present in three samples out of seven observed disk populations) could simply be due to low-mass stars (and hence lower-mass disks) being more common and hence more likely to be observed at closer distances to the ionizing star(s) than more massive stars (with more massive disks). However, we note that increasing the number of star–disk systems leads to fewer significant correlations (simply because increasing the number of objects causes a reduction in the number of random correlations caused by low-N sampling), suggesting that the reported trends in the observational data may also be hampered by small-number statistics.

Observations of protoplanetary disk locations in star-forming regions are usually only available in two dimensions, and so the position of a host star and its protoplanetary disk is projected from three dimensions into two on the sky. This is an important consideration when discussing the occurrence rate of protoplanetary disks and their individual masses as a function of distance from photoionizing massive stars.

In this section we determine how important projection effects may be on a synthetic star cluster where we have not modeled any previous dynamical evolution and we have assumed that the cluster is dense enough that significant photoevaporation has already occurred. We have assumed a simplified geometry, where the cluster has a uniform number density profile from the center to the outskirts (Figure 2(a)). In this panel, stars that have had their disks affected by photoevaporation are shown by the red crosses. It is immediately obvious from this plot that the 2D projection of fore- and background stars that are not within the photoevaporative zone (<0.3 pc from the massive star(s)) means that some of these objects would mistakenly be classed as being close to the photoionizing massive star(s).

Zoom In Zoom Out Reset image size

We present two scenarios for the severity of the photoevaporation. In both cases, we draw stellar masses randomly from the Maschberger (2013) IMF (Equation (1)) and we sample this distribution in the mass range 0.1–50 M_⊙ and randomly place these stars within a uniform density sphere. We assume that the photoionizing massive star is at the center of the cluster. This so-called "mass segregation" is observed in many star-forming regions that host massive stars, and can be either primordial (Bonnell & Davies 1998) or occur from subsequent dynamical evolution (Allison et al. 2010). Disk masses are set to be 10% of the stellar mass.

In the first scenario, photoevaporation has destroyed all of the disks within a 0.3 pc radius of the most massive star. (Note that we assume dynamics are unimportant for the moment.) We show the disk mass as a function of 3D distance from the photoionizing massive star in Figure 2(b), where the stars with disks that have been destroyed by photoevaporation are shown by the red points. The 2D projection places fore- and background stars close to the massive star(s), so in this projection we would see a mixture of stars with and without disks in the photoevaporative zone (Figure 2(c)).

In the second scenario, we change the original disk masses by a factor that is linearly proportional to the distance from the ionizing star. Stars beyond 0.3 pc retain their original disk masses. In Figure 2(d), which is a 3D projection, we show the original disk masses within 0.3 pc by the open circles, and draw a line to the new disk masses which are shown by the red plus symbols. The 2D projection is shown in Figure 2(e), and it is clearly impossible to distinguish disks that have been affected by photoevaporation from those that have not in this plot.

We have repeated this thought experiment for star-forming regions with different geometries (including very substructured fractal distributions, and much more centrally concentrated clusters than the one we present here) and unsurprisingly find the same result. We therefore expect projection effects to confuse our analysis of photoevaporation, both in observed star-forming regions and in simulated star-forming regions where we allow the regions to evolve dynamically but are projecting the positions of the stars into a 2D field of view.

In this section we describe the setup and implementation of our N-body simulations of the dynamical evolution of star-forming regions, as well as the post-processing disk mass loss due to photoevaporation. Our simulations are very similar to those presented in Parker et al. (2021) and for full details we refer the interested reader to that paper, as well as Parker et al. (2014b), which describes the initial conditions for the positions and velocities of the stars in detail.

4.1. Simulation Setup

4.2. Disk Properties

4.3. Photoevaporation

Massive stars emit significant power in the form of FUV luminosity, where the individual photon energies are h ν <13.6 eV, and EUV luminosity, where the individual photon energies are h ν > 13.6 eV. Several authors (e.g., Störzer & Hollenbach 1999; Haworth et al. 2018a; Hollenbach et al. 2000) have shown that FUV radiation is more destructive to protoplanetary disks than EUV radiation.

To implement disk mass loss from photoevaporation due to FUV radiation, we use the FRIED grid from Haworth et al. (2018a), which takes the disk radius, disk mass, and FUV luminosity flux received by the star,

$$\begin{eqnarray}&&{F}_{\mathrm{FUV}}=\displaystyle \frac{{L}_{\mathrm{FUV}}}{4\pi {d}^{2}},\end{eqnarray} \tag{ 2 }$$

expressed in terms of the background FUV flux in the interstellar medium, G₀ = 1.6 × 10⁻³ erg cm⁻² s⁻¹ (Habing 1968), and computes the mass loss from the disk due to the FUV radiation. In Equation (2), d is the distance from the FUV-emitting star, and L_FUV is the FUV luminosity of that star. We take L_FUV as a function of stellar mass from Figure 1 in Armitage (2000), who used the stellar atmosphere models of Buser & Kurucz (1992).

To a lesser extent, the disk also loses mass due to EUV radiation, and we subtract mass from the disk, ${\dot{M}}_{\mathrm{EUV}}$, according to the prescription in Johnstone et al. (1998):

$$\begin{eqnarray}&&{\dot{M}}_{\mathrm{EUV}}\simeq 8\times {10}^{-12}{r}_{d}^{3/2}\sqrt{\displaystyle \frac{{{\rm{\Phi }}}_{i}}{{d}^{2}}}\,{{M}_{\odot }\mathrm{yr}}^{-1}.\end{eqnarray} \tag{ 3 }$$

Here, Φ_i is the ionizing EUV photon luminosity from each massive star in units of 10⁴⁹ s⁻¹ and is dependent on the stellar mass according to the observations of Vacca et al. (1996) and Sternberg et al. (2003); the disk radius r_d is expressed in units of au and the distance to the massive star d is in pc.

We subtract mass from the disks according to the FUV-induced mass-loss rate in the FRIED grid and the EUV-induced mass-loss rate from Equation (3). Models of mass loss in disks usually assume the mass is removed from the edge of the disk (where the surface density is lowest) and we would expect the radius of the disk to decrease in this scenario. We employ a very simple way of reducing the radius by assuming the surface density of the disk at 1 au, Σ_{1 au}, from the host star remains constant during mass loss (see also Haworth et al. 2018a). If

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{1\mathrm{au}}=\displaystyle \frac{{M}_{d}}{2\pi {r}_{d}1\,\mathrm{au}},\end{eqnarray} \tag{ 4 }$$

where M_d is the disk mass, and r_d is the radius of the disk, then if the surface density at 1 au remains constant, a reduction in mass due to photoevaporation will result in the disk radius decreasing by a factor equal to the disk mass decrease, so that

$$\begin{eqnarray}&&{r}_{d}({t}_{k})=\displaystyle \frac{{M}_{d}({t}_{k})}{{M}_{d}({t}_{k-1})}{r}_{d}({t}_{k-1}),\end{eqnarray} \tag{ 5 }$$

where r_d (t_k ) and M_d (t_k ) are the disk radius and disk mass after photoevaporation in a given time interval, and r_d (t_k−1) and M_d (t_k−1) are disk radius and disk mass before photoevaporation. In the case where the disk radius decreases due to mass loss, the smaller radius acts to slow down the rate of photoevaporation. However, if the G₀ field is high enough, a disk can still be completely destroyed (i.e., lose all of its mass).

We do not include any other internal evolution of the disks, e.g., viscous evolution or inner truncation from internal photoevaporation from the host star, nor do we allow mass loss due to accretion from the disk onto the host star. In practice, both viscous evolution and accretion onto the star will act to deplete the disk (viscous spreading increases the disk radius and decreases the disk surface density, making it more susceptible to external photoevaporation). In Parker et al. (2021) we showed that viscous spreading and accretion onto the star both contribute to many more disks being destroyed, and on much more rapid timescales (within 1 Myr when both accretion and viscous spreading are enabled, compared to ∼5 Myr without).

For further details of the evolution of disks in these N-body simulations, we refer the interested reader to Parker et al. (2021), which also includes a discussion of the adopted timestep (10⁻³ Myr) in the disk evolution calculations.

In this section we first focus on simulations with high initial stellar densities (such as those suggested for the ONC). We then look at lower-density simulations (commensurate with the observed present-day density in the σ Orionis cluster), before comparing our results with distributions of disk-bearing and diskless stars in the Orion star-forming region.

5.1. Photoevaporation in High-density Simulations

5.1.1. Dynamical Evolution of the Star-forming Regions

The dynamical evolution of the star-forming regions is very similar to that described in Parker et al. (2014b, 2021) and Nicholson et al. (2019). The most massive stars are initially placed randomly in the substructure. Due to the high degree of spatial and kinematic substructure, as well as a subvirial energy ratio, the regions violently relax to form smooth, centrally concentrated star clusters after ∼1 Myr. During this process, the massive stars migrate to the cluster center via dynamical mass segregation (Allison et al. 2010; Parker et al. 2014b). In the following analysis we calculate the distances from the most massive star, but this is also usually close to the location of the cluster center (defined as the mean position of all stars).

We will first describe the results for one realization of the initial conditions, before discussing different realizations of the same simulations which, while formally statistically identical, contain different numbers of massive stars (with slightly different masses due to randomly sampling the stellar IMF), as well as (random) differences in the initial positions and velocities of the stars. Slight differences in the numbers of massive stars will lead to slightly different FUV fields, which in turn can translate into more, or fewer, disks being destroyed (see Parker et al. 2021).

5.1.2. 10% Disk Masses, 2D Projection

We start by looking at the effects of photoevaporation where the initial disk mass is set to be 10% of the host star's mass. The initial distribution (at t = 0 Myr) is shown in Figure 3(a), where we plot the disk mass as a function of the 2D distance of the host star from the most massive star. (This distance is calculated using the x and y coordinates only.)

Zoom In Zoom Out Reset image size

As the region evolves, the disks experience mass loss due to photoevaporation, so that after 4 Myr of evolution the majority of the disks are destroyed (leading to fewer black data points from the simulations as time increases). From an initial population of 1462 disks in this simulation, after 0.5 Myr there are only 301 disks, after 1 Myr there are 181 disks, and after 4 Myr there are 53 disks. While this mass loss and consequent destruction of disks may seem drastic, it is typical of many of the FRIED models presented by Haworth et al. (2018a). As an example, radiation fields of 10⁴ G₀ are typical in our very dense regions, and in the FRIED models such high G₀ values garner mass-loss rates of ∼10⁻⁶ M_⊙ yr⁻¹, meaning that in 0.1 Myr a disk could lose 0.1 M_⊙ of material. While not all mass-loss rates are so drastic, this serves to illustrate why the disk masses are so strongly depleted in dense star-forming regions (see also Parker et al. 2021).

We see very little correlation between the mass of the disks and the distance of the host star from the ionizing source at any snapshot in time. Interestingly, and as detailed in Section 2, observations suggest that, in some star-forming regions, there is a correlation between disk mass and projected distance from the massive stars. The orange circles are data from the Mann et al. (2014) ALMA study of disks in the ONC, and these authors claim a strong correlation with disk masses increasing the further they are from the massive stars (where the disk masses are inferred from the dust content of the disks using sub-millimeter continuum observations, assuming a gas-to-dust ratio of 100). Ansdell et al. (2017) report similar results from the σ Orionis star-forming region (shown by the blue triangles in Figures 3(a)–(d)). Other authors (Guarcello et al. 2007, 2010) report similar results in more distant star-forming regions, but Mann et al. (2015) find no correlation between mass of the disks and the distance of the host star from the ionizing source in NGC 2024 (raspberry squares).

We also plot histograms of the distance (again in two dimensions only) from the massive star of low-mass stars that host a disk (open histogram) and those that have had their disks destroyed (red histogram) in Figures 3(e)–(h). These panels indicate that, on average, a star closer to the ionizing massive star is less likely to host a disk, but because stars move in and out of the cluster center, diskless stars can also be located at large distances from the ionizing star(s).

In this particular simulation, the five most massive stars are 23, 18, 17, 15, and 15 M_⊙. The number of massive stars, and their exact masses, will result in a specific radiation field for a given stellar density, and this field may increase or decrease in a different realization of the simulation if the mass function is stochastically sampled. Additionally, the stars have different initial positions and velocities as set by a random number generator. The combination of these factors could mean that some simulations might display a correlation between the disk masses and distance to the main ionizing sources, whereas others would not.

In Figure 4 we show the rolling average of disk mass as a function of distance from the most massive stars for each simulation (for clarity, we only show the first 10 realizations of the 20 simulations, but we will report on all the simulations in the following text), in incremental bins containing 10 stars. Where the Spearman rank test reports a significant (p-value <0.1) positive correlation, ⁴ we show the rolling average by a colored line. Where there is a significant negative correlation (p-value <0.1), we show the rolling average by a black line. Where there is no significant correlation, we show the rolling average by a light gray line.

Zoom In Zoom Out Reset image size

First, we note that the random sampling of the IMF results in a correlation of disk mass with distance to the ionizing star in one simulation (the magenta line in Figure 4(a)) out of 20 (and also one simulation displays a negative correlation) before any photoevaporation has occurred. Then, at 0.5 Myr, 5/20 simulations display a positive correlation (and one is negative), at 1 Myr 2/20 display a positive correlation and 2/20 display a negative correlation. Finally, at 4 Myr 1/20 shows a positive correlation and 3/20 display a negative correlation.

We make five further points related to Figure 4. First, the majority of simulations (always more than 70%) do not show a correlation. Second, those that do show a correlation tend to be transient–the same simulation displays no correlation, or a negative correlation at earlier or later times. Third, if we can obtain a significant negative correlation of decreasing disk mass with distance from the ionizing source–which cannot be due to anything other than stochastic dynamics–it is not clear why one would attach any physical meaning to a significant positive correlation. Fourth, we note that the results do not depend on the exact numbers of massive stars, or their individual masses. For example in Figure 4, the simulation that displays a significant positive correlation (the red line) contains stars with masses 19, 19, 19, 12, and 11 M_⊙. One of the simulations that displays a negative correlation contains stars with masses 22, 21, 15, 12, and 11 M_⊙ (i.e., very similar), whereas another simulation that displays a negative correlation has stars with masses 34, 18, 11, 8, and 7 M_⊙ (which would likely produce more EUV radiation, but less FUV radiation). Finally, we note that even though some simulations display a positive correlation, the rolling average fluctuates so much it is difficult to discern a strong correlation compared to other simulations where the Spearman test does not return a significant correlation.

In the Appendix we show further scatter plots of the disk mass and distance to the most massive star, and indicate the masses of the five most massive stars in each. We conclude that the results do not depend on the exact numbers of massive stars drawn randomly from the IMF.

5.1.3. 10% Disk Masses, 3D Projection

In Figure 5 we show the same simulation as in Figure 3, but this time we use the full 3D information to calculate the positions of low-mass stars (including those with disks) with respect to the ionizing massive star. Blinking between these two figures highlights the problem with projecting distances in two dimensions.

Zoom In Zoom Out Reset image size

After 4 Myr of dynamical evolution and mass loss from disks due to photoevaporation, all of the disks within 0.5 pc of the cluster center have lost the gas content of their disks (shown by the absence of black points near the ionizing star in Figure 5(d)). In the 2D projection, this boundary is blurred (Figure 3), with disks appearing to be much closer to the ionizing stars in the 2D projection (though none is as close to the ionizing star as some of the observed disks in the ONC).

The histograms of destroyed disks (red) and surviving disks (black, open) show a similar trend to the 2D data (Figures 5(e)–(h)), where stars closest to the ionizing sources are less likely to have disks, but diskless stars can also be found at large distances from the massive stars.

For the remainder of our analysis, we will only consider 2D projections of the data, as this is the information observational studies are currently party to. However, we must bear in mind the caveat that some fore- and background stars will be present in our projections.

5.1.4. 1% Disk Masses, 2D Projection

When we follow the evolution of protoplanetary disks with initial masses of only 1% of that of the host star, we see that disks are rapidly destroyed in our dense star-forming regions (see Nicholson et al. 2019). In Figure 6 we show the positions of disk-hosting stars with respect to the ionizing massive stars, and very few disks survive beyond ages of 0.5–1 Myr (fewer than 10 from an initial population of ∼1460 disks, depending on the initial conditions). This suggests that if protoplanetary disks are to form gas giant planets in young star-forming regions where photoevaporation is active, they require much higher initial masses (of order 10% of the host star's mass).

Zoom In Zoom Out Reset image size

5.2. Photoevaporation in Low-density Simulations

5.2.1. Dynamical Evolution of the Star-forming Regions

Our low-density simulations take longer to relax and form a star cluster, and this process has not reached completion after the first 4 Myr of dynamical evolution. As a result the most massive star is not necessarily in the center of the star-forming region when we determine the distance of each low-mass disk-bearing star to the ionizing source. These lower-density regions still contain several (typically one to five) O-type stars, as well as up to ∼10 B-type stars, but have a much lower stellar density (∼100 stars pc⁻³). This results in a background radiation field of between 100 and 1000 G₀, which is a factor of at least 10 lower than in the high-density simulations. The masses of the five most massive stars in the simulation we use to illustrate the destruction of disks in low-density environments are 19, 19, 19, 12, and 11 M_⊙. In the Appendix we show that the results–as for the high-density regions–are not dependent on the exact realization of the IMF.

5.2.2. 10% Disk Masses, 2D Projection

When the disks are set to be 10% of the host star's mass in our low-density simulations, photoevaporation is still very effective at destroying disks (Figure 7), likely due to the still very high FUV radiation field (500 G₀ in this particular simulation). From an initial population of 1465 disks, after 0.5 Myr 1198 disks remain, reducing to 526 disks after 4 Myr. As in the higher-density simulations, we see very few instances of a positive correlation between disk mass and distance from the ionizing stars in these simulations.

Zoom In Zoom Out Reset image size

In Figure 7 we also show the histograms of the distances to the most massive stars at 0, 0.5, 1 and 4 Myr. These low-density regions take longer to dynamically relax than the higher-density regions, and as such the histograms of remaining (open) and destroyed (red) disks differ from those in the high-density regions until an age of ∼4 Myr, when the region has dynamically mixed. In these regions the destroyed disks tend to follow the same spatial distribution as the surviving disks until the region has collapsed to form a centrally concentrated cluster. At this point (after ∼4 Myr), the stars with destroyed disks tend to be closer to the ionizing stars than those with surviving disks.

We plot the average disk mass as a function of distance from the ionizing star in Figure 8. As in Figure 4, we show the rolling average of disk mass as a function of distance from the most massive stars for each simulation and, again, we only show the first 10 realizations of the 20 simulations for the sake of clarity. Because these simulations are lower density, the corresponding radiation fields will be smaller (by a similar factor to the decrease in density, i.e., a factor of ∼100; Parker et al. 2021). This means that more disks survive for longer, and so we increase the size of our incremental bins to contain 50 stars. As before, where the Spearman rank test reports a significant (p-value <0.1) positive correlation, we show the rolling average by a colored line. Where there is a significant negative correlation (p-value <0.1), we show the rolling average by a black line. Where there is no significant correlation, we show the rolling average by a light gray line.

Zoom In Zoom Out Reset image size

Our simulations are setup such that we can keep the mass distributions of stars constant, but vary the initial radii of the star-forming regions after the masses have been drawn from the IMF. For this reason, we see the same positive correlation at 0 Myr (before any photoevaporation has taken place) as in the high-density simulations (magenta line). Interestingly, this positive correlation remains throughout the simulation (see Figures 8(b)–(d)); however, some correlations that develop later in the simulation that could be attributed to distance from the ionizing star are short-lived (e.g., the orange line present at 0.5 and 1 Myr), again suggesting that this behavior is rather stochastic.

In the lower-density simulations, at 0.5 Myr, 4/20 simulations display a positive correlation and 3/20 display a negative correlation; at 1 Myr 5/20 display a positive correlation and 2/20 display a negative correlation. Finally, at 4 Myr 7/20 show a positive correlation and none displays a negative correlation. If we discount the simulation with the birth correlation (i.e., not attributable to photoevaporation) then again more than 70% of the simulations do not display a correlation of increasing disk mass with increasing distance from the ionizing star.

5.2.3. 1% Disk Masses, 2D Projection

If we reduce the initial disk masses to be only 1% of the host star (Figure 9), then after 4 Myr a higher fraction of the disks have lost mass or been destroyed due to photoevaporation, with only 11 disks from the initial population surviving to 4 Myr. In our simulations, the remaining disk distributions are inconsistent with the masses of disks, and the positions of their host stars with respect to the ionizing sources measured in Mann et al. (2014) for the ONC and Mann et al. (2015) for NGC 2024, but are consistent with the masses and positions in σ Orionis (Ansdell et al. 2017).

Zoom In Zoom Out Reset image size

Interestingly, in both our low-density simulations the disk-hosting stars are never as close to the ionizing stars as closest observed examples in both the ONC and NGC 2024 at the start of the simulation (compare the black points in Figures 7(a) and 9(a) with the leftmost colored points). This suggests that both of these star-forming regions must have been more dense in the past, and opens the promising avenue of the properties of disk-hosting stars being used as a diagnostic to predict the initial conditions of their host star-forming regions.

5.3. Comparison with Disks in Star-forming Regions

Many authors have reported a correlation between the disk mass and the projected position of the host star in relation to ionizing massive star(s) in young star-forming regions (Guarcello et al. 2007, 2010, 2014; Fang et al. 2012; Mann et al. 2014; Ansdell et al. 2017; Eisner et al. 2018; van Terwisga et al. 2019), while some authors have found no correlation (Mann et al. 2015). Our simulations span three orders of magnitude in initial stellar densities, and include disks that are both 10%, and 1% of the host star's mass.

In no simulations do we reproduce the strength of the observed trend of disk mass and distance to the ionizing star in the Mann et al. (2014) data for the ONC. In some of our simulations there is a correlation between the disk mass and the distance to the most massive star, but this appears to be coincidental and not dependent on the dynamical evolution of the star-forming region, or the strength of the FUV radiation field. A striking result that we will follow up in a future paper is that once the FUV field exceeds ∼100 G₀, then the numbers of ionizing stars present in the simulation do not strongly correlate with the number of disks that are destroyed. Therefore, the biggest factor in how many disks are destroyed is the initial stellar density of the star-forming region. Other changes to the initial conditions (different amounts of spatial and kinematic substructure, global virial ratio) have only a minimal effect on disk destruction rates (see also Nicholson et al. 2019).

One feature of our simulations is that in dense star-forming regions like the ONC, we would expect photoevaporation to be extremely efficient at depleting the gas content of protoplanetary disks. Most observational studies that measure gas and dust masses in disks select targets based on known infrared excesses (i.e., the presence of a disk). However, in our simulations we quickly create a population of diskless stars and, furthermore, the diskless stars move due to dynamical evolution and so are not always the closest objects to the ionizing sources.

Recently, Yao et al. (2018) conducted a study of young stars in nearby regions and determined whether individual objects were diskless or disk-bearing. We show histograms of their data for Orion in Figure 10, where stars are binned depending on their projected distance from θ¹ Ori C. As in our simulated data, stars that still host disks are shown by the black (open) histogram, and stars without disks are included in the red (solid) histogram. The data of Yao et al. (2018) are incomplete for the innermost regions of the ONC, but if anything there appear to be more diskless stars further from the ionizing stars than stars with disks. This is inconsistent with the results of our simulations, where, on average, stars that retain their disks are at similar distances from the ionizing stars than those that lose their disks (compare the open and filled histograms, respectively, in panels (e)–(h) of Figures 3 and 5).

Zoom In Zoom Out Reset image size

Intriguingly, however, the observations of Yao et al. (2018) are consistent with our dynamical picture in which stars can lose their disks due to external photoevaporation, and then move significant distances within the star-forming region. The big caveat here is that our simulations do not include other disk depletion mechanisms (internal photoevaporation, truncation through encounters, or rapid planet formation) that may have occurred within the observed disks. While internal photoevaporation would be unlikely to provide a correlation between disk mass and distance to ionizing stars, truncation of protoplanetary disks due to stellar encounters is efficient in very dense (>10⁴ stars pc⁻³) regions (Vincke et al. 2015; Portegies Zwart 2016; Winter et al. 2018). This could explain the low disk masses in the central regions of the ONC (although we reiterate that stars would move in and out of the central regions). On the other hand, the density of σ Ori is probably too low for disk truncation to explain the correlation in that region.

Finally, we note that our simulations assume only one epoch of star formation in a region, whereas recent work by Winter et al. (2019) assumes multiple star formation episodes, which means—if prevalent in nearby star-forming regions—that any observed disks could have been subject to differing amounts of photoevaporation.

Alternatively, if star formation is gradual then disks could have different ages and be in different stages of evolution. Longmore et al. (2014) proposed a "conveyor belt" scenario for star formation, where clusters like the ONC could be gradually fed low-mass stars via filaments; this would be observed as an age spread in the whole region with the younger stars being more central, and potentially with higher numbers of protoplanetary disks in the central regions (see also Hillenbrand 1997; Beccari et al. 2017; Winter et al. 2019).

Several authors have found a trend that the masses of protoplanetary disks in star-forming regions increase the further the disk is (projected in two dimensions) from the massive star (Guarcello et al. 2007, 2014; Fang et al. 2012; Mann et al. 2014; Ansdell et al. 2017; Eisner et al. 2018; van Terwisga et al. 2019), which is explained as disk mass loss due to photoevaporation from massive stars being more efficient the closer the disks are to these ionizing sources. Furthermore, this dependence on distance appears to be extremely strong. For example, in the ONC Mann et al. (2014) report that disk masses increase from 10⁻⁴ M_⊙ in close proximity (<0.1 pc) to θ¹ Ori C, up to values around 0.1 M_⊙ at distances of up to 1 pc. Similar trends are reported in σ Ori by Ansdell et al. (2017).

Recently, the advent of ALMA has enabled both the dust and gas components of protoplanetary disks to be probed, and the respective masses of the gas and dust content to be estimated. Mann et al. (2014) do not detect any CO in any of the disks in their sample, suggesting that the mass of the gas is very low, or non-existent. Ansdell et al. (2017) find that only six disks in σ Ori contain significant amounts of gas, and all of them are more than 1 pc away from the most massive star in projection (these disks are shown by the purple triangles in our plots in Figures 3, 5, 6, 7, and 9).

If the trend in increasing disk mass with distance from the massive star is not explained by photoevaporation of the gas content, can dust be photoevaporated in significant amounts? In the FRIED grid models we use from Haworth et al. (2018b), the fraction of the dust mass in a disk that is photoevaporated is low. Even in an extremely strong (10³–10⁴ G₀) FUV field, the dust mass only decreases by a factor of a few within 10 Myr, assuming that the amount of mass in small dust grains, which could be entrained in the evaporative wind, is low. (In contrast, the gas component of the disk is reduced by several orders of magnitude, and is usually completely evaporated in such strong FUV fields.) The observed result that disk masses increase with increasing distance from the massive stars in several star-forming regions (Guarcello et al. 2007; Mann et al. 2014; Ansdell et al. 2017) would therefore only be consistent with external photoevaporation if these disks had a significant gas component by mass.

If photoevaporation did significantly reduce the dust content of protoplanetary disks, it is not clear that we would expect a trend in disk mass increasing the further the disk is from an ionizing massive star. Dynamical evolution of star-forming regions depends on the initial stellar density of a star-forming region (Parker et al. 2014b; Parker 2014), but even low-density ( ≥10 M_⊙ pc⁻³) star-forming regions evolve on timescales similar to the photoevaporation timescales in the models of Störzer & Hollenbach (1999), Facchini et al. (2016), and Haworth et al. (2018b, 2018a); see also Nicholson et al. (2019).

Therefore, we would not expect a low-mass disk-hosting star to remain close to the ionizing massive star for a sustained length of time, and the disk-hosting star would spend a considerable amount of time in different locations. This is evidenced in our simulations, where the histogram of stars with disks, compared to those whose disks have been destroyed, appears to be random, centered on the average interstellar distance between stars. In other words, the distributions are random because both sets of stars—disk-hosting and diskless—have moved around in the star-forming region. The recent survey by Yao et al. (2018) suggests there is very little difference in the distributions of disk-hosting and diskless stars in Orion, although their data are incomplete in the immediate vicinity of θ¹ Ori C.

It is therefore unclear what is producing the observed trend of increasing disk mass with increasing distance from the massive stars in some star-forming regions. Mass segregation, where the most massive stars are preferentially found at the center, is observed in many star-forming regions, including the ONC. This could result in a radially dependent spatial distribution of stellar masses, which could in turn result in a radial dependence of the disk masses. However, this would predict the opposite effect to what is observed, namely that the more massive stars (and hence disks) would be found closer to the center of the star-forming region.

Finally, we reiterate that the timescales for photoevaporating the gas content of disks are extremely short (see also Concha-Ramírez et al. 2019; Nicholson et al. 2019), depleting or destroying altogether the material from which gas giant planets form. This argues for either extremely rapid (≪1 Myr) planet formation, or that giant planets exclusively form in star-forming regions like Taurus (e.g., Güdel et al. 2007; Luhman et al. 2009), which contain no massive (O- and B-type) stars (unlike the ONC). The latter scenario is in tension with mounting evidence that the Sun experienced enrichment from short-lived radioisotopes in its natal star-forming region (e.g., Gounelle & Meynet 2012; Parker et al. 2014a; Young 2016; Lichtenberg et al. 2016, 2019; Boss 2017; Portegies Zwart 2019), which requires the presence of one or more massive stars (Adams 2010; Lugaro et al. 2018). Alternatively, giant planet formation may need to exclusively occur within several astronomical units of the host star as Nicholson et al. (2019) demonstrate that disks with small radii (≤10 au) can survive in very strong radiation fields for several megayears.

We present N-body simulations of the dynamical evolution of star-forming regions in which protoplanetary disks lose mass due to external FUV and EUV radiation from the most massive stars. We have compared the positions of disk-hosting stars and the disk masses to recent observations of disks in nearby star-forming regions. In particular, we have calculated the disk mass loss due to photoevaporation as a function of projected distance from the ionizing massive star(s). Our conclusions are as follows.

(i) In agreement with other authors (e.g., Scally & Clarke 2001; Adams et al. 2004; Winter et al. 2018; Concha-Ramírez et al. 2019; Nicholson et al. 2019), we find that external FUV and EUV radiation has a highly detrimental effect on protoplanetary disks; in many simulations very few disks with initial radii of 100 au remain after only a few megayears.

(ii) In the observational data, the disks in the ONC sample published by Mann et al. (2014), those in the more extended OMC published by van Terwisga et al. (2020), and those in σ Orionis published by Ansdell et al. (2017) display a significant correlation of increasing disk mass with increasing distance from the ionizing massive star. We note that if the initial disk dust masses are proportional to the host stars' masses, then such a trend may simply be due to the imprint on the disk masses of randomly sampling the stellar IMF.

(iii) Projection effects add significant noise to the spatial distribution of disk-hosting stars. Stars with disks that appear to be close to the ionizing sources may be fore- or background members of the star-forming region.

(iv) After accounting for projection effects, we find no evidence that the observed trend of increasing disk mass with increasing distance from ionizing stars is due to external photoevaporation, unless the adopted photoevaporation models are very wrong. The reason for this is twofold. First, the gas content is uniformly destroyed on rapid timescales in the simulations. Second, the observed disks contain little or no gas, whereas even a strong FUV field will not result in significant mass loss from the dust component of the disk.

(v) Even if the photoevaporation models are reasonably accurate, in a dense star-forming region there is no reason to expect a dependence of the disk properties on the distance to the ionizing stars. Our simulations suggest that stars lose their disks, and then move around the star-forming region so that they may reside in far-flung locations from the ionizing stars.

A first step to resolving some of the issues we have identified would be to improve on existing observations of disks in star-forming regions, and to obtain observations of disk populations in other star-forming region to improve the statistical sample and any comparisons between regions.

R.J.P. and O.P. acknowledge support from the Royal Society in the form of Dorothy Hodgkin Fellowships. H.L.A. was supported by the 2018 Sheffield university Undergraduate Research Experience (SURE) scheme. We thank Yuhan Yao for providing their observational data ahead of final publication. We thank the anonymous referee for a prompt and helpful report.

In this section of the Appendix we show the non-dependence of stochastic sampling on the IMF on our results. In each realization of a star-forming region, we draw masses randomly from the IMF in the mass range 0.1–50 M_⊙. Due to the modest numbers of stars in our simulations (N = 1500), this means that the high-mass end of the IMF is not fully sampled, which in turn results in different numbers of massive stars (>17.5 M_⊙) in each simulation.

The different numbers of massive stars result in different strength radiation fields, and as we also change the initial positions and velocities of stars in the different realizations of the simulations, the effect of changing the radiation field is difficult to predict analytically. For this reason, we present the results of four additional simulations of star-forming regions with high initial stellar densities (10⁴ stars pc⁻³) and these are shown in Figure 11. For brevity, we only show the results at 1 Myr, and we indicate the masses of the five most massive stars in the panel captions. These plots should be directly compared with Figures 3(c) and (g), which show the results for a simulation where the most massive stars are 23, 18, 17, 15, and 15 M_⊙.

Zoom In Zoom Out Reset image size

The different numbers (and masses) of these ionizing stars do change the number of surviving disks after 1 Myr, but our conclusions are unchanged. There is still no correlation between the mass of the disk and the distance to the most massive ionizing star, and the process of dynamical mass segregation means that all of the ionizing stars are in central locations. Furthermore, the stars that lose their disks move around the star-forming region following disk destruction, and so the locations of stars with and without disks present a wide distribution in distance from the most massive star.

In Figure 12 we present the results of four additional simulations for our low-density simulations. Here, we compare different simulations after 4 Myr of dynamical evolution and disk photoevaporation, and so these plots should be directly compared with Figures 7(d) and (h), whose five most massive stars were 19, 19, 19, 12, and 11 M_⊙. As with the higher-density simulations, the numbers of surviving disks differs depending on the mass function (and stochastic differences in the evolution of star-forming regions), but there is no dependence of the mass of a disk on the distance of its host star to the most massive star(s).

Zoom In Zoom Out Reset image size