* [1] E. Abbe, J. Fan, K. Wang, and Y. Zhong.  Entrywise eigenvector analysis of random matrices with low expected rank.  Annals of Statistics, 48(3):1452–1474, 2020.
  * [2] E. Abbe and C. Sandon.  Community detection in general stochastic block models: Fundamental limits and efficient algorithms for recovery.  In 2015 IEEE 56th Annual Symposium on Foundations of Computer Science, pages 670–688. IEEE, 2015.
  * [3] N. Alon and N. Kahale.  A spectral technique for coloring random 3-colorable graphs.  SIAM Journal on Computing, 26(6):1733–1748, 1997.
  * [4] N. Alon, M. Krivelevich, and B. Sudakov.  Finding a large hidden clique in a random graph.  In Proceedings of the Ninth Annual ACM-SIAM Symposium on Discrete Algorithms, SODA ’98, page 594?598, USA, 1998. Society for Industrial and Applied Mathematics.
  * [5] A. Blum and J. Spencer.  Coloring random and semi-random k-colorable graphs.  Journal of Algorithms, 19(2):204–234, 1995.
  * [6] R. B. Boppana.  Eigenvalues and graph bisection: An average-case analysis.  In 28th Annual Symposium on Foundations of Computer Science (SFCS 1987), pages 280–285, 1987.
  * [7] A. Coja-Oghlan.  A spectral heuristic for bisecting random graphs.  Random Structures & Algorithms, 29(3):351–398, 2006.
  * [8] S. Dhara, J. Gaudio, E. Mossel, and C. Sandon.  Spectral algorithms optimally recover (censored) planted dense subgraphs.  arXiv preprint arXiv:2203.11847, 2022.
  * [9] S. Dhara, J. Gaudio, E. Mossel, and C. Sandon.  Spectral recovery of binary censored block models.  Proceedings of the 2022 Annual ACM-SIAM Symposium on Discrete Algorithms (SODA), 2022.
  * [10] D. Gross.  Recovering low-rank matrices from few coefficients in any basis.  IEEE Transactions on Information Theory, 57(3):1548–1566, 2011\.
  * [11] B. Hajek, Y. Wu, and J. Xu.  Achieving exact cluster recovery threshold via semidefinite programming.  IEEE Transactions on Information Theory, 62(5):2788–2797, 2016\.
  * [12] S. Janson, T. Łuczak, and A. Rucinski.  Random Graphs. Wiley, New York, 2000.
  * [13] M. Krivelevich and B. Sudakov.  Pseudo-random graphs.  In More sets, graphs and numbers, pages 199–262. Springer, 2006\.
  * [14] N. Ryu and S.-Y. Chung.  Community detection with colored edges.  In 2016 IEEE International Symposium on Information Theory (ISIT), pages 1899–1903. IEEE, 2016.
  * [15] W. Schudy and M. Sviridenko.  Concentration and moment inequalities for polynomials of independent random variables.  In Proceedings of the twenty-third annual ACM-SIAM symposium on Discrete Algorithms, pages 437–446. SIAM, 2012.
  * [16] V. Vu.  A simple SVD algorithm for finding hidden partitions.  Combinatorics, Probability and Computing, 27(1):124?140, 2018.
  * [17] S. Yun and A. Proutière.  Accurate community detection in the stochastic block model via spectral algorithms.  arXiv, abs/1412.7335, 2014.

##  Appendix A Proof of Poisson approximation

###### Proof.

_(Proof of Lemma 2.5)._ Observe that $(D_{a,r},D_{b,r})$ are independent over
$r$ as they depend on disjoint sets of independent random variables. Thus,

$\begin{split}{exact-deg=prob}&\mathbb{P}(D=d)\\\
&=\prod_{r\in[n]\setminus\\{i\\}}\frac{|S_{r}|!}{d_{1,r}!d_{2,r}!(|S_{r}|-d_{1,r}-d_{2,r})!}(\alpha\psi_{r})^{d_{1,r}}(\alpha(1-\psi_{r}))^{d_{2,r}}(1-\alpha)^{|S_{r}|-d_{1,r}-d_{2,r}}\\\
&\quad\times\frac{(|S_{i}|-|V|)!}{d_{1,i}!d_{2,i}!(|S_{i}|-|V|-d_{1,i}-d_{2,i})!}(\alpha\psi_{i})^{d_{1,i}}(\alpha(1-\psi_{i}))^{d_{2,i}}(1-\alpha)^{|S_{i}|-|V|-d_{1,i}-d_{2,i}}.\end{split}$
(A.1)

We use Stirling’s approximation and the fact that $1-\mathrm{e}^{-x}\asymp x$
as $x\to 0$. Thus, using the assumptions on $d$, for each $r\neq i$, we have

$\displaystyle\frac{|S_{r}|!}{(|S_{r}|-d_{1,r}-d_{2,r})!}$
$\displaystyle\asymp\frac{\mathrm{e}^{-|S_{r}|}|S_{r}|^{|S_{r}|+\frac{1}{2}}}{\mathrm{e}^{-|S_{r}|+d_{1,r}+d_{2,r}}(|S_{r}|-d_{1,r}-d_{2,r})^{|S_{r}|-d_{1,r}-d_{2,r}+\frac{1}{2}}}$
$\displaystyle=\mathrm{e}^{-d_{1,r}-d_{2,r}}\frac{(|S_{r}|-d_{1,r}-d_{2,r})^{d_{1,r}+d_{2,r}}}{(1-\frac{d_{1,r}+d_{2,r}}{|S_{r}|})^{|S_{r}|+\frac{1}{2}}}$
$\displaystyle\asymp\frac{(|S_{r}|-d_{1,r}-d_{2,r})^{d_{1,r}+d_{2,r}}}{(1-\frac{d_{1,r}+d_{2,r}}{|S_{r}|})^{\frac{1}{2}}}$
$\displaystyle\asymp\big{(}n\rho_{r}\big{(}1+O(\log^{-2}n)\big{)}\big{)}^{d_{1,r}+d_{2,r}}\asymp(n\rho_{r})^{d_{1,r}+d_{2,r}},$

where in the last step we have used $d_{1,r},d_{2,r}=o(\log^{3/2}n)$. Also,

$(1-\alpha)^{|S_{r}|-d_{1,r}-d_{2,r}}\asymp\mathrm{e}^{-\alpha(|S_{r}|-d_{1,r}-d_{2,r})}\asymp\mathrm{e}^{-t\rho_{r}\log(n)/2}.$

Thus, the $r$-th product term in (LABEL:exact-deg=prob) is asymptotically
equal to

$\displaystyle\mathrm{e}^{-t\rho_{r}\log n}\frac{(\rho_{r}\psi_{r}t\log
n)^{d_{1,r}}(\rho_{r}(1-\psi_{r})t\log n)^{d_{2,r}}}{d_{1,r}!d_{2,r}!}.$

The identical approximation holds for $r=i$ as well using the fact that
$|S_{i}\setminus V|=n\rho_{i}(1+O(\log^{-2}n))$. Thus the proof follows from
(LABEL:exact-deg=prob). ∎

##  Appendix B Proof of Lemma 3.3

We first state a special case of [15, Theorem 1.3].

######  Lemma B.1.

Let $X_{1},\dots,X_{n}$ be independent random variables, and let
$X=\sum_{i=1}^{n}X_{i}$. Let $L>0$. Suppose that for each $i\in[n]$ and
$k\in\mathbb{N}$,

$\mathbb{E}\left[|X_{i}|^{k}\right]\leq k\cdot
L\cdot\mathbb{E}\left[|X_{i}|^{k-1}\right].$ (B.1)

Let $\mu_{0}=\sum_{i=1}^{n}\mathbb{E}[|X_{i}|]$. Then, for any $\lambda>0$,

$\mathbb{P}\left(\left|X-\mathbb{E}[X]\right|\geq\lambda\right)\leq\mathrm{e}^{2}\max\left\\{\exp\left(-\frac{\lambda^{2}}{\mu_{0}LR}\right),\exp\left(-\frac{\lambda}{LR}\right)\right\\},$
(B.2)

where $R$ is an absolute constant.

###### Proof of Lemma 3.3.

We first verify (B.1), for $L=\max\\{|x|:x\in S\\}$. For $k\in\mathbb{N}$,

$\displaystyle\mathbb{E}\left[|X_{i}|^{k}\right]$ $\displaystyle=\sum_{x\in
S}|x|^{k}\mathbb{P}(X_{i}=x)$ $\displaystyle\leq L\sum_{x\in
S}|x|^{k-1}\mathbb{P}(X_{i}=x)$
$\displaystyle=L\mathbb{E}\left[|X_{i}|^{k-1}\right]$ $\displaystyle\leq
kL\mathbb{E}\left[|X_{i}|^{k-1}\right].$

We apply Lemma B.1, noting that $\mu_{0}=\mathbb{E}[Y]$. Set
$\lambda=\delta|\mathbb{E}[X]|$. Then $\lambda<\mu_{0}$, so that the first
bound in (B.2) applies, giving the claim. ∎

##  Appendix C Proof of eigenvalue properties

###### Proof of Lemma 5.1.

The only $2\times 2$ matrices whose eigenvalues are not distinct are the
multiples of the identity matrix. Indeed,

$\displaystyle\det\begin{pmatrix}a&b\\\
b&c\end{pmatrix}=0\implies\lambda^{2}-(a+c)\lambda+ac-b^{2}=0,$

which has same roots in $\lambda$ if and only if $(a-c)^{2}+b^{2}=0$.
Therefore, $A^{\prime}$ in (5.2) has identical eigenvalues if and only if
$y=q/(1-q)$.

Next, note that $\det(A^{\prime}(y))$ is a quadratic function in $y$, which
has at most two roots unless $\det(A^{\prime}(y))$ is the zero polynomial. To
rule out the latter possibility, note that

$\det\bigg{(}A^{\prime}\Big{(}\frac{q}{1-q}\Big{)}\bigg{)}=\rho(1-\rho)\Big{(}p_{1}-\frac{q}{1-q}(1-p_{1})\Big{)}\Big{(}p_{2}-\frac{q}{1-q}(1-p_{2})\Big{)},$

which is nonzero if $p_{1}$ and $p_{2}$ are both different from $q$. When
$p_{2}=q$, then $\det(A^{\prime}(0))=\rho(1-\rho)(p_{1}q-q^{2})$, which cannot
be zero due to our assumption that $0<p_{1}.p_{2},q<1$ cannot be all be the
same. Hence, $\det(A^{\prime}(y))$ cannot be a zero polynomial, and thus there
are at most two values of $y$ such that $\det(A^{\prime}(y))=0$. Combined with
the condition for having distinct roots, there is a set $\mathcal{Y}$ with
$|\mathcal{Y}|\leq 3$ such that, for $y\notin\mathcal{Y}$, $A^{\prime}(y)$ has
two distinct and nonzero eigenvalues. This proves Lemma 5.1 Part (1).

