External photoevaporation of protoplanetary discs in Cygnus OB2: linking discs to star formation dynamical history

TL;DR

This study links the properties of protoplanetary discs in Cygnus OB2 to its star formation history by modeling FUV-driven disc erosion, revealing the timing of gas expulsion and the impact of external photoevaporation.

Contribution

It introduces a novel method to constrain star-forming environment history using disc properties and models FUV-induced mass loss to infer gas expulsion timing and initial conditions.

Findings

Gas was expelled approximately 0.5 Myr ago under certain assumptions.

Disc mass scales steeply with stellar mass, M_disc ∝ m_star^β, with β ≥ 2.7.

Cygnus OB2 likely formed from massive clumps of about 10^4 solar masses.

Abstract

Many stars form in regions of enhanced stellar density, wherein the influence of stellar neighbours can have a strong influence on a protoplanetary disc (PPD) population. In particular, far ultraviolet (FUV) flux from massive stars drives thermal winds from the outer edge of PPDs, accelerating disc destruction. In this work, we present a novel technique for constraining the dynamical history of a star forming environment using PPD properties in a strongly FUV irradiated environment. Applying recent models for FUV induced mass loss rates to the PPD population of Cygnus OB2, we constrain how long ago primordial gas was expelled from the region; $ 0.5$ Myr ago if the Shakura & Sunyaev $\alpha$-viscosity parameter is $\alpha = 10^{-2}$ (corresponding to a viscous timescale of $\tau_\mathrm{visc} \approx 0.5$ Myr for a disc of scale radius $40$ au around a $1\, M_\odot$ star). This value of $\alpha$ is effectively an upper limit, since it assumes efficient extinction of FUV photons throughout the embedded phase. With this gas expulsion timescale we are able to produce a full dynamical model that fits kinematic and morphological data as well as disc fractions. We suggest Cygnus OB2 was originally composed of distinct massive clumps or filaments, each with a stellar mass $\sim 10^4 \, M_\odot$. Finally we predict that in regions of efficient FUV induced mass loss, disc mass $M_\mathrm{disc}$ as a function of stellar host mass $m_\mathrm{star}$ follows a power law with $M_\mathrm{disc} \propto m_\mathrm{star}^\beta$, where $\beta \gtrsim 2.7$ (depending on disc initial conditions and FUV exposure). This is steeper than observed correlations in regions of moderate FUV flux ($1 < \beta <1.9$), and offers a promising diagnostic to establish the influence of external photoevaporation in a given region.

Figures

Tables

Table 1. Table 1: Parameters of the ‘best-fit’ model, used to reproduce the properties of the observed stellar population of Cyg OB2.

Type	$M_{\mathrm{stars}}$ ($M_{\odot }$)	$a_{\mathrm{stars}}$ (pc)	$\gamma$	$P_0$	$D_0$	$M_{\mathrm{gas},0}$ ($M_{\odot }$)	$a_{\mathrm{gas}}$ (pc)	${\tau }_{\mathrm{delay}}$ (Myr)	${\tau }_{\exp }$ (Myr)
FILA	$8\cdot {10}^4$	$7$	$5.8$	$1$	$2.5$	$7.9\cdot {10}^6$	$10$	$1$	$1.5$

Equations

$$\rho ={\rho }_0{\left(1+\frac{r^2}{a_{\mathrm{stars}}^2}\right)}^{-\frac{\gamma +1}{2}}$$

$${\rho }_{\mathrm{gas}}=\frac{3M_{\mathrm{gas}}}{4\pi a_{\mathrm{gas}}}{\left(1+\frac{r^2}{a_{\mathrm{gas}}^2}\right)}^{-5/2}$$

$$\varphi (r)=-\frac{G{M}_{\mathrm{gas}}}{a_{\mathrm{gas}}}{\left(1+\frac{r^2}{a_{\mathrm{gas}}^2}\right)}^{-1/2}.$$

$$U_{\mathrm{tot}}=\sum _i^N{m}_i\left(\varphi (r_i)+\sum _{j\ne i}^N\frac{G{m}_j}{2r_{ij}}\right).$$

$$Q_{\mathrm{vir},0}\equiv \frac{{\sum }_i^N{m}_i{v}_i^2}{2U_{\mathrm{tot}}}=0.5$$

$$f(v)\propto v^2\exp (-v^2)$$

$${\mathbf{v}}_{g+1}={\mathbf{v}}_g+\delta {\mathbf{v}}_{g+1}$$

$$f(\delta v)\propto \delta v^2\exp \left(-2^g\delta v^2\right)$$

$$\frac{N_{\mathrm{c},j}}{4\pi r_{\mathrm{c},j}^2\Delta r_{\mathrm{c},j}}=\frac{\rho (r_{\mathrm{c},j})}{⟨m_{\mathrm{star}}⟩}$$

