1 Introduction

Planets form in the discs of dust and gas known as ‘planet-forming’ or ‘protoplanetary’ discs. These discs are found around most young stellar objects (YSOs) in star-forming regions up to ages of \(\sim 3{-}10\) Myr (e.g. [1, 2]). Since 2015, Atacama Large Millimetre/sub-Millimetre Array (ALMA) has instigated a revolution in our understanding of these objects. In particular, the high-resolution sub-mm continuum images, tracing the dust content of protoplanetary discs, not only offer an insight into the dust masses and radii (e.g. [3,4,5,6,7]), but also a wealth of internal substructures which serve as a window into planet formation processes (e.g. [8,9,10]). The nearby (distance \(D\lesssim 150\) pc) discs in the low-mass star-forming regions such as Taurus and Lupus can be sufficiently resolved to expose rings, gaps and spirals in dust, as well as kinematic signatures in gas [11,12,13,14,15], which may either be signatures of planets or the processes that govern their formation.

Prior to ALMA, some of the first images of discs were obtained around stars in the Orion Nebula cluster (ONC) at a distance of \(D\sim 400\) pc away, much further than the famous HL Tau [8, 16]. These discs exhibit non-thermal emission in the radio [17,18,19,20], with resolved ionisation fronts and disc silhouettes seen at optical wavelengths [21,22,23,24]. Identified as ‘solar-system sized condensations’ by Churchwell et al. [17], they were estimated to be undergoing mass loss at a rate of \(\dot{M} \sim 10^{-6}\,M_\odot\) yr\(^{-1}\). O’dell et al. [23] confirmed them to be irradiated protoplanetary discs, dubbing them with the contraction ‘proplyds’. While this term has now been dropped for the broader class of protoplanetary disc, it is still used in relation to discs with a bright, resolved ionisation front owing to the external UV radiation field. These objects are of great importance in unravelling what now appears to be an important process in understanding planet formation: external photoevaporation.

The process of external photoevaporation is distinguished from the internal photoevaporation (or more commonly just ‘photoevaporation’) of protoplanetary discs via the source that is responsible for driving the outflows. Both processes involve heating of the gas by photons, leading to thermal winds that deplete the discs. However, internal photoevaporation is driven by some combination of FUV (far-ultraviolet), EUV (extreme-ultraviolet) and X-ray photons originating from the central host star (e.g. [25,26,27]), depleting the disc from inside-out (e.g. [28,29,30,31,32]). External photoevaporation, by contrast, is driven by an OB star external to the host star-disc system. The high FUV and EUV luminosities of OB stars, combined with the heating of outer disc material that experiences weaker gravity, can result in extremely vigorous outside-in depletion [33]. Indeed, as inferred by Churchwell et al. [17], the brightest proplyds in the ONC exhibit wind-driven mass loss rates of up to \(\dot{M}_{\mathrm{{ext}}} \sim 10^{-6}\,M_\odot\) yr\(^{-1}\) (see also, e.g. [34, 35]).

At present, external photoevaporation remains an under-studied component of the planet formation puzzle. It has long been thought that the majority of star formation occurs in clusters or associations, which include OB stars [36], and that mass loss due to external photoevaporation can be more rapid than the \(\dot{M}_{\mathrm {acc/int}} \sim 10^{-8}\,M_\odot\) yr\(^{-1}\) typical of stellar accretion (e.g. [37]) and internal photoevaporation (e.g. [32]). Thus, it is reasonable to ask why external photoevaporation has, until now, not featured heavily in models for planet formation. While there is no one answer to this question, two important historical considerations are compelling in this context:

  • The first is a selection bias. Low mass star-forming regions (SFRs) are more numerous, so they are also the closest to us (although nearby star formation is probably also affected by the expansion of a local supernova bubble —[38]). Surveys with, for example, ALMA, have rightly prioritised bright and nearby protoplanetary discs in these low mass SFRs. However, when accounting for the relative number of stars, these regions are not typical birth environments for stars and their planetary systems (see Sect. 5.2). This has motivated recent studies of more typical birthplaces that experience strong external irradiation from neighbouring massive stars (e.g. [4, 6, 39]).

  • The second is the well-known ‘proplyd lifetime problem’ in the ONC (e.g. [40]). The presence or absence of a near-infrared (NIR) excess, indicating inner material, was one of the earliest available probes for studying disc physics. The high fraction (\(\sim 80\) percent) of stars with a NIR excess in the ONC [41] was an apparent paradox that appeared to undermine the efficacy of external photoevaporation in influencing disc populations. We discuss this problem in Sect. 5.4.

Given the apparently high fraction of planet-forming discs that undergo at least some degree of external photoevaporation, this process has particular relevance for modern efforts to connect observed disc populations to exoplanets (e.g. [42]). Now is therefore an opportune time to take stock of the findings from the past three decades.

In this review, we first summarise the theory of external photoevaporation in Sect. 2, including both analytic estimates that provide an intuition for the problem and state-of-the-art simulations and microphysics. We address the observational signatures of disc winds in Sect. 3, including inferences of individual mass loss rates. We consider the role of external photoevaporation for disc evolution and planet formation in Sect. 4, including evidence from disc surveys. In Sect. 5 we contextualise these studies in terms of the physics and demographics of star-forming environments. Finally, we summarise the current understandings and open questions in Sect. 6.

2 Theory of externally driven disc winds

2.1 The most basic picture

We begin by reviewing our theoretical understanding of external photoevaporation. At the very basic level, one can determine whether material will be unbound from a potential by comparing the mean thermal velocity of particles (i.e. the sound speed) in the gas with the local escape velocity. In an isothermal system, equating the sound speed and escape velocity yields the gravitational radius, beyond which material is unbound

$$\begin{aligned} R_{\mathrm{{g}}} = \frac{GM_*}{c_{\mathrm{{s}}}^2} \end{aligned}$$
(1)

where \(M_*\) is the point source potential mass, and \(c_{\mathrm{{s}}}\) is the sound speed. Therefore, if an isothermal disc were to extend beyond \(R_{\mathrm{{g}}}\), then material would be lost in a wind. Consider now an isothermal disc that is entirely interior to the gravitational radius and is hence fully bound (upper panel of Fig. 1). If this disc is then externally irradiated, the temperature and hence sound speed increase, driving down the gravitational radius to smaller values, potentially moving interior to the edge of the disc and unbinding its outer parts (lower panel of Fig. 1).

Fig. 1
figure 1

A schematic of the basic picture of external photoevaporation in terms of the gravitational radius. The gravitational radius is that beyond which mean thermal motions (propagating at the sound speed) exceed the escape velocity and are unbound. In this basic picture, a disc smaller than the gravitational radius will hence not lose mass. External UV irradiation heats the disc, leading to faster mean thermal motions (a higher sound speed), driving the gravitational radius to smaller radii and unbinding material in the disc

The details of external photoevaporation do get substantially more complicated than the above picture, for example with pressure gradients helping to launch winds interior to \(R_{\mathrm{{g}}}\). However, this picture provides a neat basic insight into how external photoevaporation can instigate mass loss in otherwise bound circumstellar planet-forming discs.

2.2 EUV- and FUV-driven flows

Fig. 2
figure 2

Schematic diagrams of the flow structure in EUV (Fig. 2a) and FUV (Fig. 2b) externally driven winds, adapted from Johnstone et al. [33]. In the EUV-driven wind, the flow from the disc edge at radius \(R_{\mathrm {d}}\) travels at a subsonic velocity through the thin photodissociation region (PDR) of thickness \(xR_{\mathrm {d}}\) with \(x\lesssim 1.5\) before reaching the ionisation front (IF). Mass loss is therefore determined by the thermal pressure at the IF. For \(x\gtrsim 1.5\), the wind is launched from the disc edge at a supersonic velocity, producing a shock front that is reached before the IF. In this case, the mass loss rate is determined by the thermal conditions in the PDR

2.2.1 Flux units

Before considering the physics of EUV- and FUV-driven winds, a note on the units canonically used to measure UV fluxes is necessary. While the ionising flux is usually measured in photon counts per square centimetre per second, FUV flux is normally expressed in terms of the Habing unit, written \(G_0\) [43]. This is the flux integral over \(912-2400\mathring{\hbox {A}}\), normalised to the value in the solar neighbourhood, i.e.

$$\begin{aligned} \left( \frac{F_{\mathrm {FUV}}}{1\,G_0}\right) = \int _{912\mathring{\hbox {A}}}^{2400\mathring{\hbox {A}}} \frac{F_\lambda d\lambda }{1.6\times 10^{-3}\,\text {erg s}^{-1}\,\text {cm}^{-2}}. \end{aligned}$$
(2)

Another measure of the FUV field strength is the Draine unit [44], which is a factor 1.71 larger than the Habing unit. Hence, \(10^3 \, G_0 \approx 585\) Draines. We highlight both because two similar units that vary by a factor of order unity can and does lead to confusion. For reference, the UV environments that discs are exposed to in star-forming regions ranges from \(\ll 1\) (i.e. embedded discs) to \(\sim 10^{7}\)G\(_0\) (discussed further in Sect. 5). For the sake of clarity, we will hereafter consider low FUV environments to be those with \(F_{\mathrm {FUV}}\lesssim 100\, G_0\), intermediate environments greater than this up to \(F_{\mathrm {FUV}} \lesssim 5000 \, G_0\) and high FUV environments for any FUV fluxes \(F_{\mathrm {FUV}}\gtrsim 5000\, G_0\).

2.2.2 Flow geometry

The basic physical picture of an externally photoevaporating protoplanetary disc was laid out by Johnstone et al. [33] and has remained largely the same since. We summarise that picture here because it is useful for what follows, but refer interested readers to the more detailed discussion in that original work as well as that of Störzer and Hollenbach [40].

The heating mechanism that launches the thermal wind may be driven by ionising EUV photons (energies \(h\nu >13.6\) eV), heating gas to temperatures \(T\sim 10^4\) K, or photodissociation FUV photons (\(6\,\mathrm {eV}< h\nu < 13.6\) eV), yielding temperatures of \(T\sim 100{-}1000\) K [45, 46]. The EUV photons penetrate down to the ionisation front, at radius \(R_{\mathrm {IF}} = (1+x)R_{\mathrm {d}}\) from the disc-hosting star, with disc outer radius \(R_{\mathrm {d}}\). The atomic gas outside of \(R_{\mathrm {IF}}\) is optically thin for the FUV photons, which penetrate down to \(R_{\mathrm {d}}\) producing a neutral PDR of thickness \(x R_{\mathrm {d}}\). Whichever photons drive the photoevaporative wind, the flow geometry far from the disc surface is approximately spherically symmetric, since it is accelerated by the radial pressure gradient (cf. the Parker wind). We will start with this approximation of spherical geometry, which guides the following analytic intuition.

2.2.3 EUV-driven winds

For an EUV-driven wind, the gravitational radius is \(R_{\mathrm {g}}\sim 10\) au for a solar mass star (and scales linearly with stellar mass), such that we are generally in a regime in the outer disc where \(R_{\mathrm {d}}\gg R_{\mathrm {g}}\), and a wind can be launched. If the EUV flux is sufficiently strong, the ionisation front (IF) sits close to the disc surface making the PDR thin (\(x \lesssim 1.5\); see Sect. 2.2.4). In this case, the basic geometry of the system is shown in Fig. 2a. The thermal pressure at the disk surface is determined by the ionisation rate, with the flow proceeding subsonically through the PDR. If we assume isothermal conditions in the PDR, then the density \(n_{\mathrm {I}}\) is constant: \(n_{\mathrm {I}} = n_0 = N_{\mathrm {D}} /xR_{\mathrm {d}}\), where \(N_{\mathrm {D}}\) is the column density of the PDR, which is the column density required to produce an optical depth \(\tau _{\mathrm {FUV}} \sim 1\). This column density is dominated by the base of the wind \(N_{\mathrm {D}}\approx n_0 R_{\mathrm {d}}\), but is dependent on the microphysics and dust content (see Sect. 2.3). For our purposes, we will simply adopt \(N_{\mathrm {D}} \sim 10^{21}\) cm\(^{-2}\) (although see the direct calculations of [40], for example).

Since density in the PDR is constant, the velocity in the flow \(v_{\mathrm {I}} \propto r^{-2}\) to maintain a constant mass flux. In order to conserve mass and momentum flux, the velocity at the ionisation front must be \(v_{\mathrm {IF}} = c_{\mathrm {s,I}}^2/2c_{\mathrm {s,II}} \sim 0.5\) km s\(^{-1}\), where \(c_{\mathrm {s,I}} \approx 3\) km s\(^{-1}\) is the sound speed in the PDR and \(c_{\mathrm {s,II}} \approx 10\) km s\(^{-1}\) beyond the IF. We can write the mass loss rate:

$$\begin{aligned} \dot{M}_{\mathrm {EUV}} = 4\pi {\mathcal {F}} (1+x)^2 R_{\mathrm {d}}^2 n_{\mathrm {I}} m_{\mathrm {I}} \frac{c_{\mathrm {s,I}}^2}{2c_{\mathrm {s,II}}} , \end{aligned}$$
(3)

where \({\mathcal {F}}\) is a geometric factor and \(m_{\mathrm {I}}\) is the mean molecular mass in the PDR, and as before x is the relative thickness of the PDR with respect to the disc radius \(R_{\mathrm {d}}\). Since \(n_{\mathrm {I}} \propto 1/x\) this mass loss rate appears to diverge in the limit of a thin PDR. However, x must satisfy the recombination condition such that:

$$\begin{aligned} \frac{f_{\mathrm {r}} \varPhi }{4\pi d^2} = \int _{R_{\mathrm {IF}}}^\infty \alpha _{\mathrm {B}} n_{\mathrm {II}}^2 \,\mathrm {d}r = R_{\mathrm {IF}}^4\left( \frac{m_{\mathrm {I}} v_{\mathrm{{IF}}}}{m_{\mathrm {II}} v_{\mathrm {II}}}\right) ^2 \int _{R_{\mathrm {IF}}}^\infty \alpha _{\mathrm {B}} \frac{n_{\mathrm{{I}}}^2}{r^4} \,\mathrm {d}r \end{aligned}$$
(4)

where \(\varPhi\) is the EUV counts of the ionising source at distance d, \(f_{\mathrm {r}}\) is the fraction of photons unattenuated by the interstellar medium (ISM), and \(\alpha _{\mathrm {B}} = 2.6 \times 10^{-13}\) cm\(^{3}\) s\(^{-1}\) is the recombination coefficient for hydrogen with temperature \(10^4\) K (e.g. [47]). Substituting \(R_{\mathrm {IF}} = (1+x)R_{\mathrm {d}}\) into Eq. 4 and adopting typical values, we can write a defining equation for x that can be solved numerically. However, in the limit of \(x\ll 1\) we can also simply estimate:

$$\begin{aligned} \dot{M}_{\mathrm {EUV}} = 5.8 \times 10^{-8} \epsilon _{\mathrm {EUV}} {\mathcal {F}} \left( \frac{d}{0.1\, \mathrm {pc}}\right) ^{-1} \left( \frac{ f_{\mathrm {r}} \varPhi }{10^{49} \, \mathrm {s}^{-1}} \right) ^{1/2} \left( \frac{ R_{\mathrm {d}}}{100 \, \mathrm {au}}\right) ^{3/2} \, M_\odot \, \mathrm {yr}^{-1}, \end{aligned}$$
(5)

which is the solution for an ionised globule with no PDR physics [48]. Here, \(\epsilon _{\mathrm {EUV}} \sim 1\) is a correction factor that absorbs uncertainties in the PDR physics. We notice that the EUV-driven wind is super-linearly dependent on \(R_{\mathrm {d}}\), but only scales with the square root of the EUV photon count.

2.2.4 FUV-driven winds

We will here proceed under the assumption that the outer radius and FUV flux are sufficient to produce a supersonic, FUV-heated wind. This means that \(R_{\mathrm {d}}\gg R_{\mathrm {g}}\), where \(R_{\mathrm {g}}\) is the gravitational radius in the PDR. In this case, where the wind mass loss is determined by FUV heating, the neutral wind must launch at (constant) supersonic velocity \(v_{\mathrm {I}} = v_0\gtrsim c_{\mathrm {s,I}}\), where \(v_0\) is the launching velocity from the disc surface with number density \(n_0\), while \(v_{\mathrm {I}}\) is the wind velocity in the PDR as before. To conserve mass, the density in the PDR drops as \(n_{\mathrm {I}} \propto r^{-2}\). The wind travels faster than at the IF (\(v_{\mathrm {IF}} = c_{\mathrm {s,I}}^2/2c_{\mathrm {s,II}}\), as before) and therefore must eventually meet a shock front at radius \(R_{\mathrm {shock}}\). Assuming this shock is isothermal, the density \(n_{\mathrm {I}}\) increases by a factor \(\sim M^2\), where \(M = v_0/c_{\mathrm {s,I}}\) is the Mach number. If the region between the shock and ionisation fronts is isothermal, then \(n_{\mathrm {I}}\) is then constant and \(v_{\mathrm {I}} \propto r^{-2}\) to conserve mass. This geometry is shown in Fig. 2b.