To prove Lemma 5.1 Parts (2),(3), consider the case where (5.3) holds. Also,
take
$y=y(p,q)=\log\big{(}\frac{1-q}{1-p}\big{)}/\log\big{(}\frac{p}{q}\big{)}$.
Due to symmetry of $y(p,q)$, (5.13) and (5.14) imply

$y(p,q)>\frac{q}{1-q}\text{ if }p>q\quad\text{and}\quad
y(p,q)<\frac{q}{1-q}\text{ if }p<q$

Also, $y(p,q)>0$ for all $0<p,q<1$. Thus, if (5.3), then
$y(p,q)\neq\frac{q}{1-q}$. Hence, the off diagonal entries of
$A^{\prime}(y(p,q))$ are nonzero and the eigenvalues of $A^{\prime}(y(p,q))$
are distinct. If $p_{2}=q$ then
$\det(A^{\prime}(y(p,q)))=\rho(1-\rho)(p-q)(1+y)(q-y(1-q))\neq 0$, so its
eigenvalues are both nonzero. This proves Part (2).

To prove Part (3), take $p_{1}=p_{2}=p$. Note that $\det(A^{\prime}(y))=0$ is
$0$ if and only if $|p-y(1-p)|=|q-y(1-q)|$. Since $p\neq q$, the latter holds
if and only if $y=\frac{p+q}{2-p-q}$ or $p+q=1$. Let $x=(p+q)/2$ and
$\epsilon=(p-q)/2$, and observe that

$\displaystyle
y(p,q)=\log\Big{(}\frac{1-q}{1-p_{1}}\Big{)}/\log\Big{(}\frac{p}{q}\Big{)}$
$\displaystyle=\log\left(\frac{1-x+\epsilon}{1-x-\epsilon}\right)/\log\left(\frac{x+\epsilon}{x-\epsilon}\right)$
$\displaystyle=\frac{\log(1+\epsilon/(1-x))-\log(1-\epsilon/(1-x))}{\log(1+\epsilon/x)-\log(1-\epsilon/x)}$
$\displaystyle=\frac{\sum_{i=0}^{\infty}2(\epsilon/(1-x))^{2i+1}/(2i+1)}{\sum_{i=0}^{\infty}2(\epsilon/x)^{2i+1}/(2i+1)}$
$\displaystyle=\frac{x}{1-x}\cdot\frac{\sum_{i=0}^{\infty}(\epsilon/(1-x))^{2i}/(2i+1)}{\sum_{i=0}^{\infty}(\epsilon/x)^{2i}/(2i+1)}.$

The ratio of infinite sums is greater than $1$ if $x<1/2$ and less than $1$ if
$x>1/2$, so $y=\frac{p+q}{2-p-q}\iff 2x=p+q=1$. This proves Part (3). ∎

###### Proof of Lemma 6.1.

First, observe that $P^{\scriptscriptstyle(y)}$ has a zero eigenvalue if and
only if its determinant is zero. Since $\det(P^{\scriptscriptstyle(y)})$ is a
polynomial in $y$ of degree at most $k$, either it is identically zero (i.e.,
(1) holds) or there exists a subset $\mathcal{Y}_{1}\subset\mathbb{R}$ with
$|\mathcal{Y}_{1}|\leq k$ such that $\det(P^{\scriptscriptstyle(y)})\neq 0$
for all $y\notin\mathcal{Y}_{1}$.

Next, observe that for any given $y$, the eigenvalues of
$P^{\scriptscriptstyle(y)}$ are the roots of the characteristic polynomial
$\chi^{\scriptscriptstyle(y)}(\lambda):=\det(P^{(y)}-\lambda I)$. Let
$f^{\scriptscriptstyle(y)}(\lambda)$ be the polynomial with leading
coefficient $1$ in $\lambda$ that is the greatest common divisor of
$\chi^{\scriptscriptstyle(y)}(\lambda)$ and
$(\chi^{\scriptscriptstyle(y)})^{\prime}(\lambda)=\frac{\mathrm{d}}{\mathrm{d}\lambda}\chi^{\scriptscriptstyle(y)}(\lambda)$.
Then $\chi^{\scriptscriptstyle(y)}(\lambda)$ has repeated roots in $\lambda$
if and only if $f^{\scriptscriptstyle(y)}(\lambda)$ is not a constant function
in $\lambda$. Now, consider $\chi^{\scriptscriptstyle(y)}(\lambda)$ and
$\frac{\mathrm{d}}{\mathrm{d}\lambda}\chi^{\scriptscriptstyle(y)}(\lambda)$ as
elements of $\mathbb{R}^{\scriptscriptstyle(y)}[\lambda]$, the ring of
polynomials in $\lambda$ with coefficients that are rational functions of $y$.
Then, there exist
$f^{\star\scriptscriptstyle(y)},g_{1}^{\scriptscriptstyle(y)},g^{\scriptscriptstyle(y)}_{2},h_{1}^{\scriptscriptstyle(y)},h_{2}^{\scriptscriptstyle(y)}\in\mathbb{R}^{\scriptscriptstyle(y)}[\lambda]$
such that the leading coefficient of $f^{\star\scriptscriptstyle(y)}$ is $1$,
and

$\displaystyle
f^{\star\scriptscriptstyle(y)}=g_{1}^{\scriptscriptstyle(y)}\chi^{\scriptscriptstyle(y)}+g_{2}^{\scriptscriptstyle(y)}(\chi^{\scriptscriptstyle(y)})^{\prime},\quad\chi^{\scriptscriptstyle(y)}=h_{1}^{\scriptscriptstyle(y)}f^{\star\scriptscriptstyle(y)},\quad(\chi^{\scriptscriptstyle(y)})^{\prime}=h_{2}^{\scriptscriptstyle(y)}f^{\star\scriptscriptstyle(y)}.$

Thus, for any $y$, $f^{\star\scriptscriptstyle(y)}$ will evaluate to
$f^{\scriptscriptstyle(y)}$, unless the denominator of at least one
coefficient of $f^{\star\scriptscriptstyle(y)}$,
$g_{1}^{\scriptscriptstyle(y)}$, $g_{2}^{\scriptscriptstyle(y)}$,
$h_{1}^{\scriptscriptstyle(y)}$, or $h_{2}^{\scriptscriptstyle(y)}$ evaluates
to $0$. Since the coefficients are rational functions in $y$, this can only
happen for $y\in\mathcal{Y}_{2}$, where $\mathcal{Y}_{2}$ is a finite set.
Therefore, if $f^{\star\scriptscriptstyle(y)}(\lambda)$ is a constant in
$\lambda$, then for all $y\notin\mathcal{Y}_{2}$, the eigenvalues of
$P^{\scriptscriptstyle(y)}$ are all distinct. Taking
$\mathcal{Y}=\mathcal{Y}_{1}\cup\mathcal{Y}_{2}$, we have shown the (3) holds
in this case.

Next, suppose that $f^{\star\scriptscriptstyle(y)}(\lambda)$ is not constant
in $\lambda$. Then, $P^{\scriptscriptstyle(y)}$ must have a repeated
eigenvalue for all $y\notin\mathcal{Y}_{2}$. The eigenvalues of
$P^{\scriptscriptstyle(y)}$ change continuously as functions of $y$. Thus, if
there was any $y$ for which its eigenvalues were all distinct they would have
to be distinct for all values of $y$ sufficiently close to that one.
Therefore, $P^{\scriptscriptstyle(y)}$ must have repeated eigenvalues for all
values of $y$ and (2) holds in this case. This completes the proof of Lemma
6.1. ∎


# Photoevaporation versus enrichment in the cradle of the Sun

Miti Patel1,2, Cheyenne K. M. Polius1⋆, Matthew Ridsdill-Smith1,3⋆, Tim
Lichtenberg4 and Richard J. Parker1
1Department of Physics and Astronomy, The University of Sheffield, Hicks
Building, Hounsfield Road, Sheffield, S3 7RH, UK
2School of Physics & Astronomy, University of Leicester, University Road,
Leicester, LE1 7RH, UK
3National Astronomical Research Institute of Thailand, 260 Moo 4, T. Donkaew,
A. Maerim, Chiangmai 50180, Thailand
4Kapteyn Astronomical Institute, University of Groningen, PO Box 800, 9700 AV
Groningen, Netherlands Authors made an equal contribution and each should be
considered joint lead authorE-mail<EMAIL_ADDRESS>Society
Dorothy Hodgkin Fellow

###### Abstract

The presence of short-lived radioisotopes (SLRs) 26Al and 60Fe in the Solar
system places constraints on the initial conditions of our planetary system.
Most theories posit that the origin of 26Al and 60Fe is in the interiors of
massive stars, and they are either delivered directly to the protosolar disc
from the winds and supernovae of the massive stars, or indirectly via a
sequential star formation event. However, massive stars that produce SLRs also
emit photoionising far and extreme ultraviolet radiation, which can destroy
the gas component of protoplanetary discs, possibly precluding the formation
of gas giant planets like Jupiter and Saturn. Here, we perfom $N$-body
simulations of star-forming regions and determine whether discs that are
enriched in SLRs can retain enough gas to form Jovian planets. We find that
discs are enriched and survive the photoionising radiation only when the dust
radius of the disc is fixed and not allowed to move inwards due to the
photoevaporation, or outwards due to viscous spreading. Even in this optimal
scenario, not enough discs survive until the supernovae of the massive stars
and so have zero or very little enrichment in 60Fe. We therefore suggest that
the delivery of SLRs to the Solar system may not come from the winds and
supernovae of massive stars.

###### keywords:

methods: numerical – protoplanetary discs – photodissociation region (PDR) –
open clusters and associations: general.

††pagerange: Photoevaporation versus enrichment in the cradle of the
Sun–Photoevaporation versus enrichment in the cradle of the Sun††pubyear: 2023

##  1 Introduction

Understanding the origin of our Solar system is one of the crucial unsolved
problems in astrophysics, as constraining the initial conditions of our
planetary system will also enable greater understanding of the many thousands
of extrasolar planets.

It is clear that planetary systems form very early on in a star’s life (Haisch
et al., 2001; Najita & Kenyon, 2014; Richert et al., 2018; Alves et al., 2020;
Segura-Cox et al., 2020), when the star is still embedded in its natal
molecular cloud in the company of tens to thousands of stellar siblings (Lada
& Lada, 2003).

