# Warm millimetre dust in protoplanetary discs near massive stars

Thomas J. Haworth1
1 Astronomy Unit, School of Physics and Astronomy, Queen Mary University of
London, London E1 4NS, UK E-mail<EMAIL_ADDRESS>
(Accepted ???. Received ???; in original form ???)

###### Abstract

Dust plays a key role in the formation of planets and its emission also
provides one of our most accessible views of protoplanetary discs. If set by
radiative equilibrium with the central star, the temperature of dust in the
disc plateaus at around $10-20$ K in the outer regions. However sufficiently
nearby massive stars can heat the outer disc to substantially higher
temperatures. In this paper we study the radiative equilibrium temperature of
discs in the presence of massive external sources and gauge the effect that it
has on millimetre dust mass estimates. Since millimetre grains are not
entrained in any wind we focus on geometrically simple 2D-axisymmetric disc
models using radiative transfer calculations with both the host star and an
external source. Recent surveys have searched for evidence of massive stars
influencing disc evolution using disc properties as a function of projected
separation. In assuming a disc temperature of $20$ K for a disc a distance $D$
from a strong radiation source, disc masses are overestimated by a factor that
scales with $D^{-1/2}$ interior to the separation that external heating
becomes important. This could significantly alter dust mass estimates of discs
in close proximity to $\theta^{1}$C in the Orion Nebular Cluster. We also make
an initial assessment of the effect upon snow lines. Within a parsec of an O
star like $\theta^{1}$C a CO snow line no longer exists, though the water snow
line is virtually unaffected except for very close separations of $\leq
0.01\,$pc.

###### keywords:

accretion, accretion discs – circumstellar matter – protoplanetary discs –
radiative transfer – planets and satellites: formation

††pagerange: Warm millimetre dust in protoplanetary discs near massive
stars–References††pubyear: 2020

##  1 introduction

Although we have known about the impact of environment on circumstellar discs
essentially for as long as we have been able to directly image them (O’dell &
Wen, 1994), there has recently been a resurgence of interest in the topic.

Most stars form in clustered environments (e.g. Lada & Lada, 2003; Krumholz et
al., 2019). The three main ways that a cluster environment affects discs is
through external photoevaporation, dynamical (gravitational) encounters and
through compositional inheritence and enrichment (e.g. of short lived
radionuclides like Aluminium 26, Lichtenberg et al., 2019; Reiter, 2020,
though we do not focus on this here). Recent developments with regard to the
first two points include a better understanding of the effect of dynamics and
external photoevaporation in different types of stellar cluster (Scally &
Clarke, 2001; Winter et al., 2018; Concha-Ramírez et al., 2019; Nicholson et
al., 2019) which were facilitated by improvements in external photoevaporation
models (Haworth et al., 2018b; Haworth & Clarke, 2019), observations of
interacting systems (e.g. Rodriguez et al., 2018; Kurtovic et al., 2018) and
an ever increasing catalogue of photoevaporating discs (e.g. Henney & O’Dell,
1999; Mann et al., 2014; Kim et al., 2016; Haworth et al., 2021). The recent
linking of possibly distinct exoplanet populations in Gaia phase space
overdensities also provides a tantalizing hint of the impact of environment
(at the formation stage, or over a longer period of time) on the resulting
planets themselves (Winter et al., 2020; Kruijssen et al., 2020; Longmore et
al., 2021).

Another powerful diagnostic comes from finding systematic statistical
variation of disc properties in a region. For example in Hubble Space
Telescope observations, disc sizes were found to vary with distance from the O
star $\theta^{1}$C by Tsamis et al. (2013) (see also Henney & Arthur, 1998;
Henney & O’Dell, 1999). The unprecedented sensitivity and resolution that ALMA
provides now also allows us to analyse disc dust mass and radius statistics
throughout star forming regions (e.g. Mann et al., 2014; Ansdell et al., 2017;
Eisner et al., 2018; Boyden & Eisner, 2020; Ansdell et al., 2020). In
particular, trends in disc properties as a function of projected separation
from the strongest UV sources is often interpreted as evidence for external
photoevaporation (Mann et al., 2014; Ansdell et al., 2017; Eisner et al.,
2018). However, the dynamical evolution of the cluster can drastically
complicate this picture (Parker et al. in preparation). Nevertheless trends do
remain that require an explanation.

The most accessible (i.e. least demanding in terms of observing time) insight
into the disc properties with ALMA comes from estimating the dust mass from
the measured continuum flux $F_{\nu}$. This is done by integrating the formal
solution to the equation of radiative transfer in the limit of being optically
thin, constant absorption and emission coefficient, no background source and
dust emitting as a blackbody (Planck function $B_{\nu}$)

$I_{\nu}=B_{\nu}\left(1-\exp(-\tau_{\nu})\right)\approx B_{\nu}\tau_{\nu}$ (1)

over the beam solid angle $\Omega_{B}$ to convert the emergent intensity
$I_{\nu}$ to a flux $F_{\nu}$

$F_{\nu}\approx B_{\nu}\tau_{\nu}\Omega_{B}=B_{\nu}\kappa_{\nu}\rho
l\Omega_{B}$ (2)

where $\kappa_{\nu}$ is the opacity and $\rho$ and $l$ are the density and
thickness of the emitter along the line of sight. The solid angle is then
related to the mass

$M_{d}=\rho l\Omega_{B}D^{2}$ (3)

hence

$M_{d}=\frac{D^{2}F_{\nu}}{\kappa_{\nu}B_{\nu}(T)}$ (4)

is the well known expression used to estimate dust masses from a measured flux
$F_{\nu}$ (Hildebrand, 1983). This relies upon an assumed distance $D$,
opacity at the frequency of observation $\kappa_{\nu}$ and dust temperature
$T$. In many prior studies of discs the dust temperature and mass was
estimated from SED fitting (e.g. Beckwith et al., 1990; Andre & Montmerle,
1994; Andrews & Williams, 2005, 2007). However often in ALMA mm continuum
surveys this dust temperature is assumed to be 20 K, which probably applies
well to nearby discs to the Sun which are very extended (hundreds of au in
size) and expected to be relatively cold over the bulk of the disc mass.
Tazzari et al. (2017) validated this assumption further with UV plane fitting
of discs in Lupus, finding no obvious dependence on disc temperature.

So for nearby discs, the host star heating only affects the very inner regions
and it does seem reasonable to estimate the disc mass using a cooler
temperature like 20 K that well represents the temperature of the bulk of the
disc mass. However when studying discs in relatively close proximity to
massive stars ($\sim 1\,$pc in the case of $\theta^{1}$C, as we will
demonstrate here), warming of the outer disc by external sources could mean
that disc dust masses in the denser parts of massive clusters are actually
being overestimated in surveys when a temperature of 20 K is assumed for all
discs. The implications here could be important for understanding the
statistical properties of discs in clusters. Furthermore, this external
heating could also affect the locations (or even existence) of snow lines in
the disc, as well as the nature of grain evolution and even the ability to
introduce pressure bumps, which could all affect the planet formation process
in these high radiation environments.