Solving mass and momentum conservation requirements for the flow, we have:

$$\begin{aligned} R_{\mathrm {shock}} = \sqrt{ \frac{v_0}{2c_{\mathrm {s, II}}}} R_{\mathrm {IF}}, \end{aligned}$$
(6)

which immediately puts a minimum distance below which EUV mass loss dominates over FUV. If \(R_{\mathrm {shock}}< R_{\mathrm {d}}\) then the shock front is inside the disc, the flow is subsonic at the base, and we are back in the EUV-driven wind regime. Given that \(\sqrt{{v_0}/{2c_{\mathrm {s, II}}}} \sim 0.4{-}0.6\) (with \(v_0 \sim 3{-}6\) km s\(^{-1}\)), then we require \(R_{\mathrm {IF}}\gtrsim 2.5 R_{\mathrm {d}}\) for an FUV-driven wind to be launched. Coupled with Eq. 4, this gives the minimum distance required for the launching of FUV-driven winds. Conversely, a maximum distance exists from the requirement that the gas can be sufficiently heated to escape the host star, meaning that FUV-driven winds only occur at intermediate separations from a UV source.

Under the assumption that the FUV flux is sufficient to launch a wind and \(R_{\mathrm {IF}}\gtrsim 2.5 R_{\mathrm {d}}\), then the overall mass loss in the flow does not care about what is going outside of the base of the FUV-launched wind. The mass loss rate in the FUV-dominated case is simply:

$$\begin{aligned} \dot{M}_{\mathrm {FUV}} = 4\pi {\mathcal {F}} R_{\mathrm {d}}^2 n_0 v_0 m_{\mathrm {I}} \approx 1.9 \times 10^{-7} \epsilon _{\mathrm {FUV}} {\mathcal {F}} \left( \frac{R_{\mathrm {d}}}{100\,\mathrm {au}} \right) \, M_\odot \, \mathrm {yr}^{-1}, \end{aligned}$$
(7)

where \({\mathcal {F}}\) is the geometric correction and \(m_{\mathrm {I}}\) is the mean molecular mass in the PDR as before. All of the difficult physics is contained within a convenient correction factor \(\epsilon _{\mathrm {FUV}}\). In reality, this expression is only helpful in the limit of a very extended disc, and computing the mass loss rate in general requires a more detailed treatment that we will discuss in Sect. 2.3.

Nonetheless, we gain some insight from the estimate from Eq. 7. First, we see that the mass loss rate is not dependent on the FUV flux \(F_{\mathrm {FUV}}\). In reality, there is some dependence due to the increased temperature in the PDR with increasing \(F_{\mathrm {FUV}}\), but this dependence is weak for \(F_{\mathrm {FUV}}\gtrsim 10^4\, G_0\) [46]. From Eq. 5, we also see that the mass loss rate \(\dot{M}_{\mathrm {FUV}}\) scales less steeply with \(R_{\mathrm {d}}\) than \(\dot{M}_{\mathrm {EUV}}\) does. This means that once a disc has been sufficiently truncated by these external winds, the FUV dictates the mass loss rate. Since the time-scale for this depletion can be short (see Sect. 4.1), we expect FUV flux to dominate mass loss over the disc lifetime for reasonable EUV fluxes.

While this picture is useful to gain some insight into the physics of external photoevaporation, accurately computing the mass loss rate in the wind requires more detailed numerical modelling of the PDR physics. We consider efforts in this direction to date as follows.

2.3 Microphysics of external photoevaporation

One of the biggest challenges surrounding external photoevaporation is that it depends upon a wide array of complicated microphysics. The wind is launched primarily by the FUV radiation field, and determining the temperature in this launching region, which is critical, requires solving photodissociation region (PDR) microphysics, which in itself can consist of many hundreds of species and reactions and complicated thermal processes (e.g. [45, 46]). In particular, as we will discuss below, the line cooling is difficult to estimate in 3D, meaning most PDR codes are limited to 1D (e.g. [49]). The dust and PAH abundance in externally driven wind also play key roles in determining the mass loss rate, but may differ substantially from the abundances in the ISM (e.g. [50, 51]). In addition to the complicated PDR-dynamics, EUV photons establish an ionisation front downstream in the wind, which affects the observational characteristics. Here we introduce some of the key aspects of the microphysics of external photoevaporation in more detail.

2.3.1 The theory of dust grain entrainment in external photoevaporative winds

We begin by considering the dust microphysics of external photoevaporation. First it is necessary to provide some context by discussing briefly how dust evolves in the disc itself. Canonically, in the ISM the dust-to-gas mass ratio is \(10^{-2}\) and grains typically follow a size distribution of the form

$$\begin{aligned} \frac{\mathrm {d} n(a_{\mathrm {s}})}{\mathrm {d} a_{\mathrm {s}}} \propto a_{\mathrm {s}}^{-q} \end{aligned}$$
(8)

[52] with grain sizes spanning \(a_{\mathrm {s}}\sim 10^{-3}-1\) \(\mu\)m [53] and \(q\approx 3.5\). In protoplanetary discs, the dust grains grow to larger sizes, which eventually (when the Stokes number is of order unity) dynamically decouples them from the gas, leading to radial drift inwards to the inner disc. This growth proceeds more quickly in the inner disc (e.g. [54]) and so there is a growth/drift front that proceeds from the inner disc outwards. It is not yet clear how satisfactory our basic models of this process are, particularly in terms of the timescale on which it operates, since if left unhindered by pressure bumps in the disc it quickly results in most of the larger drifting dust being deposited onto the central star (see e.g. [54,55,56]). However, for our purposes the key point is that the abundance of smaller (\(\sim \mu\)m) grains in the disc ends up depleted relative to the ISM due to grain growth.

The nature of dust in the external photoevaporative wind is important for three key reasons

  1. 1.

    The dust in the wind sets the extinction in the wind and hence has a significant impact on the mass loss rate.

  2. 2.

    The extraction of dust in the wind could have implications for the mass reservoir in solids for terrestrial planet formation and/or the cores of gas giants.

  3. 3.

    The entrainment of dust in winds could provide observational diagnostics of the external photoevaporation process (we will discuss this further in Sect. 3).

and so the key questions are what size, and how much, dust is entrained in an external photoevaporative wind. This problem was addressed in the semi-analytic work of [50]. They solved the flow structure semi-analytically (we discuss semi-analytic modelling of external photoevaporative winds further in Sect. 2.4) and calculated the maximum entrained grain size. The efficiency of dust entrainment in the wind is dependent on the balance of the Epstein [57] drag exerted on the dust grain with density \(\rho _{\mathrm {s}}\) by the outflowing gas of velocity \(v_{\mathrm {th}}\) versus the gravity from the host star. Thus, the condition for a dust grain to be lost to external photoevaporation is [50]:

$$\begin{aligned} a_{\mathrm {s}} < a_{\mathrm {ent}} =\frac{1}{4\pi {\mathcal {F}} } \frac{v_{\mathrm {th}}\dot{M}_{\mathrm {ext}}}{G M_* \rho _{\mathrm {s}}} , \end{aligned}$$
(9)

where \(4\pi {\mathcal {F}}\) is the solid angle subtended by the wind (see [58]).

The main outcome of the above is that only small grains are entrained in the wind and the mean cross section is reduced. Therefore, when grain growth proceeds to the disc outer edge, the dust-to-gas mass ratio, mean cross section, and hence, extinction in the wind drops substantially. This makes external photoevaporation more effective than previously considered when the dust in the wind was treated as ISM-like. This lower cross section in the wind is now accounted for in numerical models of external photoevaporation, assuming some constant low value (e.g. the FRIED grid of mass loss rates [59], discussed more in 2.5). However, what is still missed in models is that the cross section in the wind is actually a function of the mass loss rate [50] and so needs to be solved iteratively with the dynamics.

[60] also recently introduced dynamically decoupled dust into 1D isothermal models of externally irradiated discs (discussed more in Sect. 2.5.1), finding that it does indeed lead to a radial decrease in the maximum grain size, which could be searched for observationally. This gradient in grain sizes has been observed by Throop et al. [61] and Miotello et al. [62], who studied the large proplyd 114-426 on the near side of the ONC, calculating the attenuation of the background radiation through translucent parts of the disc. They found a maximum grain size decreasing from 0.7 \(\mu\)m to 0.2 \(\mu\)m, moving away from the disc outer edge over a distance of around 250 au, consistent with [60].

2.3.2 Photodissociation region physics for external photoevaporation

The FUV-excited photodissociation region (PDR) microphysics determines the composition, temperature and therefore the dynamics of the inner parts of external photoevaporative winds. As discussed above, this FUV/PDR part of the inner wind can determine the mass loss rate from the disc. This is not a review on PDR physics (for further information, see, e.g. [45, 46, 63]) but given its importance for setting the temperature, and therefore, the dynamics, we provide a brief overview of some relevant processes.

We focus primarily on the main heating and cooling contributions. These are summarised as a function of extinction for an external FUV field of 300 G\(_0\) in Fig. 3, which is taken from [50]. Note that, as we will discuss below, the exact form of these plots depends on the FUV field strength and the assumed composition, e.g. the metallicity, dust grain properties and polycyclic aromatic hydrocarbon (PAH) abundance.

Fig. 3
figure 3

A summary of the key heating and cooling mechanisms in a medium irradiated by a \(300\,G_0\) FUV radiation field. PAH-driven photoelectric heating dominates until high \(A_V\), where cosmic rays take over. The key cooling mechanism in the wind is the escape of line photons. From [50]

The heating mechanism that is anticipated to be most important for external photoevaporation is photoelectric heating (see the left-hand panel of Fig. 3) that occurs when PAHs lose electrons following photon absorption, increasing the gas kinetic energy [63]. The impact this can have on the mass loss rate is illustrated in Fig. 4, which shows the results of numerical models of an externally photoevaporating 100 au disc around a 1 M\(_\odot\) star in a \(1000 \, G_0\) FUV environment as a function of metallicity. Each coloured set of points connected by a line represents a different PAH-to-dust ratio. Reducing the PAH-to-dust ratio has a much larger impact on the mass loss rate than changing the overall metallicity. These models are previously unpublished extensions of the 1D PDR-dynamical calculations of [59], which are discussed further in 2.5. When the metallicity is reduced, the PAH abundance and heating are also lowered, but so is the line cooling. Changes in metallicity therefore only lead to relatively small changes to the mass loss rate as the heating and cooling changes compensate. Conversely, changing only the PAH-to-dust ratio can lead to dramatic changes in the mass loss rate.

A key issue for the study of external photoevaporation is that the PAH abundance in the outer parts of discs and in winds is very poorly constrained. For the proplyd HST 10 in the ONC, [51] inferred a PAH abundance relative to gas around a factor 50 lower than the ISM and a factor 90 lower than NGC 7023 [64]. Note that the models in Fig. 4 use a dust-to-gas mass ratio of \(3\times 10^{-4}\) so an \(f_PAH\) of 0.1 in Fig. 4 corresponds to a PAH-to-dust ratio of 1/330. PAH detections around T Tauri stars are generally relatively rare [65, 66], which leads us to expect that the PAH abundance is depleted in discs irrespective of external photoevaporation. This lower PAH abundance would mean less heating due to external photoevaporation, resulting in lower external photoevaporative mass loss rates. Conversely, [67] demonstrated that PAH emission from the inner disc could be suppressed when PAHs aggregate into clumps, which also crucially would not suppress the heating contribution from PAHs (Lange, priv. comm.). However, it is unclear if that same model for PAH clustering applies at larger radii in the disc, let alone in the wind itself and so this is to be addressed in future work (Lange, priv. comm.).

Fig. 4
figure 4

External photoevaporative mass loss rate as a function of metallicity (\(Z/Z_\odot\)) for a 100 au disc around a 1 M\(_\odot\) star irradiated by a \(1000 \, G_0\) radiation field. Each coloured line represents a different value of the base PAH-to-dust mass ratio scaling \(f_{\mathrm {PAH}}\). These are extensions of the FRIED PDR-dynamical models of [59]. When the overall metallicity is scaled, there are changes to both the heating and cooling contributions that broadly cancel out. Conversely, varying the PAH-to-dust ratio (which is very uncertain) can lead to large changes in the mass loss rate. Note that these calculations have a floor value of \(10^{-11}\) M\(_\odot\) yr\(^{-1}\)

Given its potential role as the dominant heating mechanism, determining the PAH abundance in the outer regions of discs is vital for understanding the magnitude of mass loss rates and so should be considered a top priority in the study of external photoevaporation. The James Webb Space Telescope (JWST) should be able to constrain abundances by searching for features such as the 15-20 \(\mu\)m emission lines (e.g. [68,69,70,71,72,73]). [74] also demonstrated with synthetic observations that the upcoming Twinkle [75] and Ariel [76, 77] missions should be able to detect the PAH 3.3 \(\mu\)m feature towards discs, at least out to 140 pc. Even if detections with Twinkle/Ariel would not succeed in high UV environments because of the larger distances to those targets, containing PAH abundances of the outer disc regions in lower UV nearby regions would also provide valuable constraints. These will be an important step for calibrating models, refining the mass loss rate estimates and hence our understanding of external photoevaporation.

Often the dominant cooling in PDRs is the escape of line photons from species such as CO, C, O and C\(^+\). Evaluating this is the most challenging component of PDR calculations, since to estimate the degree of line cooling, line radiative transfer in principle needs to sample all directions (\(4\pi\) steradians) from every single point in the calculation. For this reason, most PDR studies to date (even without dynamics) have been 1D, where it is assumed that exciting UV radiation and cooling radiation can only propagate along a single path, with all other trajectories infinitely optically thick (e.g. [49, 78,79,80]). Most dynamical models of external photoevaporation with detailed PDR microphysics have also therefore been 1D. For example, [58] used the [79] 1D PDR code to pre-tabulate PDR temperatures as inputs for 1D semi-analytic dynamical models (we will discuss these in more detail in Sect. 2.4). Note that 2D models of other features of discs have circumvented this issue by assuming a dominant cooling direction, for example vertically through the disc (e.g. [81]), or radially in the case of internal photoevaporation calculations [82]. This approach is not applicable in multidimensional simulations of an externally driven wind, where there is no obvious or universally applicable dominant cooling direction.

The 3d-pdr code developed by [83] and based on the ucl-pdr code [78] was the first code (and to our knowledge remains the only code) able to treat PDR models in 3D. It utilises a healpix scheme [84] to estimate the line cooling in 3D without assuming preferred escape directions. healpix breaks the sky into samples of equal solid angle at various levels of refinement. For applications to external photoevaporation, 3d-pdr was coupled with the Monte Carlo radiative transfer and hydrodynamics code torus [85] in the torus-3dpdr code [86], making 2D and 3D calculations possible in principle, which we will discuss more in Sect. 2.5. However, the computational expense of doing 3D ray tracing from every cell in a simulation iteratively with a hydrodynamics calculation is prohibitively expensive. Finding ways to emulate the correct temperature without solving the full PDR chemistry may offer a way to alleviate this problem (e.g. [87]).

2.4 1D Semi-analytic models of the external photoevaporative wind flow structure

In Sect. 2.3, we discussed the importance of PDR microphysics for determining the temperature structure and hence the flow structure of externally irradiated discs. We also noted that PDR calculations are computationally expensive and are usually limited to 1D geometries. Until recently, calculations of the mass loss rate that utilise full PDR physics have also been confined to 1D and solved semi-analytically. Here we briefly review those approaches.

First we describe the 1D approach to models of external photoevaporation and the justification for such a geometry. 1D models essentially follow the structure radially outwards along the disc mid-plane into the wind, as illustrated in Fig. 5. The grid is spherical, but the assumption is that the flow only applies over the solid angle subtended by the disc outer edge at \(R_{\mathrm{{d}}}\). That fraction of \(4\pi\) steradians is

$$\begin{aligned} {\mathcal {F}} = \frac{H_{\mathrm{{d}}}}{\sqrt{H_{\mathrm{{d}}}^2+R_{\mathrm{{d}}}^2}} \end{aligned}$$
(10)

for a disc scale height \(H_{\mathrm {d}}\) (again see Fig. 5). The mass loss rate at a point in the flow at R with velocity \(\dot{R}\) and density \(\rho\) is then

$$\begin{aligned} {\dot{M}} = 4\pi R^2{\mathcal {F}} \rho \dot{R}. \end{aligned}$$
(11)
Fig. 5
figure 5

A schematic of the 1D semi-analytic model structure and how it is used to estimate total mass loss rate estimates. The flow solution is solved along the disc mid-plane with appropriate boundary conditions, e.g. at the disc outer edge and at some critical point in the wind. This mid-plane flow is then assumed to apply over the entire solid angle subtended by the disc outer edge