Therefore, the star-forming environment likely governs the early evolution of
planetary systems. Our Solar system is no exception. The decay products of
short-lived radioisotopes 26Al and 60Fe in the chondritic meteorites suggest
that the Sun either formed from material enriched in these isotopes (Gounelle
& Meynet, 2012; Gounelle, 2015; Young, 2014; Boss, 2013, 2017), or the
isotopes were delivered to the Solar system through e.g. the supernovae and/or
winds of nearby massive stars (Ouellette et al., 2010; Lichtenberg et al.,
2016; Fatuzzo & Adams, 2015, 2022; Parker et al., 2023).

However, whilst massive stars have been shown to enrich the early Solar system
(and probably many other planetary systems), they have also been shown to have
a very destructive influence on protoplanetary discs, as their far and extreme
ultraviolet radiation acts to evaporate the gas content from protoplanetary
discs (Scally & Clarke, 2001; Adams et al., 2004; Nicholson et al., 2019).
Photoevaporation of protoplanetary discs is a very efficient process, and
occurs even in star-forming regions with relatively low stellar densities
(Parker et al., 2021).

This creates an interesting tension; could the Solar system be enriched in
short-lived radioactive isotopes, but retain enough of its gaseous
protoplanetary disc such that the gas giant planets (Jupiter and Saturn) would
still be able to form? Lichtenberg et al. (2019) found that the heating from
26Al determines the long-term water budget of a planet, with Solar system-
levels of 26Al resulting in liquid water oceans, and lower levels leading to
ice worlds. Gilmour & Middleton (2009) even argue that intelligent
civilisations require heating from 26Al in order to form and evolve.

In this paper, we follow the enrichment of the Solar system through stellar
winds and supernovae, but also track the depletion of gas in the protosolar
disc due to external photoevaporation from the same massive stars that provide
the enriching SLRs.

Many studies have shown that SLRs can be delivered to the protosolar disc via
supernovae (e.g. Ouellette et al., 2007, 2010; Parker et al., 2014a;
Lichtenberg et al., 2016; Nicholson & Parker, 2017), and recently Parker et
al. (2023) demonstrated that enrichment from the winds of massive stars can
enhance (and often dominate) the delivery of SLRs to protoplanetary discs.

However, none of the above studies (with the exception of Lichtenberg et al.,
2016, who implemented viscous spreading) followed the evolution of the disc,
especially under the influence of external photoevaporation, which numerous
studies have shown to have a highly detrimental effect on gaseous discs in
star-forming regions (Scally & Clarke, 2002; Adams et al., 2004; Haworth et
al., 2018b; Winter et al., 2018; Nicholson et al., 2019; Concha-Ramírez et
al., 2021; Parker et al., 2021). Portegies Zwart (2019) tracked the enrichment
of stars from winds and supernovae, but concentrated on interactions with
other low-mass stars that would sculpt the Edgeworth-Kuiper Belt and/or Oort
Cloud, rather than photoevaporation of the protoplanetary discs.

The paper is organised as follows. In Section 2 we describe the set-up of the
$N$-body simulations and our post-processing analysis to calculate the
evolution and enrichment of the protoplanetary discs. In Section 3 we present
our results, and we provide a discussion in Section 4. We conclude in Section
5.

##  2 Method

In this Section we describe the set-up of our numerical $N$-body simulations,
and the post-processing analysis used to calculate the evolution of
protoplanetary discs and the enrichment from massive star winds and
supernovae.

###  2.1 Star-forming regions

Observations of young star-forming regions (Cartwright & Whitworth, 2004;
Peretto et al., 2006; Sánchez & Alfaro, 2009; Hacar et al., 2013; Foster et
al., 2015; Buckner et al., 2019) show that pre-main sequence stars form in
filamentary structures, resulting in a spatially and kinematically
substructured distribution of the stars. Similarly, hydrodynamic simulations
that follow the conversion of gas to stars also produce spatial and kinematic
substructure (Bate, 2012; Dale & Bonnell, 2012), although Parker & Dale (2015)
show that the relation between the gas and stars is non-trivial.

To mimic the spatially and kinematically substructured distributions in young
star-forming regions, we set up our $N$-body simulations as self-similar,
fractal distributions (Goodwin et al., 2004). We use the box-fractal method,
which has the advantage that the amount of substructure is defined by just one
number, the fractal dimension, $D$.

The spatial distribution is set up by first defining a cube with side of
length $N_{\rm div}=2$, and a ‘parent’ star is placed at the centre of this
cube. The cube is then divided into $N^{3}_{\rm div}$ sub-cubes and a ‘child’
particle is placed at the centre of each sub-cube. The probability that a
child particle becomes a parent particle is $N_{\rm div}^{(D-3)}$, and child
particles that do not become parents are removed, along with their parent
stars.

Each child’s sub-cube is itself then divided into $N^{3}_{\rm div}$ sub-sub-
cubes, and the process is repeated until there is a generation of particles
with significantly more stars than required.

All of the parental generation particles are removed, so that only the final
generation is retained. The distribution is then pruned so that the star
particles lie within the boundary of a sphere, rather than the cube we used to
create the fractal. If there are more stars remaining than the desired number
of stars, $N_{\star}$ then stars are removed until we reach $N_{\star}$.

The stellar velocities follow the spatial substructure such that stars that
are close together have similar velocities, whereas more distant stars can
have very different velocities, similar to the observed Larson (1981)
relations.

During the construction of the fractal, the parent particles are assigned a
velocity drawn randomly from a Gaussian distribution of mean zero. The child
particles inherit their parent’s velocities, plus a small random offset that
decreases with the number of generations from the parent particle.

The velocities are then scaled to a global virial ratio, $\alpha_{\rm
vir}=T/|\Omega|$, where $T$ and $|\Omega|$ are the total kinetic and potential
energies of the stars, respectively, and virial equilibrium is defined as
$\alpha_{\rm vir}=0.5$.

In all of our simulations we adopt a fractal dimension $D=2.0$, which results
in a moderately substructured distribution.

Similarly, all of our simulations are set up to be subvirial, i.e.
$\alpha_{\rm vir}=0.3$, which results in a global collapse of the star-forming
region, erasing the substructure and forming a smooth, centrally concentrated
star cluster (Parker et al., 2014b).

In our simulations, we wish to maximise the number of massive stars that could
enrich protoplanetary discs with SLRs, but also produce FUV and EUV radiation
fields that could cause the evaporation of gas from protoplanetary discs.
Star-forming regions follow a mass distribution of the form (Lada & Lada,
2003)

$\frac{dN}{dM}\propto M_{\rm SF}^{-2},$ (1)

where $M_{\rm SF}$ is the mass of a star-forming region. Eq. 1 implies that
low-mass star-forming regions are much more common than high-mass star-forming
regions. However, for a typical Initial Mass Function we expect there to be
more massive stars when the mass of the star-forming region is higher
(Elmegreen, 2006; Weidner & Kroupa, 2006; Parker & Goodwin, 2007; Nicholson &
Parker, 2017).

We choose a total number of stars, $N_{\star}=1500$, such that our star-
forming regions contain enough massive stars for photoevaporation and
enrichment of the discs, but not so massive that the star-forming regions
would be extremely rare in the Galaxy.

We draw $N_{\star}=1500$ stars from a Maschberger (2013) IMF, which has a
probability density distribution of the form

$p(m)\propto\left(\frac{m}{\mu}\right)^{-\alpha}\left(1+\left(\frac{m}{\mu}\right)^{1-\alpha}\right)^{-\beta}.$
(2)

In Eq. 2, $\mu=0.2$ M⊙ is the average stellar mass, $\alpha=2.3$ is the high
mass Salpeter (1955) slope, and $\beta=1.4$ is used to describe the slope of
the low-mass end of the IMF. We sample the IMF in the range 0.1 – 50 M⊙.

Our choice for the upper limit to the IMF is informed by observations of star-
forming regions with a similar number of stars to those we model here, such as
the ONC, where the most massive star, $\theta^{1}$ Ori C, is $\sim$40 M⊙,
(Lehmann et al., 2010). Some authors posit that the most massive star that can
form in a region is fundamentally related to the mass of the region (Weidner &
Kroupa, 2006), whereas others suggest the only upper limit to the mass of a
star that can form is the total mass of the region itself (Elmegreen, 2006).
In the latter scenario, the most massive star can in theory be as massive as
the region (although this is unlikely in random sampling), and can also be
significantly lower than the most massive stellar mass from the Weidner &
Kroupa (2006) relation. The most massive stars known in the Universe are
$\sim$300 M⊙ (Crowther et al., 2010), but we restrict our IMF to 50 M⊙ as a
compromise between the constrained most massive star–region mass relation
(Weidner & Kroupa, 2006) and the expected outcome of randomly sampling the IMF
(Parker & Goodwin, 2007).

We keep the number of stars $N_{\star}$ constant, we keep the fractal
dimension $D=2.0$ constant, and vary the initial radii of the star-forming
region. We adopt radii of 1 pc, 2.5 pc and 5.5 pc, which correspond to initial
median stellar densities of $1000$ M⊙ pc-3 (“high density”), 100 M⊙ pc-3
(“moderate density”) and 10 M⊙ pc-3 (“low density”), respectively.

The star-forming regions are evolved as a pure $N$-body simulation for 10 Myr
using the kira integrator within the Starlab package (Portegies Zwart et al.,
1999, 2001). We assume the stars all form at the same time, and we do not
include stellar evolution in the simulations (although we use the models for
Limongi & Chieffi (2018) to determine the SLR yields, see Section 2.3). For
simplicity, we do not include primordial stellar binaries in the simulations,
although these are likely to be ubiquitous in star-forming regions (Raghavan
et al., 2010; Chen et al., 2013; Duchêne & Kraus, 2013; Pineda et al., 2015;
Ward-Duong et al., 2015).

###  2.2 Disc evolution

We follow the enrichment and destruction of the protoplanetary discs via a
post-processing analysis, i.e. the discs are not included in the $N$-body
simulations. We use the $N$-body simulation output (every 0.01 Myr) to track
the relative positions of the massive stars compared to the low-mass (disc-
hosting) stars, and then use the distances to the massive stars to determine
the amount of enrichment in SLRs, and the amount of mass lost from the discs
due to EUV and FUV radiation from massive stars (this has a time-step of
$10^{-4}$ Myr to accruately capture the disc evolution, Parker et al., 2021).

Each star in the mass range 0.1 – 3.0 M⊙ is assigned a disc, and we set the
initial disc masses to be 10 per cent of the host star’s mass,

$M_{\rm disc}=0.1\,m_{\star}.$ (3)

The initial disc radii are the same in each simulation, but we run different
sets of the same simulation’s initial conditions but with disc radii of either
$r_{\rm disc}=10$ or 100 au.