External irradiation has been modelled in strong radiation scenarios looking
at the effect on the gas dynamics/proplyds (e.g. Henney & O’Dell, 1999;
Richling & Yorke, 2000; García-Arredondo et al., 2001; Haworth & Clarke,
2019). However these calculations are computationally expensive, usually focus
on the dynamics, and when it comes to observables usually focus on the
photodissociation region (PDR) and photoionised gas diagnostics. Walsh et al.
(2013) compared the temperature and chemistry of an isolated disc with an
identical disc system irradiated by a strong external field, corresponding to
being within $<0.1\,$pc of an O star. They noted substantial heating and
impact upon the chemistry (e.g. the CO snow line disappears). More
intermediate/low FUV radiation fields have also been included in hydrostatic
disc models, but with a view to how it affects factors like ionisation and CO
photodesorption (e.g. Cleeves et al., 2013; Cleeves, 2016) rather than the
dust radiative equilibrium temperature. Champion et al. (2017) undertook 1D
PDR models of proplyds using the meudon code, however in those they compute
the dust temperature in the surface layers and the disc temperature is assumed
to be uniformly 19.5 K. Robberto et al. (2002) also took a similar approach
semi-analytically, considering the proplyd to be a spherical system with non-
spherical external irradiation and internal heating by the host star, but
again the focus was on the layers external to the circumstellar disc and their
infrared emission.

Sellek et al. (2020) studied the dynamical evolution of dust in discs with an
external photoevaporative wind. They accounted for the fact that only small
grains are entrained away from the disc and into a wind/proplyd envelope, as
expected from Facchini et al. (2016) and observed by Miotello et al. (2012),
Biganzoli et al. (2017). However Sellek et al. (2020) assumed a temperature
profile for the disc that is dominated by the central source. Ndugu et al.
(2018) included a parameterised heating of the outer disc in planet population
synthesis models, finding that heating due to the cluster environment is
important for suppressing large populations of cold Jupiters, which are not
observed (particularly at low metallicity). Their model included the effect of
heating on snow line locations, but the locations themselves in different
environments were not discussed.

Here we focus on understanding how external heating by massive stars affects
millimetre continuum mass estimates. Furthermore we aim to understand how this
in turn affects the statistical variation of disc properties near massive
stars and the implications that has for understanding how massive stars impact
disc evolution.

##  2 Simple optically thin radiative equilibrium model

We take a quick first look at the possible impact of an external source by
considering an extension to the classic optically thin radiative equilibrium
expression using the bolometric luminosity of the host star and external
source, similar to the approach taken by Tsamis et al. (2013) in estimating
the dust temperature of ONC proplyds.

Assume that there is a circumstellar disc that is heated to the grain
radiative equilibrium temperature by the central host star (luminosity
$L_{H})$ and an external source (luminosity $L_{ext}$). Also assume that the
disc is optically thin. The heating rate for a grain of radius $r_{d}$ and
albedo $\alpha$ at radial distance $R$ from the host star with the external
source at distance $D$ is

$H=A_{H}4\pi r_{d}^{2}\left(\frac{L_{H}}{4\pi R^{2}}+\frac{L_{ext}}{4\pi
D^{2}}\right)(1-\alpha)$ (5)

where $A_{H}$ is the fraction of the spherical grain surface that is
irradiated. In the case of a planar radiation field travelling through a
medium with zero scattering $A_{H}=1/4$ because only one side of the grain is
irradiated. However, for a more general radiation field, or a planar field
incident upon a scattering medium $A_{H}=1$. We spend the time clarifying this
because Monte Carlo radiative transfer codes (Lucy, 1999) generally compute
the averaged properties in each volume of a calculation, effectively setting
$A_{H}=1$. We assume that the heating only has a radial dependence on the
contribution from the host star and that the external heating is uniform.
Assuming both sources are blackbodies the heating term becomes

$H=A_{H}4\pi
r_{d}^{2}\sigma\left(\frac{R_{H}^{2}}{R^{2}}T_{H}^{4}+\frac{R_{ext}^{2}}{D^{2}}T_{ext}^{4}\right)(1-\alpha).$
(6)

The cooling rate of spherical grains in the disc at temperatre $T_{d}$ is

$C=4\pi r_{d}^{2}\sigma T_{d}^{4}$ (7)

Equating the heating and cooling in radiative equilibrium

$4\pi r_{d}^{2}\sigma T_{d}^{4}=A_{H}4\pi
r_{d}^{2}\sigma\left(\frac{R_{H}^{2}}{R^{2}}T_{H}^{4}+\frac{R_{ext}^{2}}{D^{2}}T_{ext}^{4}\right)\left(1-\alpha\right)$
(8)

gives a radial temperatue profile

$T_{d}=A_{H}^{1/4}\left[\left(\frac{R_{H}}{R}\right)^{2}T_{H}^{4}+\left(\frac{R_{ext}}{D}\right)^{2}T_{ext}^{4}\right]^{1/4}\left(1-\alpha\right)^{1/4}$
(9)

This is shown for a Trappist-1 type star at various distances from a star
similar to $\theta^{1}$C ($R_{ext}=10\,R_{\odot},T_{ext}=39000\,$K) in Figure
1.

For an ensemble of external sources this can be generalized to

$T_{d}=A_{H}^{1/4}\left[\left(\frac{R_{H}}{R}\right)^{2}T_{H}^{4}+\sum_{k}\left(\frac{R_{ext,k}}{D_{k}}\right)^{2}T_{ext,k}^{4}\right]^{1/4}\left(1-\alpha\right)^{1/4}$
(10)

Note that in practice most discs will not be optically thin and the albedo is
a function of wavelength and the grain size distribution.

Figure 1: Analytic disc radiative equilibrium radial temperature profiles for
a host star with $T_{H}=2550$, $R_{H}=0.121$ (Trappist-1) at various distances
from a $\theta^{1}$C analogue.

The mass estimate given by equation 4 in the Rayleigh-Jeans regime scales
inversely with the Temperature, so if assuming a temperature of 20 K the
overestimate of the disc mass is by a factor $T_{d}/20$. If the external
source dominates the bulk of the disc temperature this would hence lead to an
overestimate that scales as $D^{-1/2}$.

##  3 Radiative Transfer Models

To explore the dust temperature of externally irradiated discs in more detail
we ran Monte Carlo radiative transfer calculations using the torus code
(Harries et al., 2019). The approach used is based on that of Lucy (1999). The
code uses an adapative mesh refinement grid-based approach, with a 2D
cylindrical ($R,z$) grid geometry and the host star at the origin. The total
source luminosity (host star plus any external source) is discretised into
photon packets that undertake a physically motivated random walk over the grid
until they escape. This contributes to the energy density in each cell which
sets the dust temperature.

In these calculations, photon packets from the host star are emitted with
random direction and a frequency randomly sampled from the stellar spectrum.
When an external photon source is included, it is done so beyond the upper
bound of the grid. External photons are hence introduced in a plane parallel
fashion from the upper bound of the grid, with frequency again randomly
sampled from the external source’s spectrum. We choose this approach rather
than an isotropic external source, since many irradiated discs close to
massive stars are clearly not symmetrically irradiated. This way, we can ask
questions like whether the external radiation field affects just the upper
atmosphere, also the disc mid plane, or even the far side of the disc from the
external radiation source. In this paper we only consider a plane parallel
external field propagating parallel to the $z$ axis on the cylindrical $R,z$
grid, since there is a higher probability of the radiation impinging primarily
upon the disc surface than the disc edge for a random disc orientation. Given
that 3D models would be required for arbitrary external radiation orientations
(otherwise you end up with what is effectively a cylindrical external source)
we argue that this situation is the prudent choice.