$$\xi (m)\propto \{\begin{array}{l}{m}^{-1.3}\mathrm{for}0.08M_{\odot }\le m<0.5M_{\odot }\\ m^{-2.3}\mathrm{for}0.5M_{\odot }\le m<1.0M_{\odot }\\ m^{-2.4}\mathrm{for}1.0M_{\odot }\le m<100M_{\odot }\\ 0\mathrm{else}\end{array}$$

$${\dot{M}}_{\mathrm{acc},0}=\frac{3\alpha M_{\mathrm{disc},0}{H}_1^2{\Omega }_1}{2R_1^2}$$

$${\Sigma }_0(r<R_{\mathrm{disc},0})\approx \frac{M_{\mathrm{disc},0}}{2\pi R_1^2}{\left(\frac{r}{R_1}\right)}^{-1}{e}^{-r/R_1}$$

$$R_1=40{\left(\frac{m_{\mathrm{star}}}{1M_{\odot }}\right)}^{1/2}\mathrm{au}.$$

$${\tau }_{\mathrm{visc}}\approx \frac{5.4\cdot 10^3}{\alpha }{\left(\frac{R_1}{40\mathrm{au}}\right)}^{3/2}{\left(\frac{m_{\mathrm{star}}}{1M_{\odot }}\right)}^{-1/2}\mathrm{yr}$$

$$\mathcal{E}=\frac{T^{+}}{T^{-}+T^{+}}$$

$$M_{\mathrm{disc},0}/M_{\odot }=e^{\mu +\sigma X}$$

Full text

††pubyear: 2018††pagerange: External photoevaporation of protoplanetary discs in Cygnus OB2: linking discs to star formation dynamical history–External photoevaporation of protoplanetary discs in Cygnus OB2: linking discs to star formation dynamical history

1 Introduction

The majority of stars form and spend their early stages of evolution in regions of enhanced stellar density (Lada & Lada, 2003). This might be as part of a gravitationally bound cluster, or if star formation efficiency is sufficiently low this environment can become a short-lived association (Lada & Lada, 2003; Murray, 2011; Pfalzner & Kaczmarek, 2013). This latter case is expected to be common (although see Kruijssen, 2012), therefore understanding the effect of clustered environments on the early stages of stellar evolution is important even for apparently isolated stars.

In particular, feedback from the stellar environment has been shown to have a significant impact on planet formation in young clusters and associations. Protoplanetary disc (PPD) lifetimes are thought to be $\sim 3-10$ Myr (e.g Haisch et al., 2001; Williams & Cieza, 2011; Ribas et al., 2015), during which period a star typically remains in its formation environment. In this case, PPD destruction can be induced due to close star-disc encounters (Clarke & Pringle, 1993; Ostriker, 1994; Pfalzner et al., 2005; Breslau et al., 2014; Winter et al., 2018a; Vincke & Pfalzner, 2018; Cuello et al., 2018; Concha-Ramírez et al., 2019) or far-ultraviolet (FUV) driven external photoevaporation (and, to a lesser extent, extreme ultraviolet - EUV, Johnstone et al., 1998; Störzer & Hollenbach, 1999; Armitage, 2000; Adams, 2010; Facchini et al., 2016; Haworth et al., 2018a; Winter et al., 2018b). Recent studies suggest that, statistically speaking, the dominant disc dispersal mechanism is external photoevaporation (Scally & Clarke, 2001; Winter et al., 2018b), while encounters can set PPD initial conditions during the evolution of stellar multiple systems (Winter et al., 2018c; Bate, 2018).

FUV radiation has a significant influence on disc evolution in a wide range of environments. Mass loss from the outer edge is driven when the thermal energy exceeds the gravitational potential (e.g. Johnstone et al., 1998). The most obvious examples of such a process are proplyds, visible for example in the Orion Nebula Cluster (O’Dell & Wen, 1994), where the finite penetration of ionising radiation into a disc’s neutral wind creates an offset (cometary) ionisation front. The strong neutral winds in the Orion Nebula Cluster are driven by FUV flux of $F_{\mathrm{FUV}}\sim 3\times 10^{4}\,G_{0}$ 111$G_{0}$ is the Habing unit (Habing, 1968). It is a measure of the FUV field as a multiple of the solar neighbourhood value: $1.6\times 10^{-3}$ erg cm*-2* s*-1* .. Since then proplyds have been found for discs which experience a factor $\gtrsim 10$ lower flux (Kim et al., 2016). Even when a proplyd is not visible, a disc can exhibit significant FUV induced mass loss. Haworth et al. (2017) find that in the very extended disc around IM Lup, with $F_{\mathrm{FUV}}\sim 4\,G_{0}$, photoevaporation drives substantial mass loss. A large fraction of stars form in regions where FUV flux is considerably greater than this (Fatuzzo & Adams, 2008; Winter et al., 2018b).

The properties of young massive stellar clusters/associations and the giant molecular clouds (GMCs) from which they form are diverse, and the link between them is not well characterised (see Longmore et al., 2014, for a review). During formation the early cluster may undergo cold collapse (Tobin et al., 2009; Kuznetsova et al., 2018), or after the expulsion of gas the stellar population may become supervirial (Goodwin, 2009; Pfalzner & Kaczmarek, 2013), dependent on the density and velocity dispersion of the primordial GMC. This in turn influences the evolution of mass segregation (e.g. Bonnell & Davies, 1998) and substructure (e.g Goodwin & Whitworth, 2004) within the cluster. Because the environment of a star has an influence on the associated PPD, the dynamical history of a cluster is closely linked with the properties of a disc population.

While many authors have attempted to account for the local environment in considering PPD evolution (e.g. Scally & Clarke, 2001; Cleeves et al., 2016; Guarcello et al., 2016), none have inverted this method and used observed disc populations to put constraints on the dynamical history of a star forming environment. This work is partly motivated by this goal.

Cygnus OB2 (Cyg OB2) is a young massive OB association in the Cygnus X region, and has been used as an empirical test of feedback mechanisms on PPD evolution. It contains many massive stars up to $\sim 100\,M_{\odot}$ (e.g. Massey & Thompson, 1991; Wright et al., 2015) which contribute to strong FUV radiation fields. Guarcello et al. (2016) analysed the disc fraction within Cygnus OB2 as a function of FUV flux, and found that surviving discs were less common at small projected separations from massive stars. Other authors, such as Wright et al. (2016), have made observations which indicate a complex dynamical substructure within the association. Collating this evidence, we here aim to apply $N$-body simulations and those combining viscous disc evolution and photoevaporation to replicate observations of Cyg OB2. We will reproduce the present day stellar kinematics and a dynamical history consistent with the observed disc fraction distribution. In this way we can shed light on both the history and the likely future of the PPD population and the stellar components.

In the remainder of this work we first review the observational constraints on the properties of Cyg OB2 in Section 2. We describe our numerical method and models in Section 3. In Section 4 we compare our models with the observational data. We draw conclusions in Section 5.

2 Properties of Cygnus OB2

2.1 Stellar population

Cyg OB2 is a young association at a distance $\sim 1.33$ kpc from the Sun (Kiminki et al., 2015). The majority of members formed $3-5$ Myr ago (Wright et al., 2010), although some stars have ages as young as $\sim 2$ Myr (Hanson, 2003) and as old as $\sim 7$ Myr (Drew et al., 2008). Estimates of the total stellar mass of Cyg OB2 have varied between $2$–$10\times 10^{4}\,M_{\odot}$ (Knödlseder, 2000; Hanson, 2003; Drew et al., 2008; Wright et al., 2010), although Wright et al. (2015) find a slightly lower mass of $\sim 1.6\times 10^{4}\,M_{\odot}$ within a radius of $13$ pc of the apparent centre. This includes a population of $\sim 169$ OB stars, the most massive of which is $\sim 100M_{\odot}$ with an age of $\sim 2$ Myr. Cyg OB2 does not exhibit evidence of mass segregation, whereby the most massive stars occupy regions with the greatest gravitational potential (Wright et al., 2014b).

The initial mass function (IMF) at the high mass end is not agreed upon in the literature, with many authors arriving at different conclusions (Massey & Thompson, 1991; Massey et al., 1995; Knödlseder, 2000; Kiminki et al., 2007; Wright et al., 2010; Comerón & Pasquali, 2012). Wright et al. (2015) take into account massive stars which have evolved to their end state. In this way they find that the observed stellar masses follow an IMF $\zeta(m)\propto m^{-2.39\pm 0.19}$ at high masses, which is approximately consistent with the ‘universal’ IMF of Kroupa (2001) (or indeed a Salpeter, 1955, IMF). However, inferring the high mass IMF is problematic since the occurance of supernovae in Cyg OB2 remains a point of debate (see Butt et al. 2003 and discussion in Wright et al. 2015).

Wright et al. (2016) describe the spatial density distribution in Cyg OB2 with an Elson, Fall and Freeman profile (EFF profile hereafter - Elson et al., 1987):

[TABLE]

where $a_{\mathrm{stars}}=7.5$ pc and $\gamma=5.8$. Normalising for the total mass of the cluster this makes the central mass density $\rho_{0}\approx 22\,M_{\odot}$ pc*-3*. However, this profile was derived for a small central region of $\sim 8$ pc$\times 8$ pc, which does not enclose the estimated core radius. Hence, when we consider the distribution of the stellar mass in our models (see Section 4.3.2) we will focus on reproducing the mass enclosed within a projected radius of $13$ pc (Wright et al., 2015).

2.2 Velocity dispersion

Wright et al. (2016) presented a high-precision proper motion study of the X-ray sources in Cyg OB2. They found that the region is gravitationally unbound and exhibits an anisotropic velocity dispersion with proper motion components $\sigma_{\alpha}=13.0^{+0.8}_{-0.7}$ km s*-1* and $\sigma_{\delta}=9.1^{+0.5}_{-0.5}$ km s*-1*. The radial velocity dispersion has also been measured to be $\sigma_{r}\sim 10$ km s*-1*, although systematic overestimates due to the binary fraction introduce uncertainties (Kiminki et al., 2007, 2008).

Interestingly Wright et al. (2016) found little evidence for expansion (or contraction) in the velocity field when considering the large scale variations of the proper motion distribution. This finding is independent of the definition of the centre of the association, and was argued by dividing proper motions into radial and azimuthal components. The ratio of kinetic energy between the radial and azimuthal directions was found to be approximately $60:40$. In the radial (projected) velocities, no bias is found towards or away from the centre. Azimuthally there is some preference for a ratio $66:34$ in favour of kinetic energy in the direction of negative PA. Because Cyg OB2 is not bound, this is interpreted by Wright et al. (2016) as a remnant of the angular momentum of the primordial GMC.

2.3 Substructure

Interpretation of the internal substructure in Cyg OB2 is not straightforward. Knödlseder (2000) concluded that it is a spherically symmetric region with a diameter of $\sim 2^{\circ}$ ($46$ pc). Since then, a number of authors have suggested a more complex morphology:

• •

Bica et al. (2003) suggested that Cyg OB2 is home to two open clusters which can be seen towards the centre. However, Guarcello et al. (2013) note that the two apparent clusters are divided by a bright nebula. This makes it unclear if the two are physically separate or merely appear so due to the higher extinction in the intervening region.

• •

Wright et al. (2010) found evidence of populations within Cyg OB2 with distinct ages, $3.5^{+0.8}_{-1.0}$ Myr and $5.25^{+1.5}_{-1.0}$ Myr for central and northwestern regions respectively. Ostensibly, this suggests multiple star forming events. However, the ages in physically separated regions exhibited a wide spread such that they are almost consistent with being coeval. Further the authors acknowledge a number of sources of uncertainty, including variability in pre-MS stars (Herbst et al., 1994); binarity (Preibisch & Zinnecker, 1999); variable accretion (Baraffe et al., 2009); or non-uniform extinction (Guarcello et al., 2012).

• •

A large number of A stars were identified south of the apparent centre by Drew et al. (2008). This population appears distinct spatially and non-coeval with the OB population (although the estimated ages $\sim 5$–$7$ Myr are nearly consistent with coevality within uncertainties). As the authors note, it is also possible that these stars are actually behind Cyg OB2 along the line-of-sight and wrongly associated due to projection (Schneider et al., 2006).

• •

Guarcello et al. (2013) used a critical side length criterion in the mimimum spanning tree of the disc-bearing population to suggest that Cyg OB2 has a clumpy substructure. However, as stated by the authors, the definitions of these clumps are dependent on the definition of the critical side length. Additionally the non-uniform extinction due to foreground gas complicates this argument as in the case of the two open clusters identified by Bica et al. (2003).

• •

Perhaps the best evidence for underlying structure in Cyg OB2 has been the proper motion study of Wright et al. (2016). On small scales they found evidence for kinematic substructure, which is the correlation of proper motion vectors with position. Applying a Moran’s I statistic (Moran, 1950) they found correlation with a significance of $9.7\sigma$ and $12.5\sigma$ in RA and Dec velocity components respectively.

2.4 PPD population

Guarcello et al. (2016) studied the correlation between the fraction of surviving PPDs as a function of the local FUV and EUV intensity. They use a sample of 7924 X-ray sources (Wright et al., 2014a), for which the presence (absence) of a PPD is inferred by the presence (absence) of an infrared excess in the photometric data compiled from numerous surveys (see Guarcello et al., 2013). Subsequently they estimate the local flux as a function of projected separation from each O star (see Guarcello et al., 2007). The disc fractions as a function of $F_{\mathrm{FUV}}$ are divided into 6 bins between $\sim 10^{3}\,G_{0}$ and $\sim 5\times 10^{4}\,G_{0}$. Over this space the disc fraction drops monotonically from $\sim 40\%$ to $18\%$ with FUV intensity (see figure 3 in Guarcello et al., 2016). These observations are discussed in the context of modelling in Section 4.2.

2.5 Observational summary & modelling challenges

Cyg OB2 is a well studied young association, and as such a large number of physical characteristics serve as constraints and measures for the success of any modelling attempts. Some such metrics are as follows:

• •

Wright et al. (2016) find that the velocity dispersion in Cyg OB2 is anisotropic. This suggests that the stars share the systematic large-scale velocity field of the primordial GMC. We will find that such observations cannot be reproduced by simple models without underlying substructure, and we consider the initial properties required in Section 4.3.1.

• •

Cyg OB2 presently has a central mass of $\sim 1.6\times 10^{4}\,M_{\odot}$ within $13$ pc of the centre. Any dynamical model should match this central density after evolving for the age of the association ($\sim 3-5$ Myr). The initial stellar mass required to maintain this central density is explored in Section 4.3.2.

• •

Although Cyg OB2 is gravitationally unbound with a large velocity dispersion, there is no bias (inwards or outwards) in the radial component of kinetic energy. This apparently suggests a lack of recent rapid expansion, despite the high velocities. When we consider our final $N$-body model, we will explore what expansion metric we would ‘observe’ (Section 4.3.4).

• •

Statistical measures of the proper motion distribution suggest kinematic substructure which is probably indicative of stars travelling together as small virialised groups (Wright et al., 2016). We consider evidence for kinematic substructure for our final $N$-body model in Section 4.3.4.

• •

The disc fraction as a function of (projected) FUV field strength (Guarcello et al., 2016) provides a constraint on the length of time for which external photoevaporation has been an efficient mechanism in Cyg OB2. This allows us to put constraints on the gas expulsion timescale. We perform these calculations to constrain the appropriate $N$-body model in Section 4.2, then revisit the PPD properties obtained from our final model in Section 4.4.

3 Numerical Method

The goal of our models is to reproduce the observed dynamical properties of the stellar population and the observed fractions of surviving discs. The latter is achieved by tracking the FUV flux experienced by PPDs evolving within a given $N$-body model. In this section we discuss the modelling procedure applied to both the stellar dynamics (Section 3.1) and the disc evolution (Section 3.2).

3.1 Kinematic modelling

The dynamical evolution of the stellar population is calculated using Nbody6++GPU (Wang et al., 2015). This is an MPI/GPU accelerated version of Nbody6 (Spurzem, 1999; Aarseth, 2003), and has built-in routines to deal with the evolution of a stellar cluster within a (gas) potential. As the stellar components of Cyg OB2 are presently highly supervirial, the latter feature is necessary for an initially virialised state.

3.1.1 Gas potential

The complex nature of Cyg OB2 means that some simplifying assumptions are required for dynamical modelling. We first assume that the cluster was initially in virial equilibrium due to the contribution to the gravitational potential of the gas in the primordial association. This is achieved by invoking a Plummer potential (for numerical convenience) corresponding to a gas density profile:

[TABLE]

where $M_{\mathrm{gas}}$ is the total gas mass, $r$ is the radial distance from the centre of the association, and $a_{\mathrm{gas}}$ is the scale parameter. During initial tests, we have varied $a_{\mathrm{gas}}$ (and the stellar scale parameter $a_{\mathrm{stars}}$), and although it has a mild effect on the kinematic properties of our models, the exact value (within order unity) is not crucially important. This is particularly true in the stochastically defined initial conditions of substructured models (see Section 3.1.2 below), where kinematics are more dependent on the specific realisation. Physically we expect $a_{\mathrm{gas}}\gtrsim a_{\mathrm{star}}$, and we fix $a_{\mathrm{star}}=7$ pc (as there is no clear evidence of past expansion) and $a_{\mathrm{gas}}=10$ pc.

For a Plummer density profile the specific gravitational potential is:

[TABLE]

The total potential in a given cluster make up of $N$ stars of mass $m_{i}$ initially at distance $r_{i}$ from the centre and $r_{ij}$ from a star of mass $m_{j}$ is:

[TABLE]

Then using equations 2 and 3 we require that initially

[TABLE]

where $v_{i}$ is the magnitude of the initial velocity of the $i^{\mathrm{th}}$ star. The total gas mass is chosen to maintain initial virial equilibrium for a cluster with a given velocity dispersion.

Our prescription of gas removal introduces an expulsion timescale $\tau_{\mathrm{exp}}$ over which time the potential is removed. The gas mass is reduced linearly such that $\dot{M}_{\mathrm{gas}}=M_{\mathrm{gas,}0}/\tau_{\mathrm{exp}}$. We vary $\tau_{\mathrm{exp}}$ to investigate how this affects the disc population due to extinction (see Section 4.2). We fix the time at which gas expulsion is initiated to be $\tau_{\mathrm{delay}}=1$ Myr, consistent with the age of the most massive stars in Cyg OB2 ($\sim 2$ Myr) if the cluster age is $3$ Myr (the period for which we evolve the whole system). We further define $\tau_{\mathrm{gas}}\equiv\tau_{\mathrm{delay}}+\tau_{\mathrm{exp}}$. For a discussion of the influence of gas expulsion on the dynamical state of a young cluster, see Baumgardt & Kroupa (2007).

Clearly, while computationally necessary here, a Plummer potential is not a realistic reflection of the physical conditions of the primordial gas distribution during the embedded phase. Initially, gas density distributions would be physically expected to trace the stellar density (since the stars form from the gas). Subsequently we would expect gas expulsion to occur as expanding bubbles from the most massive stars in the region (e.g. Dale & Bonnell, 2011; Dale et al., 2014; Ali et al., 2018). In terms of the cluster evolution, we therefore expect intra-clump or -filament potential to reduce faster than inter-clump potential. The influence of a geometrically complex potential on the evolution of the stellar population is certainly of interest for accurately reproducing observed kinematics and spatial distributions. However, we assume that the Plummer potential imposed here is sufficient for reproducing the global distribution of stellar positions and velocities. This is justified since we are primarily interested in preventing the rapid escape of the high velocity stars from the central regions of Cyg OB2. This is achieved by our spherically symmetric potential.

3.1.2 Stellar initial conditions

We define three different types of stellar initial condition, which we label uniform (UNIF), fractal (FRAC), and filamentary (FILA). A uniform cluster exhibits no underlying substructure, while fractal and filamentary clusters have enhanced local number densities, and stellar positions are correlated with velocities. While the surface density of FRAC model is ‘clumpy’, and individual clumps can be spatially isolated, a FILA model is defined such that the (radially binned) surface density follows an EFF profile. A filamentary model captures the morphology seen in both observations and simulations of star forming regions (e.g. Bonnell et al., 2008; Molinari et al., 2010; André et al., 2014), however we include both FRAC and FILA models in this work to compare their properties by the metrics of interest. We discuss the generation of each set of initial conditions below.

The simplest initial conditions are UNIF models for which stellar positions are simply drawn from an EFF profile (equation 1, Elson et al., 1987). We fix $\gamma=5.8$, consistent with the present day value, and fix $a_{\mathrm{stars}}=7$ pc. The magnitude of the velocities are then chosen from a Maxwell-Boltzmann distribution

[TABLE]

with a random direction, then normalised as discussed in Section 3.1.1. We note that drawing from this distribution (particularly at the high velocity dispersions we require - see Section 4.3.1) will produce a significant number of unbound stars. It is possible that this distribution is truncated at large velocities, and indeed star formation models suggest stars might have a smaller velocity dispersion than the priomordial gas (Offner et al., 2009, see Section 4.3.3).

A FRAC model is generated using the recipe fully described by Goodwin & Whitworth (2004) which we briefly review here (see also Scally & Clarke, 2002; Craig & Krumholz, 2013). First we define a cube with side length $2$ (in arbitrary units) centred at the origin. We then divide it up into $(2P_{0})^{3}$ sub-cubes, where $P_{0}$ is an integer (initially chosen to be unity) which dictates the number of largest scale sub-clusters. The centre of each of these represent the potential positions for the first generation of stars, $g=1$. All positions are offset by a vector with magnitude uniformly drawn between [math] and $2^{-(g+1)}/P_{0}$, and random (isotropically drawn) direction. ‘Parent’ positions $\bm{r}_{g}$ have $8$ possible sites for ‘child’ positions $\bm{r}_{g+1}$ which are placed with a probability $2^{D_{0}-3}$, where $D_{0}\leq 3$ is the fractal dimension (which we fix at $2.5$). Only existing children can parent future generations. We repeat this process until the number of positions greatly exceeds the number of stars in the model, at which point the members are randomly allocated to positions. The side length of the original cube is then redefined such that the initial half-mass radius matches that of an equivalent UNIF model with parameters $a_{\mathrm{stars}}=7$ pc, $\gamma=5.8$.

The velocities for each generation of stars $g=1$ is chosen in the same way as for the UNIF model. Velocities for subsequent generations $\bm{v}_{g+1}$ are correlated to the velocity of the parent star $\bm{v}_{g}$:

[TABLE]

where $\delta\bm{v}_{g}$ is a velocity with a random direction and magnitude $\delta v_{g}$. The latter is drawn from the modified Maxwell-Boltzmann distribution

[TABLE]

such that child stars have velocities correlated to their parents.

A FILA model is a hybrid between FRAC and UNIF. To construct the initial conditions we first produce a FRAC model and then force the stellar positions into an EFF profile. This is achieved by dividing the radial positions into bins, with the $j^{\mathrm{th}}$ bin containing $N_{\mathrm{c},j}$ stars approximately at radius $r_{\mathrm{c},j}$, and rescaling the size of each bin ($\Delta r_{\mathrm{c},j}$) from inside out to produce the appropriate number density with respect to imposed EFF profile:

[TABLE]

where $\langle m_{\mathrm{star}}\rangle$ is the average stellar mass and $\rho(r)$ is defined in equation 1.

The results of this process are filament-like structures, as shown in figure 1, in which FRAC and FILA models are compared. In a FRAC model the density profile has hard edges and a clump-like substructure (as in figure 1a), while in the FILA model stellar density drops off smoothly with radius and produces filament-like substructure (as in figure 1b). This is because the initial clumps from which the stellar density is composed become ‘stretched’ radially when we impose the EFF profile, to produce several elongated distributions of stars. Both types of model demonstrate spatial and kinematic asymmetry with respect to the gas potential.

We draw stellar masses from a Kroupa (2001) IMF:

[TABLE]

where at high mass end $>1\,M_{\odot}$ we use the observed mass function in Cyg OB2 (Wright et al., 2015). In our models stars are not primordially mass segregated.

Reproducing observations requires estimating an appropriate field of view, which is in turn dependent on the definition of the cluster centre. Observationally some discrepancy exists in this definition between different works, although authors generally agree within a few pc (see Wright et al., 2016, for a discussion). For our purposes an approximate estimate of this centre is sufficient as we find that all of the metrics we consider are only weakly dependent on our choice. For UNIF models, the cluster centre remains the centre of mass of the original set up. For FILA and FRAC models, where considerable anistropy exists in the stellar kinematic distribution, it is necessary to estimate the centre of mass for each snapshot in time. We choose an efficient (approximate) algorithm, in which we find the point which maximises the mass within a given projected radius (chosen to be $R_{\mathrm{cent}}=10$ pc) by sampling recursively over a grid of points. Providing the grid is sufficiently highly resolved, the centre we obtain is insensitive to the exact value of $R_{\mathrm{cent}}$ and number of iterations. All subsequent results should be understood in this context.

3.2 Disc evolution model

In this section we discuss the prescription we apply to calculate PPD evolution. Each disc is exposed to an FUV flux resulting from tracking the contributions from the stellar components within a given $N$-body model.

3.2.1 FUV flux and mass loss rate

Flux in the FUV energy range ( $6$ eV$<h\nu<13.6$ eV222Photons with $h\nu>13.6$ eV are considered extreme ultraviolet (EUV).) is the most efficient contributer to external photoevaporation of PPDs in an intermediate regime of distances from massive stars. As discussed in Winter et al. (2018b), extreme ultraviolet radiation can also contribute to the thermal wind, however it only dominates the mass loss rate in regions where FUV flux is low ($F_{\mathrm{FUV}}\lesssim 50\,G_{0}$) or very high ($F_{\mathrm{FUV}}\gtrsim 10^{5}\,G_{0}$). Both of these thresholds lie outside the region of interest in Cyg OB2. Additionally, in these cases disc destruction timescales are extremely long and short respectively, and the additional mass loss contributed by EUV flux in young massive clusters is therefore not significant to our consideration of the survival rates of discs on a timescale of Myr.

To calculate the FUV luminosity for a star of a given stellar mass we apply the precription discussed by Armitage (2000) and Winter et al. (2018b). Model grids of Schaller et al. (1992) contain the total luminosities and effective temperatures of stars of a given mass, for which we use the results for metallicity $Z=0.02$ and the output age closest to $1$ Myr. We then obtain the wavelength dependent luminosity from the atmosphere models of Castelli & Kurucz (2004). In this way we obtain an FUV luminosity for all stars with a mass $m$ in the range $1\,M_{\odot}<m<100\,M_{\odot}$.

For the FUV induced mass loss rates $\dot{M}_{\mathrm{wind}}$ we use the recent grid of models calculated by Haworth et al. (2018b). The grid covers a wide range of parameter space in outer disc radius ($1$ – $400$ au), disc masses ($\sim 10^{-8}$ – $0.1\,M_{\odot}$), FUV field strengths ($10$ – $10^{4}$ $G_{0}$) and stellar masses ($0.05$ – $1.9\,M_{\odot}$). These mass loss rates are interpolated linearly and applied to a viscously evolving disc to establish the expected disc properties in a given cluster environment.

3.2.2 Viscous disc evolution

To calculate the state of PPDs evolving within the cluster we must take into account vicous expansion as well as the photoevaporative mass loss. The viscous disc evolution is modelled using the method of Clarke (2007, and subsequently , and ). In such a paramerisation, viscosity is assumed to scale linearly with radius $r$ within the disc, corresponding to a constant $\alpha$-viscosity parameter (Shakura & Sunyaev, 1973) and a temperature which scales with $r^{-1/2}$. The accretion rate at the inner edge of the disc is initially

[TABLE]

where $M_{\mathrm{disc,}0}$ is the initial disc mass, $H_{1}$ and $\Omega_{1}$ are the scale height and Keplerian frequency at the disc scale radius $R_{1}$ (see Hartmann et al., 1998). Providing that the outer radius of the disc $R_{\mathrm{disc},0}$ is significantly greater than $R_{1}$, the initial surface density can be written

[TABLE]

and $\Sigma_{0}(r>R_{\mathrm{disc},0})=0$. We choose $R_{\mathrm{disc},0}=2.5R_{1}$, such that the integral of equation 7 over the disc area yields $92\%$ of $M_{\mathrm{disc},0}$. We take the scale radius to be

[TABLE]

Unless otherwise stated we will assume that $M_{\mathrm{disc},0}$ is uniformly distributed between $0.01$ – $0.1\,m_{\mathrm{star}}$ (Andrews et al., 2013); then equation 8 means that the distribution of initial surface densities at $R_{1}$ remains independent of stellar host mass. The scale height $H$ is proportional to the radius throughout the disc, and we choose $H_{1}/R_{1}=0.05$. The maximum initial accretion rate (equation 6) as a function of stellar mass is indicated in figure 2. The associated viscous timescale is

[TABLE]

or $\tau_{\mathrm{visc}}\approx 0.5$ Myr for a solar mass star with $\alpha=10^{-2}$, $R_{1}=40$ au. Equations 8 and 9 yield $\tau_{\mathrm{visc}}\propto m_{\mathrm{star}}^{1/4}$.

During the course of this work we will explore the effect of altering disc initial conditions on their final properties, however we always consider a distribution of initial disc masses. Allowing a range of initial disc masses accounts for the observed range of stellar ages (since internal processes deplete the disc over time) and variable disc initial conditions. Additionally, for a given FUV flux environment some fraction of discs survive, and the findings of Guarcello et al. (2016) indicate that this survival fraction reduces monotonically with increasing $F_{\mathrm{FUV}}$. We find that, depending on initial disc properties, dynamical mixing between regions of different FUV flux and the range of $F_{\mathrm{FUV}}$ within a single bin alone are insufficient to produce the observed survival fractions (i.e. observed survival fractions between $F_{\mathrm{FUV}}$ bins do not jump rapidly from $\sim 0\%$ to $\sim 100\%$ at a certain threshold). The chosen initial conditions and the variation between discs are therefore important in reproducing observations. Whether or not this dispersion is inherited from a tight correlation between stellar mass and PPD initial conditions is explored by considering host mass independent disc models in Section 4.4.3.

Using the initial conditions described above we calculate the evolution of the disc over a 1D grid, spaced evenly in $r^{1/2}$. The cell at the inner edge has a zero torque boundary condition, and the cell at the outer edge experiences mass depletion due to the photoevaporation induced wind (i.e. for a given mass loss rate $\dot{M}_{\mathrm{wind}}$, the disc material is removed from the outer edge). The edge cell is then redefined depending on whether there is net mass loss or mass accumulation (see Section 4.2 in Winter et al., 2018b).

4 Results and discussion

4.1 Modelling approach

We aim to produce an $N$-body model, including initial gas potential, with a self-consistent treatment for the photoevaporation of the PPD population. We simplify the modelling procedure by the following approach:

First we estimate the gas expulsion timescale $\tau_{\mathrm{gas}}\equiv\tau_{\mathrm{delay}}+\tau_{\mathrm{exp}}$ by considering the observed surviving disc fractions as a function of FUV field strength (Section 4.2). The presence of primordial gas influences the models in two ways: it imposes a gravitational potential on the stellar population, and suppresses photoevaporation by extincting FUV photons. Using the latter effect we can calibrate the period of efficient exposure to the observed survival rate of PPDs. There exists a degeneracy in this calculation with the assumed disc viscosity, and we explore the influence of varying both parameters (see Section 4.2). 2. 2.

Having established the rate of gas expulsion, we apply the appropriate time-dependent potential to establish the dynamical evolution of the stellar population consistent with kinematic and spatial data. In Section 4.3.1 we deduce the initial conditions required to reproduce the observed anistropic velocity dispersion. In Section 4.3.2 we vary the initial stellar mass required to reproduce the observed central density (within $13$ pc of the apparent centre). 3. 3.

Finally we combine our findings into a ‘best-fitting’ model, and explore the evolution of PPDs over time in Section 4.4. This allows us both to test whether external photoevaporation is a viable mechanism for disc depletion in Cyg OB2 and to make predictions regarding the disc properties for future observations.

4.2 Disc fractions and gas expulsion

We first estimate the influence of the timescale for gas expulsion on the PPD population survival fractions. A number of physical effects must be taken into account when considering the statistical distribution of disc fractions. Most obviously, projection effects can lead to a distribution of physical distances between stars for every apparent (projected) separation, and therefore a range of $F_{\mathrm{FUV}}$. Binarity and the initial PPD properties will also influence the total disc fraction. However, to first order, the steepness of the drop-off of the disc fraction with $F_{\mathrm{FUV}}$ indicates the length of time for which external photoevaporation has been an efficient mechanism for disc destruction in a given enivironment. Initial tests suggest that a massive gas mass ($\gtrsim 10^{6}M_{\odot}$) is required to maintain virial equilibrium in the primordial Cyg OB2. Due to extinction of FUV photons, this gas mass is sufficient to dramatically reduce photoevaporation efficiency.

The current PPD population therefore allows us to constrain when the gas component of Cyg OB2 was expelled. The relatively short period of expulsion ($\lesssim 2$ Myr) means that we are free to consider the influence of FUV photons on disc evolution from the time $\tau_{\mathrm{gas}}$ at which gas is completely removed. Practically this means that we can apply a simplified UNIF model with the present day mass and ‘switch on’ photoevaporation at different times. We then compare the disc fractions as a function of FUV flux after $3$ Myr of evolution, and thus estimate the gas expulsion timescale. The FUV flux in those bins is calculated in the same way as in Guarcello et al. (2016), using the projected distance between stars - this we call ‘projected’ $F_{\mathrm{FUV}}$ (as opposed to ‘real’ $F_{\mathrm{FUV}}$, as experienced by a given disc). As an estimate of the influence of dynamical mixing between projected FUV flux bins, we start with a stellar velocity dispersion $17$ km/s and hold the stars in virial equilibrium with an external potential.

The rate at which irradiated discs are destroyed is also dependent on the $\alpha$-viscosity parameter (Section 3.2.2). The chosen $\alpha$ dictates the rate at which a disc viscously expands into a region where photoevaporation rates are high, as well as dictating mass loss through accretion. This adds a degeneracy to our approach which, given uncertainties in $\alpha$, introduces similar uncertainties in $\tau_{\mathrm{gas}}$. We investigate this degeneracy by calculating surviving PPD fractions for a number of different values for $\tau_{\mathrm{gas}}$ ($1$, $2$, $2.5$ and $2.75$ Myr) and $\alpha$ ($5\times 10^{-3}$, $10^{-2}$, $2\times 10^{-2}$).

We calculate disc evolution for a subset of discs with host stars in the mass range $0.5$ – $2$ $M_{\odot}$. Stars less massive than this are not present in the sample used by Guarcello et al. (2016), while the disc mass loss rates calculated by Haworth et al. (2018b) do not apply for higher mass stars. A disc is considered destroyed if it has a mass $<10^{-5}\,M_{\odot}$, which is estimated by Guarcello et al. (2013) as a limit below which SEDs are more difficult to interpret. In fact discs of such low mass ($\sim 10^{-5}\,M_{\odot}$) are destroyed quickly in regions of strong FUV fields, so our results are insensitive to the exact value of this threshold.

The results of this preliminary modelling procedure compared with the observational findings of Guarcello et al. (2016) are summarised in figure 3. We expect disc fractions at lower $F_{\mathrm{FUV}}$ are overestimated as we do not consider other disc dispersal processes (such as internal photoevaporation, see Section 4.4). We therefore focus on matching disc fractions in regions of higher $F_{\mathrm{FUV}}$. We find that $\tau_{\mathrm{gas}}=2.5$ Myr and $\alpha=10^{-2}$ give a good fit to the data for projected $F_{\mathrm{FUV}}\gtrsim 3000\,G_{0}$. Since we have assumed $100\%$ extinction of FUV photons before gas is completely expelled, this $\alpha$ is an upper limit. As discussed above, reducing the expulsion timescale to $\tau_{\mathrm{gas}}=1$ Myr while decreasing the viscosity such that $\alpha=5\times 10^{-3}$ also yields the correct disc fractions. However, if gas ejection was initiated at the time of formation of the most massive stars ($\sim 2$ Myr) this would suggest instantaneous expulsion, which is not physical. In this case early supernovae may be responsible for driving gas mass loss (see discussion in Section 2.1). However, for short $\tau_{\mathrm{gas}}$ Cyg OB2 must have had a extremely large initial stellar mass to maintain the present day central density ($\gg 10^{5}\,M_{\odot}$, see Section 4.3.2 and figure 7), which is neither supported by observations nor computationally practicable for the range of models we wish to explore.

While the simplified model presented in this section is not an accurate representation of the dynamical conditions in the region, it represents a first-order approximation on which we can base our choice of kinematic parameters in Section 4.3. We will again consider the disc fractions for a more realistic model in Section 4.4.

4.3 Stellar population properties

In this section we consider the initial conditions for our $N$-body models of the stellar population of Cyg OB2. We proceed by first aiming to reproduce the observed velocity dispersion in the region by varying the initial velocity distribution and substructure (Section 4.3.1). Subsequently we match the observed central density by varying the initial stellar mass (Section 4.3.2).

4.3.1 Velocity dispersion, anisotropy and substructure

In terms of the stellar population, the first observable quantity we aim to reproduce is the central velocity dispersion and its anisotropy ($\langle v^{2}\rangle^{1/2}\approx 17$ km/s, and the proper motion dispersions $\sigma_{\alpha,\delta}\gtrsim\sigma_{r}$ the line of sight component – see Wright et al., 2016). This is because it is not strongly sensitive to the initial stellar mass of the cluster. We demonstrate this by considering UNIF cluster models (where cluster conditions are relatively non-stochastic) for a range of initial 3D velocity dispersions and stellar masses; we plot the resuts in figure 4. As the potential is dominated by the gas in the cluster, and the stellar component is itself highly supervirial, the final central velocity dispersion is insensitive the total stellar mass.

Figure 4 demonstrates that UNIF models fail to reproduce the observed velocity dispersion. This is because in a non-substructured model we find supressed proper motion velocity dispersions $\sigma_{\alpha,\delta}$ with respect to the radial component $\sigma_{r}$, since stars with high velocities in the plane of the sky preferentially leave the central field of view over time (velocity sorting). To remedy this, we need to incorporate substructure into our model.

To demonstrate the influence of substructure on the components of the velocity dispersion within a finite field of view, we run FILA and FRAC models with large scale substructure ($P_{0}=1$), and rotate a snapshot at $3$ Myr such that the radial velocity dispersion is approximately minimised. The results of this process are shown in figure 5, where the decomposed velocity components are again compared with the initial three dimensional velocity dispersion. The anisotropies in the velocity dispersion are reproduced in both FILA and FRAC cluster models.

The degree of anisotropy in the velocity dispersion is dependent on the scale of substructure. Figure 6 shows the results for initially smaller filaments/clumps (using $P_{0}=2$). While there still exist stochastic fluctuations in the components of the velocity dispersion, the degree of anisotropy is not sufficient to yield proper motion velocity dispersions greater than the radial component. This suggests that the initial region consisted of large clumps or filaments of mass $\sim 10^{4}$ $M_{\odot}$.

We estimate that the required initial velocity dispersion in Cyg OB2 is $\langle v_{0}^{2}\rangle^{1/2}\sim 50$ km/s, although the stochasticity of the substructured models makes precise estimates impractical. In future models we will assume this is the initial three dimensional velocity dispersion. Other required parameters, such as the initial stellar mass, are chosen to be consistent with this property.

4.3.2 Central mass

We wish to alter the initial stellar mass such that the central density at the end of the simulation is consistent with observations. A cluster model with a central mass of $1.6\times 10^{4}\,M_{\odot}$ after $3$ Myr of evolution is required, where $\langle v_{0}^{2}\rangle^{1/2}=50$ km/s as discussed in Section 4.3.1. To find the appropriate initial mass we run UNIF, FRAC and FILA models at stellar masses $M_{\mathrm{stars}}=2\times 10^{4}\,M_{\odot}$, $4\times 10^{4}\,M_{\odot}$ and $8\times 10^{4}\,M_{\odot}$. The results are shown in figure 7. We find that $M_{\mathrm{stars}}=8\times 10^{4}\,M_{\odot}$ (figure 7c) is sufficient to yield the required central mass after $3$ Myr of evolution.

Neither the initial mass nor the substructure has a significant effect on the fraction of the stellar mass which remains within a $13$ pc projected radius (approximately a quarter in each case). As in the case of our velocity investigation, this is due to the fact that that the stellar component of the potential energy is much smaller than the total kinetic energy.

The central mass in all of our models undergoes a similar temporal evolution. In figure 7 we see an initial rapid mass loss as stars with the highest energies escape the potential. This is because we do not truncate the velocity dispersion such that escapers are initially forbidden; it is not clear whether or not this is physically realistic, and the effect of choosing such initial conditions is discussed in Section 4.3.3 in the context of primordial gas mass. However, since stars which escape the central regions are not considered in our PPD models and the gravitational potential is dominated by the gas component, whether or not these early escapers are initially included in the model is of secondary importance. After this initial decline in mass, some high energy stars remain bound, and therefore return to the central regions, causing a modest oscillation in the mass. This is again a consequence of using a Boltzmann velocity distribution without truncating the high velocity end. The magnitude of this oscillation is more variable for substructured regions where velocities (and therefore kinetic energies) are correlated. Gas expulsion (starting at $\tau_{\mathrm{delay}}=1$ Myr and continuing over $\tau_{\mathrm{exp}}=1.5$ Myr) results in a decrease in gravitational potential, and as the number of stars with energies sufficient to escape increases, the central mass decreases.

4.3.3 Summary of best fitting model

By considering gas expulsion, stellar mass, velocity dispersion and initial substructure, we have found a model for the evolution of Cyg OB2, summarised by the parameters in Table 1, which fits observations of the stellar population. While we refer to this as our ‘best-fit’ model, this is to be understood as the result our process of deduction in terms of the appropriate parameters, and not as the optimisation of a statistical metric or parameter space exploration. In our models we find gas expulsion was completed $\sim 0.5$ Myr ago (for a PPD viscous timescale of $0.5$ Myr), the initial velocity dispersion was $\sim 50$ km/s and the initial cluster mass was $\sim 8\times 10^{4}\,M_{\odot}$. We further suggest that the largest scales of initial coherent clumps within the primordial Cyg OB2 had a mass of $\sim 10^{4}\,M_{\odot}$. No significant dynamical differences can be found between FRAC and FILA type models at the present time, and we hereafter use a FILA model in our analysis of the disc population (Section 4.4).

The gas mass required for initial virial equilibrium is $M_{\mathrm{gas}}\sim 8\cdot 10^{6}\,M_{\odot}$ (with a scale parameter $a_{\mathrm{gas}}=10$ pc). This would make the primordial GMC massive compared to known Milky Way molecular clouds, although the census is not complete (see Longmore et al., 2014, for a review). It is also possible that the initial velocity dispersion is overestimated due to the truncation of equation 4 at high velocities (simulations suggest stars have a subvirial initial velocity dispersion with respect to the primordial gas - see Offner et al., 2009; Kruijssen, 2012). While this would not influence the central velocity dispersion (high velocity stars leave the centre in any case) it would reduce the number of escapers early on in the cluster evolution and therefore the required initial stellar mass. It would also reduce the gas mass necessary for virial equilibrium as $M_{\mathrm{gas}}\propto\langle v^{2}\rangle$ if $M_{\mathrm{gas}}\gg M_{\mathrm{stars}}$ (equation 3). If our gas mass estimate is accurate then this makes the star formation efficiency $\sim 1\%$, although this is probably a lower limit (an upper limit on $M_{\mathrm{gas}}$).

4.3.4 Substructure and expansion observables

Alternate kinematic constraints not considered in the previous analysis include the measures of substructure and the absence of expansion signatures in the stellar kinematics (Wright et al., 2016). We find that determining these metrics is problematic for a given cluster model. This is because the values obtained differ stochastically depending on initial conditions, the time of ‘observation’ and the subset of stars used in taking a measurement. To obtain an accurate probability of finding the observed values for these metrics, a large number of models would need to be tested, which would be computationally expensive. However, for our chosen model the velocity field is illustrated in figure 8. We find that the correlation between positions and proper motions is clear (figure 8a), while figure 8b does not show clear evidence of expansion. However, the stellar components are in fact expanding globally since they are unbound.

To illustrate this point more fully, we define the expansion parameter:

[TABLE]

where $T^{+/-}$ is defined as the total kinetic energy of stars directed in the postive/negative projected radial direction (in the plane of the sky). Thus $\mathcal{E}\rightarrow 1$ or [math] if the velocity dispersion indicates rapid expansion or contraction respectively. A value of $\mathcal{E}\approx 0.5$ would usually be taken as evidence that a stellar population is not expanding. In Figure 9 we show the cumulative distribution of the measured expansion parameter for random subsets of $800$ stars in the Wright et al. (2016) field of view (the central $8$ pc $\times$ $8$ pc). We find that a wide range of values for $\mathcal{E}$ is possible at any given time. Depending on the time at which the velocities are observed and the chosen subsample, our model is found to be consistent with an observed value $\mathcal{E}\approx 0.5$. This is because the $\mathcal{E}$ distribution varies considerably and non-monotonically in time even for a single model. However, at any given time the stellar population is expanding (filaments are moving away from each other). Alternative geometrical signatures may be more successful at gauging such expansion. Due to the stochasticity and wide range of possible $\mathcal{E}$ values for our model, we conclude that $\mathcal{E}$ alone is not a sufficient metric to draw conclusions on the expansion of a substructured association. For further discussion on the kinematic indicators of expansion in OB associations, see Baumgardt & Kroupa (2007) and Ward & Kruijssen (2018).

In the remainder of this work we will first revisit the disc population in our model, checking that the population is consistent with the known disc fractions in the region. Subsequently we will explore predictions for the disc mass and radius distribution relevant for future observations in the region.

4.4 Disc properties

In what follows we will consider a FILA model with the properties described in Table 1. We fix the $\alpha$-viscosity of the disc population with the derived value $\alpha=10^{-2}$ (Section 4.2). As in Section 4.2 all discs are assumed to have an initial mass that is uniformly distributed between $1$ – $10\%$ of their host star mass (see Section 3.2.2), and scale radius as in equation 8. We consider the same range of stellar masses $0.5$ – $2$ $M_{\odot}$ and calculate the evolution of a subset of $5000$ discs.

4.4.1 Disc fractions in best-fit model

The disc fraction distribution in our best-fit model after $\sim 3$ Myr is summarised by figure 10, in which we show the surviving disc fraction as a function of projected FUV flux (figure 10a) with the spatial distribution and projected $F_{\mathrm{FUV}}$ in figure 10b. We find that after $3$ Myr the number of surviving discs is slightly overestimated in our model, particularly at the lower FUV fluxes. There are two reasons why we expect this to be the case. Firstly, extinction in the FUV may decrease before gas is fully removed from the cluster. If gas is expelled due to the flux from massive stars, then the central, highly irradiated environments would become less dense and allow efficient photoevaporation at earlier times. Additionally, clumpy gas distributions may have a similar effect in reducing extinction during this period. Secondly, internal photoevaporation due to the stellar host depletes the gas content even when a PPD evolves in isolation. This speeds up destruction timescales, particularly in regions of lower $F_{\mathrm{FUV}}$, where the mass loss rates induced by internal and external photoevaporation become comparable.

The influence of the above considerations is uncertain and we therefore do not attempt to model them here. Our model does however reproduce the correct disc fractions within a reasonable period of time, particularly at $F_{\mathrm{FUV}}\gtrsim 3000\,G_{0}$. We conclude that the observations of Guarcello et al. (2016) can be explained by external photoevaporation of the PPD population in Cyg OB2.

4.4.2 Disc mass and FUV flux environment

To make predictions about the disc population in Cyg OB2 for comparison with future observations we consider PPD mass as a function of FUV flux (figure 11). We find that, if used in isolation, the local flux experienced by a given star is a poor predictor of the disc mass. Correlation with real or projected FUV flux only becomes clear when the host mass, to which the final PPD mass is closely correlated (as demonstrated by the colour gradient in figure 11 and discussed in Section 4.4.3), is also taken into account. In particular, we expect that low mass stars ($<1\,M_{\odot}$) in the centre of Cyg OB2 host exclusively low mass discs ($\lesssim 10^{-3}\,M_{\odot}$, if any).

Without taking into account stellar mass, can we find differences in disc properties between PPDs in apparently high versus low FUV flux environments? To answer this, we must consider the sensitivity limit for future observations. Considering ALMA band 6 sensitivity, reasonable integration times ($\sim 30$ minutes) for a survey sample suggest flux densities down to $F_{\nu}(850\,\mu\mathrm{m})\sim 40$ $\mu$Jy can be detected. At the distance of Cyg OB2 this means that dust masses can be established down to a few $M_{\earth}$ (Andrews & Williams, 2005). The corresponding total disc masses are $M_{\mathrm{disc}}\sim 10^{-3}\,M_{\odot}$ if the gas to dust ratio is $\Sigma_{\mathrm{gas}}/\Sigma_{\mathrm{dust}}=10^{2}$. Given that the PPDs are likely to be gas depleted by external photoevaporation, the latter assumption is probably not accurate for many discs (Ansdell et al., 2016), and this complicates the interpretation of observations. Nonetheless, we show the cumulative PPD mass fraction after $3$ Myr for two FUV flux bins in figure 12. While disc masses are indeed suppressed at higher projected FUV flux, considering the sensitivity limit of ALMA makes finding differences between the two populations impractical. A sample of several $100$s of PPDs would be required to find a difference between the total population in high and low FUV flux environments (including non-detections). We find that similar sample sizes would be required to find differences in disc outer radius distributions in the two environments.

4.4.3 Stellar mass independent disc initial conditions

In obtaining the results in Section 4.4.2 we have already assumed that the initial disc mass is dependent on the stellar mass of the host. We wish to confirm our results are not sensitive to this assumption. We therefore recalculate the PPD evolution in our model under the assumption that the initial disc mass is uncorrelated to the stellar mass. To this end we draw the initial disc masses from a log-normal distribution:

[TABLE]

where the random variable $X\sim\mathcal{N}(0,1)$ is drawn from a standard normal distribution, uncorrelated with $m_{\mathrm{star}}$. The values of $\mu=-3.25$ and $\sigma=0.7$ are chosen such that the mean and dispersion of $M_{\mathrm{disc},0}$ match those chosen for the stellar mass dependent initial conditions. The scale radius $R_{1}$ is again defined according to equation 8. This distribution of initial conditions recovers similar disc fractions as a function of projected $F_{\mathrm{FUV}}$ (figure 13). The resulting disc mass distribution is shown in figure 14. Our findings do not differ greatly from those in figure 11, except for a predictable weaker correlation between final disc mass and host mass.

4.4.4 Disc mass dependence on host mass

The relationship between final disc mass and stellar mass is considered in figure 15a (figure 15b) for disc masses initially correlated (uncorrelated) with the host star mass. We obtain power law indices ($M_{\mathrm{disc}}\propto m_{\mathrm{star}}^{\beta}$) $\beta=3.93\pm 0.11$ and $2.75\pm 0.12$ respectively. In both cases this is a substantially more superlinear relationship than in local PPD populations that have not been significantly photoevaporated ($1<\beta<1.9$ Andrews et al., 2013; Pascucci et al., 2016). Physically this is because $\dot{M}_{\mathrm{wind}}$ is strongly dependent on stellar mass such that discs around low mass stars are depleted much faster (due to a shallower potential) than those around high mass stars for a fixed flux and disc radius (Haworth et al., 2018b). In regions where a disc population has undergone significant photoevaporation we then expect lower mass stars to host lower mass PPDs. Even if it proves impossible to obtain sufficiently large samples of mm-based mass determinations in Cyg OB2 in order to conduct the comparison made in figure 12, an alternative way to test the role of external photoevaporation would be to examine any evidence for steep disc mass–stellar mass relationship. Indeed, Ansdell et al. (2017) find a steepening of the relationship with age across different regions that could be due to external photoevaporation. To the contrary, Eisner et al. (2018) find a shallow disc mass–host mass relationship in the strongly irradiated discs of the Orion Nebula Cluster; this could be the result of the youth or complex formation history of the stellar population (Beccari et al., 2017; Kroupa et al., 2018). Disentangling these effects requires a detailed modelling of those regions, such as that presented in this paper.

To estimate whether detecting large $\beta$ values is possible in Cyg OB2, we extract a subset of $20$ discs from our model, choosing them such that they lie in a region of projected FUV flux $3000\,G_{0}<F_{\mathrm{FUV}}<8000\,G_{0}$. A lower limit of $3000\,G_{0}$ is appropriate since below this figures 10a and 13 suggest that alternative processes to external photoevaporation have a comparable influence on disc evolution. We also select stellar masses to be approximately uniformly distributed between $0.5$ – $2\,M_{\odot}$. We assume error of $30\%$ in $M_{\mathrm{disc}}$, with $M_{\mathrm{disc}}<10^{-3}\,M_{\odot}$ constituting a non-detection, and a $10\%$ error in $M_{\mathrm{star}}$. We then attempt to fit a power law using the Linmix package (Kelly, 2007) for a number of different PPD subsets. An illustrative example of such an exploration is shown in figure 16. In all cases we find that the fit effectively discounts power law relations with $\beta\leq 1.9$. Values of $\beta\leq 1.9$ are similarly discounted in the model with stellar host mass independent disc initial conditions. Our findings suggest that, if a sample is carefully selected, evidence of external photoevaporation should be detectable in future observations of Cyg OB2.

One might ask whether it is possible to find signatures of underlying substructure within a population of PPDs. Given that the dynamical history is linked to the irradiation by FUV photons, it may be possible to detect distinct groups in the mass distribution of a disc sample. As discussed previously, in an externally photoevaporated population, disc properties are more strongly correlated with host mass than FUV flux. Therefore considering disc properties as a function of host mass is the best chance of finding distinct groups within a sample. This possibility is explored in figure 15, where the largest scale fractal membership is indicated for each PPD. We find that the mass distribution does not demonstrate a clear segregation between groups. The subgroups are in themselves massive ($\sim 10^{4}\,M_{\odot}$) and the high mass end of the IMF is therefore well sampled in each filament. It is possible that for smaller scale substructure (with filaments of mass $<10^{3}\,M_{\odot}$, see Winter et al. 2018b) a less well sampled IMF might mean that different filaments have quite different local $F_{\mathrm{FUV}}$. In this case we should expect to see a distinct PPD mass and outer radius distributions between filaments.

4.4.5 Disc radii

An alternative observable property, the outer disc radius $R_{\mathrm{out}}$, may prove to be a better probe of substructure in future observations. This is because, while $M_{\mathrm{disc}}$ is depleted over the entire lifetime of the disc, the outer radius in regions of strong FUV flux is set by a balance between photoevaporative mass loss and viscous expansion on short timescales. We have so far neglected discussion of disc radius simply because resolving disc radii down to $\sim 10$s au is challenging in a region such as Cyg OB2; the highest resolution of ALMA is $\sim 0.02^{\prime\prime}$, corresponding to $\sim 30$ au at $1.33$ kpc distance. However, because it may be possible to resolve disc radii to sufficient accuracies in closer young stellar environments, we illustrate the distribution of disc radii in figure 17. There are some indications of different radius distributions between different fractal associations, but no clear segregation between them. The majority of discs have $R_{\mathrm{disc}}<60$ au which would require $<0.05^{\prime\prime}$ resolution at the distance of Cyg OB2. We note that outer disc radius measurements can be subject to large uncertainties, and that the dust emission is frequently less extended than the gas (e.g. Birnstiel et al., 2010; Guilloteau et al., 2011; de Gregorio-Monsalvo et al., 2013). This discrepancy may be less problematic in regions where disc radii are externally suppressed (for example close stellar encounters may result in wavelength independent PPD outer radius measurements, as in the case of the disc around HV Tau C - Monin & Bouvier, 2000; Winter et al., 2018c). This assumption requires further study into gas-dust physics in externally photoevaporating discs (see Cleeves et al. 2016 for discussion of gas and dust evolution in an FUV irradiated PPD). Overall we see the same trend as in disc masses that the outer radius is correlated with host mass. This is again because $\dot{M}_{\mathrm{wind}}$ decreases with $m_{\mathrm{star}}$, such that the radius at which viscous expansion is in equilibrium with FUV induced mass loss is more extended for more massive host stars.

5 Conclusions

Using $N$-body simulations and viscous disc evolution models we have successfully reproduced the properties of Cygnus OB2, including the stellar kinematics and the surviving PPD fractions as a function of projected FUV flux. Our modelling supports the following scenario:

• •

If the viscous evolution of PPDs is well described by a Shakura & Sunyaev $\alpha$-parameter $\alpha=10^{-2}$ ($\tau_{\mathrm{visc}}\approx 0.5$ Myr for a solar mass star), then expulsion of the primordial gas content in the region must have concluded $\sim 0.5$ Myr ago. Approximately $0.5$ Myr of exposure to strong FUV fields is required to reproduce the current disc survival rates as a function of (projected) flux. This value of $\alpha$ will be an upper limit if FUV extinction due to the primordial gas was not efficient (because of a clumpy spatial distribution, for example).

• •

The initial three dimensional velocity dispersion must have been $\sim 50$ km/s in order to be consistent with the present day central velocity dispersion. This is large with respect to observed stellar populations, even for OB associations (e.g. Ward & Kruijssen 2018 and references therein). However, velocity gradients $\sim 10$ kms*-1pc-1* are found across giant molecular clouds, particularly in the case of cloud-cloud collision (Wu et al., 2015; Bisbas et al., 2018; Pols et al., 2018).

• •

Anisotropy in the present day velocity dispersion requires the presence of primordial subtructure. We find that Cyg OB2 must have been comprised of large scale filaments (or fractal clumps), with a mass $\sim 10^{4}$ $M_{\odot}$.

• •

In such a model, a gas content of $\sim 8\times 10^{6}\,M_{\odot}$ is required to maintain initial virial equilibrium. This is massive compared to the known distribution of molecular cloud masses in the Milky Way (Longmore et al., 2014).

• •

The total stellar mass required to sustain a sufficient central mass after $3$ Myr is $\sim 8\times 10^{4}$ $M_{\odot}$. This suggests a star formation efficiency of $\sim 1\%$.

• •

The apparent lack of expansion measured by Wright et al. (2016) need not be interpreted as evidence of no physical expansion of the stellar population. We find that our (expanding) model is equally consistent with an expansion parameter $\mathcal{E}\approx 0.5$. Such a measurement can only be interpreted probabilistically for an association with large scale substructure.

The primary caveat of these conclusions is that the gas expulsion timescale derived by disc fractions is degenerate with $\alpha$. Our choice ($\alpha=10^{-2}$, although this is in turn dependent on choice of scale radius $R_{1}$ - see equation 9) yields a good fit to the disc fractions and also allows us to reproduce kinematic observations with gas expulsion timescale $\tau_{\mathrm{gas}}=2.5$ Myr. In future, observations of PPD populations in regions with a strong FUV flux but with more modest velocity dispersions (and primordial gas density) might offer further constraints for $\alpha$.

Finally we make predictions for future observations of the disc mass and radius distribution in Cyg OB2. We find that samples of PPDs in highly FUV irradiated environment have significantly reduced masses and outer radii than in regions of more modest flux. However, taking into account sensitivity limits, statistically observing these differences requires sample sizes of $\gtrsim 100$s of discs in both high and low FUV flux bins. This is similarly true for disc outer radii, and makes probing the difference between disc property distributions as a function of FUV environment directly impractical at present.

In Cygnus OB2 we expect a strong correlation between stellar mass and disc mass ($M_{\mathrm{disc}}\propto m_{\mathrm{star}}^{\beta}$, with $\beta>2$) which is a consequence of the fact that, for a disc with given radius and in a given FUV environment, the mass loss rate is higher in the shallower potential of low mass stars. This effect should be clear in a sufficiently large sample of irradiated discs in any environment. In Cygnus OB2, we demonstrate that - taking into account the finite sensitivity of ALMA - it would be sufficient to target around $20$ stars of known mass in order to demonstrate a value of $\beta$ that is steeper than the canonical value $\sim 1$ – $1.9$ which is seen in non-irradiated disc samples. We conclude that, pending empirical confirmation, discs around low mass stars born in environments of strong FUV flux are likely to have a significantly depleted mass budget for planet formation. Indeed, unless dust within PPDs can rapidly grow to size scales where it is immune to photoevaporative stripping (e.g. Youdin & Goodman, 2005), planet formation may be completely supressed for such host stars (Haworth et al., 2018a).

Acknowledgements

We would like to thank Tom Haworth for providing access to the FRIED grid of mass loss rates ahead of publication. AJW thanks Sverre Aarseth and Long Wang for helpful discussion regarding $N$-body simulations. We thank the anonymous referee for a considered report which improved the clarity of the manuscript. This work has been supported by the DISCSIM project, grant agreement 341137 funded by the European Research Council under ERC-2013-ADG. AJW thanks the Science and Technology Facilities Council (STFC) for their studentship.