We first calculate how the discs respond to the FUV and EUV radiation fields
from massive stars. Haworth et al. (2018a) show that very little dust is lost
from discs due to the FUV and EUV fields, but significant amounts of gas are
destroyed, depending on the strength of the radiation field.

To estimate the FUV field incident on each low-mass star, we determine the
flux of FUV radiation $F_{\rm FUV}$ using the FUV luminosity $L_{\rm FUV}$
from Armitage (2000), which in turn is derived from the stellar atmosphere
models of Buser & Kurucz (1992) and Schaller et al. (1992):

$F_{\rm FUV}=\frac{L_{\rm FUV}}{4\pi d^{2}}.$ (4)

We use the distance to each massive star, $d$, to calculate the flux and then
sum the fluxes, as all simulations contain more than one massive star (defined
as $>$5 M⊙ for the purposes of generating an FUV radiation field). We convert
the FUV flux from cgs units into the dimensionless Habing (1968) unit,
$G_{0}=1.8\times 10^{-3}$ erg s-1 cm-2, which is the ambient background FUV
flux in the interstellar medium.

To calculate the mass lost due to the FUV radiation field, we use the FRIED
grid from Haworth et al. (2018b), which is a grid of mass-loss rates for
different combinations of stellar mass, $G_{0}$, disc mass, disc radius and
disc surface density. For each low mass star, we use the stellar mass, disc
mass, disc radius and $G_{0}$ field to establish via linear interpolation the
most appropriate mass-loss rate due to FUV radiation, $\dot{M}_{\rm FUV}$.

We calculate the mass-loss due to the EUV radiation field using Eq. 5 from
Johnstone et al. (1998), which depends on the distance to the massive star $d$
and the disc radius $r_{\rm disc}$:

$\dot{M}_{\rm EUV}\simeq 8\times 10^{-12}r^{3/2}_{\rm
disc}\sqrt{\frac{\Phi_{i}}{d^{2}}}\,\,{\rm M_{\odot}\,yr}^{-1}.$ (5)

$\Phi_{i}$ is the ionizing EUV photon luminosity from each massive star in
units of $10^{49}$ s-1 and is dependent on the stellar mass according to the
observations of Vacca et al. (1996) and Sternberg et al. (2003).

We subtract mass from the discs every $10^{-4}$ Myr assuming the FUV and EUV-
induced mass-loss rates from the FRIED models and Eq. 5. Models of mass-loss
from discs usually assume material is lost from the edge of the disc (where
the surface density is lowest) and so we would expect the radius to decrease
in this scenario.

Previous work (e.g. Haworth et al., 2018b; Haworth & Clarke, 2019; Parker et
al., 2021) assumes that if the surface density in the inner regions of the
disc (e.g. 1 au from the star) remains constant, then the radius is reduced
linearly in the presence of photoevaporative mass-loss, thus:

$\Sigma_{\rm 1\,au}=\frac{M_{\rm disc}}{2\pi r_{\rm disc}[{\rm 1\,au}]},$ (6)

where $M_{\rm disc}$ is the total disc mass, and $r_{\rm disc}$ is the disc
radius, as before, and $\Sigma_{\rm 1\,au}$ is the surface density at 1 au. We
then recalculate the disc radius after mass-loss, $r_{\rm disc}(t_{k})$ using

$r_{\rm disc}(t_{k})=\frac{M_{\rm disc}(t_{k})}{M_{\rm disc}(t_{k-1})}r_{\rm
disc}(t_{k-1}).$ (7)

The radius after mass-loss is therefore a function of radius before mass-loss,
$r_{\rm disc}(t_{k-1})$ multiplied by the ratio of the disc mass after
($M_{\rm disc}(t_{k})$) and before ($M_{\rm disc}(t_{k-1})$) mass-loss. The
subscript k refers to the (arbitrary) time after mass-loss, and the subscript
k-1 refers to the time before mass-loss.

As well as the inward evolution caused by mass-loss, the discs also evolve due
to viscous evolution, which causes the radius to increase. We implement
viscous expansion in our discs by implementing the diffusion equation (Lynden-
Bell & Pringle, 1974; Pringle, 1981) with the parameterisation given in
Hartmann et al. (1998) and Hartmann (2009).

In this formulation, the surface density $\Sigma$ at a given radius $R$ is

$\Sigma=1.4\times 10^{3}\frac{{\rm
e}^{-R/(R_{1}t_{d})}}{(R/R_{1})t_{d}^{3/2}}\left(\frac{M_{\rm disc}(0)}{\rm
0.1M_{\odot}}\right)\left(\frac{R_{1}}{\rm 10\,au}\right)^{-2}{\rm
g\,cm}^{-2},$ (8)

where $M_{\rm disc}(0)$ is the disc mass before viscous evolution and $R_{1}$
is a radial scaling factor, which we set as $R_{1}=10$ au. $t_{d}$ is a non-
dimensional time, such that at a given physical time $t$

$t_{d}=1+\frac{t}{t_{s}},$ (9)

and the viscous timescale, $t_{s}$ is given by

$t_{s}=8\times 10^{4}\left(\frac{R_{1}}{\rm
10\,au}\right)\left(\frac{\alpha}{10^{-2}}\right)^{-1}\left(\frac{M_{\star}}{\rm
0.5\,M_{\odot}}\right)^{1/2}\left(\frac{T_{100\,{\rm au}}}{\rm
10\,K}\right)^{-1}{\rm yr}.$ (10)

In Eq. 10, $\alpha$ is the disc viscosity (Shakura & Sunyaev, 1973) and
$T_{\rm 100\,au}$ is the temperature of the disc at a distance of 100 au from
the host star. We assume a disc temperature profile of the form

$T(R)=T_{\rm 1\,au}R^{-q},$ (11)

where $T_{\rm 1\,au}$ is the temperature at 1 au from the host star and is
derived from the stellar luminosity for pre-main sequence stars by Luhman et
al. (2003); Luhman (2004) and Kirk & Myers (2011). We adopt $q=0.5$ and
$\alpha=0.01$ (Hartmann et al., 1998). Recent observations (see e.g. Manara et
al., 2022, and references therein) suggest that $\alpha<<10^{-2}$, which would
result in little viscous expansion of the discs in the timeframe of the
simulations. We explore the consequences of this by running sets of
simulations in which we keep the dust radius fixed.

Given the mass of the star, we calculate the temperature at 100 au and then
calculate the viscous timescale $t_{s}$. We then use this to calculate the
surface density $\Sigma$ as a function of radius $R$. We set the truncation
radius, $R_{\rm trunc}$ as the radius $R$ where the surface density falls
below $10^{-6}$g cm-2. After viscous expansion, the position of the truncation
radius changes, and we use the truncation radius before ($R_{\rm
trunc}(t_{n-1})$) and after ($R_{\rm trunc}(t_{n})$) to calculate the change
to the disc radius after viscous expansion:

$r_{\rm disc}(t_{n})=r_{\rm disc}(t_{n-1})\frac{R_{\rm trunc}(t_{n})}{R_{\rm
trunc}(t_{n-1})}.$ (12)

Note that the subscripts n and n-1 differ from the subscripts used in Eq. 7 to
denote inward evolution of the disc radius due to photoevaporation because we
cannot model these two physical processes at exactly the same time. We
calculate the change in disc radius due to photoevaporation _before_ any
change in radius due to viscous evolution.

Formally the mass-loss due to photoevaporation is the gas component of the
disc, but simulations have shown that the dust radius evolves in step with the
gas radius (Sellek et al., 2020). To test the effect of this assumption, we
will perform analysis where the dust radius is kept constant, and then where
the dust radius moves with the gas radius.

###  2.3 Disc enrichment

To estimate the amount of enrichment received by a protoplanetary disc we need
to determine the amount of 26Al and 60Fe deposited from the winds and
supernovae of massive stars, and how this compares to the amount of the
respective stable isotopes, 27Al and 56Fe.

We assume the usual gas-to-dust ratio of 100:1 for Solar metallicity stars, so
the total dust mass, $m_{\rm dust}$, in the disc is

$m_{\rm dust}=0.01M_{\rm disc}.$ (13)

The amount of 26Al is expressed in terms of the stable version of the element,
27Al, and we determine the mass of 27Al in the disc, $m_{{}^{27}{\rm Al}}$,
from the fraction of Al in chondrites, $f_{\rm Al,CI}=8500\times 10^{-6}$
(Lodders, 2003)

$m_{{}^{27}{\rm Al}}=8500\times 10^{-6}m_{\rm dust}.$ (14)

Similarly, the amount of 60Fe is expressed in terms of its stable version,
56Fe, and we determine the mass of 56Fe in the disc, $m_{{}^{56}{\rm Fe}}$,
from the fraction of Fe in chondrites, $f_{\rm Fe,CI}=1828\times 10^{-4}$
(Lodders, 2003)

$m_{{}^{56}{\rm Fe}}=1828\times 10^{-4}m_{\rm dust}.$ (15)

For each $<$3 M⊙ star at each snapshot in time, we calculate the geometric
cross section of the disc for capturing material from a supernova,

$\eta_{\rm SN}=\frac{\pi r_{\rm disc}^{2}}{4\pi d^{2}}{\rm cos}\,\theta,$ (16)

where $d$ is the distance to the massive star(s) at the instant of a supernova
and $\theta$ is the inclination of the disc. Following Lichtenberg et al.
(2016), for each star we adopt $\theta=60^{\circ}$ as the likely average
inclination between the disc and the ejecta for a random distribution.

The supernovae times for each massive star are calculated by summing the main,
and post-main sequence timescales provided in table 5 from Limongi & Chieffi
(2018). We assume the massive stars are Solar metallicity, and are rotating at
300 km s-1 (de Mink et al., 2013) – see Parker et al. (2023) for a comparison
with models where the massive stars are not rotating. For our adopted IMF mass
range Limongi & Chieffi (2018) provide data for 13, 15, 20, 25, 40 and 60 M⊙
stars, and we perform a linear interpolation for stars whose mass lies in
between the values in the data table. In the rotating models, stars with
initial mass 25 M⊙ explode as supernovae after 9.9 Myr, whereas in the non-
rotating models the supernovae of 25 M⊙ stars occur at 7.7 Myr. In the models
of Limongi & Chieffi (2018), stars more massive than 25 M⊙ collapse directly
to a black hole and do not explode as supernovae.