###  3.1 Disc density distribution and dust properties

Our calculations will solve for the dust radiative equilibrium temperature,
but density and grain distributions need to be specified.

We construct the disc using a truncated power law with a surface density
profile

$\Sigma(R)=\Sigma_{1\textrm{au}}\left(\frac{R}{\textrm{au}}\right)^{-1}$ (11)

and scale height of the form

$H(R)=H_{1\textrm{au}}\left(\frac{R}{\textrm{au}}\right)^{m}$ (12)

where we choose $m=1$. The inner disc radius is 0.1 au in each case. For such
a surface density profile the mass encapsulated scales linearly with the
radius, and the surface density normalization is given by

$\Sigma_{1\textrm{au}}=\frac{M_{d}}{2\pi R_{d}\times(1\,\textrm{au})}$ (13)

where $R_{d}$ is the disc outer radius, beyond which the density is set to a
negligibly low value ($10^{-35}$ g cm-3). The density structure on the
cylindrical grid is then described by

$\rho(R,z)=\frac{\Sigma(R)}{\sqrt{2\pi}H(R)}\exp\left(-\frac{z^{2}}{2H(R)^{2}}\right).$
(14)

A summary of the parameters of discs for the models used in this paper is
given in Table 1. We include a mixture of larger discs (100 au) and smaller
ones that are more representative of continuum disc sizes in the ONC (10-50
au, Eisner et al., 2018).

For the exploration in this paper we use a single grain population. These are
Draine & Lee (1984) silicates with a minimum and maximum grain size of 0.1
$\mu$m and 2 mm respectively. The power law of the $dn(a)/da\propto a^{-q}$
distribution is $q=3.3$.

Figure 2: MIST v1.2 temperature and radius evolutionary tracks for stars of
different masses. We use the values at 1 Myr (the vertical dotted line).

We do not solve for the hydrostatic equilibrium structure of the disc in these
models. Rather, a parametric disc structure is imposed and the radiative
equilibrium temperature calculated this caveat is discussed in section 3.3. We
consider a few different values for the scale height at 1 au in the solar mass
case, from the canonical $H/R=0.1$ down to $H/R=0.025$ to represent a more
settled case. If optically thin the dust radiative equilibrium temperature of
an isolated disc scales with the host star luminosity as $T_{d}\propto
L_{*}^{1/4}$. At lower stellar masses we hence scale the imposed $H/R=0.1$ by
a factor

$f_{H}=\left(\frac{L_{*}}{L_{\odot}}\right)^{1/8}\left(\frac{M_{*}}{M_{\odot}}\right)^{-1/2}.$
(15)

We note that we do not attempt to parameterise an external envelope in these
models since the bulk of the dust remains within the disc (e.g. Sellek et al.,
2020). Dust in an envelope will have been liberated from the disc in a wind so
is just a redistribution of a small fraction of the material in our model,
where we are focusing on the disc itself. This does however mean that we
cannot model things like infrared emission, to which the proplyd envelope
makes an important contribution (e.g. Champion et al., 2017).

We place each model star-disc at 1, 0.5, 0.1, 0.05 and 0.01 pc from a
$\theta^{1}$C analogue, with radius 10$\,R_{\odot}$ and effective temperature
of 39000 K. We also include an isolated version of each model disc with no
external source.

Stellar | $T_{*}$ | $R_{*}$ | $R_{d}$ | $\Sigma_{1\textrm{au}}$ | Disc dust | $H_{1\textrm{au}}$
---|---|---|---|---|---|---
Mass ($M_{\odot}$) | (K) | ($R_{\odot}$) | (au) | g cm-2 | mass ($M_{\oplus}$) | (au)
1 | 4468 | 2.35 | 100 | 0.425 | 10.0 | 0.1
1 | 4468 | 2.35 | 100 | 0.425 | 10.0 | 0.05
1 | 4468 | 2.35 | 100 | 0.425 | 10.0 | 0.025
1 | 4468 | 2.35 | 50 | 0.851 | 10.0 | 0.1
1 | 4468 | 2.35 | 50 | 0.851 | 10.0 | 0.05
1 | 4468 | 2.35 | 50 | 0.851 | 10.0 | 0.01
1 | 4468 | 2.35 | 25 | 1.70 | 10.0 | 0.1
1 | 4468 | 2.35 | 10 | 4.25 | 10.0 | 0.1
0.5 | 3784 | 1.86 | 50 | 0.851 | 10.0 | 0.12
0.5 | 3784 | 1.86 | 25 | 1.70 | 10.0 | 0.12
0.5 | 3784 | 1.86 | 10 | 4.25 | 10.0 | 0.12
0.5 | 3784 | 1.86 | 10 | 2.12 | 5.0 | 0.12
0.5 | 3784 | 1.86 | 10 | 0.425 | 1.0 | 0.12
0.1 | 2970 | 0.98 | 50 | 0.851 | 10.0 | 0.21
0.1 | 2970 | 0.98 | 25 | 1.70 | 10.0 | 0.21
0.1 | 2970 | 0.98 | 10 | 4.25 | 10.0 | 0.21
0.1 | 2970 | 0.98 | 10 | 2.12 | 5.0 | 0.21
0.1 | 2970 | 0.98 | 10 | 0.425 | 1.0 | 0.21
Table 1: Summary of the star-disc properties of models in this paper. Each of
these models are placed at different distances from an external UV source.

###  3.2 Defining the sources

The stellar luminosity provides an additional complication because the pre-
main sequence luminosity is a different function of time for different stellar
masses and needs to be computed with a model such as those by Siess et al.
(2000) or MESA (Paxton et al., 2011; Paxton et al., 2013, 2015; Choi et al.,
2016; Dotter, 2016). In particular, before the zero age main sequence the
luminosity is a shallower function of stellar mass. For this initial study we
adopt stellar temperature and radii taken from Mesa Isochrones and Stellar
Tracks (MIST, v1.2) web interpolator
111http://waps.cfa.harvard.edu/MIST/interp_tracks.html at a time of 1Myr,
assuming Fe/H = 0 (the default, using a protosolar Fe/H=0.0142 changes both
the stellar temperature and radius by less than 0.2 per cent) and initial
v/v${}_{\textrm{crit}}$ = 0.4. The adopted stellar temperatures and radii are
shown in Figure 2 in the wider context of the stellar evolution. Table 1 also
specifies the stellar parameters used.

An argument could be made about the short lifetimes of discs in close
proximity to massive stars (the proplyd lifetime problem) but this problem is
coupled to uncertainties on the disc properties, such as mass, age and time
actually spent in the high UV environment. We therefore deem this choice of
PMS stellar properties at a single age to be pragmatic and a much better
representation of the probable reality than zero age main sequence
luminosities.