The 1D geometry is justified based on the expectation that it is there that material is predominantly lost from the disc outer edge. That expectation arises because:

  1. 1.

    Material towards the disc outer edge is the least gravitationally bound.

  2. 2.

    The vertical scale height is much smaller than the radial, which results in a higher density at the radial sonic point than the vertical one [58].

This is demonstrated analytically in the case of compact discs in the appendix of [58], who show that the ratio of mass lost from the disc surface to disc edge is

$$\begin{aligned} \frac{\dot{M}_{\mathrm {surface}}}{\dot{M}_{\mathrm {edge}}} \approx \left( \frac{R_{\mathrm{{d}}}}{R_{\mathrm{{g}}}}\right) ^{1/2}. \end{aligned}$$
(12)

where \(\dot{M}_{\mathrm {surface}}\) and \(\dot{M}_{\mathrm {edge}}\) are the mass loss rates from the disc upper layers and mass loss rates from the outer edge, respectively. As before, \(R_{\mathrm{{d}}}\) and \(R_{\mathrm{{g}}}\) are the disc outer radius and gravitational radius, respectively. That is for larger, more strongly heated discs there is a more significant contribution from the disc surface. This has also been tested and validated in 2D radiation hydrodynamic simulations by [88] (that we will discuss fully in Sect. 2.5) who showed that, at least in a \(10^{3}\, G_0\) environment, the majority of the mass loss comes from the disc outer edge and the rest from the outer 20 percent of the disc surface. Mass loss rates in 2D and analogous 1D models were also determined to be similar to within a factor two, with the 1D mass loss rates being the lower values. Mass loss rates computed in one dimension are therefore expected to be somewhat conservative but reasonable approximations.

[58] took a semi-analytic approach to solving for the flow structure by using pre-tabulated PDR temperatures from the code of [79] and using those in the flow equations. They found that the flow structure is analogous to a [89] wind, but non-isothermal and with centrifugal effects. At each point in the flow, the pre-tabulated PDR temperatures are interpolated as a function of local density incident FUV and extinction. The boundary conditions used were the conditions at the disc outer edge and the sonic point in the flow. They demonstrated both that FUV-driven winds are dominant for setting the mass loss rate and that winds could be driven interior to the gravitational radius, down to \(\sim 0.1-0.2\,R_{\mathrm {g}}\) (see also, e.g. [90, 91]).

[50] took a similar approach to 1D semi-analytic models with pre-tabulated PDR temperatures from [83]. As already discussed above, their main focus was on dust entrainment and the impact of grain growth in the disc on the dust properties in the wind. They found that the entrainment of only small grains, coupled with grain growth in the disc, reduces the extinction in the wind and can enhance the mass loss rate. In addition, they used a different approach to the outer boundary condition, finding a critical point in the modified Parker wind solution and taking into account deviations from isothermality at that point. They then integrated from that critical point, inwards to the disk. Thanks to this different approach [50] were able to compute solutions over a wider parameter space than before, particularly down to low FUV field strengths.

Semi-analytic models have offered a powerful and efficient tool for estimating mass loss rates in different regimes. However, there are still regions of parameter space where solutions are not possible and semi-analytic models are still only limited to 1D. To alleviate those issues, we require radiation hydrodynamic simulations.

2.5 Radiation hydrodynamic models of external photoevaporation

Above we discussed semi-analytic calculations of the external wind flow structure. Those have the advantage that they are quick to calculate. However, they are limited by being restricted to 1D solutions and by solutions not always being calculable. This leaves a demand for radiation hydrodynamic calculations capable of solving the necessary radiative transfer and PDR temperature structure in conjunction with hydrodynamics. Such calculations can solve for the flow structure in 2D/3D and in any scenario.

The radiation hydrodynamics of external photoevaporation is one of the more challenging problems in numerical astrophysics because the key heating responsible for launching the wind is described by PDR physics. That is, we are required to iteratively solve a chemical network that is sensitive to the temperature with a thermal balance that is sensitive to the chemistry. To make matters worse, the cooling in PDRs is non-local, being dependent on the escape probability of line photons into \(4\,\pi\) steradians from any point in the simulation (see Sect. 2.3.2). In other scenarios, this does not cause significant issues if there is a clear dominant escape direction. For example, within a protoplanetary disc the main cooling can quite reasonably be assumed to occur vertically through the disc, since other trajectories have a longer path length (and column) through the disc (“1+1D” models, e.g. [81, 92]). Similarly, for internal winds there are models with radiation hydrodynamics and PDR chemistry, where the line cooling is evaluated along single paths radially on a spherical grid [82, 93,94,95]. In the complicated structure of an external photoevaporative wind, however, this sort of geometric argument cannot be applied and so multiple samplings of the sky (\(4\pi\) steradians) are required from every point in a simulation.

Although 3D cooling is ideally required, simulations have been performed using approximations to the cooling in high UV radiation fields in particular, where the PDR is small. For example, [96] ran 2D axisymmetric simulations of discs irradiated face-on. In their calculations, the optical depth and cooling are estimated using a single ray from the cell to the irradiating UV source (this same path is used for calculating the exciting UV and cooling radiation). They also employed a more simple PDR microphysics model compared to work at the time such as [33], which enabled the move from 1D to 2D-axisymmetry. [96] studied the mass loss of proplyds as well as the observational characteristics using intensity maps derived from their dynamical models, some of which are illustrated in Fig. 6. They found, in the first geometrically realistic EUV+FUV irradiation models, that rapid disc dispersal gives morphologies in various lines similar to those observed in the ONC.

The torus-3dpdr code [86] is a key recent development in the direct radiation hydrodynamic modelling of external photoevaporation. It was constructed by merging components of the first fully 3D photodissociation region code 3d-pdr [83] with the torus Monte Carlo radiative transfer and hydrodynamics code [85]. 3d-pdr (and hence torus-3dpdr) address the 3D line cooling issue using a healpix scheme, which breaks the sky into regions of equal solid angle. torus-3dpdr has been used to run a range of 1D studies of external photoevaporation. It has been shown to be consistent with semi-analytic calculations [97]. It was used to study external photoevaporation in the case of very low mass stars, with a focus on Trappist-1 [98]. The approach was to run a grid of models to provide the mass loss rate as a function of the UV field strength, disc mass/radius for a 0.08 M\(_\odot\) star and interpolate over that grid in conjunction with a disc viscous evolutionary code based on that of Clarke [99] to evolve the disc. The usefulness of such a grid led to the FRIED (FUV Radiation Induced Evaporation of Discs) grid of mass loss rates that has since been employed in a wide range of disc evolutionary calculations by various groups (e.g. [100,101,102], others). Note, given the discussion on the importance of PAHs above that the FRIED models use a dust-to-gas ratio a factor 33 lower than the canonical \(10^{-2}\) for the ISM and a PAH-to-dust mass ratio a further factor of 10 lower than ISM-like. This is conservative (i.e. a PAH-to-gas ratio of 1/330 that in the ISM) compared to the factor 50 or so PAH depletion measured by [51]. So models predict that the PDR heating is still capable of driving significant mass loss, even when the PAH abundance is heavily depleted compared to the ISM.

These applications were all still 1D though and so there is a growing theoretical framework that is based on that geometric simplification. [88] ran 2D-axisymmetric external photoevaporation calculations (with 3D line cooling, utilising the symmetry of the problem) and found that 1D calculations are, if anything, conservative since the 2D mass loss rates were slightly higher. True 3D calculations with 3D line cooling and the disc not irradiated isotropically or face-on, are yet to be performed. Though in principle torus-3dpdr is capable of this, in practice the 3D ray tracing of the healpix scheme makes such calculations computationally expensive.

Fig. 6
figure 6

A gallery of intensity maps (\(\log _{10}\)erg s\(^{-1}\) cm\(^{-2}\) strad\(^{-1}\)) resulting from the 2D-axisymmetric radiation hydrodynamic simulations of proplyds by [96]. The simulations utilised a simplified microphysics approach which enabled them to model proplyds in 2D-axisymmetry, with UV radiation incident from the top. These multi-wavelength synthetic line emission maps were compared and found to be consistent with the properties of ONC proplyds

2.5.1 The disc-wind transition

Most of the models discussed so far impose a disc as a boundary condition from which a wind solution is derived. In reality, the wind is not launched from some arbitrarily imposed point in an irradiated disc. [60] recently implemented a model without an imposed inner boundary and a smooth transition from disc to wind using a slim disc approach [103]. In this first approach, they assume an isothermal disc, but demonstrate that while the transition from disc to wind is narrow, it is not negligibly thin. They also introduced dust of different sizes into their model, which predicts a radial gradient in grain size in the outer disc/inner wind. Although the fixed inner boundary models are valid for computing steady state mass loss rates for disc evolutionary models, a worthwhile future development will be to include detailed microphysics in a slim-disc approach like that of [60].

2.6 Summary and open questions for the theory of externally driven disc winds

  1. 1.

    Numerical models of external photoevaporation require some of the most challenging radiation hydrodynamics models in astrophysics. This is primarily because it is necessary to include 3D PDR microphysics, including line cooling with no obvious dominant cooling direction.

  2. 2.

    Limited by the above, 1D models of external photoevaporation are now well established and are used to estimate disc mass loss rates. But 2D and 3D simulations are still limited.

Some of the many open questions and necessary improvements to current models are

  1. 1.

    What is the PAH abundance in external photoevaporative winds and the outer regions of discs? This is key to setting the wind temperature and mass loss rate.

  2. 2.

    Including mass loss rate dependent dust-to-gas mass ratios and maximum grain sizes (and hence extinction) in numerical models of external photoevaporation. At present a single representative cross section is assumed, irrespective of the mass loss rate.

  3. 3.

    Can accurate temperatures from PDR microphysics be computed at vastly reduced computational expense (e.g. via emulators)?

  4. 4.

    What is the interplay between internal and external winds?

  5. 5.

    3D simulations of external photoevaporation with full PDR-dynamics

  6. 6.

    Non-isothermal slim-disc models of externally photoevaporating discs.

3 Observational properties of externally photoevaporating discs

Here we discuss observations to date of individual externally photoevaporating discs. We discuss the diversity in their properties such as UV environment and age, and summarise key diagnostics.

3.1 Defining the term proplyd

The term “proplyd” was originally used to describe any disc directly imaged in Orion with HST in the mid-90s (e.g. [23, 33]) as a portmanteau of “protoplanetary disc”. Since then, use of the term has adapted to only refer to cometary objects resulting from the external photoevaporation of compact systems. However, this use of the term is ambiguous, since a cometary morphology can result from both externally irradiated protoplanetary discs and externally irradiated globules, which may be host no embedded star/disc, such as many of the compact globulettes as illustrated in Fig. 7 [104, 105]. We therefore propose to define a proplyd as follows

Proplyd:

A circumstellar disc with an externally driven photoevaporative wind composed of a photodissociation region and an exterior ionisation front with a cometary morphology.

We have chosen this definition such that it makes no distinction as to whether EUV or FUV photons drive the wind, but specifically define that for an object to be a proplyd then the wind must be launched from a circumstellar disc. In the absence of a disc, it is a globule or globulette (also sometimes referred to as an evaporating gas globules, or EGG, e.g. [106]). To further clarify, a globule or globulette with an embedded YSO (identified through a jet for example) with cometary morphology would also not be defined as a proplyd since it is the ambient material being stripped rather than the disc.

Fig. 7
figure 7

The left-hand panel is a gallery of proplyds in the ONC – evaporating discs around YSO’s (Credit: NASA/ESA and L. Ricci). The right-hand panel is a gallery of globulettes in Carina from [105]. Globulettes have radii from hundreds to thousands of au and can also take on a cometary morphology when externally irradiated, though in many cases do not contain any YSOs, which by our definition would mean that they are not proplyds (note that [105] never referred to them as such, rather using the term globulette)

3.2 Where and what kind of proplyds have been detected?

Proplyds were originally discovered in the ONC and there are now over 150 known in the region (e.g. [23, 33, 107, 108]). The ONC is around 1-3 Myr old with the primary UV source the O6V star \(\theta ^1\)C. In general, the surface brightness of proplyds if the local EUV exposure is lower than proplyds can be harder to detect. For example, the surface brightness of \(\mathrm {H}\alpha\) is [109]:

$$\begin{aligned} \langle S(\mathrm {H}\alpha ) \rangle \approx 7 \times 10^{11} \frac{\alpha _{\mathrm {H}\alpha }^{\mathrm {eff}}}{\alpha _B}\frac{\varPhi }{10^{49} \, \mathrm {s}^{-1}} \left( \frac{d}{0.1 \, \mathrm {pc}} \right) ^{-2} \, \mathrm {s}^{-1} \, \mathrm {cm}^{-2} \, \mathrm {s.r}^{-1} \end{aligned}$$
(13)

where \(\alpha _{\mathrm {H}\alpha }^{\mathrm {eff}} = 1.2\times 10^{-13}\) cm\(^{-3}\) s\(^{-1}\) and \(\alpha _B = 2.6\times 10^{-13}\) cm\(^{-3}\) s\(^{-1}\) at temperature \(T=10^4\) K are hydrogen recombination coefficients (e.g. [47]). Nonetheless, in recent years there have been important detections of proplyds in lower EUV environments. Proplyds can be found in the ONC at separations of up to \(\sim 1\) pc from \(\theta ^1\)C [20, 110, 111]. Meanwhile, [112] found proplyds in the vicinity of the B1V star 42 Ori in NGC 1977, demonstrating that B stars can drive external photoevaporative winds. [113] also presented the discovery of proplyds in NGC 2024. In that region it appears that both an O8V and B star are driving external photoevaporation. The main significance of proplyds there is the \(\sim 0.2-0.5\) Myr age of a subset of the region where proplyds have been discovered. This is important since it implies that external photoevaporation can even be in competition with our earliest stage evidence for planet formation [114, 115].

In these regions, proplyds have been detected with host star masses from around solar down to almost planetary masses (\(<15\,M_{\text {jup}}\) [112, 116,117,118]). The UV fields incident upon these proplyds ranges from \(>10^5 \, G_0\) down to possibly around \(100 \, G_0\) [62]. Mass loss rates are estimated to regularly be greater than \(10^{-7}\) M\(_\odot\) yr\(^{-1}\) and sometimes greater than \(10^{-6}\) M\(_\odot\) yr\(^{-1}\) (e.g. [34, 35, 113, 119]). Examples of binary proplyds have also been discovered [120]. For sufficiently close separation binaries, the disc winds merge to form a so-called interproplyd shell, which was studied by [121]. These nearby regions, the ONC, NGC 1977 and NGC 2024, show the clearest evidence for external photoevaporation.

Due to resolution and sensitivity issues, unambiguous evidence for external photoevaporation is more difficult to obtain in more distant star-forming regions than the \(D \sim 400\) pc of Orion. [122] identified candidate proplyds in Trumpler 14, at a distance of \(D\sim 2.8\,\)kpc, and [106] subsequently detected discs towards those candidates with ALMA. Although many of those candidates are large evaporating globules (in some cases, with embedded discs detected), some are much smaller and so could be bona fide proplyds.

There are other regions where “proplyd-like” objects have been discovered, including Cygnus OB2 [123, 124], W5 [125], NGC 3603 [126], NGC 2244, IC 1396 and NGC 2264 [127]. However, our evaluation of those systems so far is that they are all much larger than ONC proplyds and likely evaporating globules. Given the high UV environments of those regions and the identification of evaporating globules, we do expect external photoevaporation to be significant in the region. However, it remains unclear for many of these objects whether the winds are launched from an (embedded) star-disc system. Future higher-resolution observations (e.g. extremely large class telescopes should resolve ONC-like proplyds out to Carina) and/or new diagnostics of external photoevaporation that do not require spatially resolving the proplyd (for example line ratios, or tracers that show a wind definitely emanates from a disc) are required in these regions.

We provide a further discussion of proplyd demographics and particularly the demographics of discs in irradiated star-forming regions in Sect. 4.3.

3.3 Estimating the mass loss rates from proplyds

As discussed above, proplyds have a cometary morphology with a “cusp” pointing towards the UV source responsible for driving the wind and an elongated tail on the far side, pointing away from the UV source. The leading hemisphere that is directly irradiated by the UV source is referred to as the cometary cusp. On the far side of the proplyd is a trailing cometary tail.

The extent of the cometary cusp is set by the point beyond which all of the incident ionising flux is required to keep the gas ionised under ionisation equilibrium. A higher mass loss rate and hence denser flow increase the recombination rate in the wind and would move the ionisation front to larger radii. Conversely, increasing the ionising flux will reduce the ionisation front to smaller radii. As a result, the ionisation front radius \(R_{\textsc {if}}\) (i.e. the radius of the cometary cusp) is related to the ionising flux incident upon the proplyd and the mass loss rate \(\dot{M}_{\mathrm {ext}}\) [33]. This is independent of the actual wind driving mechanism, being enforced simply by photoionisation equilibrium downstream of the launching region of the flow. This provides a means to estimate the mass loss rate from the disc:

$$\begin{aligned} \left( \frac{\dot{M}_{\mathrm {ext}}}{10^{-8}\,M_\odot \,\mathrm {yr}^{-1}}\right) = \left( \frac{1}{1200}\right) ^{3/2}\left( \frac{R_{\textsc {if}}}{\text {au}}\right) ^{3/2}\left( \frac{d}{\text {pc}}\right) ^{-1}\left( \frac{\varPhi }{10^{45}\,\mathrm {s}^{-1}}\right) ^{1/2} \end{aligned}$$
(14)

where \(\varPhi\) is the ionising photons per second emitted by the source at distance d responsible for setting the ionisation front. This has been applied to estimating mass loss rates in NGC 2024 [113] and NGC 1977 [117]. Note that this neglects extinction between the UV source and the proplyd and could underestimate the true separation between the UV source and proplyd. Both of these effects would reduce the true ionising flux incident upon the proplyd and so Eq. 14 provides an upper limit on the mass-loss rate.

Generally, the mass loss rate could alternatively be inferred if one knows the density and velocity through a surface in the wind enclosing the disc. However, the ionisation front is very sharp, making it an ideal surface through which to estimate the mass loss rate.

Other more sophisticated model fits to proplyds have been made, such as [119], who use photoionisation, hydrodynamics and radiative transfer calculations to model the proplyd LV 2, but Eq. 14 provides a quick estimate that gives mass loss rates comparable to those more complex estimates.

3.4 Multi-wavelength anatomy of proplyds

Here we provide a brief overview of some key observational tracers of proplyds. We also highlight possible alternative tracers that might prove useful to identify external photoevaporation when proplyds cannot be inferred based on their cometary morphology (i.e. in weaker UV environments and distant massive clusters). A schematic of the overall anatomy and location of various tracers of a proplyd are shown in Fig. 8.

Fig. 8
figure 8

The anatomy of a proplyd, highlighting the observables in different parts of the system. The three main zones are the disc, PDR and ionised gas

3.4.1 Ionised gas tracers in absorption and emission

Proplyds are observed in two ways using ionised gas emission lines. Proplyds are typically found in H ii regions, which are associated with widespread ionised emission lines. Proplyds on the side of the H ii region closest to the observer can therefore manifest in the absorption of those ionised gas emission lines. In addition, the ionisation front at the rim of the proplyd cusp is an additional source of ionised emission lines that can be directly detected. Ionised gas tracers in emission probe the region close to or outside of the Hi-H ii ionisation front (e.g. [35]). These ionised gas tracers detected in emission are valuable for estimating mass loss rates using the procedure discussed in 3.3. Prominent ionised gas emission lines include H\(\,\alpha\), Paschen\(\,\alpha\), [O iii], [Ne ii] and [N ii]. Combinations of these photoionised gas tracers have also been used to constrain models of proplyds. For example by [119], who compared simulated images and observations in H\(\,\alpha\), [N ii] (6583Å) and [O iii] (5007Å) to model the proplyd LV 2.

3.4.2 Disc tracers

CO Thanks to its brightness, CO is one of the most common line observations towards protoplanetary discs. [50] pointed out that because the specific angular momentum in the wind scales as \(R^{-2}\) rather than the Keplerian \(R^{-3/2}\) that non-Keplerian rotation is a signature of external photoevaporation. Although this deviation grows with distance into the wind, it does not get a chance to do so to detectable levels before CO is dissociated [97, 128]. CO is therefore not expected to kinematically provide a good probe of external winds for proplyds (though it may be useful for more extended discs with slow external winds in low UV environments, as we will discuss in 3.5). However, this does not preclude CO line ratios or intensities showing evidence of external heating in spatially resolved observations.

Dust continuum External photoevaporation influences the dust by i) entraining small dust grains in the wind [50] and ii) by heating the dust in the disc. Directly observing evidence for grain entrainment would provide a key test of theoretical models both of external photoevaporation and dust-gas dynamics. In addition to the prediction that small grains are entrained, Owen and Altaf [60] predict a radial gradient in the grain size distribution.

Evidence for such a radial gradient in grain size was inferred in the ONC 114-426 disc by [62] (see also [61]). This disc is the largest in the central ONC, on the near side of the H ii region. Although the UV field incident upon it is expected to be of order \(10^2\, G_0\), it is clearly evaporating, with the wind resulting in an extended diffuse foot-like structure (though no clear cometary proplyd morphology). They mapped the absorption properties of the translucent outer parts of the disc, finding evidence for a radially decreasing maximum grain size, which would be consistent with theoretical expectations. Revisiting 114-426 to obtain further constraints on grain entrainment would be valuable, as well as searching for the phenomenon in other systems. JWST will offer the capability to similarly study the dust in the outer parts of discs in silhouette in the ONC, comparing JWST Paschen \(\alpha\) absorption with HST H\(\alpha\) absorption (as part of PID: GTO 1256 [129]).

The dust in discs is also influenced by the radiative heating from the environment. If a proplyd is sufficiently close to a very luminous external source the grain heating can be comparable to, or in some parts of the disc exceed, the heating from the disc’s central star. If this is not accounted for when estimating the dust mass in a proplyd (i.e. if one assumes some constant characteristic temperature, typically \(T=20\,\)K) the mass ends up increasingly overestimated in closer proximity to the luminous external source [130] and so may suppress spatial gradients in disc masses at distances within around 0.1 pc of an O star like \(\theta ^1\)C (e.g. [6, 131]).

3.4.3 Photodissociation region tracers

Photodissociation region (PDR) tracers are valuable because they trace the inner regions of the flow. In particular, as discussed in 2.3.2, it is the FUV/dissociation region that determines the mass loss rate. PDR tracers in the wind are also valuable because they are what we will rely on to identify externally photoevaporating discs that are non-proplyds (discussed further in Sect. 3.5). PDR tracers of external photoevaporation have received a lot of attention in recent years and so are explored somewhat more thoroughly than photoionised gas tracers here.

[OI] 6300Å [132] modelled the [OI] 6300Å  emission from proplyds motivated by [133] observations of the ONC proplyd 1822-413. They found that in the case of external photoevaporation the line is emitted following the photodissociation of OH, with the resulting excited oxygen having a roughly 50 per cent chance of de-exciting through the emission of a [OI] 6300Å  photon. For the density/temperature believed to be representative in the wind, this dissociation is expected to dominate over OH reacting out by some other pathway. This model approximately reproduced the flux of 1822-413. Ballabio et al. (in preparation) generalised this, studying how the [OI] 6300Å  line strength varies with UV field strength and star disc parameters. They found that the line is a poor kinematic tracer of external winds because the velocity is too low to distinguish it from [OI] emission from internal winds (though spectro-astrometry, e.g. [134], and future instrumentation may solve this issue by spatially distinguishing the internal and external winds). However, the [OI] luminosity increases significantly with UV field strength. The ratio of [OI] 6300Å  luminosity to accretion luminosity is therefore expected to be unusually high in strong UV environments. The ratio usually has a spread of about \(\pm 1\) orders of magnitude in low UV environments [135], but the models predict the ratio increases by an additional order of magnitude above the usual upper limit in a \(10^5\)G\(_0\) environment. This could make [OI] a valuable tool for identifying external photoevaporation at large distances where proplyds cannot be spatially resolved.

There are some observational challenges associated with utilising this diagnostic though. One is that the targeted YSO’s emission has to be distinguished from emission from the star-forming region. For example, [OI] emission from a background PDR. Furthermore, estimating the accretion luminosity of proplyds appears to have only been attempted in a handful of cases, with tracers like H \(\alpha\) also possibly originating from the proplyd wind.

C I As we discussed above, the most commonly used disc gas tracer, CO, is dissociated in the wind. This means that it is ineffective for detecting deviations from Keplerian rotation that are expected in external photoevaporative winds. C I primarily resides in a layer outside of the CO abundant zone until it is ionised at larger distances in the wind, so could therefore trace the deviation from Keplerian rotation. Haworth and Owen [128] proposed that C I offers a means of probing the inner wind kinematically. A key utility of this would be its possible use as an identifier of winds in systems where there is no obvious proplyd morphology. Haworth et al. [117] used APEX to try and identify the C I 1-0 line in NGC 1977 proplyds, which are known evaporating discs in an intermediate (\(F_{\mathrm {FUV}} \sim 3000 \, G_0\)) UV environment. However, they obtained no detections, which they explain in terms of those proplyd discs being heavily depleted of mass. However, an alternative explanation would be if the disc were just depleted in carbon. Distinguishing these requires independent constraints on the mass in those discs. Overall, the utility of C I remains to be proven and should be tested on higher mass evaporating discs. Based on the expected flux levels (see also [136]), it seems unlikely though that C I will be suitable for mass surveys searching for the more subtle externally driven winds when there is no ionisation front. Though it could be used for targeted studies of extended systems that are suspected to be intermediate UV environments.

Far-infrared lines: C II 158 \(\mu\)m and [OI] 63\(\mu\)m The [C II] 158 \(\mu\)m and [O i] 63 \(\mu\)m lines are both bright tracers of the PDR. They have both been observed with Herschel [137] towards a small number of proplyds by [138], who compared the line fluxes with uniform density 1D PDR models with the Meudon code [80] to constrain parameters such as the mean flow density, which is supported by ALMA observations. They suggested that the proplyd PDR self-regulates to maintain the H-H\(_2\) transition close to the disc surface and maintain a flow at \(\sim 1000\) K in the supercritical regime (\(R_{\mathrm {d}} > R_{\mathrm {g}}\)). Their models also pointed towards a number of heating contributions being comparably (or more) important than PAH heating. However, those calculations also assumed uniform density and ISM-like dust, whereas we now know it is depleted in the wind. Overall, this highlights the need for further detail in PDR-dynamical models.

Clearly these PDR tracers do have enormous utility for understanding the conditions in the inner part of external photoevaporative winds. The main limitation to using these far-infrared lines now is the lack of facilities to observe them, with Herschel out of commission and SOFIA (e.g. [139]) due to end soon. We are unaware of any short-term concepts to alleviate this, but in the longer term there are at least two relevant far-infrared probe class mission concepts being prepared, FIRSST and PRIMA, which would address this shortfall. However, these missions would not launch until the 2040s.

3.5 External photoevaporation of discs without an observed ionisation front

Proplyds are most easily identified because of their cometary morphology. However in weaker UV fields it is still possible to drive a significant wind that is essentially all launched only from close to the disc outer edge (where material is most weakly bound). Recent years have seen the discovery of possible external winds from very extended discs in very weak UV environments, down to \(F_{\mathrm {FUV}}\lesssim 10\,G_0\). The first example of this was the extremely large disc IM Lup, which has CO emission out to \(\sim 1000\,\)au. IM Lup had previously been demonstrated to have an unusual break in the surface density profile at around 400 au in submillimetre array (SMA) observations by [140]. [141] then observed IM Lup in the continuum and CO isotopologues with ALMA, similarly finding that the CO intensity profile could not be simultaneously reproduced at all radii by sophisticated chemical models, inferring a diffuse outer halo of CO. Using Hipparcos to map the 3D locations of the main UV sources within 150 pc of IM Lup and geometrically diluting their UV with various assumptions on the extinction, [141] had determined that the UV field incident upon IM Lup is only \(F_{\mathrm {FUV}}\sim 4\,G_0\), so not expected to be sufficient to drive an external wind. [142] demonstrated using 1D radiation hydrodynamic models that the CO surface density profile as a disc/halo could be explained by a slow external photoevaporative wind launched from around 400 au by an extremely low FUV environment, \(F_{\mathrm {FUV}} \sim 4 \, G_0\). This is possible because the disc is so extended that the outer parts are only weakly gravitationally bound, so even modest heating can drive a slow molecular wind. However, 2D models are required to give a more robust geometric comparison between simulations and observations to verify an outer wind in IM Lup.

Another candidate external wind in a \(F_{\mathrm {FUV}}<10\,G_0\) UV environment was identified in HD 163296 by [14] and [143]. They developed a framework in which the 3D CO emitting surface of the disc is traced, which can then be translated into a map of the velocity as a function of radius and height in the disc as illustrated in Fig. 9. Their main focus was the meridional flows identified at smaller radii in the disc, but they serendipitously discovered evidence for an outer wind launched from \(\sim 350-400\) au (see the right-hand blue box of Fig.  9). This is yet to be interpreted with any numerical models or synthetic observations, which will be necessary to support the interpretation that it is external photoevaporation that is responsible.

[143] also carried out a similar analysis of the disc MCW 480 but found no evidence of an outer wind despite having a similar radial extent and environment. Whether this is a consequence of the face-on orientation (\(\sim 33^\circ\) compared to \(\sim 45^\circ\) for HD 163296) or because there is no outer wind remains to be determined. Indeed, the HD 163296 kinematic feature that appears to be an outflow may also be due to some other mechanism. Furthermore, a similar approach is yet to be applied to IM Lup to search for kinematic traces of an outer wind.

Fig. 9
figure 9

The azimuthally averaged velocity (vectors) at the height of CO emission as a function of radius in HD 163296 by [143]. In addition to meridional flows, there is detection of a possible outer wind at \(\sim 320-400\) au, highlighted by the blue dashed box on the right

Looking ahead, determining whether external irradiation can really launch winds from discs down to FUV fluxes \(F_{\mathrm {FUV}}<10\,G_0\) is important for understanding just how pervasive external photoevaporation is.

3.6 Summary and open questions for observational properties of externally photoevaporating discs

  1. 1.

    External photoevaporation has been directly observed (e.g. proplyds) for almost 30 years. The vast majority of proplyd observations are in the ONC.

  2. 2.

    In recent years direct evidence for external photoevaporation is being identified in more regions, and B stars are now also known to facilitate it. However, the range of environments in which it has been observed is still very limited.

Some of the many open questions are:

  1. 1.

    What are robust diagnostics and signposts of external winds in weak UV environments?

  2. 2.

    What are diagnostics of external photoevaporation in distant clusters where a cometary proplyd morphology is spatially unresolved?

  3. 3.

    How widespread and significant are external winds from discs in weak UV environments?

4 Impact on disc evolution and planet formation

In this section, we consider how a protoplanetary disc evolves when exposed to strong external UV fields, and the consequences for planet formation. We will only briefly describe some relevant processes for planet formation in isolated star-disc systems. This is a vast topic with several existing review articles on both protoplanetary discs (e.g. [144,145,146,147]) and the consequences for forming planets (e.g. [148,149,150,151]) to which we refer the reader. Adams [152] reviewed general influences on planet formation (see also [153]), with a focus on the Solar System. Here we focus on a detailed look on how external photoevaporation affects the formation of planetary systems.

4.1 Gas evolution

4.1.1 Governing equations

The gas dominates over dust by mass in the interstellar medium (ISM) by a factor \(\rho _{\mathrm {g}}/\rho _{\mathrm {d}}\sim 100\) [154, 155]. This ratio is usually assumed to be similar in (young) protoplanetary discs, although CO emission suggests that the gas may be somewhat depleted with respect to the dust [5]. Nonetheless, gas is a necessary ingredient in instigating the growth of planetesimals by streaming instability [156] and represents the mass budget for assembling the giant planets [157, 158]. Thus, how the gas evolves in the protoplanetary disc is one of the first considerations in planet formation physics.

Despite its importance, the gas evolution of the disc remains uncertain. The observed accretion rates of \(\dot{M}_{\mathrm {acc}} \sim 10^{-10} {-}10^{-7} \, M_\odot\) yr\(^{-1}\) onto young stars [37, 159, 160], depending on the stellar mass [161, 162], imply radial angular momentum transport. For several decades, this angular momentum transport has widely been assumed to be driven by an effective viscosity, mediated by turbulence that may originate from the magnetorotational instability [163]. In the absence of perturbations, the surface density \(\varSigma _{\mathrm {g}}\) of the gaseous disc as a function of cylindrical radius r then follows [164]:

$$\begin{aligned} {{\dot{\varSigma }}_{\mathrm {g}}} = \frac{1}{r} \partial _r \left[ 3 r^{1/2} \partial _r \left( \nu \varSigma _{\mathrm {g}} r^{1/2} \right) \right] - {\dot{\varSigma }}_{\mathrm{{int}}} -{\dot{\varSigma }}_{\mathrm{{ext}}}. \end{aligned}$$
(15)

The loss rates \({\dot{\varSigma }}_{\mathrm{{int}}}\) and \({\dot{\varSigma }}_{\mathrm{{ext}}}\) are the surface density change due to the internally and externally driven winds, respectively. The kinematic viscosity \(\nu\) is usually parametrised by an \(\alpha\) parameter [165] such that \(\nu (r) = \alpha c_{\mathrm{{s}}}^2/\varOmega _{\mathrm {K}}\), for sound speed \(c_{\mathrm {s}}\) and Keplerian frequency \(\varOmega _{\mathrm {K}}\). In a disc with a mid-plane temperature \(T\propto r^{-1/2}\), this yields \(\nu \propto r\) [166].