Table 1: A summary of the different initial conditions of our simulated star-forming regions. The columns show the initial radius of the star-forming region, $r_{F}$, the initial median local stellar density, $\tilde{\rho}$, the disc evolution physics and the stellar evolution model. The final column indicates whether the simulation is shown in a Figure in Section 3.  $r_{F}$ | $\tilde{\rho}$ | $r_{{\rm disc},\,i}$ | Disc evolution | Stellar evolution | Fig.
---|---|---|---|---|---
1 pc |  $1000$ M⊙ pc-3 | 100 au | Photoevap. & dust fixed | LC18, Rotating (300 km s-1) |
1 pc |  $1000$ M⊙ pc-3 | 100 au | Photoevap. & dust fixed | LC18, non-rotating (0 km s-1) |
2.5 pc |  $100$ M⊙ pc-3 | 100 au | Photoevap. & dust fixed | LC18, Rotating (300 km s-1) | Figs. 1 & 2
2.5 pc |  $100$ M⊙ pc-3 | 100 au | Photoevap. & dust fixed | LC18, non-rotating (0 km s-1) |
1 pc |  $1000$ M⊙ pc-3 | 100 au | Photoevap. & dust follows gas | LC18, Rotating (300 km s-1) |
2.5 pc |  $100$ M⊙ pc-3 | 100 au | Photoevap. & dust follows gas | LC18, Rotating (300 km s-1) | Fig. 3
1 pc |  $1000$ M⊙ pc-3 | 10 au | Photoevap. & viscous | LC18, Rotating (300 km s-1) |
1 pc |  $1000$ M⊙ pc-3 | 10 au | Photoevap. & viscous | LC18, non-rotating (0 km s-1) |
1 pc |  $1000$ M⊙ pc-3 | 100 au | Photoevap. & viscous | LC18, Rotating (300 km s-1) |
1 pc |  $1000$ M⊙ pc-3 | 100 au | Photoevap. & viscous | LC18, Non-rotating (0 km s-1) |
2.5 pc |  $100$ M⊙ pc-3 | 100 au | Photoevap. & viscous | LC18, Rotating (300 km s-1) |
2.5 pc |  $100$ M⊙ pc-3 | 10 au | Photoevap. & viscous | LC18, Rotating (300 km s-1) | Fig. 4
2.5 pc |  $100$ M⊙ pc-3 | 10 au | Photoevap. & viscous | LC18, Non-rotating (0 km s-1) |
5.5 pc |  $10$ M⊙ pc-3 | 100 au | Photoevap. & viscous | LC18, Rotating (300 km s-1) |

Because we adopt the rotating models (in line with observations of massive
stars), the majority of the 26Al originates in the winds of the massive stars
_before_ they explode as supernovae. In their tables 8 and 9, Limongi &
Chieffi (2018) provide the total yields for different isotopes (table 8) and
then the yields solely from the winds of the massive stars (table 9). We use
this information to determine the yield from supernovae only, and the yield
from the wind over the lifetime of the massive star(s).

We assume that the isotope production in the wind is linear with time, and we
use this to obtain a yield per time, which we then use to calculate the amount
of 26Al and 60Fe produced by each massive star at a given time in the
simulation. In reality, the mass-loss from the stars is unlikely to be linear,
and we are probably overestimating the production of 26Al at early times.
Similarly, we do not account for the radioactive decay of the 26Al and
therefore our results can be thought of as an upper limit to the amount of
enrichment that can be obtained from the massive stars’ winds.

The cross section for capture of wind material is calculated as the volume of
material swept out by a low-mass star as it traverses a distance $\Delta
r_{\star}$ through a wind bubble of radius $r_{\rm bub}$

$\eta_{\rm wind}=\frac{3}{4}\frac{\pi r_{\rm disc}^{2}\Delta r_{\star}}{\pi
r_{\rm bub}^{3}}.$ (17)

We implement two different regimes for the density of the wind bubbles. First,
we assume a very compact bubble around each massive star with a radius $r_{\rm
bub}=0.1$ pc. Secondly, we assume the bubble(s) disperse rapidly, and the
bubble has a radius $r_{\rm bub}=2r_{1/2}$, where $r_{1/2}$ is the half-mass
radius of the star-forming region (i.e. the radius within which half the total
stellar mass in the region is enclosed).

We use the total mass of 26Al and 60Fe captured/swept up by each low mass star
and divide this by the mass in stable isotopes, $m_{{}^{27}{\rm Al}}$ and
$m_{{}^{56}{\rm Fe}}$, to determine the 26Al and 60Fe ratios:

$Z_{\rm Al}=\frac{m_{{}^{26}{\rm Al}}}{m_{{}^{27}{\rm Al}}},$ (18)

and

$Z_{\rm Fe}=\frac{m_{{}^{60}{\rm Fe}}}{m_{{}^{56}{\rm Fe}}}.$ (19)

For each low-mass star, we calculate $Z_{\rm Al}$ three times; for supernovae
only, for supernovae and local winds ($r_{\rm bub}=0.1$ pc), and then finally
for supernovae and dispersed winds ($r_{\rm bub}=2r_{1/2}$). The contribution
of 60Fe comes from supernovae only.

Limongi & Chieffi (2018) provide different yields of 26Al and 60Fe depending
on whether the massive stars are rotating at 0, 150 or 300 km s-1. For the
non-rotating models (0 km s-1) massive stars explode as supernovae earlier,
but less material is emitted in their winds. In contrast, the fast rotating
stars have longer lifetimes and explode as supernovae later, but more material
is emitted in their winds.

A summary of the different simulations (and disc physics) is given in Table 1.

##  3 Results

In this Section we will present three sets of simulations, each with a
different treatment of the protoplanetary disc physics. We will then discuss
the effects of varying other initial conditions, such as the stellar density
and the stellar evolution models used to calculate the amount of 26Al and 60Fe
delivered to the protoplanetary discs.

###  3.1 Photoevaporation with fixed dust radius, $r_{\rm disc,i}=100$ au

The models in Parker et al. (2023) assumed the disc radius remained constant
throughout the duration of the simulation. In this first set of models from
our initial conditions, we adjust the gas radius of the disc according to Eqn.
7, and we also record the reduction in the disc’s (gas) mass. However, the
dust radius (and mass) is assumed to remain constant, and we therefore
calculate the $Z_{\rm Al}=^{26}$Al/27Al and $Z_{\rm Fe}=^{60}$Fe/56Fe ratios
assuming a dust radius $r_{\rm disc}=100$ au.

Figure 1: The abundances of 26Al (panel a) and 60Fe (panel b) after 10 Myr in
our moderately dense ($\tilde{\rho}\sim 100$ M⊙ pc-3) simulations where the
disc loses gas mass, and the gas radius moves inwards, but the dust radius is
fixed to the initial value ($r_{\rm disc}=100$ au). Each line represents the
distribution of low-mass ($\leq$3 M⊙) disc-hosting stars in each separate
simulated star-forming region. Only three simulations contain supernovae that
explode before 10 Myr (the orange, dark blue and cyan lines in both panels).
The 26Al enrichment shown in panel (a) is from supernovae and localised
stellar winds. In seven simulations, a significant fraction (shown where the
coloured lines meet the y-axis in panel (a)) do not experience any enrichment.
The vertical dotted line in panel (a) shows the measured 26Al abundance in the
Solar system Thrane et al. (2006). The vertical dot-dashed line in panel (b)
shows the (lower) estimate for the measured 60Fe abundance in the Solar system
(Tang & Dauphas, 2012), whereas the vertical dashed line shows the (higher)
estimate for the measured 60Fe abundance in the Solar system (Mishra et al.,
2016).

We first show the expected enrichment distributions for these simulations in
Fig. 1. In panel (a) we show the $Z_{\rm Al}$ ratios in these simulations,
where we assume the enrichment comes from both supernovae and stellar winds,
where the winds are restricted to local ($<0.1$ pc) bubbles around the massive
star(s). Each line is the $Z_{\rm Al}$ distribution in a separate run of the
same simulation initial conditions, and comprises all of the low-mass ($\leq
3$ M⊙) stars in each simulation. Due to the stochastic sampling of the IMF,
the total number of low-mass (disc-hosting) stars varies between simulations,
but is usually higher than 1450 out of the total number of 1500 stars in each
simulation. The vertical dotted line is the measured value for the Solar
system ($Z_{\rm Al,SS}=5.85\times 10^{-5}$, Thrane et al., 2006)

Supernovae explode before 10 Myr in only three out of ten simulations. This
may seem low, but the Limongi & Chieffi (2018) models we adopt assume stellar
rotation and that massive stars $>$25 M⊙ collapse directly to form a black
hole. The stellar rotation delays the supernovae, and stars less massive than
25 M⊙ do not explode until very late on in the simulations.

In panel (b) we show the enrichment distributions for $Z_{\rm Fe}$, where the
enrichment comes from supernovae only. Here, almost no systems attain the
levels of enrichment commensurate with the high estimate for the Solar system
($10^{-6}$, Mishra et al., 2016; Cook et al., 2021, the right vertical dashed
line in Fig. 1); and around 25 per cent of systems experience enrichment
commensurate with the low estimate for the Solar system ($10^{-8}$, Tang &
Dauphas, 2012; Trappitsch et al., 2018, the left vertical dashed line in Fig.
1).

We now determine that – for the systems that are enriched – how many of their
discs survive the photoevaporation from the same massive stars that are
enriching the discs.

The star-forming regions have moderate initial stellar densities
($\tilde{\rho}\sim 100$ M⊙ pc-3) and radiation fields of $10^{3}-10^{4}G_{0}$
(Parker et al., 2021). The fraction of discs that lose all of their gas in
these simulations is $\sim$60 per cent after 2.5 Myr and 90 per cent after 10
Myr.

In Fig. 2 we show the disc gas properties (radius and mass) against the Al and
Fe abundances. In all panels we show the values at 2.5 Myr by the blue points,
and at 10 Myr by the black points. Panels (a)–(c) show the disc gas mass as a
function of the 26Al abundance. Panel (a) shows the abundances from supernovae
only, and therefore there are no data at 2.5 Myr because the massive stars
have yet to explode.

Panels (b) and (c) show the gas mass as a function of abundance from both
supernovae and stellar winds, either assuming the winds are dispersed in the
star-forming regions (panel b) or concentrated in local bubbles (panel c). In
both cases, there are more data points at 2.5 Myr because fewer discs have
been destroyed by the EUV/FUV fields.

Interestingly, the enrichment from supernovae alone (panel a), and supernovae
and dispersed winds (panel b) does not reach the levels measured in the Solar
system ($Z_{\rm Al}\leq 10^{-5}$, compared to the Solar system value of
$Z_{\rm Al,SS}=5.85\times 10^{-5}$).