There is a further consideration required, which is how to treat the external
radiation field entering the 2D cylindrical grid. In section 3 we argued for
treating this as a plane parallel field entering from the upper bound of the
model grid. However, on a cylindrical grid the cylindrical area of a cell
centred on $R_{i}$ with thickness $\Delta R$ increases as

$\Delta A=2\pi\left[\left(R_{i}+\frac{\Delta
R}{2}\right)^{2}-\left(R_{i}-\frac{\Delta R}{2}\right)^{2}\right]$ (16)

so randomly introducing photon packets at the upper grid boundary with uniform
probability in $R$ would not correspond to a uniform planar field, but rather
one with a density per unit area that decreases with $R$. To address this, the
square root of the randomly sampled $R/R_{grid}$ is used, where $R_{grid}$ is
the outermost radial coordinate of the grid.

The energy carried per photon packet $\Delta E$ is always the same, with
different frequencies simply corresponding to different numbers of photons in
the packet. For $N$ packets and total source luminosity $L$ the energy (per
unit time) carried by a packet is $\delta E=L/N$. For the purpose of computing
the photon packet energy, the luminosity of the external source, which is
intrinsically $L_{ext,0}$, is scaled by

$L_{ext}=L_{ext,0}\frac{\pi R_{grid}^{2}}{4\pi D_{grid}^{2}}$ (17)

where $R_{grid}$ is the radial grid size and $D_{grid}$ is the distance of the
external source from the grid.

Figure 3: Points represent analytic optically thin disc radiative equilibrium
mid-plane radial temperature profiles for a host star with $T_{H}=2550$,
$R_{H}=0.121$ (Trappist-1) at various distances from a star similar to
$\theta^{1}$C. The solid lines are results from torus radiative transfer
models. The left hand panel is an extremely optically thin disc
($M_{d}=10^{-17}\,$M⊕) and the right hand panel is for a $R_{d}=$50 au,
$M_{d}=10\,$M⊕ disc.

###  3.3 Obvious caveats

Before proceeding to look at calculation results we note some of the obvious
caveats with our approach to the modelling. This is not designed to be a
detailed study of grain/disc evolution and planet formation, but rather a
first look at the effect of external irradiation on the thermal structure of
the circumstellar millimetre dust disc. To this end we consider only a single
snapshot of any given disc. We consider a distribution of grain sizes that is
the same everywhere (i.e. there is no midplane larger grain population and
diffuse atmospheric population due to settling and/or entrainment of smaller
grains in a wind). We also do not solve for hydrostatic equilibrium as these
are essentially models of the dust, and the combination of
growth/settling/drift with a photoevaporative wind that entrains small grains
means that a hydrostatic model is not necessarily better than a parametric
dust structure. However, we do explore the impact of the the scale height of
the dust which can be considered a rough proxy for settling. For the systems
in close proximity to a strong external radiation field we are assuming that
the embedded dust is still geometrically disc-like, despite there being a
larger scale proplyd structure. For this reason we do not explore synthetic
SEDs in this paper, since the warm envelope makes an important contribution to
the SED (e.g. Robberto et al., 2002; Champion et al., 2017). We will address
the dynamic evolution of the dust in future work.

Another caveat is that although we consider a realistic pre main sequence
stellar luminosity, we are not including the effects of accretion luminosity
which may occasionally increase the heating from the central source.

###  3.4 Benchmarking in the optically thin regime

We begin by comparing torus against equation 9 in the extremely optically thin
regime. We still use a disc model, as described above, only with a dust mass
of $10^{-17}$ M⊕. In this first look we use host star properties similar to
those of Trappist-1, i.e. $T_{H}=2550\,$K and $R_{H}=0.121$ R⊙ (though note
that in our later calculations we use appropriate pre main sequence stellar
properties from evolutionary models as described above) and the external
source is similar to $\theta^{1}$C ($T_{ext}$=39000 K, $R_{ext}=10\,$R⊙).

A comparison of equation 9 and the mid-plane temperature of this extremely
optically thin disc is shown in the left hand panel of Figure 3. The right
hand panel is the same setup, only with the disc mass increased to a more
typical value of $10$ M⊕. For the analytic approximation in this case we
assume a negligible albedo.

Overall the optically thin simulations give good agreement with the analytic
solution, particularly beyond 10 au. When a more realistic disc mass is used
the agreement is still good in the outer disc (dependent upon the distance of
the disc from the radiation source). As expected, in the more realistic disc
mass case the optically thick inner disc drops to a lower temperature value
than equation 9 predicts.

We also study the deviation from the simple analytic solution in our other
models in section 4.1, finding that the analytic approximation is typically
good in the regions of the disc where the external field sets the temperature,
but typically overestimates the disc temperature by a factor of around 4 where
the host star dominates. Nevertheless equation 9 could still be used to place
upper limits on the temperature and therefore corresponding lower limits on
the disc mass.

###  3.5 Testing against ONC proplyd disc temperature estimates

Tsamis et al. (2013) estimated the disc dust temperature of three proplyds in
the ONC: LVC 2, HS1 and HST10. They also provided disc mass/radius estimates
and the projected separation from $\theta^{1}$ C. We used their
disc/separation parameters, employing the host star parameters and canonical
scale height at 1 au that we use for a $0.5$ M⊙ star (see table 1) to run
comparison models of these systems.

Figure 4: Comparison of our mid-plane dust temperature structure with the disc
dust temperatures inferred by Tsamis et al. (2013) for the ONC proplyds LV2,
HST 1 and HST 10. The agreement in the outer disc where the external radiation
dominates is to within 20 per cent.

The Tsamis et al. (2013) dust temperature estimates were made by assuming that
25 per cent of the $\theta^{1}$ C bolometric luminosity reaches the disc
surface. A comparison of their inferred dust temperatures with our model mid-
plane dust temperature structure is given in Figure 4. In the outer disc where
the external field dominates the temperature the agreement is to within at
worst 20 per cent.

##  4 Results and discussion

We begin by presenting the thermal structure of externally irradiated dust
discs, before turning our attention to observables and effect upon disc dust
mass estimates.

###  4.1 Thermal structure of internally and externally irradiated discs

Figure 5 shows the dust temperature structure of one of our model discs at
different distances from a $\theta^{1}$C-like external source. This scenario
is a 0.5 M⊙ host star with a 25 au, 10 M⊕ dust disc and the result is
representative of the typical behaviour of our models. The upper left panel is
an isolated disc, so there is no external radiation source, and illustrates
the motivation for the usual assumption of a disc temperature of 20 K when
computing disc mass estimates. The other panels from left to right, top to
bottom, move the disc closer to an external radiation field, which is
impinging from the upper boundary. This results in the upper half of the disc
being warmer than the lower side, though the lower side is still heated
relative to the isolated disc case. Since we only consider a static snapshot
we do not assess the dynamical impact of this asymmetric heating here.

Figure 5: The dust radiative equilibrium temperature structure of a 10 M⊕ dust
mass, 25 au disc around a 0.5 M⊙ star at various distances from a 39000 K 10
R⊙ star. The top left panel is isolated (no external source) and in the other
panels the label in the lower right denotes the distance from the external
source which irradiates the grid from the top part of each panel, as
illustrated in the top-middle panel.