In the following discussion, we will assume angular momentum transport is viscous. However, recently several empirical studies of discs have suggested a low \(\alpha \sim 10^{-4} {-} 10^{-3}\) [167,168,169]. This is difficult to reconcile with observed accretion rates. Alternative candidates, particularly magnetohydrodynamic (MHD) winds, have been suggested to drive angular momentum transport (e.g. [170]). In this case, angular momentum is not conserved but extracted from the gas disc, with consequences for the disc evolution [171, 172]. In the following, we will assume a standard viscous \(\alpha\) disc model, with the caveat that future simulations may offer different predictions by coupling the externally driven photoevaporative wind with MHD-mediated angular momentum removal.

4.1.2 Implementing wind driven mass loss

In order to integrate Eq. 15, we must define the form of \({\dot{\varSigma }}_{\mathrm{{int}}}\) and \({\dot{\varSigma }}_{\mathrm{{ext}}}\). The internal wind may be driven by MHD effects (e.g. [170, 173]) or thermally due to a combination of EUV (e.g. [174, 175]), X-ray (e.g. [32, 176, 177]) and FUV (e.g. [178]) photons, or probably a combination of the two (e.g. [179,180,181]). We do not focus on the internally driven wind in this review, but note that the driving mechanism influences the radial profile of \({\dot{\varSigma }}_{\mathrm{{int}}}\) (see [182], for a review).

Several authors have included the external wind in models of (viscous) disc evolution (e.g. [29, 99, 101, 183,184,185]). In general, these studies follow a method similar to that of Clarke [99] in removing mass from the outer edge, because winds are driven far more efficiently where the disc is weakly gravitationally bound to the host star. We discuss the theoretical mass loss rates in Sect. 2; in brief, the analytic expressions by Johnstone et al. [33] are applied to compute the mass loss rate in the EUV-driven wind, while early studies applied the expressions by Adams et al. [58] for an FUV-driven wind. The latter has now been improved upon using more detailed models by Haworth et al. [59], such that one can interpolate over the FRIED grid to find an instantaneous mass loss rate. For typical EUV fluxes, once the disc is (rapidly) sheared down to a smaller size, any severe externally driven mass loss is expected to be driven by FUV rather than EUV photons (see Sect. 2.2). For this reason, the EUV mass loss is often neglected in studies of disc evolution.

Since the mass loss rate is sensitive to the outer radius \(R_{\mathrm {out}}\), care must be taken when implementing a numerical scheme that the value of \(R_{\mathrm {out}}\) is sensibly chosen. In practice, a sharp outer radius quickly develops for a disc with initial mass loss rate higher than the rate of viscous mass flux (accretion). For a viscous disc with \(\nu \propto r\), the surface density \(\varSigma \propto r^{-1}\) in the steady state, which is the same profile as adopted for the numerical models in the FRIED grid. One can then interpolate using the total disc mass and outer radius. The latter is evolved each time-step by considering the rate of wind-driven depletion versus viscous re-expansion (e.g. [99, 186]). However, the physically correct way to define the outer radius is to find the optically thin-thick transition, since the flow in the optically thin region is set by the wind launched from inside this radius. Mass loss scales linearly with surface density in the optically thin limit [50], such that one can define \(R_{\mathrm {out}}\) to be the value of r that gives the maximal mass loss rate for the corresponding \(\varSigma (r)\) in the disc evolution model (see discussion by [101]). Under the assumption of a viscous disc with \(\nu \propto r\) both approaches yield similar outcomes, but the approach of Sellek et al. [101] should be adopted generally. For example, this prescription would be particularly important in the case of a disc model incorporating angular momentum removal via MHD winds [172].

Fig. 10
figure 10

Evolution of the externally driven mass loss rates \(\dot{M}_{\mathrm {ext}}\) (black, from the FRIED grid), accretion rates \(\dot{M}_{\mathrm {acc}}\) (blue) and radii containing 90 percent of the disc mass (red) for a disc evolving under Eq. 15 with a viscous \(\alpha =3\times 10^{-3}\), no internal photoevaporation and varying the external FUV flux \(F_{\mathrm {FUV}}= 50 \,G_0\) (dotted), \(500\,G_0\) (dashed) and \(5000 \, G_0\) (solid). We stop integrating when the mass loss rate or inner surface density reaches the floor in the FRIED grid (\(\dot{M}_{\mathrm {ext}} = 10^{-10}\,M_\odot\) yr\(^{-1}\) and \(\varSigma _{\mathrm {g}}(1\,\mathrm {au}) = 0.1\) g cm\(^{-2}\), respectively)

4.1.3 Viscous disc evolution with external depletion of gas

In Fig. 10, we show examples for the evolution of the disc radius that contains 90 percent of the mass, \(R_{90}\), and corresponding externally driven mass loss rates (from [59]) due to the numerical integration of Eq. 15 (see also Figs. 3 and 4 of [101], for example). To illustrate the variation in mass loss and radius, we choose an initial scale radius \(R_{\mathrm {s}} = 100\) au and mass \(M_{\mathrm {disc}}=0.1\,M_\odot\), truncated outside of 200 au, and with a viscous \(\alpha =3\times 10^{-3}\). For simplicity, we ignore internal winds (\({\dot{\varSigma }}_{\mathrm {int}} = 0\) everywhere), in order to highlight the main consequences of the externally driven photoevaporative wind in isolation.

From such simple models, we gain some insights into how we expect disc evolution to be affected by externally driven mass loss. In the first instance, the efficiency of these winds at large r leads to extremely rapid shrinking of the disc. In a viscous disc evolution model, this shrinking continues until the outwards mass flux due to angular momentum transport balances the mass loss due to the wind such that the accretion rate \(\dot{M}_{\mathrm {acc}} \sim \dot{M}_{\mathrm {ext}}\) [187, 188]. This offers a potential discriminant for disc evolution models: while \(\dot{M}_{\mathrm {acc}} \sim \dot{M}_{\mathrm {ext}}\) for a viscously evolving disc in an irradiated environment, this need not be the case if angular momentum is removed from the disc via MHD winds or similar. In either case, the initial rapid mass loss rates of \(\dot{M}_{\mathrm {ext}} \sim 10^{-7} \, M_\odot\) yr\(^{-1}\) are only sustained for relatively short timescales of a few \(10^5\) yr. Because the mass loss rate is related to the spatial extent of proplyds (Eq. 14), this implies easily resolvable proplyds should be short-lived and therefore rare.

This rapid shrinking of the disc has consequences for the apparent viscous depletion timescale of the disc, leading Rosotti et al. [185] to expound the usefulness of the dimensionless accretion parameter:

$$\begin{aligned} \eta \equiv \frac{\tau _{\mathrm {age} }\dot{M}_{\mathrm {acc}}}{M_{\mathrm {disc}}}, \end{aligned}$$
(16)

where \(\tau _{\mathrm {age}}\), \(M_{\mathrm {disc}}\) are the age and mass of the disc, respectively, while \(\dot{M}_{\mathrm {acc}}\) is the stellar accretion rate. For disc evolution driven by viscosity, and where the disc can reach a quasi-steady-state, we expect \(\eta \approx 1\). Indeed, \(\eta \approx 1\) for discs in low mass local SFRs when using the dust mass \(M_{\mathrm {dust}}\) (or, more precisely, sub-mm flux) as a proxy for the total disc mass \(M_{\mathrm {disc}} = 100\cdot M_{\mathrm {dust}}\). While numerous processes can interrupt accretion and yield \(\eta <1\), only an outside-in depletion process can yield \(\eta >1\). Rosotti et al. [185] showed that this applies to a number of discs in the ONC, hinting that external photoevaporation has sculpted this population.

With the inclusion of internal disc winds, a number of disc evolution scenarios become possible for an externally irradiated disc. The internal wind drives mass loss outside of some launching radius \(R_{\mathrm {launch}} \approx 0.2 R_{\mathrm {g}}\), which is inside of the gravitational radius due to hydrodynamic effects [90]. Once (viscous) mass flux through the disc becomes sufficiently small (\(\dot{M}_{\mathrm {acc}} \lesssim \dot{M}_{\mathrm {int}}\), the mass loss in the internal wind), the internal wind depletes material at \(r\approx R_{\mathrm {launch}}\) faster than it is replenished, leading to gap opening. Subsequently, the inner disc is rapidly drained and inside-out disc depletion proceeds [28, 189]. Coleman and Haworth [29] discussed how the balance of internal and external photoevaporation can alter this evolutionary picture. In the limit of a vigorous externally driven wind, the outer disc may be depleted down to \(R_{\mathrm {launch}}\). In this case, the internal wind no longer drives inside-out depletion, and the disc is dispersed outside-in. In the intermediate FUV case (\(F_{\mathrm {FUV}}=500\, G_0\)), external disc depletion erodes the outer disc without reaching \(R_{\mathrm {launch}}\). In models for disc material where outer disc material is eventually transported inwards, then outer disc depletion should still shorten the lifetime to some degree. In Fig. 10, we see that only the disc exposed to \(F_{\mathrm {FUV}}=5000\,G_0\) (above which we consider ‘high FUV environments’ by our definitions in Sect. 2.2.1) is sheared down to \(R_{\mathrm {out}}< 20\) au before the inner surface density becomes small (lower than the FRIED grid range). Thus, if angular momentum transport is inefficient beyond this radius, then sustained exposure to high FUV fluxes \(F_{\mathrm {FUV}}\gtrsim 5000\, G_0\) should be required to shorten the disc lifetime. In this case, we may expect inside-out depletion for discs exposed to more moderate \(F_{\mathrm {FUV}}\), although the external depletion still reduces their overall mass.

4.2 Solid evolution

We now consider how the evolution of solids within the gas disc may be influenced by externally driven winds. The growth of ISM-like dust grains to larger aggregates and eventually planets is the result of numerous interrelated physical processes, covering a huge range of scales. We do not review these processes in detail here, but refer the reader to the recent reviews by Lesur et al. [146, with a focus on dust-gas dynamics] and Drazkowska et al. [149, with a focus on planet formation]. Due to the complexity of the topic, and the fact that the primary empirical evidence comes from local, low mass SFRs without OB stars, most studies to date have focused on dust growth in isolated protoplanetary discs. We here consider the results of the few investigations focused on dust evolution specifically in irradiated discs.

Fig. 11
figure 11

The effect of external FUV irradiation on the dust content of a protoplanetary disc, adapted from Sellek et al. [101]. Figure 11a shows the maximum grain size, relative to an isolated disc (\(F_{\mathrm {FUV}}=0\, G_0\)) as a function of time and radius in an irradiated disc with viscous \(\alpha =10^{-3}\), exposed to FUV flux \(F_{\mathrm {FUV}}=1000\, G_0\). Figure 11b shows the total fraction of dust removed by the wind over the lifetime of a disc as a function of the viscous \(\alpha\) and \(F_{\mathrm {FUV}}\)

4.2.1 Dust evolution

Sellek et al. [101] investigated the drift and depletion of dust in a viscously evolving protoplanetary disc. Dust is subject to radial drift, wherein dust moves towards pressure maxima (i.e. inwards in the absence of local pressure traps – [190]), as well as grain growth dependent on the local sticking, bouncing and fragmentation properties (see [54] and references therein). The drift velocity is determined by the Stokes number, which in this context is the ratio of the stopping time of the dust grain to the largest eddy timescale \(\sim \varOmega _{\mathrm {K}}^{-1}\). Near the mid-plane and in the Epstein regime, this can be approximated:

$$\begin{aligned} {\mathrm{St}} \approx \frac{\pi }{2} \frac{ \rho _{\mathrm {s}} a_{\mathrm {s}}}{\varSigma _{\mathrm {g}}}, \end{aligned}$$
(17)

where \(a_{\mathrm {s}}\) and \(\rho _{\mathrm {s}}\) are the grain size and density, respectively. The draining of large dust grains can be understood in terms of the balance of viscous mass flux and radial drift [54], such that the equilibrium value of \({\mathrm{St}}\) is:

$$\begin{aligned} {\mathrm{St}}_{\mathrm {eq}} = 3\alpha \left| \frac{\frac{3}{2} + \frac{\mathrm {d}\ln \varSigma _{\mathrm {g}}}{\mathrm {d}\ln r} }{\frac{7}{4} - \frac{\mathrm {d}\ln \varSigma _{\mathrm {g}}}{\mathrm {d}\ln r} } \right| \end{aligned}$$
(18)

for the standard \(\alpha\) disc model with \(\nu \propto r\). For \({\mathrm{St}}>{\mathrm{St}}_{\mathrm {eq}}\), dust drifts inwards regardless of the viscosity. At the outer edge of the disc \(|{\mathrm {d}\ln \varSigma _{\mathrm {g}}}/{\mathrm {d}\ln r}|\) becomes large, so that \({\mathrm{St}}_{\mathrm {eq}} \rightarrow 3\alpha\). Hence dust in the outer disc that grows above \({\mathrm{St}}>3\alpha\) will migrate rapidly inwards. External evaporation acts to increase \(|{\mathrm {d}\ln \varSigma _{\mathrm {g}}}/{\mathrm {d}\ln r}|\) in the outer disc, clearing this region of large dust grains. Figure 11a, adapted from Sellek et al. [101], shows how the external wind can rapidly evacuate the outer disc of large grains, dependent on the value of \(\alpha\). Given that this occurs on short timescales compared to the disc lifetime, it has consequences for planet formation and observational properties of discs, possibly explaining why the sub-mm flux-radius relationship seen in low mass SFRs [191] does not hold in the ONC [6].

The clearing of large grains from the outer disc also has consequences for the quantity of solid material that can be lost in the wind. Sellek et al. [101] use the entrainment constraints given by Eq. 9 to estimate an entrainment fraction \(f_{\mathrm {ent}}\) for a given grain size distribution. The fraction of dust removed in their global externally irradiated disc model is shown in Fig. 11b. When viscosity is large (\(\alpha \gtrsim 10^{-2}\)), this fraction can exceed 50 percent for \(F_{\mathrm {FUV}}\gtrsim 10^3\, G_0\). This trend of higher dust depletion with higher \(\alpha\) is both due to the faster replenishment of disc material in the outer regions where the wind is launched, and because inward drift of large grains is less efficient (higher \({\mathrm{St}}_{\mathrm {eq}}\)). However, in general depletion of dust is not efficient, with typically less than half of the dust mass being lost. However, in the models of Sellek et al. [101] this does not result in an enhancement of the dust-to-gas ratio due to the enhanced loss of the large dust grains via rapid inwards drift.

A big caveat of the above discussion is that it does not consider the role of local pressure traps in halting the radial drift (e.g. [192, 193]). This local accumulation of solids can also lead to a mutual aerodynamic interaction that leads to a local unstable density growth that can seed planet formation [156]. If sufficient quantities of dust remain when the gas component is depleted by external photoevaporation, then this may serve to initiate planetesimal formation in the outer disc [194], perhaps leading to low mass or rocky planets rather than gas giants. Future studies may consider how the efficiency of dust trapping in the outer disc affects this picture.

4.2.2 Planet formation and evolution

As discussed in the introduction, the influence of external photoevaporation is still rarely included in models for planet formation, despite its apparent prevalence. However, Ndugu et al. [195] studied how increases in the disc temperature due to external irradiation alter the formation process. By increasing the disc scale height (thus decreasing mid-plane density), this heating reduces the efficiency of pebble accretion and giant planet formation in the outer disc. In this framework, giant planets that do form at high temperatures more frequently orbit with short periods (hot/warm Jupiters) because their planet cores need to form early or in the inner disc. Temperature also has an influence after the formation of a low mass planet core. Low mass planets that have not opened up a gap in the gas disc undergo type I migration, which is due to torques associated with a number of local resonances [196]. These torques are sensitive to thermal diffusion, such that they also depend on local temperature and associated opacity. This may lead to complex migration behaviour for the low mass planets [197], which would also be influenced by increasing the disc temperature due to external irradiation [130].