Figure 2: Properties of discs as a function of the amount of enrichment in a
moderately dense star-forming region ($\tilde{\rho}\sim 100$ M⊙ pc${}^{-3})$
where the discs all have initial radii $r_{\rm disc,i}=100$ au. The disc gas
radii are evolved inwards due to photoevaporation (and the discs lose gas-
mass), but the dust radii (and masses) remain constant. Panels (a)–(c) show
the gas mass as a function of the amount of 26Al enrichment (expressed as the
26Al/27Al ratio). Panel (a) shows enrichment from supernovae only, panel (b)
shows enrichment from supernovae and dispersed wind ejecta and panel (c) shows
enrichment from supernovae and localised wind ejecta. Panels (d)–(f) show the
gas radii as a function of the amount of 26Al enrichment, where panel (d)
shows enrichment from supernovae only, panel (e) shows enrichment from
supernovae and dispersed wind ejecta and panel (f) shows enrichment from
supernovae and localised wind ejecta. Panel (g) shows the disc gas mass as a
function of the 60Fe enrichment and panel (h) shows the disc radius as a
function of the 60Fe enrichment. Vertical dotted lines in panels (a)–(f) show
the measured 26Al abundance in the Solar system Thrane et al. (2006). The
vertical dotted lines in panels (g) and (h) show the range of quoted values
for the measured 60Fe abundance in the Solar system (Tang & Dauphas, 2012;
Mishra et al., 2016). The horizontal lines in panels (a)–(c) and (g) indicate
1 MJup, which we take to be the minimum mass required to form Jupiter and
Saturn in our Solar system. In all panels, black points indicate the values at
10 Myr, whereas the blue points indicate the values much earlier, at 2.5 Myr.
Stars that are roughly Solar-mass (0.5 - 1.5 M⊙) are highlighted by a red
circle.

However, for winds entrained in localised ($<0.1$ pc) bubbles, several stars
per simulation reach Solar system-like abundances. Note that there are no
systems with low 26Al/27Al ratios after 2.5 Myr in the simulations where the
winds are entrained in localised bubbles (the blue points in Fig. 2 and 2);
this is because the yields from these localised bubbles are much higher than
when the winds are dispersed thoughout the star-forming region. The handful of
points at 10 Myr with low 26Al/27Al ratios are stars that are enriched by
supernovae only, and previously did not encounter the small-volume wind
bubbles. Furthermore, around half of these objects still retain enough gas
($>$1 MJup) to form gas giant planets (systems lying above the horizontal line
in panels (a)–(c) and on or around the vertical line). A tiny subset of these
systems are Solar-mass stars (indicated by the red circles).

Similarly, for the gas radii (Figs. 2–2), around a third of the systems that
are enriched to Solar system levels retain reasonably large ($>1$ au) disc
radii.

In these simulations, the enrichment in 60Fe comes from supernovae alone and
the gas mass-enrichment plot is shown in Fig. 2 and the gas radii-enrichment
plot is shown in Fig. 2. The very few surviving gas-rich discs do not attain
enrichment levels even at the lower-end of the measured values for our Solar
system (Tang & Dauphas, 2012, the lefthand vertical line in panels (g) and
(h)).

###  3.2 Photoevaporation where dust radius follows gas radius inwards,
$r_{\rm disc,i}=100$ au

We now run the same set of simulations (moderate stellar density –
$\tilde{\rho}\sim 100$ M⊙ pc-3), initial disc radii $r_{\rm disc,i}=100$ au
with external photoevaporation, but this time the dust radius follows the gas
radius (e.g. Sellek et al., 2020). For example, when we adjust the radius due
to photoevaporative mass-loss using Eq. 7, we also reduce the radius of the
disc when calculating the cross sections for sweeping up supernovae and wind
ejecta (Eqs. 16 and 17).

We show the disc gas properties (radius and mass) as a function of the amount
of 26Al enrichment in Fig. 3. In this figure we do not show the results for
enrichment from supernovae alone, as no discs that are in close proximity to
massive stars as they explode retain any gas content. Similarly, we do not
show plots for 60Fe enrichment because the gas content of discs enriched in
60Fe has been destroyed.

The lefthand panels (Fig. 3 and 3) show the enrichment versus disc mass (panel
a) and enrichment versus disc radius (panel c) for dispersed winds and
supernovae. As before, more discs survive at 2.5 Myr (the blue points) than at
10 Myr (the black points); however, the enrichment levels at 2.5 Myr are much
lower than at 10 Myr. However, even at 10 Myr, the enrichment levels from
winds that have dispersed around the star-forming region are always $Z_{\rm
Al}<10^{-7}$, more than two orders of magnitude lower than the measured value
in the Solar system.

The enrichment levels of discs are higher when we assume the winds are
retained in localised bubbles (Fig. 3 and 3), but fewer stars’ discs survive
due to their closer proximity to the photoionising radiation. Even then, the
enrichment levels are not close to the lower-limit or “canonical value" of the
$Z_{\rm Al}$ ratio for the Solar system ($\sim 10^{-5}$, Jacobsen et al.,
2008; Kita et al., 2013).

Figure 3: Properties of discs as a function of the amount of enrichment in a
moderately dense star-forming region ($\tilde{\rho}\sim 100$ M⊙ pc${}^{-3})$
where the discs all have initial radii $r_{\rm disc,i}=100$ au. The disc gas
radii are evolved inwards due to photoevaporation (and the discs lose gas-
mass), and the dust radii also follow the gas radii inwards. The dust masses
remain constant. Panels (a)–(b) show the gas mass as a function of the amount
of 26Al enrichment (expressed as the 26Al/27Al ratio). Panel (a) shows
enrichment from supernovae and dispersed wind ejecta and panel (b) shows
enrichment from supernovae and localised wind ejecta. Panels (c)–(d) show the
gas radii as a function of the amount of 26Al enrichment, where panel (c)
shows enrichment from supernovae and dispersed wind ejecta and panel (d) shows
enrichment from supernovae and localised wind ejecta. The vertical dotted
lines in panels (a)–(d) show the measured 26Al abundance in the Solar system
Thrane et al. (2006). The horizontal lines in panels (a) and (b) indicate 1
MJup, which we take to be the minimum mass required to form Jupiter and Saturn
in our Solar system. In all panels, black points indicate the values at 10
Myr, whereas the blue points indicate the values much earlier, at 2.5 Myr.
Stars that are roughly Solar-mass (0.5 - 1.5 M⊙) are highlighted by a red
circle.

###  3.3 Photoevaporation and viscous evolution of disc, $r_{\rm disc,i}=10$
au

Many theories of the evolution of protoplanetary discs posit that a
significant amount of expansion due to viscosity occurs (e.g. Shakura &
Sunyaev, 1973; Lynden-Bell & Pringle, 1974; Hartmann et al., 1998; Concha-
Ramírez et al., 2019), although recent observations suggest that the amount of
viscous expansion may be quite modest (e.g. Manara et al., 2022).

In similar simulations to those we present here, Parker et al. (2021) show
that viscous discs that expand are depleted much faster than static discs,
because the movement of the disc radius inwards preserves – to some extent –
the surface density of the disc, whereas a reduction in mass from
photoevaporation, followed by the gas radius moving back outwards, reduces the
surface density of the disc. The reduction in surface density therefore makes
the disc even more susceptible to photoevaporation at the next timestep.

In one set of simulations, we implemented photoevaporation, disc mass-loss and
viscous evolution with an initial disc radius of $r_{\rm disc,i}=100$ au, but
all of the discs were destroyed before any significant enrichment had
occurred. Therefore, in the set of simulations we describe here, the disc
radii are initialised to $r_{\rm disc,i}=10$ au, which was adoped by
Lichtenberg et al. (2016) in their simulations.

As with the previous simulations where the dust radius is allowed to move
inwards with the gas radius, the discs are either destroyed, or have such
small radii when the supernovae explode, that we do not get any significant
enrichment of discs from supernovae alone. Again, this means that the discs in
these simulations are not enriched in 60Fe and the enrichment in 26Al comes
from the stellar winds.

We show the results for the simulations where we assume the winds are
dispersed in the lefthand panels (a and c) in Fig. 4. A significant number of
discs experience low levels of enrichment ($Z_{\rm Al}<10^{-8}$) after 2.5
Myr, but only a handful of systems that are enriched retain any gas after 10
Myr, and even then the enrichment is $Z_{\rm Al}<10^{-7}$, two orders of
magnitude lower than in the Solar system.

If we assume that the winds from the massive stars stay in localised bubbles
(Fig. 4 and 4), then the enrichment levels increase to $Z_{\rm Al}\leq
10^{-6}$, but this is still an order of magnitude lower than the Solar system
value. Furthermore, in these models, no discs are enriched and retain gas
content at 10 Myr (note the absence of black points in Figs. 4 and 4).

Figure 4: Properties of discs as a function of the amount of enrichment in a
moderately dense star-forming region ($\tilde{\rho}\sim 100$ M⊙ pc${}^{-3})$
where the discs all have initial radii $r_{\rm disc,i}=10$ au. The disc gas
radii are evolved inwards due to photoevaporation (and the discs lose gas-
mass), and outwards due to viscous spreading. The dust radii follow the gas
radii inwards due to photoevaporation (though the dust masses remain constant)
and outwards due to viscous spreading. Panels (a)–(b) show the gas mass as a
function of the amount of 26Al enrichment (expressed as the 26Al/27Al ratio).
Panel (a) shows enrichment from supernovae and dispersed wind ejecta and panel
(b) shows enrichment from supernovae and localised wind ejecta. Panels (c)–(d)
show the gas radii as a function of the amount of 26Al enrichment, where panel
(c) shows enrichment from supernovae and dispersed wind ejecta and panel (d)
shows enrichment from supernovae and localised wind ejecta. The vertical
dotted lines in panels (a)–(d) show the measured 26Al abundance in the Solar
system Thrane et al. (2006). The horizontal lines in panels (a) and (b)
indicate 1 MJup, which we take to be the minimum mass required to form Jupiter
and Saturn in our Solar system. In all panels, black points indicate the
values at 10 Myr, whereas the blue points indicate the values much earlier, at
2.5 Myr. Stars that are roughly Solar-mass (0.5 - 1.5 M⊙) are highlighted by a
red circle.

###  3.4 Stellar density

In the results presented in the previous subsections, we adopt initial stellar
densities of $\sim$100 M⊙ pc-3. The initial stellar density of observed star-
forming regions is difficult to determine, but analysis of the current spatial
(Parker et al., 2014b; Blaylock-Squibbs et al., 2022) and kinematic
(Schoettler et al., 2020) distributions of star-forming regions can place
strong constraints on the initial densities.

Recent work (Parker & Schoettler, 2022) has shown that some star-forming
regions may be as dense as $\sim 10^{4}$ M⊙ pc-3, whereas some very low
density regions (e.g. Taurus, Cyg OB2) exhibit spatial and kinematic
substructure, which means they were unlikely to have ever been more dense than
they are observed today (Wright et al., 2014; Wright et al., 2016).