Figure 6: Mid-plane temperature profiles of externally irradiated disc models.
The upper panel has a 0.5 M⊙ host star, 25 au disc outer radius and 10 M⊕ of
dust. The lower panel has a 1 M⊙ host star, 100 au disc outer radius and 10 M⊕
of dust. The solid lines are the radiative transfer calculation results and
the dotted lines are the simple analytic radiative equilibrium approximations.

Figure 6 shows the mid-plane temperature of two of our star-disc models. In
each case we include the isolated result, as well as at distances of 1, 0.1
and 0.01 pc from the external source. The models are M${}_{*}=0.5\,$M⊙,
R${}_{d}=25$ au, $M_{d}=10\,$M⊕ (upper panel) and M${}_{*}=1\,$M⊙,
R${}_{d}=100$ au, $M_{d}=10\,$M⊕ (lower panel), and are representative of the
typical behaviour we see in all of the models. In addition to the model
results (solid lines) we also include the analytic approximation given by
equation 9 (dashed lines). For the discs nearest the external source (0.01 pc)
the analytic approximation is reasonable but there is more significant
deviation (a factor $\sim$4) in the inner parts of the disc at larger
separations when the external source doesn’t set the disc temperature.

A key point to make here is that the disc temperature is in excess of 20 K
once a disc is at a separation of less than a parsec from the $\theta^{1}$ C
like source. We will assess how this affects disc mass estimates in section
4.2, but note now that an upper limit on the temperature (which could be
provided by the analytic approximation) corresponds to a lower limit on the
inferred mass using equation 4, which is a valuable addition to an estimate
made assuming a temperature of 20 K.

####  4.1.1 Dependence of temperature structure on dust scale height

The main calculations in this paper assume a single mixed grain population
with small (0.1 $\mu$m) grains right the way through to 2 mm sized grains. In
reality the larger grains settle towards the mid-plane and some of the smaller
grains may be lifted upwards in a photoevaporative wind. To make a simple
assessment of the possible impact of these processes we ran calculations with
smaller scale heights, but otherwise identical parameters (stellar properties,
disc mass/radius).

A comparison of the mid-plane temperature structure of a canonical scale
height model and others with a factor two and four smaller scale heights is
given in Figure 7. The temperature of the outer parts of the disc heated by
the external source doesn’t change, but the radius at which the external
source dominates (in the mid-plane at least) moves inwards. This is because
the mid-plane is denser and so becomes optically thick to the host star
radiation more quickly. However, the vertical column at any given radius (and
hence to the external radiation) is the same.

There are two important points to make from this comparison. The first is that
although we are considering a single density distribution for each model,
there is only a finite amount of dust in the model so even if small grains are
elevated away from the mid-plane we do not expect it to significantly affect
the importance of external heating on the disc dust temperature (e.g. if the
system were a proplyd). The second point is that what lifting dust away from
the mid-plane may do is increase the ability of the host star to contribute to
warming the disc (and vice-versa for settling).

Figure 7: Mid-plane temperature profiles of externally irradiated disc models
with different scale heights at a distance of 0.1 pc from the external source.
A smaller scale height results in lower inner-disc temperatures, with the
external radiation field dominating the mid-plane temperature at smaller
distances from the host star. This is because the mid-plane is more optically
thick to the host-star radiation, but the overall vertical column is
unchanged.  Figure 8: The radius at which certain temperatures are reached in
the disc mid-plane. For example the CO snow line is at 22 au until about 1 pc
distance from a $\theta^{1}$C analogue, but at closer distances no CO snow
line would exist. The water snow line, of possible importance to inner planet
formation, is relatively unaffected by environment except for at very close
separations. This is for a 100 au disc with a $10\,M_{\oplus}$ dust mass and
$H/R=0.1$

####  4.1.2 Comments on the effect of external irradiation upon snow lines

Another important effect of external irradiation is the impact upon the
locations (or even existence of) snow lines (e.g. as illustrated in a single
value of the external UV field in Walsh et al., 2013). The idea that planetary
composition could link to formation locations in the disc through snow lines
(Öberg et al., 2011) has been complicated in recent years due to factors like
time dependent chemistry and pebble drift (Booth & Ilee, 2019), thermal
instability at snow lines (Owen, 2020) and the time evolution of the host star
luminosity (Miley et al., 2021). External irradiation adds yet another factor.

Here we do not study the time evolution of snow lines like e.g. Booth & Ilee
(2019) and Miley et al. (2021), but in Figure 8 we show the radius in a 100
au, 10 M⊕ disc around a 1 M⊙ star at which the temperature drops below certain
values as a function of distance from the external radiation source. The
rightmost set of points represent the thermal structure of an isolated disc
(placed at 100 pc in Figure 8). Then moving left from this rightmost point one
can follow how the location of particular mid-plane temperatures migrates
outwards.

A temperature of 20 K roughly corresponds to the temperature of the CO snow
line, 50 K is roughly the temperature of the CO2 snow line and 150 K is
roughly the temperature of the water snow line. The CO snow line moves
outwards significantly compared to an isolated disc, doing so by a factor two
even at separations of about a parsec from the O star. At closer separations,
the CO snow line ceases to exist (Walsh et al., 2013, inferred no CO snow line
for their single externally irradiated disc at what would have been
$D<0.1\,$pc). The water snow line has been proposed as an important part of
the planet forming mechanism for ensembles of planets on small orbits such as
in Trappist-1 (Ormel et al., 2017; Schoonenberg et al., 2019). Naively based
on our results we would anticipate that this would be rather resilient to the
effect of the external irradiation as the water snow line is only
significantly affected at separations of $\leq 10^{-2}$ pc. Of course even if
the water snow line location isn’t affected, the mass reservoir for planet
formation in Trappist-1 like systems could still be significantly affected by
external photoevaporation (Haworth et al., 2018a).

###  4.2 Dust mass estimates

As dicussed above, trends in disc dust masses near massive stars are assessed
using equation 4, restated here

$M_{d}=\frac{D^{2}F_{\nu}}{\kappa_{\nu}B_{\nu}(T)}$

where $D$ is the source distance, $F_{\nu}$ the measured flux, $\kappa_{\nu}$
the opacity and $B_{\nu}(T)$ the Planck function. In the Rayleigh-Jeans regime
there is a linear scaling of the Planck function with the temperature (Ansdell
et al., 2017; Eisner et al., 2018, assumed $T=20$ K).

The significant additional heating from external sources considered here could
therefore conceivably lead to overestimated disc dust masses which might
suppress any imprint of external photoevaporation on the disc mass. To test
this we produce 850 $\mu$m synthetic observations from our calculations, again
using the Monte Carlo radiative transfer approach in torus. We assume a
distance of 414 pc in our synthetic observations. In the Monte Carlo radiative
transfer calculation the opacity is as described in section 3.1. We also have
to assume an opacity to compute the dust mass. We follow Eisner et al. (2018)
and use of $\kappa_{\nu}=3.1$ cm2 g-1 at 850 $\mu$m for this.