Perhaps more directly, where there is sufficient mass loss from the disc due to an externally driven wind, this also reduces the time and mass budget available for (giant) planet formation and migration. Internal winds have long been suggested to play an important role in stopping inward planet migration [198]. Armitage et al. [199] and Veras and Armitage [200] investigated a how type II migration of giant planets (within a gap) can be severely altered by the loss of the disc material in a photoevaporative wind, even leading to outwards migration if the outer disc material is removed. Since then, a number of authors have investigated how giant planet migration can be halted by internally driven disc winds (e.g. [30, 201, 202]). Planet population synthesis models have now started to implement prescriptions for mass loss due to an external wind (such as in the Bern synthesis models – [203]). However, this is presently based on a single (typical) estimated \({\dot{\varSigma }}_{\mathrm {ext}}\) that is constant throughout time and with radius r, rather than the more physical outside-in, radially dependent depletion models discussed in Sect. 4.1. Recently, Winter et al. [204] looked at how growth and migration are altered by external FUV flux exposure, and we show the outcomes of some of these models in Fig. 12. The outside-in mass loss prescription leads to more dramatic consequences than those of Emsenhuber et al. [203], curtailing growth and migration early. As a result, low FUV fluxes \(F_{\mathrm {FUV}}\lesssim 100\, G_0\) produce planets that are massive, \(M_{\mathrm {p}}\gg 100\,M_\oplus\) and on short orbital periods (\(P_{\mathrm {orb}}\ll 10^4\) days), similar to those that are most frequently discovered by radial velocity (RV) surveys (e.g. [205]). The typical FUV fluxes in the solar neighbourhood are \(\langle F_{\mathrm {FUV}}\rangle \sim 10^3 {-}10^4\, G_0\) (see Sect. 5.2), which yield typical planet masses and orbital periods closer to those of the massive planets in the Solar System. Such relatively low mass planets fall close to, or below, typical RV detection limits, and the anti-correlation between planet mass and semi-major axis may therefore contribute to the inferred dearth of detected planets with periods \(P_{\mathrm {orb}}\gtrsim 10^3\) days [206]. Testing the role of external photoevaporation for populations of planets further requires coupling these prescriptions with population synthesis efforts.

While we can try to identify the role of external photoevaporation via correlations in disc populations, more direct ways to connect environment with the present-day stellar (and exoplanet) population would be useful. For an example of how this might work in practice, Roquette et al. [207] have highlighted that premature disc dispersal may leave an impact on the rotation period distribution of stars. The rotation of a star is decelerated due to ‘disc-locking’ during the early pre-main-sequence phase accelerates rotation. Thus, by shortening the inner disc lifetime (see Sect. 4.3.2), external photoevaporation may leave an imprint on the stellar population up to several 100 Myr after formation. This may explain, for example, the increased number of fast rotating stars in the ONC compared to Taurus [208]. In future, models and rotation rates may be used to interpret the disc dispersal timescales for main sequence stars, which may complement investigations into planet populations in stellar clusters (e.g. [209,210,211]).

Fig. 12
figure 12

The final masses \(M_{\mathrm {p}}\) and orbital periods \(P_{\mathrm {orb}}\) of giant planets evolving under a growth and evolution model for planetary cores of mass \(10\, M_\oplus\) and formation time of \(\tau _{\mathrm {form}} = 1\) Myr with disc viscosity \(\alpha =3\cdot 10^{-3}\) (adapted from [204]). The circular points, connected by faint dotted green lines for fixed starting semi-major axis \(a_{\mathrm {p},0}\), represent the final location of the planet in \(P_{\mathrm {orb}}-M_{\mathrm {p}}\) space. Points are coloured by the external FUV flux experienced in that model, with the same colour is shown in a Voronoi cell to emphasise trends. Red crosses show the locations of HARPS radial velocity (RV) planet discoveries presented by Mayor et al. [205]. The 50 percent detection efficiency of the HARPS survey is shown by the red contour. We show orbital period at longer than which Fernandes et al. [206] infers a dearth of planets (dotted black line) and the limit above which planets may no longer form by core accretion (black dashed line – [212]). We also show the planets in the Solar System (green triangles)

4.3 Disc surveys of local star-forming regions

In this section so far, we focus on the theoretical influence of external photoevaporation on planet formation. However, the most important evidence for or against the influence of external photoevaporation on forming planets must be found statistically in the surveys of local protoplanetary disc populations. Such surveys and more detailed observations of individual discs offer the most direct way to probe the physics of planet formation. Here we report the evidence in the literature for the role of external photoevaporation in sculpting disc populations. We first consider the observational approaches and the challenges in inferring disc properties in Sect. 4.3.1. We then consider the evidence for variations in (inner) disc survival fractions (Sect. 4.3.2) and outer disc properties (Sect. 4.3.3) with FUV flux exposure.

4.3.1 Methods and challenges

Inner disc lifetimes Photometric censuses of young stars in varying SFRs can yield insights into the survival time of discs in regions of similar ages. Young stars exhibit luminous X-ray emission due to the magnetically confined coronal plasma [213], such that X-ray surveys with telescopes such as Chandra offer the basis for constructing a catalogue of young members of a SFR. These can be coupled with photometric surveys to infer the existence or absence of a NIR excess. Comparison of the fractions of discs with varying FUV flux either within the same SFR, or between different regions with similar ages, allows one to identify regions where discs have shorter lifetimes than the \(\sim 3\) Myr typical of low mass SFRs.

While this principle appears simple, in fact several steps required in achieving this comparison bear with them numerous pitfalls. One issue is reliable membership criteria, which were historically photometric or spectroscopic (e.g. [214, 215]), improved recently through proper motions (and to a lesser extent, parallax measurements) from Gaia DR2 [216].

Uncertainty and heterogeneity in age determination for young SFRs also represent a significant challenge, particularly in comparing across different SFRs (e.g. [217]). Michel et al. [218] point out that considering only low-mass, nearby regions, the isolated disc life-time might be a factor \(2{-}3\) longer than the canonical 3 Myr obtained by aggregating across many star-forming regions with a ride range of properties. Further, since disc life-times are shorter around high mass stars [2], care must be taken when comparing across samples of discs that may have different sensitivity limits.

Binning young stars by incident FUV flux also carries complications. The three-dimensional geometry (projection effects) and dynamical mixing in stellar aggregates may also hide correlations between historic FUV exposure and present-day disc properties (e.g. [219]). Even empirically quantifying the luminosity and spectral type of neighbour OB stars can be challenging. Massive stars are often found in multiple systems (e.g. [220]), which can lead to mistaken characterisation (e.g. [221]), while these massive stars are also expected to go through transitions in the Hertzsprung–Russell diagram in combination with rapid mass loss (see [222], for a review). Statistically, any studies attempting to measure correlations in individual regions must choose binning procedures for apparent FUV flux with care. For example, the number of stars per bin must be sufficient such that uncertainties are not prohibitively large and binning should be performed by projected UV flux rather than separation alone (i.e. controlling for the luminosity of the closest O star).

Finally, studies of NIR excess probe the presence of inner disc material, which represents the part of the disc expected to be least affected by external photoevaporation (see discussion at the end of Sect. 4.1.3). Therefore, inner disc survival fractions should be interpreted as the most conservative metric by which to measure the role of external disc depletion.

Outer disc properties

The outer regions of the protoplanetary disc are relatively weakly bound to the host star than the inner disc and are therefore more easily unbound by externally driven photoevaporative winds. Probing disc mass and radius also offers much more information than simply the presence or absence of circumstellar material. Outer disc properties are frequently inferred via probing the dust content, then assuming a canonical dust-to-gas ratio of \(10^{-2}\) to infer a total mass, although this may be significantly higher in some cases [5]. In surveys of protoplanetary discs, the dust mass is usually inferred by making the assumption that the disc is optically thin such that:

$$\begin{aligned} M_{\mathrm {dust}} \approx \frac{S_{\nu , \mathrm {dust}}D^2}{\kappa _{\nu , \mathrm {dust}} B_{\nu }(T_{\mathrm {dust}}}, \end{aligned}$$
(19)

where D is the distance to the source, \(F_{\nu , \mathrm {dust}}\) is the flux from the cool dust at frequency \(\nu\), \(T_{\mathrm {dust}}=20\) K is the dust temperature, \(B_{\nu }\) is the Planck function and \(\kappa _{\nu , \mathrm {dust}}\approx 2\) cm\(^{2}\) g\(^{-1}\) is the assumed opacity [223]. Even for studies of discs in low mass star-forming regions, the estimate from Eq. 19 comes with several assumptions, including the fixed dust temperature and opacity. In the context of discs in irradiated environments, this can even be even further complicated by the heating of the dust by neighbouring stars [130]. Even further, in some regions background free–free emission may contribute to the continuum flux and must therefore be subtracted by extrapolating from cm-wavelength observations (see [6], and references therein). However, since free–free emission may be variable, this is another source of uncertainty.

Beyond dust masses, estimates of the outer disc extents can be made using spatially resolved continuum observations by either fitting a surface density profile or defining an effective radius enclosing a fraction of the total disc emission (see discussion by [191, 224]). This comes with the significant caveat that the size-dependent inward drift of dust grains mean that the inferred radii may not trace the physical disc radius, with long integration times required to trace the small grains that remain well-coupled to the gas [225].

Probing the gas content of discs cannot be achieved by directly measuring the hydrogen. Molecular hydrogen is a symmetric rotator and thus has no electric dipole moment. Quadrupole transitions are not excited at the typical temperature for the bulk of the gas in protoplanetary discs (e.g. [226]). Instead CO isotopologues, and sometimes HD [227] are commonly used to infer disc masses [228], which requires some assumption for the carbon and oxygen abundances. Outer disc radii might be inferred from spatially resolved moment 0 maps using a similar approach as discussed for resolved dust observations (e.g. [229]), or by assuming Keplerian rotation and fitting a model to the gas kinematics (e.g. [230]). All of these methods rely on some assumptions about disc abundances and chemistry, which should be applied with caution in the case of irradiated discs. In particular the heating, winds/sub-Keplerian rotation and photodissociation of molecules close to the outer edge of the disc will all influence the observed molecular line emission [128].

Proplyd definitions Proplyds are discussed in detail in Sect. 3, and we do not revisit them in detail here. As discussed in that section, we distinguish between photoevaporating globules and proplyds, the latter of which are typically smaller than a few 100 au and are synonymous with an ionised star-disc system. For this reason, we here consider the examples in Trumpler 14 [122], Pismis 24 [231], NGC 3603 [126] and Cygnus OB2 [124] as candidate proplyds only. These range in size from \(\sim 1000\) au to several \(10^4\) au in size and are generally undetected in X-rays. We therefore consider these large photoevaporating objects as candidate proplyds or globules, rather than confirmed proplyds. This does not necessarily mean that true proplyds do not exist in these regions, as these are challenging to resolve at such large distances.

Table 1 Properties of some notable SFRs that show evidence of external depletion, except Westerlund 1 which we include since it is the subject of an upcoming JWST campaign [232]

4.3.2 Disc survival fractions

Despite the numerous difficulties, many studies have reported correlations between disc lifetimes and local FUV flux. One of the earliest was Stolte et al. [254], who found an increase from \(20\pm 3\) percent in the central 0.6 pc of NGC 3603, increasing to \(30\pm 10\) percent at separations of \(\sim 1\) pc from the centre. Later, Balog et al. [264] presented a Spitzer survey of NGC 2244 (of age \(\sim {2}\) Myr – [265, 266]) found a marginal drop off in the disc hosting stars for separations \(d< 0.5\) pc from an O star (\(27\pm 11\) percent) versus those at greater separations (\(45\pm 6\) percent). Guarcello et al. [267] obtained similar results for NGC 6611 (of age \(\sim 1\) Myr), finding \(31.1\pm 4\) percent survival in their lowest bolometric flux bin, versus \(16\pm 3\) percent in their highest (see also [268]). Similar evidence of shortened inner disc lifetimes has been reported in Arches [246, 269], Trumpler 14, 15 and 16 [248], NGC 6357 (or Pismis 24 – [231]) and the comparatively low density OB association Cygnus OB2 [252]. We highlight that correlations between location and inner disc fraction are not found ubiquitously—for example, Roccatagliata et al. [270] do not recover evidence of depleted discs in the central 0.6 pc of IC 1795. We return to discuss some of these cases in greater detail below.

Interestingly, the study by Fang et al. [231] of Pismis 24 revealed not only a correlation of disc fraction with FUV flux, but also a stellar mass-dependent effect wherein disc fractions are found to be lower for higher mass stars. While this is similar to the case of non-irradiated discs [2], it is the opposite of what might be expected from the dependence of \(\dot{M}_{\mathrm {ext}}\) on the stellar mass; lower mass stars are more vulnerable to externally driven winds due to weaker gravitational binding of disc material (e.g. [59]). This finding, if generally true for irradiated disc populations, would therefore put constraints on how other processes in discs scale with stellar host mass. For example, more rapid (viscous) angular momentum transport for discs around higher mass stars would sustain higher mass loss rates for longer, due to the mass flux balance in the outer disc (see discussion in Sect. 4.1). However, although accretion rates correlate with stellar mass [161, 162], this does not necessarily imply faster viscous evolution for discs with high mass host stars [271]; thus, the physical reason for the Fang et al. [231] findings remains unclear. Whether discs around high mass stars are generally depleted faster than those around low mass stars may be further confirmed by the upcoming JWST GTO program by Guarcello et al. [232], who aim to map out disc fractions and properties in Westerlund 1 for stars down to brown dwarf masses. This large dataset should allow to control both for stellar mass and location in the starburst cluster.

All of the regions mentioned above, in which evidence of inner disc dispersal has been inferred, share the property that they are sufficiently massive to host many O stars. A notable example for which shortened disc lifetimes are not observed is in the ONC. Despite rapid mass loss rate for the central proplyds, Hillenbrand et al. [41] inferred a high disc fraction of \(\sim 80\) percent. This may be due to a large spread in ages (typically estimated at \(\sim 1{-}3\) Myr – e.g. [242, 243, 272]), such that the resultant dynamical evolution [273, 274] leads to the central concentration of young stars [243, 272, 275], as discussed in Sect. 5.4. However, as always, interpreting the luminosity spread as a physical spread in ages is not so simple. Luminosity spreads may result from disc orientation, accretion, multiplicity and extinction, while the latter two effects may contribute to systematic gradients in the cluster. High survival rates may also reflect the inefficiency of inner disc clearing of photoevaporative winds in intermediate \(F_{\mathrm {FUV}}\) environments, which would hint at inefficient angular momentum transport at large radii (due to dead zones, for example).

Some other studies have also obtained null results when trying to find evidence of spatial gradients in the fraction of stars with NIR excess within SFRs. Studies searching for spatial correlations in NIR excess fraction are useful because they are not subject to the same large uncertainties in the stellar ages. However, for studies of this kind the question of membership criteria to the region is of utmost importance, since given a constant number of foreground/background contaminants then the outer regions should be more affected, presumably suppressing the apparent disc fraction. One example of a search for these correlations is that of Richert et al. [276], who found an absence of any correlation in the MYStIX catalogue [277] of NIR sources across several O star hosting SFRs. The methodology of Richert et al. [276] differed from other studies in that they adopted the metric of the aggregate density of disc-hosting or disc-free stars around O stars, rather than distances of these lower mass stars to their nearest O star. Similarly, Mendigutía et al. [278] used Gaia photometry and astrometry to compare relative disc fractions in a sample of SFRs, finding no spatial correlations. This study did not use FUV flux/O star proximity, but binned stars into the inner 2 pc versus the outer \(2-20\) pc for each region. Both studies have not attempted an absolute measure of disc survival fractions, but relative numbers. They highlight that environmental depletion of inner discs is not necessarily ubiquitous in high mass SFRs. Physical considerations such as age and dynamics may also play an important role in determining whether correlations in disc survival fractions can be uncovered.

In Fig. 13, we show the disc survival fraction versus SFR age for a composite sample, including a number of the massive regions discussed in this section. This figure is subject to the caveats discussed above, and particularly that observed samples of stars in massive, distant regions may be biased towards higher stellar masses than close by regions. Setting aside the caveats, Fig. 13 appears to demonstrate a shortening of disc lifetimes across numerous massive SFRs. Low mass SFRs typically have \(\tau _{\mathrm {life}}\approx 3\) Myr, while the most massive SFRs may have \(\tau _{\mathrm {life}} \lesssim 1\) Myr. This shortening of disc lifetimes has so far only been found in regions with a total FUV luminosity \(L_{\mathrm {FUV}}\gtrsim 10^{39}\) erg s\(^{-2}\) (see Table 1). In the case of many of these massive regions, the shortening of disc lifetimes can also be seen in a local gradient, for which stars that are far from an O star have a greater probability of hosting a disc. This is more difficult to explain by variations in the stellar masses of surveyed stars, since dynamical mass segregation occurs on much longer time-scales (e.g. [279]). This suggests that external photoevaporation can shorten inner disc disc lifetimes. If this is the case, it implies that if gas giants form in these environments, they must do so early (as suggested by some recent results – e.g. [114, 280]).

4.3.3 Outer disc properties

In recent times, ALMA has become the principle instrument for surveying outer disc properties in local star-forming regions. Nonetheless, prior to this revolution a number of studies had already demonstrated, by degrees, the role of external photoevaporation in sculpting the outer regions of protoplanetary discs. For example, in the ONC, Henney and Arthur [34] and Henney and O’Dell [35] used HST images and Keck spectroscopy to demonstrate the rapid mass loss rates up to \(\sim 10^{-6}\, M_\odot\) yr\(^{-1}\) of several proplyds in the core of the ONC, and their concentration around the O star, \(\theta ^1\)C. Later, Mann et al. [281] used SMA data to show a statistically significant dearth of discs with high dust masses close to \(\theta ^1\)C (see also [282]). This result has since been confirmed via ALMA observations towards the core [6] and outskirts [39] of the region. This is a clear demonstration that external photoevaporation depletes the dust content, although it remains unclear whether this is via instigating rapid radial drift Sellek et al. [101] or via entrainment in the wind [62]. Boyden and Eisner [283] also find the gas component of the disc is more compact than in other SFRs, and correlated with distance to \(\theta ^1\)C, suggestive of truncation by external photoevaporation.

Dust depletion has also been inferred in \(\sigma\) Orionis. Here the dominant UV source is the \(\sigma\) Ori multiple system, for which the most massive component has a mass \(17\,M_\odot\) [261]. From Herschel spectral energy distributions of discs in this region, Maucó et al. [240] found evidence of external depletion in the abundance of compact discs (\(R_{\mathrm {out}}<80\) au) consistent with the models of Anderson et al. [183]. The authors also found evidence of ongoing photoevaporation in one disc exhibiting forbidden [Ne ii] line emission (see also [284]). Later, Ansdell et al. [4] presented an ALMA survey that uncovered an absence of discs with high dust masses at distances \(\lesssim 0.5\) pc from \(\sigma\) Ori, consistent with models of dynamical evolution and disc depletion [285].

In NGC 2024, an early survey of the discs with the SMA presented by Mann et al. [286] revealed a more extended tail of high mass discs than those of other SFRs. The authors suggested this is indicative of a young population of discs that had not (yet) been depleted by external photoevaporation. However, van Terwisga et al. [7] presented an ALMA survey of the region and found evidence of two distinct disc populations. The eastern population is very young (\(\sim 0.2{-}0.5\) Myr old) and still embedded close to a dense molecular cloud that may shield them from irradiation. The western population is slightly older, exposed to UV photons from the nearby stars IRS 1 and IRS 2b. The western discs are lower mass than those in the east and those in similar age SFRs, which is probably due to external depletion. This conclusion is supported by the subsequent discovery of proplyd (candidates) in the region [113].

We have focused here on regions where evidence of external depletion appears to be present. It is challenging to convincingly demonstrate the converse by counter example, particularly considering the many potential issues discussed above (e.g. accurate stellar aging, projection effects and dynamical histories). As an example, in an ALMA survey of discs in \(\lambda\) Orionis, Ansdell et al. [3] reported no correlation between disc dust mass and separation from the OB star \(\lambda\) Ori. As discussed by the authors, a number of explanations for this are possible. One explanation is that, given the older age (\(\sim 5\) Myr) of the region, the discs may have all reached a similar state of depletion, with spatial correlations washed out by dynamical mixing (e.g. [219]). Meanwhile, in a survey of non-irradiated discs, van Terwisga et al. [287] demonstrated that discs in SFRs of similar ages appear to have similar dust masses, suggesting that any observed depletion in higher mass SFRs is probably the result of an environmental effect.

Fig. 13
figure 13

The fraction of surviving inner discs as a function of age of star-forming regions. Yellow points show a compilation of various star-forming regions using data from Richert et al. [288] and Mamajek [289]. We further show some specific SFRs discussed in the text and listed in Table 1. Massive and dense ‘starburst’ clusters are shown as star symbols. Where a local gradient in the disc fraction has been reported in the literature, we show the outer regions as a square symbol. We show the fraction of discs following \(f_{\mathrm {discs}} = \exp (-t_{\mathrm {age}}/\tau _{\mathrm {life}})\) for \(\tau _{\mathrm {life}} = 1\) and 3 Myr (dotted and dashed lines, respectively). An important caveat for this figure is that disc fractions in distant SFRs are based on samples that may include fewer low mass stars, which typically exhibit longer disc life-times

4.3.4 Summary and omissions

For quick reference, in Table 1 we summarise the properties of some local SFRs in which some observational evidence of external photoevaporation has been uncovered, and note some references for further reading. However, we emphasise that this is not a complete or representative list; we have chosen SFRs that exhibit strong evidence of external photoevaporation and well-studied stellar populations. Some notable exclusions include Trumpler 13/IC 1396A (the Elephant Trunk nebula) and IC 1805 (the Heart nebula). Trumpler 13 contains an O6.5V and O9V central binary [290], and the young stellar objects exhibit marginal spatial gradients in the disc survival fraction [291]. However, given the abundance of class I discs that indicate ongoing star formation, this may also originate from heterogeneous ages [292]. Similarly, in the \(2-3.5\) Myr old cluster IC 1805, with mass \(M_{\mathrm {tot}}\approx 2700\, M_\odot\) (\(L_{\mathrm {FUV}} \approx 10^{39.4}\) erg s\(^{-1}\)), any gradient in the surviving disc fraction is marginal [293]. A complex dynamical evolution, involving collapse and ejection of members [294], may confuse signatures of external disc depletion.

4.4 Missing aspects of disc evolution and planet formation in irradiated environments

Notably, we have not discussed chemistry in this section. The chemistry of discs (e.g. [295, 296]) and the planets they produce (e.g. [297,298,299]) is a complex function of the stellar metallicity and local disc temperature. For irradiated discs, temperatures increase in the outer regions and disc surface layers with respect to isolated discs, altering the chemistry in these regions, although not necessarily in the disc mid-plane (e.g. [300, 301]). However, how this alters planet formation and chemistry may dependent on disc evolution, and this remains a matter for future investigation. This issue may also soon be addressed empirically via JWST observations, with the GTO programme by Ramirez-Tannus et al. [244] aiming to probe the chemistry of discs of similar age but varying FUV flux histories in NGC 6357.

Apart from chemistry, this section has highlighted the many gaps in the understanding of the role of external photoevaporation in shaping the evolution of the disc and the formation of planets. From both theory and observations, we understand that gas disc lifetimes are shortened by external photoevaporation. This shortened lifetime of the gas disc alone should be enough to influence the formation and disc-induced migration of giant planets. Meanwhile, how the aggregation of dust produces planets in irradiated discs has only just started to be addressed (e.g. [101, 195]). When coupled with the apparent prevalence of external photoevaporation in shaping the overall disc population, the role of star formation environment must be considered as a matter of urgency for planet population synthesis efforts (e.g. [203, 302,303,304]). Meanwhile, the connection between stellar rotation periods and premature disc dispersal via external photoevaporation may in future offer a window into the birth environments of mature exoplanetary systems up to \(\sim 1\) Gyr after their formation [207].

4.5 Summary and open questions for irradiated protoplanetary disc evolution

While many aspects of disc evolution in irradiated environments remain uncertain, we can be relatively confident in the following conclusions:

  1. 1.

    External photoevaporation depletes (dust) masses and truncate the outer disc radii in regions of high FUV flux, at least for \(F_{\mathrm {FUV}} \gtrsim 10^3\, G_0\).

  2. 2.

    Inner disc lifetimes appear to be shorted for discs in regions where the total FUV luminosity \(L_{\mathrm {FUV}} \gtrsim 10^{39}\) erg s\(^{-1}\), while does not appear to be strongly affected for regions with lower \(L_{\mathrm {FUV}}\).

  3. 3.

    Dust evolution, low mass planet formation and giant planet formation all have the potential to be strongly influenced by external photoevaporation through changes in the temperature, mass budget, and time available for planet formation.

Some of the many open questions in this area include:

  1. 1.

    Do the externally driven mass loss rates in photoevaporating discs balance with accretion rates, as expected from viscous angular momentum transport?

  2. 2.

    If angular momentum transport is extracted via MHD winds rather than viscously transported, how does this influence the efficiency of external photoevaporation?

  3. 3.

    Is the lack of correlation shortened disc lifetime with FUV flux in intermediate mass environments related to a dead-zone or similar suppression of angular momentum transport?

  4. 4.

    Is the observed dust depletion in near to O stars due to entrainment in the wind or rapid inward drift?

  5. 5.

    In relation to this, how do external winds influence dust trapping and the onset of the streaming instability?

  6. 6.

    Conversely, how does the trapping of dust influence the dust depletion and dust-to-gas ratio in discs?

  7. 7.

    How does disc chemistry vary between isolated and externally FUV irradiated discs?

  8. 8.

    How does the local distribution of FUV environments influence the planet population as a whole?

5 Demographics of star-forming regions

In this section, we consider the properties of star-forming regions and the discs they host from an observational and theoretical perspective. We are interested in understanding: ‘what is the role of external photoevaporation of a typical planet-forming disc in the host star’s birth environment?’ The degree to which protoplanetary discs are influenced by external irradiation depends on exposure to EUV and FUV photons. Hence, the overall prevalence of external photoevaporation for discs (and probably the planets that form within them) depends critically on the physics of star formation and the demographics of star-forming regions.

5.1 OB star formation

OB stars have spectral types earlier than B3 and are high mass stars with luminosities \(> 10^3\, L_\odot\) and masses \(>5 \, M_\odot\). These stars are responsible for shaping a wide range of physical phenomena, from molecular cloud-scale feedback on Myr timescales [305, 306] to galactic nucleosynthesis over cosmic time [307]. The radiation feedback of these stars on their surroundings already plays a significant role during their formation stage, where stars greater than a few \(10\,M_\odot\) in mass must overcome the UV radiation pressure problem [308]. Forming OB stars may do this through monolithic turbulent collapse (e.g. [309]), early core mergers (e.g. [310]), or competitive accretion (e.g. [311]) – see Krumholz [312] for a review.

Two questions regarding OB stars are important with respect to local circumstellar disc evolution: ‘when do they form?’ and ‘how many of them are there?’. The former determines when the local discs are first exposed to strong external irradiation. The latter determines the strength and, importantly, the spatial uniformity of the UV field (we discuss this point further in Sect. 5.4).

The question of the timescale for formation of massive stars is empirically tied to the frequency of ultra-compact HII (UCHII) regions. These regions are the small (diameters \(\lesssim 0.1\) pc) and dense (HII densities \(\gtrsim 10^4\) cm\(^{-3}\)) precursors to massive stars. During the main sequence lifetime of the central massive star, the associated HII region will evolve from the embedded ultra compact state to a diffuse nebula. In the context of the surrounding circumstellar discs, the lifetime of a UCHII region represents the length of time for which the surroundings of massive stars are efficiently shielded from UV photons. These lifetimes can be inferred by comparing the number of main-sequence stars to the number of observed UCHII regions, yielding timescales of a few \(10^5\) yrs [313, 314]. The star-less phase prior to excitation of the UCHII region is short (\(\sim 10^4\) yr – [315, 316]), hence the UCHII lifetime represents the effective formation timescale for massive stars in young clusters and associations. This formation timescale is much shorter than the typical lifetime for protoplanetary discs evolving in isolation (\(\sim 3\) Myr – e.g. [1, 2]). Therefore, in regions with many OB stars we expect the local discs to be irradiated throughout this lifetime. However, this need not be the case when the number of OB stars is \(\sim 1\), in which case the time at which an OB star forms may be statistically determined by the spread of stellar ages in the SFR. We discuss this point further in Sect. 5.4.

The relative contribution of stars of a given mass to the local UV radiation field depends on the relative numbers of stars with this mass, or the initial mass function (IMF), \(\mathrm {d} n /\mathrm {d} \log m_*\). We can write the mean total luminosity \(\langle L_{\mathrm {SFR}} \rangle\) of a star-forming region (SFR) with N members:

$$\begin{aligned} \langle L_{\mathrm {SFR}} \rangle = N\langle L \rangle = N \int \frac{\mathrm {d} n }{\mathrm {d} \log m_* } \cdot L(m_*) \,\mathrm {d} \log m_* , \end{aligned}$$
(20)

where L is the luminosity of a single star. Hence, to understand the contribution of massive stars to the UV luminosity, we are interested in the shape of the IMF. The IMF in local SFRs exhibits a peak at stellar masses \(m_* \sim 0.1{-}1\, M_\odot\) and a steep power law \(\mathrm {d} n /\mathrm {d} \log m_* \propto m_*^{-\varGamma }\) with \(\varGamma \approx 1.35\) at higher masses (see [317], for a review). This power-law appears reasonable at least up to \(m_* \sim 100\, M_\odot\) [257, 318].

We combine the mass-dependent luminosity for young stars, as in the left panel of Fig. 14, with the Chabrier [319] IMF to produce the right panel of Fig. 14. Here we multiply the stellar mass-dependent luminosity by the IMF, as in the integrand on the RHS of Eq. 20, to yield the average contribution of stars of certain mass to the FUV and EUV luminosity of a SFR. Note that this is only accurate for a very massive SFR, where the IMF is well sampled. However, the figure demonstrates that for low mass SFRs—which here means where the IMF is not well sampled for stellar masses \(m_*\gtrsim 30\, M_\odot\) —both EUV and FUV luminosities are dominated by individual massive stars, rather than many low mass stars. In the following, we consider this in the context of the properties of local star-forming regions.

Fig. 14
figure 14

Left panel: stellar FUV and EUV luminosity as a function of stellar mass. The total luminosity is calculated using the effective temperature and total luminosity from the model grids of Schaller et al. [320], with metallicity \(Z=0.02\) and output time closest to \(1\,\) Myr. We apply the atmosphere models of Castelli and Kurucz [321] to give the wavelength-dependent luminosity. Right panel: The relative contribution of stars of a given mass to the total luminosity averaged across the IMF. Here we use a log-normal Chabrier [319] sub-solar IMF, and power-law with \(\varGamma =1.35\) for \(m_*>1\, M_\odot\)

5.2 Demographics of star-forming regions

We now consider the FUV flux experienced by protoplanetary discs in typical SFRs. We will focus on FUV, since as discussed in Sect. 4.1, the FUV is expected to dominate disc evolution over the lifetime of the disc. The word ‘typical’ in the context of SFRs is in fact dependent on the local properties within a galaxy, and we here refer exclusively to the solar neighbourhood (at distances \(\lesssim 2\) kpc from the sun). A number of studies have approached this problem in differing ways [188, 322,323,324]; however, all such efforts require estimating the statistical distribution \({\mathrm {d} n_{\mathrm {SFR}}}/{\mathrm {d} N}\) in the number of members N of SFRs, which dictates how many OB stars there are and therefore the local UV luminosity. For example, Fatuzzo and Adams [322] used two different approaches to this problem. One was to directly take the distribution from the list of nearby SFRs compiled by Lada and Lada [325]. While this is the most direct, it suffers from small number statistics for massive SFRs. The other approach by Fatuzzo and Adams [322] is to assume a log-uniform distribution in the number of stars existing in a SFR with a number of members N, truncated outside of a minimum \(N_{\mathrm {min}}=40\) and maximum \(N_{\mathrm {max}}=10^5\). Alternatively, one can produce a similar distribution using a smooth Schechter [326] function (see e.g. [327]):

$$\begin{aligned} \frac{\mathrm {d} n_{\mathrm {SFR}}}{\mathrm {d} N} \propto N^{-\beta } \exp \left( \frac{-N}{N_{\mathrm {max}}} \right) \exp \left( -\frac{N_{\mathrm {min}}}{N} \right) , \end{aligned}$$
(21)

where \(\beta \approx 2\) is expected from the hierarchical collapse of molecular clouds [328] and is consistent with empirical constraints (e.g. [327, 329, 330]), as well as the log-uniform distribution adopted by Fatuzzo and Adams [322, n.b. that to obtain the fraction of stars, one must multiply the mass function described by Eq. 21 by a factor N].

In order to interpret Eq. 21, we further need to estimate \(N_{\mathrm {min}}\) and \(N_{\mathrm {max}}\). The upper limit \(N_{\mathrm {max}}\) can be inferred empirically [327, 331], or by appealing to the theoretical limits. As elucidated by Toomre [332], the origin of the maximum mass of a SFR can be understood in terms of the maximum size of a region that can overcome the galactic shear (Toomre length) and local kinetic pressure (Jeans length). This length scale, the Toomre length, may then be coupled to a local surface density of a flattened disc to obtain a maximum mass for a SFR. However, in the outer regions of the Milky Way, the stellar feedback in massive SFRs can interrupt star formation and further reduce \(N_{\mathrm {max}}\). Adopting the simple model presented by Reina-Campos and Kruijssen [333] for this limit with a typical stellar mass \(m_*\sim 0.5\,M_\odot\), this yields a maximum \(N_{\mathrm {max}} \sim 7 \times 10^4\) for all SFRs (not necessarily gravitationally bound). This is broadly consistent with the most massive local young stellar clusters and associations [334].

Meanwhile, for the minimum number of members in a SFR, Lamers et al. [335] used the sample of Kharchenko et al. [336] to infer a dearth of SFRs with \(N<N_{\mathrm {min}} \sim 280\). Such a lower limit might be understood as the point at which the lower mass end of the hierarchical molecular cloud complex merges due to high star formation efficiency and long formation timescales [337]. However, obtaining constraints on \(N_{\mathrm {min}}\) in general is made difficult by the lack of completeness in surveys of low mass stellar aggregates, and it remains empirically unclear how \(N_{\mathrm {min}}\) varies with galactic environment. Note that generally SFR demographics may depend on galactic environment and therefore also vary over cosmic time (see [338], for a review).

Fig. 15
figure 15

Left panel: the FUV luminosities at an age of 1 Myr of \(6\times 10^4\) realisations drawing star-forming regions (SFRs), with the number of members drawn from Eq. 21 (cf. [322], Fig. 5b). We draw individual stellar masses from a Chabrier [319] IMF in the range \(0.08{-}100\, M_\odot\). The black points represent a subset of 1000 of the synthetic SFRs. The solid brown line follows the mean luminosity as a function of the number of members (Eq. 20), the dashed brown line is the median, while the dotted lines bound the \(14^{\mathrm {th}}{-}86^{\mathrm {th}}\) percentile range. The colour bar shows an estimate for the relative density of stars (i.e. density of SFRs multiplied by the number N of members) in logarithmic space using a Gaussian kernel density estimate with a bandwidth 0.1 dex. We estimate the numbers of stars and total luminosity of three local SFRs: Taurus (green), the ONC (purple) and Westerlund 1 (orange). Right panel: the corresponding distribution in the number of stars (red) and SFRs (black) per logarithmic bin in total FUV luminosity space

With the above considerations, we can now generate the distribution of UV luminosities experienced by stars in their stellar birth environment. We adopt \(N_{\mathrm {min}}=280\) and \(N_{\mathrm {max}}=7\times 10^4\) and couple Eq. 21 with the stellar mass dependent luminosity discussed in Sect. 5.1. The resultant FUV luminosity distribution when drawing \(6\times 10^4\) SFRs, each with N members drawn from Eq. 21, is shown in Fig. 15. For context, we estimate the FUV luminosity in Taurus using the census of members by Luhman [339], assuming the local luminosity is dominated by four B9 and three A0 stars in that sample of 438 members. We also add the ONC, whose total UV luminosity is dominated by the O7-5.5 type star \(\theta ^1\)C, with a mass of \(\sim 35 \, M_\odot\) (e.g. [340, 341]), and Westerlund 1 which has a total mass of \(\sim 6\times 10^4 \, M_\odot\) [259, 342] and thus a well sampled IMF. Regions like Taurus with only a few hundred members are common, and therefore we expect to find such regions nearby. However, each one hosts \(\sim 1/1000\) as many stars as a starburst cluster like Westerlund 1. Taurus thus represents an uncommon stellar birth environment in terms of the FUV luminosity. The most common type of birth environment for stars lies somewhere between the ONC and Westerlund 1, with FUV luminosity \(\sim 10^{40}\) erg s\(^{-1}\) (at an age of 1 Myr).

In order to understand the typical flux experienced by a circumstellar disc, we further need to consider the typical distance between it and nearby OB stars – i.e. the geometry of the SFR. To this end, Fatuzzo and Adams ([322], see also [343]) and Lee and Hopkins [323] use the Larson [344] relation that giant molecular clouds and young, embedded SFRs have a size-independent surface density (see also [345]); hence, the radius R of the SFR scales as \(R \propto \sqrt{N}\) [343, 346]. However, this relationship does not hold for very young massive clusters (see [334], for a review), with several exhibiting comparatively high densities (e.g. Westerlund 1 – [259]), while older clusters follow a shallower mass–radius relationship \(R \propto N^\alpha\) with \(\alpha \sim 0.2{-}0.3\) [347]. Instead of the Larson relation, Winter et al. [188] attempt to relate SFR demographics to galactic-scale physics by using the lognormal density distribution of turbulent, high Mach number flows combined with a theoretical star formation efficiency. This yields the stellar density distribution in a given galactic environment. In massive local SFRs, the local FUV flux \(F_{\mathrm {FUV}} \approx 1400 \, (\rho _*/10^3 M_\odot \,\mathrm {pc}^{-3})^{1/2} \, G_0\), for stellar density \(\rho _*\), due to their radial density profiles [186]. Using these two ingredients, it is possible to estimate the distribution of FUV fluxes experienced by stars in the solar neighbourhood from a mass function for SFRs. Despite the differences in the two approaches, both Fatuzzo and Adams [322] and Winter et al. [188] find that, neglecting extinction, young stars in typical stellar birth environments experience external FUV fluxes in the range \(\langle F_{\mathrm {FUV}} \rangle \sim 10^3{-}10^4 \, G_0\).

While these efforts produce some intuition as to the typical FUV radiation fields, this is not the end of the story for understanding how discs in SFRs are exposed to UV photons. In Sects. 5.3 and 5.4 we will discuss the role of interstellar extinction and the dynamical evolution of stars in SFRs.

Fig. 16
figure 16

FUV flux experienced by stars formed in the hydrodynamic simulations of Ali [348], following the collapse of a molecular cloud of mass \(10^5 \, M_\odot\) at time 0.5, 1.0 and 1.5 Myr from left to right. Top row: the gas column density and location of sink particles as points, with white points representing ionising sources. Bottom row: histogram of the number of sources binned by instantaneous FUV flux (see also [349]). Comparatively few stars are born in the early stages of the cluster formation, where attenuation of FUV photons is efficient

5.3 Extinction due to the interstellar medium

After stars begin to form, residual gas and, more importantly, the dust in the interstellar medium (ISM) attenuates FUV photons such that circumstellar discs may be shielded from external photoevaporation. The column density required for one magnitude of extinction at visual magnitudes is \(N_{\mathrm {H}}/A_{\mathrm {V} }= 1.8\times 10^{-21}\) cm\(^{-2}\) mag\(^{-1}\) [350], and this can be multiplied by a factor \(A_{\mathrm {FUV}}/A_{\mathrm {V} }\approx 2.7\) [351] to yield the corresponding FUV extinction. The main problem comes in estimating the column density between OB stars and the cluster members. Both Fatuzzo and Adams [322] and Winter et al [188] approach this problem by adopting smooth, spherically symmetric density profiles with some assumed star formation efficiency to estimate this column density. However, even during the embedded phase this approach overestimates the role of extinction because a more realistic, clumpy geometry of the gas makes the attenuation inefficient (e.g. [352]). In addition, stellar feedback from massive stars acts to drive gas outflows from the SFR even before local supernovae ignition (e.g. [353,354,355]), reducing the quantity of attenuating matter.

Due to the above concerns, understanding the influence of extinction on shielding the young disc population requires direct simulation of feedback in the molecular cloud. To approach this problem, Ali and Harries [356] simulated feedback from a single \(34\, M_\odot\) mass star in a molecular cloud of mass \(10^4 \, M_\odot\), similar to the conditions in the ONC. The authors find that many discs are efficiently shielded for the first \(\sim 0.5\) Myr of evolution, while the discs may experience short drops in UV exposure at later stages. This result is somewhat dependent on cloud metallicity, since lower metallicity increases the efficiency of feedback by lengthening the cooling time [348].

If more than one massive star forms in a SFR, this further increases feedback efficiency and geometrically reduces the efficiency of extinction. Qiao et al. [349] investigated the role of feedback in the simulations of a molecular cloud of mass \(10^5 \, M_\odot\) by Ali [348]. The resultant FUV flux experienced by the stellar population is shown in Fig.16. For such a massive region with several O stars, lower mass disc-hosting stars typically experience extremely strong irradiation by the time they reach an age of \(\sim 0.5\) Myr. Although extinction may be efficient for the stars that form early, the majority that form later are quickly exposed to \(F_{\mathrm {FUV}} \sim 10^5 \, G_0\). Hence, the embedded phase of evolution in massive SFRs can only shield discs early on. While the FUV flux experienced by discs in real SFRs also depends on the density of the regions, inefficient shielding may partially explain why disc life-times appear to be shortened in SFRs with several O stars (see Table 1). Any giant planet formation in such a region must therefore initiate early and may be strongly influenced by their environment (see Sect. 4.2.2).

5.4 The role of stellar dynamics

In some instances, the FUV exposure history of stars may strongly depend on the dynamical evolution of the SFR. This is particularly true only one dominant OB star is present. For example, the difference in \(F_{\mathrm {FUV}}\) for a disc at a separation \(d = 0.05\) pc from a single massive star (as for some of the brightest proplyds in the ONC) and those at \(d=2\) pc (a typical distance for discs in the ONC) is a factor 1600. Hence the dynamical evolution of the star and the SFR can play a major role in the historic UV exposure of any observed star-disc system.

Many studies have performed N-body simulations of SFRs to track the exposure of star-disc systems, either with aim of quantifying general trends [100, 219, 357, 358] or modelling specific regions [274, 285, 359]. These studies generally adopt the FUV driven mass loss rates given by the FRIED grid [59] with an external flux computed directly from N-body simulations. In general, these studies find that SFRs with typical stellar densities \(\rho _* \gtrsim 100\, M_\odot\) are sufficient to rapidly deplete protoplanetary discs (e.g. [100, 184, 358]). Similarly, discs are more rapidly depleted in sub-virial initial stellar populations that undergo cold collapse, as a result of the higher densities and therefore UV flux exposure [358]. However, initial substructure has a minimal effect on a external photoevaporation since UV fields are generally dominated by the most massive local stars and not necessarily nearest neighbours [219, 358]. Thus it is volume-averaged rather than star-averaged density measures in a SFR that act as a proxy for the typical external UV flux. This is not necessarily true for the bolometric flux, which may still be dominated by the closest neighbours in a highly structured SFR [323].

For an example of how the dynamics in SFRs may alter disc properties, we need only look at the closest intermediate mass SFR to the sun: the ONC. This region has historically dominated the study of externally photoevaporation protoplanetary discs since the discovery of proplyds (e.g. [22]). Surprisingly, up to \(\sim 80\) percent of stars exhibit a NIR excess [41], indicating inner-disc retention. This finding is apparently inconsistent with the \(\sim 1{-}3\) Myr age [242, 272] when accounting for reasonable stellar dynamics, mass loss rates and initial disc masses (e.g. [17, 359]). The ONC contains one O star, \(\theta ^1\)C, that dominates the local UV luminosity, tying the FUV history of the local circumstellar discs precariously to the history of this single star. Indeed, for this reason intermediate mass star-forming regions may be subject to multiple bursts of star formation due to the periodic ejection of individual massive stars [273]. This possibility seems to be supported by the clumped age distribution of stars in the ONC, hinting at multiple populations or phases of star formation [243, 360]. This possibility bears the usual caveat that luminosity spreads do not necessarily imply age spreads, as discussed in Sect. 4.3.2. However, Winter et al. [274] showed how, when considering the gravitational potential of residual gas, such episodic star formation can yield inward migration of young stars and outward migration of older stars. This yields a stellar age gradient as observed in the ONC [272]. As a result, much younger discs experiencing the strongest UV fluxes such that inner disc survival over their lifetime is feasible. Similar core-halo age gradients have been reported in other star-forming regions (e.g. [275]), highlighting the importance of interpreting disc properties through the lens of the dynamical processes in SFRs.

Unlike regions with one O star, dynamical evolution may not be such an important consideration for discs in massive SFRs. For example, tracking the UV fluxes experienced by stars in the neighbourhood of numerous OB stars, Qiao et al. [349] find that the vast majority of stars experience \(F_{\mathrm {FUV}} \sim 10^5 \, G_0\) at an age of 1 Myr. The relatively uniform flux in regions with multiple OB stars may go someway to explaining the absence of correlation between disc mass and location found in some simulations [219]. By contrast, Winter et al. [285] reproduce a disc mass-separation correlation in the relatively low mass region \(\sigma\) Orionis region [4], which is occupied by a single massive star (see Sect. 4.3). In either case, it is clear that the physics of star formation cannot be ignored when interpreting the properties of protoplanetary discs, and probably the resultant planetary systems.

5.5 Summary and open questions for the demographics of star-forming regions

Based on the previous discussion, the following conclusions represent the current understanding of the demographics of SFRs in the context of external disc irradiation:

  1. 1.

    When accounting for a standard IMF, the total FUV and EUV luminosity of a SFR is dominated by stars of mass \(\gtrsim 30 \, M_\odot\).

  2. 2.

    Although most SFRs are low mass and have few OB stars, the majority of stars form in regions with a total FUV luminosity greater than that of the ONC (\(\gtrsim 10^{39}\) erg s\(^{-1}\)).

  3. 3.

    In the absence of interstellar extinction, typical FUV fluxes experienced by star-disc systems in the solar neighbourhood are \(\langle F_{\mathrm {FUV}} \rangle \sim 10^3{-} 10^4 \, G_0\). This is enough to shorten their lifetime with respect to the typical \(\sim 3\) Myr for isolated discs.

  4. 4.

    Extinction due to residual gas in SFRs can shield some young circumstellar discs, with ages \(\lesssim 0.5\) Myr. However, at later times extinction is inefficient at shielding discs, and unattenuated FUV photons may reach protoplanetary discs.

Some of the open questions that remain in quantifying the typical external environments for planet formation:

  1. 1.

    How important are the dynamics of SFRs in determining FUV exposure for populations of discs, and how much does this vary between SFRs?

  2. 2.

    How does FUV exposure end external photoevaporation of typical protoplanetary discs vary with cosmic time?

  3. 3.

    Do particular types of planetary system preferentially form in certain galactic environments?

Fig. 17
figure 17

A cartoon for how protoplanetary disc evolution and planet formation proceeds in weakly irradiated, low mass SFRs (top) and strongly irradiated, high mass SFRs (bottom). We consider two identical, initially large discs with outer radii of \(\sim 100\) au. In the high mass SFR, any neighbouring massive stars may initially be extincted by an UCHII region for \(\sim 10^5\) yrs. However, before the star-planet system is 0.5 Myr old, extinction in the SFR will typically become inefficient and the disc is irradiated by a strong UV field. This produces a bright, extended ionisation front (IF) and rapidly truncates the disc down to a few 10 au. This may also interrupt the early giant planet formation that proceeds in the outer regions of an isolated disc in a low FUV environment. The irradiated disc is now small, which results in slow mass loss rates and a smaller IF. Given its compact size, the timescale for viscous draining of the irradiated disc becomes short compared to its isolated counterpart. This can lead to premature clearing of disc material earlier than the typical \(\sim 3{-}10\) Myr for which isolated discs survive

6 Summary and future prospects

In this review, we have discussed many aspects of the process of external photoevaporation, covering basic physics, observational signatures, consequences for planet formation and prevalence across typical star-forming regions. While numerous open questions remain, with the current understanding we can make an educated guess at how protoplanetary discs evolve in different star-forming regions. In Fig. 17, we show a cartoon of how disc evolution proceeds in low and high FUV flux environments. The early period of efficient extinction in SFRs typically lasts less than 0.5 Myr of a disc’s lifetime [349]. After this, studies of individual externally irradiated discs in proplyd form have demonstrated mass loss rates up to \(\sim 10^{-7}{-}10^{-6} \, M_\odot\) yr\(^{-1}\) [22, 35, 190]. This bright and extended proplyd state is short-lived, but the disc is rapidly eroded outside-in during this period. The influence of the external wind on the outer disc can result in dust loss via entrainment in the wind [62] or due to rapid inward migration of large grains [101]. It can also interrupt any giant planet formation that occurs on this timescale in the outer disc [114, 115, 280]. Given the shorter viscous timescale of the compact disc, this can lead to a rapid clearing of the disc material similar to expected after gap opening in internally photoevaporating discs [28]. This is corroborated by the observed shortening of inner disc lifetimes in several massive local SFRs (e.g. [246, 248, 252]), however only those in which there are several O stars and a total FUV luminosity \(L_{\mathrm {FUV}} \gtrsim 10^{39}\) erg s\(^{-1}\). In these extreme environments, where high FUV fluxes are sustained over the disc lifetimes, external photoevaporation presumably shuts off accretion onto growing planets and curtails inward migration, possibly leaving behind relatively low mass outer planets [204].

Within this framework several questions arise, requiring both theoretical and empirical future study. As an example, models for the chemical evolution of planetary discs and planets in irradiated environments, and their expected observational signatures, are urgently needed. Upcoming JWST investigations of high UV environments may shed some light on disc chemistry from an observational perspective (e.g. [244]). Inferring mass loss rates for moderately photoevaporating discs that do not exhibit bright ionisation fronts is also a crucial test of photodissociation region models. Meanwhile, perhaps the biggest problem for planet population synthesis efforts is determining how the solid content evolves differently in photoevaporating discs; when and where do solid cores form in the discs in high mass versus low mass SFRs? These are just some of the numerous questions that remain open towards the goal of understanding external photoevaporation.

In conclusion, the process of external photoevaporation appears to be an important aspect of planet formation, although the full extent of its influence remains uncertain. Future efforts in both theory and observations are required to determine how it alters the evolution of a protoplanetary disc, and the consequences for the observed (exo)planet population.