The more dense a star-forming region is, the more enrichment will occur, but
conversely, there will be more destructive photoevaporation than in lower
density regions. For several of our models, we also ran high density versions
of the same simulation and found largely similar results; some systems
experienced greater enrichment, but also lost more gas content from their
discs and the two effects largely cancel each other out. Similarly, for very
low-density regions ($\tilde{\rho}\sim 10$ M⊙ pc-3, very little
photoevaporation occurs, but there is also less (usually zero) enrichment.

###  3.5 Stellar evolution of massive stars

For the majority of our initial conditions sets we calculated the yields from
winds and supernovae of massive stars with stellar rotation (300 km s-1; the
‘recommended’ models from Limongi & Chieffi, 2018) but we also ran simulations
where the massive stars have zero rotation (also from Limongi & Chieffi,
2018).

As detailed in Parker et al. (2023), the effect of including the rotation of
the massive stars in the models (see de Mink et al., 2013) is to prolong the
lifetimes of the massive stars before they explode as supernovae, but rotation
also dredges up more of the interior of the massive stars, so more SLR
material is then entrained in the stellar winds.

In our simulations, when we include rotation in the massive stars we find that
the levels of enrichment can be higher by an order of magnitude, but almost
none of the enrichment comes from supernovae. In contrast, when the massive
stars have zero rotation, they explode as supernovae earlier, but the overall
enrichment is lower because not as much SLR material leaves the massive stars
via their winds.

##  4 Discussion

Despite exploring a wide range of parameter space, we have been unable to
enrich a protoplanetary disc with Solar system levels of 26Al and 60Fe which
also retains a significant mass ($>1$ MJup) in gas with a radius of around 10
au (which would facilitate the formation of gas giants like those in our Solar
system). We almost match these constraints for simulations in which the dust
component of the disc dynamically decouples from the gas component, and the
dust radius remains static when the gas radius moves inwards.

Current observational consensus (see Manara et al., 2022; Miotello et al.,
2022, and references therein) is that the viscosity parameter $\alpha$ is no
higher than $\sim 10^{-3}$, which would mean that the dust radius does not
increase significantly due to viscous evolution. The abundance of dust
substructures in observed protoplanetary disks (Bae et al., 2022) similarly
suggests a lack of movement of at least parts of the dust mass in a
substantial number of disks. Therefore, our simulations in which we keep the
dust radius fixed and allow the gas radius to move inwards due to
photoevaporation are likely to be the most realistic.

However, for other scenarios for the disc physics (including viscous evolution
(e.g. Concha-Ramírez et al., 2019), or when we allow the dust radius to follow
the gas radius inwards during photoevaporation, e.g. Sellek et al. (2020)),
the enrichment levels are not high enough for discs that contain enough mass
to form the gas giants in our Solar system.

In any scenario, we do not reproduce even the measured lower limit for the
60Fe content in the Solar system, because the 60Fe is delivered exclusively by
supernovae, which only explode once the majority of the discs have lost all of
their gas mass.

The natural interpretation of our results is that the Solar system’s
enrichment in SLRs cannot be from pollution of the protosolar disc from
supernovae and wind ejecta from massive stars. There are, however, several
caveats to our work which we discuss below.

First, the models of mass-loss due to photoevaporation (the FRIED grid,
Haworth et al., 2018b) could be overestimating the amount of mass lost due to
photoevaporation. Second, gas and dust leftover from star formation in the GMC
(collectively referred to as extinction) could shield the protoplanetary discs
from external photoevaporation (Qiao et al., 2023). However, simulations of
the formation of massive stars in GMCs (e.g. Dale & Bonnell, 2012; Dale et
al., 2013) show that the massive star winds blow huge cavities in their local
vicinities, which would expose protoplanetary discs to the FUV and EUV
radiation.

In addition, our disc evolution and enrichment models are a post-processing
analysis, so the discs are not reactive to the dynamics within the star-
forming regions. Whilst attempts have been made at including realistic discs
in $N$-body simulations (e.g. Rosotti et al., 2014), there is yet to be a
self-consistent simulation that can model the dynamics of the star-forming
region and the interaction between the disc and the radiation fields.

Similarly, our simulations are a purely gravitational $N$-body problem, and do
not include any physics to model the formation of stars (aside from a quasi-
realistic prescription to set up the spatial and kinematic distribution of the
stars) or the formation of the discs.

Our results diverge from the conclusions of previous work in this area (e.g.
Adams et al., 2006; Fatuzzo & Adams, 2008; Adams, 2010), who demonstrated that
our Solar system could form in a relatively populous ($N_{\star}\sim$ 1000)
star-forming region, survive photoevaporation and also obtain the measured
enrichment levels in the Solar system. However, our work makes use of newer
FUV-induced mass-loss models by Haworth et al. (2018b), which Parker et al.
(2021) demonstrate leads to significantly more disc destruction than the older
models (Johnstone et al., 1998; Störzer & Hollenbach, 1999; Scally & Clarke,
2001). Furthermore, previous research on enrichment (Adams, 2010; Parker et
al., 2014a; Lichtenberg et al., 2016) invoked supernovae from massive stars
($>$25 M$\odot$), which explode as supernovae at early ages ($<5$ Myr) before
the discs have been photoevaporated. Supernovae provide more enrichment than
stellar winds, but the most recent stellar evolution models (Limongi &
Chieffi, 2018) find that very massive stars do not explode as supernovae. The
combination of these factors are responsible for the large differences beween
our results and previous work in this area.

##  5 Conclusions

We present $N$-body simulations of the dynamical evolution of star-forming
regions to calculate the amount of enrichment in the short-lived radioisotopes
(SLRs) 26Al and 60Fe from both the stellar winds and supernovae explosions of
massive stars. Simultaneously, we calculate the effect of photoionising
radiation from massive stars on the gas masses and radii of the protoplanetary
discs to determine whether a disc can be enriched to Solar system levels, but
retain enough gas to form Jupiter and Saturn. Our conclusions are the
following:

(i) The only simulations in which we can produce Solar system levels of 26Al
enrichment are where the dust radius remains fixed, i.e. it does not move
inwards during photoevaporation (as the gas radius does) and does not move
outwards due to viscous spreading.

(ii) When we allow the dust radius to move inwards, as is thought to happen
during external photoevaporation, and/or model the viscous spreading of the
disc, no discs attain Solar system levels of enrichment and retain enough gas
to form the giant planets.

(iii) In no simulations do discs attain the Solar system levels of 60Fe, in
spite of the large uncertainty in the measured value in the Solar system. This
is because 60Fe is only produced in the supernovae explosions of massive stars
in our simulations, and photoevaporation has depleted the gas component of
protoplanetary discs by the time of the supernovae. Indeed, in the majority of
our simulations supernovae do not occur in 10 Myr due to the high stellar
rotation of the massive stars. If we switch off stellar rotation, many more
discs are enriched by the supernovae, but at levels orders of magnitude lower
than the Solar system values.

(iv) Taken together, our results suggest that enrichment of the Solar system
via the pollution of the protoplanetary disc from winds and supernovae of
massive stars in a star-forming region seems unlikely. However, we emphasise
that this result is only valid for Solar-like planetary systems that host gas
giant planets; enrichment in SLRs and the related geophysical processes (e.g.
Lichtenberg et al., 2019; Lichtenberg et al., 2022; Lichtenberg & Clement,
2022) could still occur in planetary systems that do not form gas giants.

(v) We also note that in the scenario in which the formation of the Solar
system occurs in a series of sequential or triggered star-forming events
(Gounelle & Meynet, 2012) photoevaporation from massive stars may still be
problematic for forming gas giants planets.

Given the above constraints, we encourage further research into the origin of
the SLRs in the Solar system, and suggest that their possible origin in
Asymptotic Giant Branch stars (Trigo-Rodríguez et al., 2009; Lugaro et al.,
2018; Battino et al., 2023; Parker & Schoettler, 2023) be investigated in
further detail, as this is the only scenario in which photoionising radiation
from massive stars is not present.

## Acknowledgements

We thank the anonymous referee for their comments and suggestions on the
original manuscript. RJP acknowledges support from the Royal Society in the
form of a Dorothy Hodgkin Fellowship. TL was supported by a research grant
from the Branco Weiss Foundation.

## Data availability statement

The data used to produce the plots in this paper will be shared on reasonable
request to the corresponding author.

## References

  * Adams (2010) Adams F. C., 2010, ARA&A, 48, 47
  * Adams et al. (2004) Adams F. C., Hollenbach D., Laughlin G., Gorti U., 2004, ApJ, 611, 360
  * Adams et al. (2006) Adams F. C., Proszkow E. M., Fatuzzo M., Myers P. C., 2006, ApJ, 641, 504
  * Alves et al. (2020) Alves F. O., Cleeves L. I., Girart J. M., Zhu Z., Franco G. A. P., Zurlo A., Caselli P., 2020, ApJ, 904, L6
  * Armitage (2000) Armitage P. J., 2000, A&A, 362, 968
  * Bae et al. (2022) Bae J., Isella A., Zhu Z., Martin R., Okuzumi S., Suriano S., 2022, arXiv e-prints, p. arXiv:2210.13314
  * Bate (2012) Bate M. R., 2012, MNRAS, 419, 3115
  * Battino et al. (2023) Battino U., et al., 2023, MNRAS, 520, 2436
  * Blaylock-Squibbs et al. (2022) Blaylock-Squibbs G. A., Parker R. J., Buckner A. S. M., Güdel M., 2022, MNRAS, 510, 2864
  * Boss (2013) Boss A. P., 2013, preprint,  (arXiv:1306.2220)
  * Boss (2017) Boss A. P., 2017, ApJ, 844, 113
  * Buckner et al. (2019) Buckner A. S. M., et al., 2019, A&A, 622, A184
  * Buser & Kurucz (1992) Buser R., Kurucz R. L., 1992, A&A, 264, 557
  * Cartwright & Whitworth (2004) Cartwright A., Whitworth A. P., 2004, MNRAS, 348, 589
  * Chen et al. (2013) Chen X., et al., 2013, ApJ, 768, 110
  * Concha-Ramírez et al. (2019) Concha-Ramírez F., Wilhelm M. J. C., Portegies Zwart S., Haworth T. J., 2019, MNRAS, 490, 5678
  * Concha-Ramírez et al. (2021) Concha-Ramírez F., Wilhelm M. J. C., Portegies Zwart S., van Terwisga S. E., Hacar A., 2021, MNRAS, 501, 1782
  * Cook et al. (2021) Cook D. L., Meyer B. S., Schönbächler M., 2021, ApJ, 917, 59
  * Crowther et al. (2010) Crowther P. A., Schnurr O., Hirschi R., Yusof N., Parker R. J., Goodwin S. P., Kassim H. A., 2010, MNRAS, 408, 731
  * Dale & Bonnell (2012) Dale J. E., Bonnell I. A., 2012, MNRAS, 422, 1352
  * Dale et al. (2013) Dale J. E., Ercolano B., Bonnell I. A., 2013, MNRAS, 430, 234
  * Duchêne & Kraus (2013) Duchêne G., Kraus A., 2013, ARA&A, 51, 269
  * Elmegreen (2006) Elmegreen B. G., 2006, ApJ, 648, 572
  * Fatuzzo & Adams (2008) Fatuzzo M., Adams F. C., 2008, ApJ, 675, 1361
  * Fatuzzo & Adams (2015) Fatuzzo M., Adams F. C., 2015, ApJ, 813, 55
  * Fatuzzo & Adams (2022) Fatuzzo M., Adams F. C., 2022, ApJ, 925, 56
  * Foster et al. (2015) Foster J. B., et al., 2015, ApJ, 799, 136
  * Gilmour & Middleton (2009) Gilmour J. D., Middleton C. A., 2009, Icarus, 201, 821
  * Goodwin et al. (2004) Goodwin S. P., Whitworth A. P., Ward-Thompson D., 2004, A&A, 423, 169
  * Gounelle (2015) Gounelle M., 2015, A&A, 582, A26
  * Gounelle & Meynet (2012) Gounelle M., Meynet G., 2012, A&A, 545, A4
  * Habing (1968) Habing H. J., 1968, BAIN, 19, 421
  * Hacar et al. (2013) Hacar A., Tafalla M., Kauffmann J., Kovács A., 2013, A&A, 554, A55
  * Haisch et al. (2001) Haisch Jr. K. E., Lada E. A., Lada C. J., 2001, ApJL, 553, L153
  * Hartmann (2009) Hartmann L., 2009, Accretion Processes in Star Formation: Second Edition
  * Hartmann et al. (1998) Hartmann L., Calvet N., Gullbring E., D’Alessio P., 1998, ApJ, 495, 385
  * 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
  * Jacobsen et al. (2008) Jacobsen B., Yin Q.-z., Moynier F., Amelin Y., Krot A. N., Nagashima K., Hutcheon I. D., Palme H., 2008, Earth and Planetary Science Letters, 272, 353
  * Johnstone et al. (1998) Johnstone D., Hollenbach D., Bally J., 1998, ApJ, 499, 758
  * Kirk & Myers (2011) Kirk H., Myers P. C., 2011, ApJ, 727, 64
  * Kita et al. (2013) Kita N. T., et al., 2013, Meteoritics and Planetary Science, 48, 1383
  * Lada & Lada (2003) Lada C. J., Lada E. A., 2003, ARA&A, 41, 57
  * Larson (1981) Larson R. B., 1981, MNRAS, 194, 809
  * Lehmann et al. (2010) Lehmann H., Vitrichenko E., Bychkov V., Bychkova L., Klochkova V., 2010, A&A, 514, A34
  * Lichtenberg & Clement (2022) Lichtenberg T., Clement M. S., 2022, ApJ, 938, L3
  * Lichtenberg et al. (2016) Lichtenberg T., Parker R. J., Meyer M. R., 2016, MNRAS, 462, 3979
  * 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
  * Lichtenberg et al. (2022) Lichtenberg T., Schaefer L. K., Nakajima M., Fischer R. A., 2022, arXiv e-prints, p. arXiv:2203.10023
  * Limongi & Chieffi (2018) Limongi M., Chieffi A., 2018, ApJS, 237, 13
  * Lodders (2003) Lodders K., 2003, ApJ, 591, 1220
  * Lugaro et al. (2018) Lugaro M., Ott U., Kereszturi Á., 2018, Progress in Particle and Nuclear Physics, 102, 1
  * Luhman (2004) Luhman K. L., 2004, ApJ, 602, 816
  * Luhman et al. (2003) Luhman K. L., Stauffer J. R., Muench A. A., Rieke G. H., Lada E. A., Bouvier J., Lada C. J., 2003, ApJ, 593, 1093
  * Lynden-Bell & Pringle (1974) Lynden-Bell D., Pringle J. E., 1974, MNRAS, 168, 603
  * Manara et al. (2022) Manara C. F., Ansdell M., Rosotti G. P., Hughes A. M., Armitage P. J., Lodato G., Williams J. P., 2022, arXiv e-prints, p. arXiv:2203.09930
  * Maschberger (2013) Maschberger T., 2013, MNRAS, 429, 1725
  * Miotello et al. (2022) Miotello A., Kamp I., Birnstiel T., Cleeves L. I., Kataoka A., 2022, arXiv e-prints, p. arXiv:2203.09818
  * Mishra et al. (2016) Mishra R. K., Marhas K. K., Sameer 2016, Earth and Planetary Science Letters, 436, 71
  * Najita & Kenyon (2014) Najita J. R., Kenyon S. J., 2014, MNRAS, 445, 3315
  * Nicholson & Parker (2017) Nicholson R. B., Parker R. J., 2017, MNRAS, 464, 4318
  * 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
  * Ouellette et al. (2007) Ouellette N., Desch S. J., Hester J. J., 2007, ApJ, 662, 1268
  * Ouellette et al. (2010) Ouellette N., Desch S. J., Hester J. J., 2010, ApJ, 711, 597
  * Parker & Dale (2015) Parker R. J., Dale J. E., 2015, MNRAS, 451, 3664
  * Parker & Goodwin (2007) Parker R. J., Goodwin S. P., 2007, MNRAS, 380, 1271
  * Parker & Schoettler (2022) Parker R. J., Schoettler C., 2022, MNRAS, 510, 1136
  * Parker & Schoettler (2023) Parker R. J., Schoettler C., 2023, ApJ, 952, L16
  * Parker et al. (2014a) Parker R. J., Church R. P., Davies M. B., Meyer M. R., 2014a, MNRAS, 437, 946
  * Parker et al. (2014b) Parker R. J., Wright N. J., Goodwin S. P., Meyer M. R., 2014b, MNRAS, 438, 620
  * Parker et al. (2021) Parker R. J., Nicholson R. B., Alcock H. L., 2021, MNRAS, 502, 2665
  * Parker et al. (2023) Parker R. J., Lichtenberg T., Patel M., Polius C. K. M., Ridsdill-Smith M., 2023, MNRAS, 521, 4838
  * Peretto et al. (2006) Peretto N., André P., Belloche A., 2006, A&A, 445, 979
  * Pineda et al. (2015) Pineda J. E., et al., 2015, Nature, 518, 213
  * Portegies Zwart (2019) Portegies Zwart S., 2019, A&A, 622, A69
  * Portegies Zwart et al. (1999) Portegies Zwart S. F., Makino J., McMillan S. L. W., Hut P., 1999, A&A, 348, 117
  * Portegies Zwart et al. (2001) Portegies Zwart S. F., McMillan S. L. W., Hut P., Makino J., 2001, MNRAS, 321, 199
  * Pringle (1981) Pringle J. E., 1981, ARA&A, 19, 137
  * Qiao et al. (2023) Qiao L., Coleman G. A. L., Haworth T. J., 2023, MNRAS, 522, 1939
  * Raghavan et al. (2010) Raghavan D., et al., 2010, ApJSS, 190, 1
  * Richert et al. (2018) Richert A. J. W., Getman K. V., Feigelson E. D., Kuhn M. A., Broos P. S., Povich M. S., Bate M. R., Garmire G. P., 2018, MNRAS, 477, 5191
  * Rosotti et al. (2014) Rosotti G. P., Dale J. E., de Juan Ovelar M., Hubber D. A., Kruijssen J. M. D., Ercolano B., Walch S., 2014, MNRAS, 441, 2094
  * Salpeter (1955) Salpeter E. E., 1955, ApJ, 121, 161
  * Sánchez & Alfaro (2009) Sánchez N., Alfaro E. J., 2009, ApJ, 696, 2086
  * Scally & Clarke (2001) Scally A., Clarke C., 2001, MNRAS, 325, 449
  * Scally & Clarke (2002) Scally A., Clarke C., 2002, MNRAS, 334, 156
  * Schaller et al. (1992) Schaller G., Schaerer D., Meynet G., Maeder A., 1992, A&AS, 96, 269
  * Schoettler et al. (2020) Schoettler C., de Bruijne J., Vaher E., Parker R. J., 2020, MNRAS, 495, 3104
  * Segura-Cox et al. (2020) Segura-Cox D. M., et al., 2020, Nature, 586, 228
  * Sellek et al. (2020) Sellek A. D., Booth R. A., Clarke C. J., 2020, MNRAS, 492, 1279
  * Shakura & Sunyaev (1973) Shakura N. I., Sunyaev R. A., 1973, A&A, 500, 33
  * Sternberg et al. (2003) Sternberg A., Hoffmann T. L., Pauldrach A. W. A., 2003, ApJ, 599, 1333
  * Störzer & Hollenbach (1999) Störzer H., Hollenbach D., 1999, ApJ, 515, 669
  * Tang & Dauphas (2012) Tang H., Dauphas N., 2012, Earth and Planetary Science Letters, 359, 248
  * Thrane et al. (2006) Thrane K., Bizzarro M., Baker J. A., 2006, ApJL, 646, L159
  * Trappitsch et al. (2018) Trappitsch R., et al., 2018, ApJ, 857, L15
  * Trigo-Rodríguez et al. (2009) Trigo-Rodríguez J. M., García-Hernández D. A., Lugaro M., Karakas A. I., van Raai M., García Lario P., Manchado A., 2009, Meteoritics and Planetary Science, 44, 627
  * Vacca et al. (1996) Vacca W. D., Garmany C. D., Shull J. M., 1996, ApJ, 460, 914
  * Ward-Duong et al. (2015) Ward-Duong K., et al., 2015, MNRAS, 449, 2618
  * Weidner & Kroupa (2006) Weidner C., Kroupa P., 2006, MNRAS, 365, 1333
  * Winter et al. (2018) Winter A. J., Clarke C. J., Rosotti G., Ih J., Facchini S., Haworth T. J., 2018, MNRAS, 478, 2700
  * Wright et al. (2014) Wright N. J., Parker R. J., Goodwin S. P., Drake J. J., 2014, MNRAS, 438, 639
  * Wright et al. (2016) Wright N. J., Bouy H., Drew J. E., Sarro L. M., Bertin E., Cuillandre J.-C., Barrado D., 2016, MNRAS, 460, 2593
  * Young (2014) Young E. D., 2014, Earth and Planetary Science Letters, 392, 16
  * de Mink et al. (2013) de Mink S. E., Langer N., Izzard R. G., Sana H., de Koter A., 2013, ApJ, 764, 166