It is important to note that the opacity is also a source of significant
uncertainty when estimating disc masses. Millimetre opacity estimates can vary
by an order of magnitude, particularly depending on whether grain growth is
accounted for, which corresponds to around an order of magnitude uncertainty
on the disc mass (Miyake & Nakagawa, 1993; Natta et al., 2004; Ricci et al.,
2010). Based on the discussion that follows below, in the ONC an order of
magnitude overestimate of the dust mass due to external heating is only
achieved at around 0.01 pc from $\theta^{1}C$, so the uncertainty in disc mass
due to the opacity is generally going to be larger than that due to external
heating. On the other hand, here we are mainly concerned with the effect of
heating on estimated disc masses as a function of projected separation from a
strong radiation source, i.e. the key thing is the variation in the opacity
between sources, not uncertainty in its actual value. In such a scenario,
where the discs are in the same region and similarly aged, it is reasonable to
expect that the degree of grain growth will be fairly similar and the
millimetre opacity will not be changing by an order of magnitude on a disc-to-
disc basis.

We do not account for any instrumentation effects such as finite beam size or
interferometric effects as we are mostly interested in the impact of assuming
a single temperature when computing the mass regardless of distance from the
external source, which is insensitive to such factors.

Figure 9 shows the estimated dust masses as a function of projected separation
from $\theta^{1}$ C from Mann et al. (2014) and Eisner et al. (2018).
Overplotted are a series of lines corresponding to the masses inferred from
our synthetic observations, showing how the inferred disc mass changes for any
given star-disc system as it is moved closer to the external source. As
expected the inferred disc mass increases nearer to the external source due to
the increased temperature, with an overestimate of the disc mass by a factor
$\sim 10$ at a separation of 0.01 pc. We also find that the overestimate
scales with distance from the external source as $D^{-1/2}$ as expected from
the discussion in section 2, and significant deviations from the 20 K
assumption start at about 1 pc (in this particular case). The $D^{-1/2}$
scaling interior to 1 pc is illustrated by the dotted line in Figure 9.

Ideally each individual source in the Mann et al. (2014) and Eisner et al.
(2018) surveys would be subject to bespoke modelling, but for a first
assessment of the possible impact of this bias we can simply scale the
observed masses by the overestimate factor inferred from our models
$f=\left(D/\textrm{pc}\right)^{-1/2}$, which we do in Figure 10. The upper
panel shows the original data and the lower panel accounting for the
temperature-induced overestimate of the disc mass by multiplying by the factor
$f$. The dotted and solid lines approximate the moving average in each
original and scaled data respectively. Of course there is uncertainty as to
the true separation (since all we plot in Figure 10 is the projected
separation), but our results suggest that a trend in disc mass with separation
may be being suppressed.

It is worth noting that a trend in disc dust masses as a function of projected
separation is not necessarily expected in an external photoevaporation
scenario. The dynamical evolution of a star cluster on such small spatial
scales is relatively fast, so any such signature might not survive very long.
On the other hand it may also be quickly imprinted.

Figure 9: Inferred disc 850 $\mu$m dust masses as a function of projected
separation from $\theta^{1}C$. The blue points are detections (we ignore upper
limits here) from Mann et al. (2014) and the orange points are detections from
Eisner et al. (2018). The solid coloured lines are torus radiative equilibrium
synthetic 850 $\mu$m models, with tracks representing an identical star-disc,
just at different distances from the external radiation source. As the star-
disc approaches the massive external source the mass of the dust disc is
overestimated when assuming a temperature of 20 K, since the true dust
temperature is in fact higher. The dotted line denotes a $D^{-1/2}$ scaling.

Figure 10: Inferred disc dust masses as a function of projected separation
from $\theta 1C$. The upper panel is a combination of the detections (we
ignore upper limits here) from Mann et al. (2014) (blue points) and Eisner et
al. (2018) (orange points). The lower panel is scaled by the
$(D/\textrm{pc})^{0.5}$ factor by which our models suggests the masses are
being overestimated. The dotted and solid lines approximate the moving average
of the upper and lower distributions respectively.

####  4.2.1 Effect on other surveys and observations

Here we have focused our attention on the ONC and $\theta^{1}$ C. Another
notable example of disc dust masses as a function of projected separation is
in the older $\sigma$ Orionis cluster by Ansdell et al. (2017), where an
increase in disc dust mass as a function of projected separation from $\sigma$
Ori was observed (with a stronger correlation than observed in the ONC). The
distance $D$ from a single source at which a disc would transition to being
externally heated at a radius $R$ in the disc is

$D=R\left(\frac{L_{ext}}{L_{H}}\right)^{1/2}$ (18)

where $L_{ext}$ and $L_{H}$ are the external and host star luminosities.

The correlation of disc mass with separation in $\sigma$ Orionis was at
separations of more like beyond a parsec from $\sigma$ Ori and the luminosity
of $\sigma$ Ori is also about a factor 5 weaker than $\theta^{1}$ C. We hence
expect any overestimate of disc masses as a function of separation due to
external heating plays a much weaker role in that region than in the ONC. van
Terwisga et al. (2020) also surveyed disc masses with ALMA towards NGC 2024 at
projected separations from the O8V star IRS 2b comparable to the separations
of discs considered from $\theta^{1}$ C in the ONC (down to 0.01 pc) by Eisner
et al. (2018). IRS 2b is also around a factor 5 less luminous than
$\theta^{1}$ C, however the square root dependency on this in equation 18
means that external heating should be important for disc heating/mass
estimates out to separations of around 0.5 pc in that region.

Outside of regions like those discussed above within around a parsec of
massive stars we do not expect that it is important to consider external
heating of the dust disc for mass estimates, especially given the order of
magnitude uncertainty on the opacity (see section 4.2 and Miyake & Nakagawa,
1993; Natta et al., 2004; Ricci et al., 2010). However it is also worth noting
that the bolometric luminosity has a shallower decline with stellar mass than
the FUV, so depending on the cluster density the aggregate of lower mass close
neighbours may also lead to substantial heating beyond the reaches of massive
stars. Stellar luminosities and 3D positions using gaia distances (e.g. Gaia
Collaboration et al., 2016a, b, 2018) could be combined with equation 10 (or
equation 18) to make a first assessment of this.

####  4.2.2 Effect of disc inclination and position of the external radiation
source

Estimating disc dust masses using equation 4 assumes that the emission is
optically thin, in which case there should not be a sensitivity to
inclination. If the disc were optically thick then the situation becomes
substantially more complicated. On the one hand the disc mass may be being
underestimated if it is optically thick. The disc mass inferred would then
also be sensitive to the disc inclination, since if the warmer irradiated side
were presented to the observer the emission would be brighter and (if a
temperature of 20 K were assumed) some of the mass estimate deficit from being
optically thick accounted for. Conversely if the cold side were presented the
inferred mass would be lower. For simplicity it is best to assume that the
disc is generally optically thin and consider the inclination/optical depth
effects when undertaking bespoke modelling of individual discs.

###  4.3 Possible effect on grain evolution and planet formation

So far we have focused on the effect that the increased disc dust temperature
has on continuum mass estimates and commented on the effect on snow lines. It
is also worth briefly considering the possible impact of higher disc dust
temperatures and a flatter radial temperature profile on grain evolution.

Firstly and most simply, our results indirectly support the work of Ndugu et
al. (2018), where higher outer disc temperatures were imposed in planet
population synthesis models. They found that this outer heating is important
for suppressing populations of cold Jupiters, particularly at low metallicity,
which are not observed.

We can also consider the impact of externally heated disc temperature
structures on grain growth and drift following the simple model of Birnstiel
et al. (2012). For example the drift timescale and maximum grain size before
the onset of radial drift both scale inversely with the square of the sound
speed and log pressure gradient
($\left|\textrm{d}\ln(P)/\textrm{d}\ln(R)\right|$). For a heated irradiated
disc both of these quantities will be decreased by a factor
$T_{\textrm{isolated}}/T_{\textrm{irradiated}}$ due to the change in sound
speed. Since the temperature profile in the externally heated disc is much
flatter, the reduced log pressure gradient will also act to slightly further
reduce the drift timescale and maximum grain size for drift. So overall drift
may happen to smaller grains, and more rapidly. Whether this helps to promote
or further hinder planet formation (in addition to the photoevaporative
depletion of the disc) remains to be explored.

Another factor is that the Stokes number (grain size) for fragmentation of
grains through collisions scales inversely with the temperature, so the Stokes
number/size for fragmentation will be decreased by a factor
$T_{\textrm{isolated}}/T_{\textrm{irradiated}}$ in a warmed irradiated disc.
Once again, this could act to further suppress planet formation in tandem with
the removal of material by external photoevaporation.

Finally is worth noting that if the disc temperature distribution is flatter
(as it is in the bulk of an irradiated disc) then it is easier to introduce
pressure bumps through density perturbations which could lead to dust
trapping.

It is important to note here that the total fraction planet hosts for which
the dust discs would have been exposed to external heating at some point in
their lifetime is uncertain and cannot be assessed by counting the
instantaneous fraction of stars in nearby regions within $\sim$1 pc of O
stars. Stellar clusters are dynamically evolving systems, with ongoing star
formation and stars moving in and out of high radiation parts of a cluster
(e.g. Winter et al., 2019). In addition, most stars formed at around a
redshift of 2 (Madau & Dickinson, 2014) when star forming regions may
typically have resembled more extreme/massive clusters with strong radiation
fields (e.g. Adamo et al., 2020).

##  5 Summary and conclusions

We use Monte Carlo radiative transfer models to study the radiative
equilibrium structure of protoplanetary discs irradiated by both the host star
and an external source. In particular we are interested in how external
irradiation affects trends in disc mass estimates as a function of separation
from massive stars compared to assuming a constant temperature (as is often
the case in recent ALMA continuum surveys) since these are used to infer the
impact of external photoevaporation on disc evolution. We draw the following
main conclusions from this work.

1) The majority of the dust component of discs in the vicinity ($<1\,$pc) of
massive stars can be heated by the external radiation field to well in excess
of the 20 K recently assumed in ALMA millimetre continuum surveys for
estimating disc masses in these high radiation environments. This leads to
disc masses being overestimated when assuming a 20 K disc, since warmer discs
are brighter. The overestimate scales with the separation from the external
source as roughly $D^{-1/2}$, until the disc is sufficiently far from the
external source that the host star dominates the temperature structure. In the
vicinity of $\theta^{1}$C in the ONC the external field starts to play an
important role at a separation of about a parsec and results in roughly a
factor 10 overestimate on the dust mass at separations of 0.01 pc.

2) Applying a simple $D^{-1/2}$ scaling to observed disc dust masses as a
function of projected separation from $\theta^{1}$C (at separations less than
1 pc, as motivated by our models) results in a significantly stronger
variation in disc mass with projected separation. Though of course this has
the caveat that projected separation is not necessarily the true separation.
Bespoke SED modelling of ONC discs on a case-by-case basis would be required
for a more rigorous assessment of the impact of external heating on disc mass
estimates.

3) External heating from nearby massive stars can affect mid-plane snow line
locations, which further complicates the idea of being able to relate planet
compositions to formation radii. The CO snow line is quite readily removed
completely from discs near massive stars, whereas the water snow line (which
may be important for inner planet formation) is resilient to cluster heating
except for very small separations ($\leq 0.01$ pc) from a $\theta^{1}$C type
external source. The fraction of stars subject to this kind of heating is is
yet to be determined.

4) Generally when estimating the mass of any given disc the uncertainty on the
opacity will dominate over that due to external heating. The main importance
of external heating is when searching for trends in disc properties in the
vicinity of massive stars and this paper provides the tools for making a first
assessment of this.

##  6 Data availability

The torus code is publicly available222https://bitbucket.org/tjharries/torus/.
The observational data included is available directly from Mann et al. (2014)
and Eisner et al. (2018). All plotting scripts and any other information are
available on request.

## Acknowledgements

We thank the referee for their careful assessment of this paper. We also thank
Andrew Sellek, Andrew Winter, Josh Eisner, Nick Ballering and Will Henney for
useful comments and discussions. TJH is funded by a Royal Society Dorothy
Hodgkin Fellowship. This work used the DiRAC@Durham facility managed by the
Institute for Computational Cosmology on behalf of the STFC DiRAC HPC Facility
(www.dirac.ac.uk). The equipment was funded by BEIS capital funding via STFC
capital grants ST/P002293/1, ST/R002371/1 and ST/S002502/1, Durham University
and STFC operations grant ST/R000832/1. DiRAC is part of the National
e-Infrastructure.

## References

  * Adamo et al. (2020) Adamo A., et al., 2020, Space Sci. Rev., 216, 69
  * Andre & Montmerle (1994) Andre P., Montmerle T., 1994, ApJ, 420, 837
  * Andrews & Williams (2005) Andrews S. M., Williams J. P., 2005, ApJ, 631, 1134
  * Andrews & Williams (2007) Andrews S. M., Williams J. P., 2007, ApJ, 671, 1800
  * Ansdell et al. (2017) Ansdell M., Williams J. P., Manara C. F., Miotello A., Facchini S., van der Marel N., Testi L., van Dishoeck E. F., 2017, AJ, 153, 240
  * Ansdell et al. (2020) Ansdell M., et al., 2020, AJ, 160, 248
  * Beckwith et al. (1990) Beckwith S. V. W., Sargent A. I., Chini R. S., Guesten R., 1990, AJ, 99, 924
  * Biganzoli et al. (2017) Biganzoli D., Potenza M. A. C., Robberto M., 2017, ApJ, 840, 55
  * Birnstiel et al. (2012) Birnstiel T., Klahr H., Ercolano B., 2012, A&A, 539, A148
  * Booth & Ilee (2019) Booth R. A., Ilee J. D., 2019, MNRAS, 487, 3998
  * Boyden & Eisner (2020) Boyden R. D., Eisner J. A., 2020, ApJ, 894, 74
  * Champion et al. (2017) Champion J., Berné O., Vicente S., Kamp I., Le Petit F., Gusdorf A., Joblin C., Goicoechea J. R., 2017, A&A, 604, A69
  * Choi et al. (2016) Choi J., Dotter A., Conroy C., Cantiello M., Paxton B., Johnson B. D., 2016, ApJ, 823, 102
  * Cleeves (2016) Cleeves L. I., 2016, ApJ, 816, L21
  * Cleeves et al. (2013) Cleeves L. I., Adams F. C., Bergin E. A., 2013, ApJ, 772, 5
  * Concha-Ramírez et al. (2019) Concha-Ramírez F., Wilhelm M. J. C., Portegies Zwart S., Haworth T. J., 2019, MNRAS, 490, 5678
  * Dotter (2016) Dotter A., 2016, ApJS, 222, 8
  * Draine & Lee (1984) Draine B. T., Lee H. M., 1984, ApJ, 285, 89
  * Eisner et al. (2018) Eisner J. A., et al., 2018, ApJ, 860, 77
  * Facchini et al. (2016) Facchini S., Clarke C. J., Bisbas T. G., 2016, MNRAS, 457, 3593
  * Gaia Collaboration et al. (2016a) Gaia Collaboration et al., 2016a, A&A, 595, A1
  * Gaia Collaboration et al. (2016b) Gaia Collaboration et al., 2016b, A&A, 595, A2
  * Gaia Collaboration et al. (2018) Gaia Collaboration et al., 2018, A&A, 616, A1
  * García-Arredondo et al. (2001) García-Arredondo F., Henney W. J., Arthur S. J., 2001, ApJ, 561, 830
  * Harries et al. (2019) Harries T. J., Haworth T. J., Acreman D., Ali A., Douglas T., 2019, Astronomy and Computing, 27, 63
  * Haworth & Clarke (2019) Haworth T. J., Clarke C. J., 2019, MNRAS, 485, 3895
  * Haworth et al. (2018a) Haworth T. J., Facchini S., Clarke C. J., Mohanty S., 2018a, MNRAS, 475, 5460
  * Haworth et al. (2018b) Haworth T. J., Clarke C. J., Rahman W., Winter A. J., Facchini S., 2018b, MNRAS, 481, 452
  * Haworth et al. (2021) Haworth T. J., Kim J. S., Winter A. J., Hines D. C., Clarke C. J., Sellek A. D., Ballabio G., Stapelfeldt K. R., 2021, MNRAS, 501, 3502
  * Henney & Arthur (1998) Henney W. J., Arthur S. J., 1998, AJ, 116, 322
  * Henney & O’Dell (1999) Henney W. J., O’Dell C. R., 1999, AJ, 118, 2350
  * Hildebrand (1983) Hildebrand R. H., 1983, QJRAS, 24, 267
  * Kim et al. (2016) Kim J. S., Clarke C. J., Fang M., Facchini S., 2016, ApJ, 826, L15
  * Kruijssen et al. (2020) Kruijssen J. M. D., Longmore S. N., Chevance M., 2020, ApJ, 905, L18
  * Krumholz et al. (2019) Krumholz M. R., McKee C. F., Bland-Hawthorn J., 2019, ARA&A, 57, 227
  * Kurtovic et al. (2018) Kurtovic N. T., et al., 2018, ApJ, 869, L44
  * Lada & Lada (2003) Lada C. J., Lada E. A., 2003, ARA&A, 41, 57
  * Lichtenberg et al. (2019) Lichtenberg T., Golabek G. J., Burn R., Meyer M. R., Alibert Y., Gerya T. V., Mordasini C., 2019, Nature Astronomy, 3, 307
  * Longmore et al. (2021) Longmore S. N., Chevance M., Kruijssen J. M. D., 2021, arXiv e-prints, p. arXiv:2103.01974
  * Lucy (1999) Lucy L. B., 1999, A&A, 344, 282
  * Madau & Dickinson (2014) Madau P., Dickinson M., 2014, ARA&A, 52, 415
  * Mann et al. (2014) Mann R. K., et al., 2014, ApJ, 784, 82
  * Miley et al. (2021) Miley J. M., Panić O., Booth R. A., Ilee J. D., Ida S., Kunitomo M., 2021, MNRAS, 500, 4658
  * Miotello et al. (2012) Miotello A., Robberto M., Potenza M. A. C., Ricci L., 2012, ApJ, 757, 78
  * Miyake & Nakagawa (1993) Miyake K., Nakagawa Y., 1993, Icarus, 106, 20
  * Natta et al. (2004) Natta A., Testi L., Neri R., Shepherd D. S., Wilner D. J., 2004, A&A, 416, 179
  * Ndugu et al. (2018) Ndugu N., Bitsch B., Jurua E., 2018, MNRAS, 474, 886
  * Nicholson et al. (2019) Nicholson R. B., Parker R. J., Church R. P., Davies M. B., Fearon N. M., Walton S. R. J., 2019, MNRAS, 485, 4893
  * Öberg et al. (2011) Öberg K. I., Murray-Clay R., Bergin E. A., 2011, ApJ, 743, L16
  * O’dell & Wen (1994) O’dell C. R., Wen Z., 1994, ApJ, 436, 194
  * Ormel et al. (2017) Ormel C. W., Liu B., Schoonenberg D., 2017, A&A, 604, A1
  * Owen (2020) Owen J. E., 2020, MNRAS, 495, 3160
  * Paxton et al. (2011) Paxton B., Bildsten L., Dotter A., Herwig F., Lesaffre P., Timmes F., 2011, ApJS, 192, 3
  * Paxton et al. (2013) Paxton B., et al., 2013, ApJS, 208, 4
  * Paxton et al. (2015) Paxton B., et al., 2015, ApJS, 220, 15
  * Reiter (2020) Reiter M., 2020, A&A, 644, L1
  * Ricci et al. (2010) Ricci L., Testi L., Natta A., Neri R., Cabrit S., Herczeg G. J., 2010, A&A, 512, A15
  * Richling & Yorke (2000) Richling S., Yorke H. W., 2000, ApJ, 539, 258
  * Robberto et al. (2002) Robberto M., Beckwith S. V. W., Panagia N., 2002, ApJ, 578, 897
  * Rodriguez et al. (2018) Rodriguez J. E., et al., 2018, ApJ, 859, 150
  * Scally & Clarke (2001) Scally A., Clarke C., 2001, MNRAS, 325, 449
  * Schoonenberg et al. (2019) Schoonenberg D., Liu B., Ormel C. W., Dorn C., 2019, A&A, 627, A149
  * Sellek et al. (2020) Sellek A. D., Booth R. A., Clarke C. J., 2020, MNRAS, 492, 1279
  * Siess et al. (2000) Siess L., Dufour E., Forestini M., 2000, A&A, 358, 593
  * Tazzari et al. (2017) Tazzari M., et al., 2017, A&A, 606, A88
  * Tsamis et al. (2013) Tsamis Y. G., Flores-Fajardo N., Henney W. J., Walsh J. R., Mesa-Delgado A., 2013, MNRAS, 430, 3406
  * Walsh et al. (2013) Walsh C., Millar T. J., Nomura H., 2013, ApJ, 766, L23
  * Winter et al. (2018) Winter A. J., Clarke C. J., Rosotti G., Ih J., Facchini S., Haworth T. J., 2018, MNRAS, 478, 2700
  * Winter et al. (2019) Winter A. J., Clarke C. J., Rosotti G. P., Hacar A., Alexander R., 2019, MNRAS, 490, 5478
  * Winter et al. (2020) Winter A. J., Kruijssen J. M. D., Longmore S. N., Chevance M., 2020, Nature, 586, 528
  * van Terwisga et al. (2020) van Terwisga S. E., et al., 2020, A&A, 640, A27