© ESO 2019

1. Introduction

There are several observables that put the solar system at odds with other planetary systems (Beer et al. 2004). Apart from its planetary topology, these include the curiously small disk of only ∼45 au, the morphology at the outer edge, and the tilt of the ecliptic with respect to the equatorial plane of the Sun. The high abundance of ²⁶Al with respect to the Galactic background also seems odd. Each of these observables may be the result of the early evolution of the solar system. It is controversial to think that the solar system really differs from other planetary systems (Galilei 1632; Kapteyn 1922), and naively one would expect that its characteristics are a natural consequence of its environment. Let us start with the curious abundance of ²⁶Al/²⁷Al = 4.5–5.2 × 10⁻⁵ as observed today in calcium–aluminum inclusions (CAIs) and vitreous chondrules (MacPherson et al. 1995; Jacobsen et al. 2008). These solids formed at temperatures of ≳1500 K (Hewins & Radomsky 1990), but they have a much lower abundance in ⁶⁰Fe of only ⁶⁰Fe/⁵⁶Fe = 3.8 ± 6.9 × 10⁻⁸ (Trappitsch et al. 2018). The ratio of ⁶⁰Fe/²⁶Al excludes an origin from a nearby core-collapse supernova explosion because this would result in a very low (Ouellette et al. 2010) but comparable abundances in ²⁶Al as well as in ⁶⁰Fe (Nomoto et al. 2006). An even earlier enrichment of the pre-solar nebula by the wind of a 1.6–6 M_⊙ asymptotic giant-branch star (Mowlavi & Meynet 2000; Wasserburg et al. 2006) is hard to reconcile with the timescales of star formation and disk evolution (Isella et al. 2009), and an early pre-solar enrichment through a ≳20 M_⊙ Wolf-Rayet star (Dearborn & Blake 1988; Gaidos et al. 2009; Tatischeff et al. 2010; Gounelle & Meynet 2012; Dwarkadas et al. 2017) and its subsequent supernova would lead to an anomalously high abundance in ⁶⁰Fe. These scenarios have difficulty explaining the observed short-lived radionuclides (SLRs) and neither of these explain the outer edge of the solar system’s planetesimal disk, its tilt with respect to the equatorial plane of the Sun, or the high temperatures needed for producing vitreous droplets in chondrules. The alternatives to the latter, such as electric discharge (Horányi et al. 1995) and asteroidal collisions (Arakawa & Nakamoto 2016), are controversial (Sanders & Scott 2012).

Instead of enriching the molecular cloud before the Sun formed, maybe the parental cluster hosted a Wolf-Rayet star. This star may have enriched the proto-planetary disk of the Sun directly by accretion from its copious stellar wind before it exploded in a supernova. The hosting stellar cluster has to be sufficiently massive (≳500 M_⊙) to assure a Wolf-Rayet star to be present and sufficiently dense to have its wind enrich the proto-planetary disks through accretion. Wolf-Rayet stars require a few million years to develop a massive ²⁶Al-rich wind, and in a dense environment the majority of the proto-planetary disks will by that time already have been severely truncated by stellar encounters (Punzo et al. 2014; Portegies Zwart 2016; Vincke & Pfalzner 2016) or may have evolved to a transient disk (Ribas et al. 2015). Such truncation would be a natural consequence of a dense birth-cluster (Portegies Zwart 2009), and is reconcilable with the short half-life for stars with disks in a clustered environment (Richert et al. 2018). The expectation is that disks which are affected in the early cluster evolution eventually evolve into planetary systems comparable to that of the Sun (Ronco et al. 2017).

We take these effects, the truncation of the disk due to close stellar encounters and the accretion of ²⁶Al-enriched material from a Wolf-Rayet wind, together with the effect of nearby supernovae into account in simulations of the birth cluster of the Sun. Disks in our calculations tend to be considerably truncated even before they can be enriched through accreting material from the wind of a Wolf-Rayet star. SLRs start to decay from the moment they are released from the surface of the Wolf-Rayet star. The accretion onto a circum-stellar disk will not prevent further decay of SLRs. In the solar system, the left-over by-products of ²⁶Al are found in vitreous droplets and CAIs. These form at temperatures ≳1500 K (Davis et al. 2005). Such high temperatures are sufficient to melt the dust particles and encapsulate earlier accreted ²⁶Al into vitreous droplets. Such high temperatures could be the result of a nearby supernova that has irradiated the disk. A supernova would therefore provide a natural means to embed the SLRs in vitreous droplets. The Wolf-Rayet star that initially delivered the ²⁶Al could be responsible for preserving the SLRs when it explodes at an age of 3–9 Myr (for a ≳60 − 20 M_⊙ star of solar composition; Vuissoz et al. 2004), but it could also have been another subsequent supernova.

In order to heat a circumstellar disk to ≳1500 K, the supernova has to be in close proximity. Such a nearby supernova has considerable consequences for the further evolution of the circumstellar disk. Apart from being heated, the protoplanetary disk is also harassed by the nuclear blast-wave of the supernova shell. This may lead to the truncation of the disk through ram-pressure stripping (Portegies Zwart et al. 2018), and induces a tilt to the disk due to the hydrodynamical equivalent of the Stark effect (Wijnen et al. 2017). Both processes may be responsible for shaping the outer edge of the solar system, truncating it at about 45 au and tilting the disk with respect to the equatorial plane of the Sun by i_disk = 5.6 ± 1.2°. The supernova blast-wave is insufficiently dense to copiously enrich the surviving proto-planetary disk with SLRs produced in the exploding star (Ouellette et al. 2007, 2010), which is consistent with the low abundance in observed ⁶⁰Fe. In addition, the accreted ⁶⁰Fe decays and the information of its induced abundance can only be preserved when the ⁶⁰Fe is captured in vitreous droplets, just like the accreted ²⁶Al. This requires a second heating of the disk to ≳1500 K.

The chain of events to explain the size of the disk, its tilt, and the abundances of ²⁶Al and ⁶⁰Fe, requires a windy Wolf-Rayet star and two supernovae in short succession. This seems extraordinary but appears to be a natural consequence of being born in a cluster. Here we discuss this chain of events and how it is reconcilable with the birth environment of the solar system. We quantify these processes by means of simulations in which we take into account the effects which naturally follow from the sun being born in a star cluster. These effects include the truncation of the circumstellar disk due to close stellar encounters, the accretion of ²⁶Al-enriched material from the copious wind of a nearby Wolf-Rayet star, and the effects of nearby supernovae. It turns out that these processes lead naturally to a reasonably high formation rate of systems with characteristics similar to the solar system. We further constrain the fundamental parameters of the clusters in which the solar system was born and argue that the environment from which the solar system is most likely to have emerged was a reasonably rich star cluster of moderate density.

2. The computational approach

2.1. The clustered birth envionment of the Sun

We expand the analysis of Portegies Zwart et al. (2018) by simulating the evolution of young clusters of 50 to 10⁴ stars. The integration of the equations of motion is performed in the potential of a two-armed spiral galaxy with a bar (Martínez-Barbosa et al. 2016, 2017), and we keep track of the disk truncation and mass loss due to close stellar encounters (see Sect. A.4), the Bondi-Hoyle accretion of ²⁶Al from Wolf-Rayet winds (see Sect. A.5), and the effect of supernova explosions on the protoplanetary disks (see Sect. B). Each time a supernova heats a protoplanetary disk to a temperature of > 1500 K, we preserve its present composition by stopping the decay of previously accreted SLRs. Multiple supernovae may then lead to multiple epochs of preservation with a different relative composition. When insufficiently heated, nuclear decay continues to reduce the concentration of SLRs in the disk. After the last supernova explosion occurred, at an age of ≲50 Myr, any disk that was not preserved shows only weak traces of SLRs; only sufficiently heated disks show high abundances, considerable truncation and a finite tilt angle with respect to the initial orientation of the disk.

2.2. The numerical procedure

The simulations are performed using the Astrophysical Multipurpose Software Environment (AMUSE, see Sect. A.1; Portegies Zwart et al. 2013; Pelupessy et al. 2013; Portegies Zwart 2011). The calculations are separated into three distinct parts, in which we simulate

• the effect of the supernova irradiation on a nearby protoplanetary disk (Sect. B.1),

• the effect of the supernova blast wave (Sect. B.2),

• and the consequences of encounters, accretion from stellar winds, and supernovae on protoplanetary disks in young star clusters (this section).

The parametrized effects of the supernova irradiation and blast-wave impact are integrated together with the equations of motion for the stars in the cluster. By the time a star loses mass in a wind or explodes in a supernova the effect on the other stars is taken into account. In these calculations, the stellar mass-loss parameters and supernovae are provided by the staller evolution code, whereas the masses, positions, and relative velocities of the stars with respect to the mass-losing star are provided by the N-body code. We perform a grid of calculations in cluster density, structure, mass, and virialization to find the parameters for which it is most probable to form a planetary system with characteristics (stellar mass, disk size, mass, and inclination) similar to the early solar system. In the appendix, we briefly discuss the various ingredients in the simulations. All calculations are performed to an age of 50 Myr, after which we analyze the results of searches for solar system analogs.

2.3. Initial conditions

Each simulation starts by selecting the masses, positions, and velocities of N stars. The stars in the clusters are distributed according to a Plummer (1911) sphere density distribution or using a fractal with the box-counting dimension F = 1.6 (Goodwin & Whitworth 2004). The stars are selected randomly from a broken power-law initial mass function (Kroupa 2002) between the hydrogen-burning limit and an upper mass limit according to Kroupa & Weidner (2003); in Sect. 5 we relax this assumption.

2.3.1. Selecting cluster properties

The size of the cluster is characterized with a virial radius of r_vir = 0.1 pc, 0.3, 0.5, 0.7, 1.0, 1.5, 2.0, 2.5, 3.0 and 4 pc. The number of stars in each cluster is N = 900 (the smallest number which is still expected to host a Wolf-Rayet star; see Sect. 2.3.2), 1000, 1500, 2000, 2500, 3000, 3500, 3500, 6000 and 10 000. In Sect. 5 we discuss the consequences of using a mass function with a rigid upper-mass limit of 120 M_⊙ in the randomly sampled mass function, rather than introducing a lower value for this upper-mass limit that we base on the total cluster mass. In this new set of simulations there is a ∼10% probability that a cluster of only 100 stars contains a star of ∼20 M_⊙ (see red curve in Fig. 1). In the first set of simulations, with the cluster-mass-dependent maximum stellar mass, this is not a possibility (see blue curve).

Fig. 1. Probability density function for acquiring at least one Wolf-Rayet star (solid curves) or at least one star sufficiently massive to experience a supernova (dotted curves) from the adopted broken power-law initial mass function. The red curves give the probabilities when we adopt a single fixed upper limit of 120 M⊙ to the initial mass function. For this mass function, even a very low-mass cluster could still host, by chance, a Wolf-Rayet star or a star sufficiently massive to experience a supernova explosion. The blue curve gives the probabilities of finding a massive star for a mass function with a variable upper mass limit according to Eq. (2). For the latter mass function a cluster should contain at least ∼200 stars to host a star sufficiently massive to experience a supernova, and at least ∼900 stars before it can host a Wolf-Rayet star.

For each combination of r_vir and N we perform ten simulations up to an age of 50 Myr, which is well beyond the moment of the last supernova in the cluster, and which assures that nonpreserved SLRs have decayed to an unmeasurably small abundance. In addition to the cluster mass and size we also varied the initial virial ratio using q = 0.4, 0.5 (virial equilibrium), 0.7, and 1.3.

At birth, each star receives a disk of 400 au in radius with a mass of 10% of the stellar mass in a random orientation around the parent star. This disk size is close to the upper limit of circumstellar disks observed using the Atacama Large Millimetre Array in the nearby open star clusters in Lupus (Tazzari et al. 2017), Upper-Scorpius (Barenfeld et al. 2017), Ophiuchus (Tripathi et al. 2017), Orion (Eisner et al. 2018), and Taurus (Tripathi et al. 2017). During the integration of the equations of motion, we follow the mass and size evolution of the disks, in particular, the truncation due to encounters from other passing stars (see Sect. A.4.1), their enrichment by the winds of Wolf-Rayet stars (Sect. A.5), and the influence of nearby supernovae (Appendix B).

2.3.2. Selecting masses from the initial mass function

Masses are selected from a broken power-law mass function (Kroupa 2002) between 0.08 M_⊙ and an upper limit depending on the total cluster mass according to Weidner & Kroupa (2006). We approximate this upper mass limit with

$$\begin{array}{cc}\hfill {log}_{10}\left(\frac{m_{\max }}{M_{\odot }}\right)& =-0.76+1.06{log}_{10}\left(\frac{M}{M_{\odot }}\right)\hfill \\ \hfill & -0.09{log}_{10}{\left(\frac{M}{M_{\odot }}\right)}^2·\hfill \end{array}$$(1)

We performed simulations for those clusters that, upon generating the initial mass functions have at least one star of at least 20 M_⊙.

In Fig. 1 we present the distribution for the probability of having at least one Wolf-Rayet star or at least one star sufficiently massive to experience a supernova in a cluster. For a mass function with a general upper limit (of 120 M_⊙) independent of the number of stars in the cluster, there is a finite probability of finding a Wolf-Rayet star even in a very small cluster of only ten stars, whereas for an upper-mass limit according to Eq. (2) such a cluster needs to be composed of at least approximately 900 stars. In Sect. 5.1 we discuss this upper-mass limit in more detail.

Eventually, in the calculation of the relative birth rates we correct for this bias by multiplying the cluster probability-density function with the probability of generating such a massive star (see Sect. 4). With the adopted mass function, this probability entails

$$\begin{array}{c}\hfill P_{m>20M_{\odot }}\simeq 1.0-\frac{1.2}{1+{\left(N/1350\right)}^{2.70}}·\end{array}$$(2)

We present this distribution in Fig. 1, including the probability for finding a Wolf-Rayet star in a cluster with a fixed upper limit to the mass function of 120 M_⊙.

3. Comparison with the solar system

We study the probability that a cluster produces a star with a disk similar to the solar system. For this, we assign a solar-system similarity parameter, 𝒮_sp, to each star at the end of each simulation (at t = 50 Myr). This parameter is the normalized phase-space distance in terms of stellar mass and disk size, mass, and relative inclination. The phase space distance for each of these parameters is determined by comparing the result of the simulations with Gaussian distributions around the solar system characteristics (for stellar mass 1.0 ± 0.1 M_⊙, and the disk parameters: radius r_disk = 45 ± 10 au, mass m_disk = 0.01 ± 0.003 M_⊙, and relative inclination i_disk = 5.6 ± 1.2°).

$$\begin{array}{c}\hfill {\mathcal{S}}_{\mathrm{sp}}=S_{\mathrm{stellar}-\mathrm{mass}}\times S_{\mathrm{disk}-\mathrm{mass}}\times S_{\mathrm{disk}-\mathrm{size}}\times S_{\mathrm{inclination}}.\end{array}$$(3)

For the solar system 𝒮_sp ≡ 1. The sum of the values of 𝒮_sp for all the stars in a cluster is the expectation value for solar system analogs in a particular star cluster, N_ssa.

The abundance of SLR is not part of the definition of 𝒮_sp, but in practice, it turns out that stars with a high value of 𝒮_sp also have high concentrations in ²⁶Al and sometimes also in ⁶⁰Fe (see Fig. 8). To further mediate the discussion, we directly compare the observed abundances of ²⁶Al and ⁶⁰Fe with the results of the simulations. In the latter, it is conceptually easier to consider absolute abundances, and therefore we also calculate these for the solar system. For the solar system, the observed value of ²⁶Al/²⁷Al = 4.5–5.2 × 10⁻⁵ (MacPherson et al. 1995; Jacobsen et al. 2008). The total amount of ²⁶Al in the solar system can be calculated from ²⁶Al/²⁷Al × ²⁶Al/¹H × Z_⊙ × M₂₆/M₁ = 2.9–3.4 × 10⁻⁹ M_⊙/M_⊙ (see also Gounelle 2015). Here we adopted the canonical value of ²⁷Al/¹H = 3.5 × 10⁻⁶ (Lodders 2003) and Z_⊙ = 0.71. For ⁶⁰Fe we use ⁵⁶Fe/¹H = 7.3 × 10⁻⁷ (Lodders 2003), and the observed ⁶⁰Fe/⁵⁶Fe = 3.8  ±  6.9 × 10⁻⁸ (Trappitsch et al. 2018) which results in ⁶⁰Fe/⁵⁶Fe × ²⁶Al/¹H × Z_⊙ × M₅₆/M₁ = 1.1 ± 2.0 × 10⁻¹² M_⊙/M_⊙.

We note that the abundance in ⁶⁰Fe derived by Trappitsch et al. (2018) is several orders of magnitude lower than those of Mishra & Marhas (2016, ⁶⁰Fe/⁵⁶Fe ≃7 × 10⁻⁷). Earlier work on the origin of the history of the solar system was in part motivated by these higher abundances, invoking models in which large amounts of ⁶⁰Fe could enrich the early solar system. With the newly determined lower abundances in ⁶⁰Fe some of these models show a discrepancy with the observations by over-producing ⁶⁰Fe. The model presented here is able to explain the high abundance of ²⁶Al without overproducing ⁶⁰Fe.

4. Results

4.1. The most probable birth cluster

After having performed the grid of calculations, we analyze the results, in particular with respect to the number of solar system analogs N_ssa. In Fig. 2 we present this number for the Galaxy (per Myr per cluster) by convolving the probability density function (N_ssa from the results of the simulations) with the expected number of clusters formed in the Galaxy (in mass and size; Portegies Zwart et al. 2010) and with the observed Galactic star-formation rate of 1 M_⊙ yr⁻¹ (Robitaille & Whitney 2010). The expectation value for the number of solar system equivalents born per million years is then

$$\begin{array}{c}\hfill N_{\mathrm{sse}}=N_{\mathrm{ssa}}\times f_{\mathrm{Sch}}(N)\times f_{\ln }(r_{\mathrm{vir}})\times P_{m>20M_{\odot }}.\end{array}$$(4)

Fig. 2. Expected number of solar-system equivalents formed per million years in the Galaxy for clusters with a virialized (q = 1) Plummer density distribution. The calculations are performed in a grid of ten in mass (between 900 stars and 104 stars) and radius (between 0.1 pc and 4.0 pc). The sum of the phase-space distances for each cluster is convolved with the star-cluster birth function in mass and radius, the probability of hosting a star of at least 20 M⊙ (see Sect. A.5), and assuming a star-formation rate of 1 M⊙ yr−1 (Robitaille & Whitney 2010). For the cluster-birth mass function we use a Schechter function with α = 1, β = −2.3 and mbreak = 2 × 105 M⊙ (Schechter 1976). For the cluster size distribution we adopt a log-normal distribution with rmean = 5 pc and σ = 3 pc (van den Heuvel & Portegies Zwart 2013).

Here f_Sch(N) is the Schechter (Schechter 1976) function with α = 1, β = −2.3, and m_break = 2 × 10⁵ M_⊙ and f_ln(r_vir) is the log-normal distribution with r_mean = 5 pc and σ = 3 pc (van den Heuvel & Portegies Zwart 2013).

Massive clusters tend to produce more solar-system equivalents, but when integrated over the star-cluster mass function and expected size distribution the most probable host appears to have a radius of r_vir ≃ 0.75 ± 0.25 pc and contains about 2500 ± 300 stars (see Fig. 2). This value is consistent with the earlier derived quantification of the birth cluster of the Sun (Portegies Zwart 2009; Adams 2010; Parker & Quanz 2012, assuming 𝒮_sp >  0.01 for each parameter). Clusters with these parameters produce N_ssa = 21.0 ± 5.1 solar system equivalents. In comparison, a cluster with N = 1500 and r_vir = 1.0 pc produces only N_ssa = 7.0 ± 2.8. We performed a total of 180 runs with these optimal parameters (r_vir ≃ 0.75 pc and N = 2500).

In Fig. 3 we present the distribution of disk size and disk inclination at an age of 50 Myr. The fact that this distribution is centred around the Sun should not be a surprise, because both criteria were used to determine a high value of 𝒮_sp.

Fig. 3. Probability density function for the final disk size vs. relative inclination, for those systems with a high value of 𝒮sp. The error bar indicates the current parameters of the Sun.

In Fig. 4 we present the time evolution of the average disk size, and the cumulative distribution of supernova occurrences. In the first few (∼6) million years the truncation of disks is driven by close encounters. The introduction of the supernova truncation process is clearly visible in the mean disk-size evolution. After about 14 Myr, supernovae become less dominant again in terms of disk truncation, and both processes contribute approximately equally to the size evolution of the disks.

Fig. 4. Time evolution of the mean disk size (thick orange curve) for all 180 simulations with N = 2500 and rvir = 0.75 pc. The two dotted orange curves give the standard deviation below and above the mean distribution of disk sizes. The green curve (below the orange curve) gives the evolution of mean disk size for stars with a high value of 𝒮sp. The blue curve (axis to the right, also in blue) gives the cumulative distribution of the number of supernovae in these simulations.

Close stellar encounters are more effective in truncating disks than a nearby supernova, and they more frequently (in ∼80% of the cases) truncate disks to their final value. In Fig. 5 we present the distribution of solar-system equivalents for which the disk was more severely truncated by an encounter than by a nearby supernova.

Fig. 5. Final disk radius (at 50 Myr) vs. minimum disk radius induced by a nearby supernova for those systems in which the disk was most severely truncated by the supernovae. These account for the majority (80%) of the cases, in the other ∼20% the supernova was more effective in truncating the disk than a nearby encounter. The iso-density contours are given for equal decadal percentages. The dotted curve gives the equal radius.

The majority (∼80%) of the protoplanetary disks are severely truncated by a dynamical encounter before the supernova affects the disk, but generally close stellar encounters have the greatest impact on disk size (see also Vincke et al. 2015). We illustrate this in Fig. 5, where we plot for solar-system equivalents the eventual disk radius as a function of the truncation due to the most constraining supernova.

In Fig. 6 we also present the cumulative distribution of the number of encounters as a function of time; about 60% of the truncating encounters have already occurred by the time the first Wolf-Rayet star evolves. By the time the supernova starts to preserve the proto-planetary disks, about half those disks have already experienced their strongest encounter with another star. But when supernovae effectively peter out after about 20 Myr, encounters continue to be important until the end of the simulation.

Fig. 6. Cumulative distributions of the time of the strongest encounter (orange) and the moment of the first preserving supernova (blue) for the systems equivalent to the solar-system. We only plot the distributions for the first 9 supernovae from the 180 simulations with N = 2500 and rvir = 0.75 pc. A similar distribution of supernova distance was also presented in Parker et al. (2014).

We counted as many as 19 supernovae in the clusters that are most likely to produce a solar system analog, but some clusters only experience three supernovae. In Fig. 7 we present the cumulative distribution of the distance between a supernova and the other stars in the cluster, for the first nine supernovae. We separated the curves for the different supernovae in the simulation and normalized to the first supernova. In the legend, we present the moment in time and 1σ uncertainty for the various supernovae. The first supernova typically occurs at an age of 7.0 ± 1.0 Myr, and the ninth supernova, which occurs in about 60% of the simulations, is at an age of 30.3 ± 8.7 Myr. This figure indicates that all stars are affected to some degree by at least two supernovae and that in the majority of simulations quite a large number of supernovae occur (see also Nicholson & Parker 2017). There is no correlation, however, between the number of supernovae and the expected number of Solar system equivalents.

Fig. 7. Cumulative distributions of the distance between an exploding star and the other stars in the cluster, separated in the order in which the supernovae occur and normalized to the number of runs. We only show the distance distribution for the first 9 supernovae, even though one run experienced as many as 17.

The distribution of the abundances of ²⁶Al and ⁶⁰Fe for the most likely birth cluster, presented in Fig. 8, is consistent with the observed abundances. For other cluster parameters the value of 𝒮_sp drops and the comparison with the amounts of accreted SLRs is worse. Sub-virial (q = 0.4 and 0.7) Plummer distributions tend to produce fewer solar-system equivalents by about a factor of two, but the abundance in ⁶⁰Fe is typically about an order of magnitude higher. Super virial clusters (q = 1.3) and fractal initial density distributions (using a fractal dimension of F = 1.6) tend to underproduce solar system equivalents by about an order of magnitude.

Fig. 8. Distribution of the preserved abundances in 26Al and 60Fe of proto-planetary disks for which 𝒮sp >  0.01. The equi-composition contours at constant elevation result from the superposition of 180 simulations performed with N = 2500 stars with a virial radius of rvir = 0.75 pc that resulted in 808 solar system equivalents. The error bar indicates the abundance of the solar system, with a size equivalent to the uncertainty in the measured values. The green arrow indicates the 10 Myr change in the abundance in 26Al and 60Fe due to nuclear decay. We note that the distribution is very broad, even the peak of the distribution ranges over more than two orders of magnitude in abundance.

In clusters with an initial fractal distribution, many close encounters tend to reduce the disk size of most of the stars to well below the observed solar-system disk. These disks subsequently have too small a cross-section to effectively accrete material from the Wolf-Rayet winds or supernova blast waves. Super-virial clusters tend to expand so rapidly that by the time the Wolf-Rayet stars and supernovae become effective, the stars have already moved away too far to be strongly affected by their outflows and irradiation.

4.2. The story of the three stars

The evolution of planetary systems that closely resemble the Solar system is rather typical in our simulations. This indicates that the currently observed characteristics are not extraordinary, but actually rather common. To illustrate this we present here the story of three stars, which we call green, orange, and blue. Figure 9 illustrates their evolution in terms of disk size, relative ²⁶Al abundance of the disk, and inclination; we ran 180 realizations of the system.

Upon the birth of the stars, dynamical encounters start to become effective in truncating the protoplanetary disks after about 1 Myr. The first Wolf-Rayet star starts to inject ²⁶Al-enriched material into the cluster from about 5 Myr, enriching the protoplanetary disks in its vicinity. We follow three of the stars which we eventually classify as viable solar-system analogs. Their evolutions are presented in Fig. 9.

Fig. 9. Example of the evolution of three (out of seven solar-system equivalent systems) stars in one particular simulation for N = 2500 and rvir = 0.75 pc towards a high value of 𝒮eqiv. The coloured lines (green, orange and blue) each identify the evolution of one star in time. The vertical lines identify the moments when supernova are produced in the cluster; the thickness of these lines is proportional to the amount of mass that was injected in the supernova shell.

4.2.1. The story of green

The green star starts like any other with a disk of 400 au and without any ²⁶Al. Shortly after 5 Myr the Wolf-Rayet star in the simulation starts to lose mass, which the green star accretes quite effectively resulting in an increase in its relative abundance in ²⁶Al. The explosion of the Wolf-Rayet star does not further affect the amount of accreted ²⁶Al but truncates its disk to about 300 au. The supernova was sufficiently close to preserve the disk composition, which remains constant for its remaining lifetime. The next nearby supernova occurs at an age of about 15 Myr, which further truncates the disk, but does not further increase the amount of ²⁶Al. Later, at an age of about 28 Myr, the star has a close encounter, further reducing the disk size to about 70 au. The evolution of the inclination of the disk represents a random walk and ends at a relative inclination of about 8.2°. Interestingly, by the end of the simulation, this star has no appreciable amount of ⁶⁰Fe. The first supernova preserved the disk and deposited some small amount of ²⁶Al. The accreted SLRs, however, continued to decay until the seventh supernova at about 16 Myr truncated the disk, but it did not deposit much new ⁶⁰Fe and was too far away to heat the disk appreciably.

The green star was classified as a potential solar-system candidate on the basis of its stellar mass, disk size, mass, and inclination; its ²⁶Al abundance was sufficiently copious for it to resemble the solar system, but its lack of ⁶⁰Fe makes it a less suitable candidate.

4.2.2. The story of orange

The orange star also accretes copiously from the Wolf-Rayet wind, but when this latter explodes the star is too far from the supernova for it to be heated or to have its disk truncated. The subsequent decay reduces the ²⁶Al content in the disk until the tenth and eleventh supernovae in short succession enrich the disk with ²⁶Al and ⁶⁰Fe, truncate it to about 70 au, and preserve the composition. The disk accretes copiously from the eleventh supernova, but since the irradiation of this supernova arrived before the ²⁶Al- and ⁶⁰Fe-enriched blast wave this accreted material decays again in the following few million years. Eventually, the composition settles due to the heating that occurred in the tenth and eleventh supernova. The orange star accreted enormous amounts of ⁶⁰Fe to a level of ⁶⁰Fe/²⁶Al = 0.064, which exceeds the observed value by a factor of about 42.

4.2.3. The story of blue

The blue star experienced multiple rather strong encounters before the first Wolf-Rayet star started to blow an ²⁶Al-enriched wind. By this time the disk was already truncated to about 20 au, and only a small fraction of the enriched Wolf-Rayet wind was accreted. The first five supernovae are not close enough to heat the disk to a temperature sufficiently high to preserve the SLRs, but eventually the composition of the disk is preserved at an age of about 15 Myr. The blue star accreted some ⁶⁰Fe to a final abundance of ⁶⁰Fe/²⁶Al = 8.53 × 10⁻⁶.

5. Discussion and conclusions

We performed simulations of star clusters to study the possible formation and birth environment of the solar system. We envision the early solar system before the planets formed, but using the current solar system parameters as a template. This includes the mass of the host star, the current size of the disk and the relative inclination of the ecliptic with respect to the equatorial plane of the Sun. These aspects of the solar system are prone to external influences, and they may carry information on the sun birth environment. Additionally, we anticipate that morphological changes to the young disk perpetuate to later time, and may lead to observable peculiarities today (see also, Ronco et al. 2017).

In the calculations, we follow the evolution of these circumstellar disks, while focussing on the change in disk parameters due to stellar encounters, the winds of nearby massive stars, and the effect of supernovae. The underlying direct N-body code used to calculate the dynamical evolution of the cluster in the Galactic tidal field is used to measure the encounter rate and parameters, but also the distance and relative velocity of nearby Wolf-Rayet stars and supernovae. The effects of the latter are used with parametrizations of earlier radiative-hydrodynamical calculations, which we fitted with respect to distance and relative inclination angle between the disk orientation and the line-of-sight to the exploding star. These fit results are presented in Appendices A and B. In the following section we briefly discuss some of the advances and limitations of this model.

5.1. The importance of the rigid upper limit to the stellar mass function for a specific cluster mass

We performed two sets of simulations, one in which we adopted a rigid upper limit to the initial mass fuction as a function of cluster mass (see Sect. 2.3.2), and a second series in which we relax this assumption by allowing low-mass clusters to host, by chance, a massive star.

A consequence of the first mass function is that low-mass clusters are unable to host a massive Wolf-Rayet star (see Fig. 1 in Sect. 2.3.2) or even a star sufficiently massive to experience a supernova. The possible existence of such a mass limit was already discussed in Reddish (1978) and Vanbeveren (1982), but quantified in Weidner & Kroupa (2006) which, according to Krumholz et al. (2015), is the result of the sampling bias and an underestimate of their error bars (Krumholz 2014). The Orion star-formation region does at least appear to show evidence of a deficiency of high-mass stars in low-density environments (Hsu et al. 2013), but it remains unclear if this is a global phenomenon.

The upper limit of a mass function as a function of cluster mass may have a profound influence on the dynamical (and chemical) evolution of star clusters (Kouwenhoven et al. 2014). The consequences are also profound for our study. Earlier calculations for studying the chemical enrichment of the solar system by Parker et al. (2014) take this effect into account. When neglecting cluster-mass dependency for the upper limit to the initial stellar mass function, much lower star clusters could still host stars sufficiently massive to experience a supernova or even a Wolf-Rayet star. Lichtenberg et al. (2016) and Nicholson & Parker (2017) studied the consequences for the possible proximity of an exploding star near the proto-solar system in young (< 10 Myr old) star clusters. They find that the most likely parent cluster would be about 50 − 200 M_⊙ (or about 150 to 600 stars) with a characteristic radius of about 1 pc. According to Portegies Zwart (2009), a proto-solar system would not have the opportunity to be dynamically truncated, but if sufficiently compact, such clusters could still lead to the proper amount of enrichment due to a nearby supernova, and their close proximity could sufficiently truncate the disks. We confirmed this result by performing a second series of calculations where we relaxed the rigid upper limit to the initial mass function, the consequence being that relatively low-mass sclusters can still, by chance, host a massive star. We find that the most optimum cluster size drops to about 0.2 pc with a total of about 50–100 stars per cluster, which is somewhat smaller and less massive than what is found by Lichtenberg et al. (2016) and Nicholson & Parker (2017). With the rigid upper-mass limit, low-mass clusters contribute considerably to the formation ofsolar-system equivalents. The actual contribution of the number of solar system equivalents in the Galaxy then depends quite sensitively on the star-cluster mass function. How often is a solar-mass star born together with a few other stars of which at least one evolves into a Wolf-Rayet star?

Regretfully, we do not answer these questions here, but to help to evaluate the results and possibly scale them to another –lower-mass– environment in which massive stars are common, we present Fig. 10. There we present the results of a small statistical study in which we randomly populate an initial mass function and count how many stars are sufficiently massive to experience a supernova (red) or turn into a Wolf-Rayet star (blue). We perform this experiment with a rigid upper limit to the stellar mass function of 120 M_⊙ (dotted curves), and for a mass function with an upper-mass limit that depends on the cluster mass (see Eq. (2)). The latter are presented as shaded regions. The largest discrepancy between the curves and the shaded regions are, not surprisingly, encountered in low-mass clusters.

Fig. 10. Expected number of Wolf-Rayet stars (blue) and supernovae (red) in a cluster as function of the number of stars. The red-shaded areas indicate the expectation range (one standard deviation around the mean) of stars sufficiently massive to experience a supernova for the broken power-law initial mass function with an upper mass-limit as provided by Eq. (2). The blue-shaded area gives the same statistics for the number of Wolf-Rayet stars. The dotted curves give the mean of the number of Wolf-Rayet stars (blue) and expected number of supernovae (red) for a mass-function with a fixed upper-limit of 120 M⊙.

5.2. The expected disk lifetime and survival

The processes discussed here to enrich, truncate, and tilt circumstellar disks seem to result in an appreciable number of systems with characteristics not dissimilar from the solar system. The initiation of these processes however take a while before they become effective. A Wolf-Rayet star requires at least a few million years to start developing an ²⁶Al-rich wind and the first supernova explosion typically occurs some 8–10 Myr after the birth of the cluster. In order to make these processes suitable for affecting the young solar system, we still require the circumstellar disk to be present. This may pose a serious limitation on the proposed chain of events because disks tend to dissolve on a timescale of ∼3.5 Myr (Richert et al. 2018). These lifetimes are measured in terms of timescale on which the disks of half the stars in a particular cluster drop below the detection limit. This may indicate that the disk goes from being gaseous or dust-rich to having a more rocky composition. The latter is harder to observe with infrared telescopes and therefore tends to be counted as the absence of a disk.

Other studies argue that disks may survive for much longer. Statistically, the short lifetimes of disks may be a selection effect (Pfalzner et al. 2014), and indeed, there are several environments which are much older than a few million years and still have a rich population of circumstellar disks, such as in the η Chamaeleontis Cluster with a mean age of ≳10 Myr (Lyo & Lawson 2005; Lyo et al. 2003), or the ≳30 Myr-old Tucana-Horologium Association (Mamajek et al. 2004).

At this point, it is unclear whether or not the long timescale poses a problem for the proposed chain of events. A method to shorten the timescale for this model is the introduction of multiple star clusters with slightly different ages as proposed by Parker & Dale (2016), or a Wolf-Rayet star that runs into the Sun’s birth cluster (Tatischeff et al. 2010).

The relatively short lifetimes of protoplanetary disks as anticipated from observations may be the result of the early truncation and the effects of supernovae, as we discussed here (see also Fig. 4). We find, in our most probably parent cluster that the majority of the circumstellar disks are truncated to about half their initial size on a timescale of about 10 Myr. In less smooth initial density profiles (such as the fractal distributions) or more compact clusters, this timescale can be considerably shorter. In contrast, both clusters for which disk lifetimes are observed to be long are either low in mass, as is η Chamaeleontis, or have low density, as does the Tucana-Horologium Association.

5.3. The effect of encounters and supernova on the later disk morphology

Most disks in our simulations are severely affected by dynamical encounters as well as by nearby supernovae. In particular, those disks that later turn out to acquire parameters similar to those of the solar system are strongly influenced by both processes. In many cases, dynamical encounters truncate the disk, after which it is further abrased and heated by a nearby supernova explosion. Calculations of disk harassment have been performed already and they tend to conclude that stellar encounters have a dramatic effect on the disk-edge morphology Punzo et al. (2014), Vincke & Pfalzner (2016), and that these effects may be quite common in young star clusters (Vincke et al. 2015). Both arguments are consistent with our findings. How to quantize the long-term consequences for a potential solar system is not clear however. Some of the induced fringes on the circumstellar disk may be lost in subsequent encounters or supernovae, but other effects may persist and remain recognizable also at later times (Ronco & de Elía 2018). A subsequent study, focussing on the long-term consequences of early truncation via supernova harassment of the disk would be interesting, but this was not our focus. As discussed in Sect. 4.2 in relation to Fig. 9, disks with similarities to that of the Sun tend to be truncated by dynamical encounters at an early epoch, and later again by the blast wave of a nearby supernova explosion. The former mediates the small size of the currently observed disk, but the latter is important to explain the observed tilt and the accretion and preservation of SLRs, effectively freezing the disk composition in SLRs. Recently, Cai et al. (2018) argued that the misalignment angle between a planetary system and the equatorial plane of the host star can be explained by dynamical encounters between the planetary systems in a relatively low-density star cluster. Such interaction also naturally leads to some of the processes described here, such as the truncation of the disk.

5.4. The parental star cluster

A planetary system typically experiences more than one strong encounter and more than one nearby supernova. In particular, the latter is universal to all stars in the same parental cluster. We therefore expect the composition of the SLRs of many stars to be similar to that of the Solar system, and for many stars to have truncated and inclined disks.

The star clusters which turn out to be most favorable for producing solar-system equivalents are sufficiently massive that there will be multiple stars which are sufficiently massive to experience a supernova. The typical number of supernovae is 9–15, but some clusters experience as few as 3 and some as many as 19. The first supernova tends to explode at an age of about 8 Myr. This is somewhat constructed because our clusters are initialized with the guarantee of a >  20 M_⊙ star, which for solar composition has a lifetime of ≲9 Myr, and all stars are born at the same time. In Fig. 6 we present the cumulative distributions for the time of the first conserving supernova. The time between the first and the second supernova ranges from less than 1 Myr to as long as 30 Myr. A short timescale between supernovae is essential in order to preserve ⁶⁰Fe in enriched disks (see Fig. 9). With a half-life of about 2.62 Myr for ⁶⁰Fe and even shorter for ²⁶Al, a large amount of the SLRs will have decayed if the second supernova takes too long. The clusters that produce a substantial number of solar-system equivalents experience two or more supernovae in short succession.

Interestingly, stars born in a Plummer distribution with a high value of 𝒮_sp tend to produce a wide range of abundances in ²⁶Al and ⁶⁰Fe with a reasonable match to the observed abundance in the solar system (see Fig. 8). Fractal clusters fail to reproduce the ²⁶Al and ⁶⁰Fe abundances by several orders of magnitude. The failure of fractal clusters to produce solar-system equivalents is attributed to the high degree of dynamical activity in these clusters. As a consequence, close stellar encounters in the first few million years tend to obliterate proto-planetary disks, which subsequently have too small a cross section to accrete substantial amounts of ²⁶Al from Wolf-Rayet winds and ⁶⁰Fe from the supernova blast-wave. This may change when we take the viscous evolution of the disks into account (Concha et al., in prep.), in which case truncated disks may grow back in time.

The clusters that are most likely to form a solar system analog contain N = 2500 ± 300 stars (∼900 M_⊙) with a virial radius r_vir = 0.75 ± 0.25 pc, but comparable values are obtained for somewhat-lower-mass clusters (∼550 M_⊙); the cluster radius can be as large as r_vir ≃ 1.5 pc. Increasing the virial ratio causes a dramatic drop in the expectation value of 𝒮_sp, and cool initial conditions (with a virial ratio of 0.4–0.7) tend to result in a somewhat smaller value of 𝒮_sp, but still with a distinct peak for clusters of N ∼ 1000 stars and r_vir ≃ 0.5 pc. Clusters for which the stars are initially distributed in a fractal (with dimension F = 1.6) produce fewer solar-system equivalents by about a factor of two compared to Plummer distributions. In these cases, the cluster that produces the most solar system equivalents is somewhat less massive (∼1000 stars) but is considerably larger, with a virial radius of 2–3 pc.

6. Summary

Integrating over the phase space in terms of cluster mass and size results in a Galactic birth rate of ∼30 solar-system equivalents per million years. With an expected lifespan of ∼12 Gyr we expect the Galaxy to contain some 36 000 systems with host mass, disk size, inclination angle, and abundances in ²⁶Al and ⁶⁰Fe similar to the solar system.

We argue that the Sun was born in a cluster of N = 2500 ± 300 stars (∼900 M_⊙) distributed in a smooth potential near virial equilibrium with a characteristic radius r_vir = 0.75 ± 0.25 pc. These parameters are consistent with the parameters derived for the Orion Nebula cluster by Stutz (2018). According to the calculations carried out here, such a cluster produces about 25 planetary systems with characteristics similar to the solar system. These systems have an abundance in ²⁶Al and ⁶⁰Fe that is within about two orders of magnitude of that of the current Sun’s planetesimal system, and a disk that is truncated to about 45 au and inclined with respect to the equatorial plane of the star.

The evolution of these solar-system analogs appears to be rather typical. The first few million years are characterized by repeated close encounters with other stars. This process causes the disk to be truncated to a radius of ≲250 au on average, or to ≲100 au for systems that are more similar to the solar system. A nearby Wolf-Rayet star will enrich the surviving disk with ²⁶Al isotopes, which are captured in the vitreous droplets of the disk upon the first nearby supernova that is able to heat the disk to ≳1500 K. The close proximity of this supernova again harasses the disk and also injects short-lived ⁶⁰Fe isotopes into it. This second generation of enrichment is again captured in vitreous droplets upon a subsequent nearby supernova. The entire process then takes at least two supernovae in order to explain the abundances in ²⁶Al as well as in ⁶⁰Fe. In principle, a single supernova can be sufficient, but in such cases, the ⁶⁰Fe abundance has to be primordial.

Acknowledgments

It is a pleasure to thank referee Richard Parker for useful comments, Arjen van Elteren, Maxwell Cai and Francisca Ramirez-Concha, and MSc students Fedde Fagginger and Tom Sweegers for discussions. In addition, I would like to thank Norm Murray and CITA for their hospitality during my long-term visit.

References

Appendix A: Boundary conditions and model parameters

A.1. The Astrophysical MUltipurpose Software Environment

The coupling of the various numerical methods is realized using the Astrophysical MUltipurpose Software Environment (AMUSE) Portegies Zwart et al. (2013), Portegies Zwart & McMillan (2018). AMUSE provides a homogeneous interface to a wide variety of packages which enables the study of astrophysical phenomena where complex interactions occur between different physical domains, such as stellar evolution, gravitational dynamics, hydrodynamics, and radiative transport. For this project we focus on the coupling between gravitational-dynamics solvers, stellar-, hydrodynamical, and radiative transfer solvers. The major advantage of AMUSE over other methods is the flexibility in which a wide variety of solvers can be combined to address the intricate interactions in the nonlinear multi-scale systems over all relevant scales. Also relevant is the way in which we can expand methods by incorporating additional effects. In this paper, these effects include the ablation of the circum-stellar disk by a nearby supernova, the accretion of the winds of nearby stars and the effects of close stellar encounters.

A.2. Background Galactic tidal field

Each cluster is initialized in a circular orbit at a distance of 8.5 kpc from the Galactic center. Since we mainly study relatively young star clusters we ignore the global evolution of the Galaxy but assume it to be a slowly varying potential with contributions from the bar, bulge, spiral arms, disk, and halo. We adopt the same Galactic parameters as those used in Martínez-Barbosa et al. (2016, 2017) to study the spatial distribution of solar siblings in the Galaxy.

A.3. The effects of stellar evolution, mass loss, and supernovae

Stellar evolution is taken into account using the SeBa stellar and binary-evolution code (Portegies Zwart & Verbunt 1996, 2012; Portegies Zwart & Yungelson 1998; Toonen et al. 2012).

Here we use the event-driven time-stepping scheme between the stellar evolution and the gravitational dynamics codes to assure that the stellar positions are consistent with the moment a star explodes in a supernova (see Portegies Zwart & McMillan 2018, for details).

A.4. Integrating the equations of motion of the stars in the cluster

The equations of motion of the cluster stars are integrated using the ph4 fourth-order Hermite predictor-corrector direct N-body code (for implementation details see Portegies Zwart & McMillan 2018), with a time-step parameter η = 0.01 and a softening of 100 au. The Galaxy model (see Sect. A.2), is coupled to the cluster dynamics via the sixth-order rotating-Bridge method (Portegies Zwart & McMillan 2018). This is a special high-order derived version of the classic bridge (Fujii et al. 2007) but for a rotating reference frame (Martínez-Barbosa et al. 2017). After every 0.1 Myr (equivalent to the bridge-timestep) we synchronize the gravity solver, check for energy conservation, and print diagnostics. The energy of the N-body integrator is preserved better than 1/10⁸, which is sufficient to give a reliable result (Portegies Zwart & Boekholt 2014).

While integrating the equations of motion, we check for close encounters. When two stars approach each other in a pre-determined encounter radius (initially 0.02 pc) the N-body integrator is interrupted and the system is synchronized. Such an interruption is called a stopping condition in AMUSE terminology and should be perceived as a numerical trick to stop one code in order to allow another code to resolve a particular part of the physics. The encounter is subsequently resolved and the effect on the disks of the two stars calculated (see Sect. A.4.1).

Fig. A.1. Schematic view of the the bar and spiral arms of the Galaxy at the present time. The blue bullet marks the current position of the Sun measured in an inertial system that is fixed at the center of the Galaxy. The axes Xrot and Yrot correspond to a system that corotates with the bar. We note that the spiral arms start at the edges of the bar, and the coordinates (Xrot1, Yrot1) and (Xrot, Yrot) overlap at present.

We also keep track of the mass lost by massive ≳20 M_⊙ stars during integration. This mass can be accreted by the circumstellar disks of nearby stars. As a consequence, these disks can be enriched with the SLR-rich material in the wind. This can lead to the accretion of ²⁶Al from nearby Wolf-Rayet stars (see Sect. A.5).

When a star explodes in a supernova, we generate an interrupt and subsequently calculate the effects of the irradiation and impacting SLR-loaded blastwave on all the other stars in the cluster. Short-lived radionuclides accreted by any of the disks (either by a supernova or a Wolf-Rayet star) decay with time. The half-life time for ²⁶Al is ∼1.1 Myr, and for ⁶⁰Fe this is ∼2.62 Myr (Rugel et al. 2009). As soon as a supernova heats the disk to a temperature of at least 1500 K, we freeze the abundance of SLRs. At such high temperatures, the intense radiation of the supernova melts the disk and encapsulates SLRs in vitreous droplets (see Portegies Zwart et al. 2018). Subsequent accretion may further increase the amount of SLRs, which will also decay until another supernova heats the disk again to a sufficiently high temperature. This may lead to multiple generations of enrichment. We apply the same procedure to ⁶⁰Fe which is accreted from a supernova shell.

A.4.1. The effect of encounters on disk size

The effect of a two-body encounter on the disks of both stars is solved semi-analytically. Once a two-body encounter is detected we calculate the pericenter distance, p, by solving Kepler’s equation (using the Kepler-module from the Starlab package, Portegies Zwart et al. 2001). We note that the closest approach may be well within the adopted softening radius of 100 au. The new disk-radius for a star with mass m is calculated as

$$\begin{array}{c}\hfill r_{\mathrm{disk}}^{\prime }=0.28p{\left(\frac{m}{M}\right)}^{0.32},\end{array}$$(A.1)

which was calibrated for parabolic co-planar prograde encounters (Breslau et al. 2014; Jílková et al. 2016). Here M is the mass of the other star. This equation is also applied for calculating the new disk radius of the encountering star. These new radii are adopted only if they are smaller than the pre-encounter disk radii.

In order to reduce the number of disk truncations at runtime, and therewith the number of interrupts (and synchronizations) in the N-body integration, the new encounter distance for both stars is reset to half the pericenter distance p. This prevents two stars from being detected at every integration time step while approaching pericenter, which would cause the disk to be affected repeatedly during a single encounter. This procedure therefore limits the number of encounters to the most destructive one at pericenter.

A.4.2. The effect of encounters on disk mass

The truncated disks of the encountering stars lose mass. We estimate the amount of mass lost from each disk using

$$\begin{array}{c}\hfill \mathrm{d}m=m_{\mathrm{disk}}\frac{r_{\mathrm{disk}}^{1/2}-r_{\mathrm{disk}}^{{}^{\prime }1/2}}{r_{\mathrm{disk}}^{1/2}}·\end{array}$$(A.2)

Both encountering stars may accrete some of the material lost from the disk of the other star, which we calculate with

$$\begin{array}{c}\hfill \mathrm{d}{m}_{\mathrm{acc}}=\mathrm{d}mf\frac{m}{M+m}·\end{array}$$(A.3)

Here, f ≤ 1 is a mass transfer efficiency factor (we adopted f = 1). Both equations are applied symmetrically in the two-body encounter, and as a consequence, both stars lose some mass and gain some of what the other has lost. Disks may become enriched by accreting material from earlier enriched disks in close encounters. Here we assume that accreted material has the same composition as the mean host disk.

A.5. Alumunum-26 enrichment due to Wolf-Rayet winds

The copious mass-loss of a Wolf-Rayet star is enriched in ²⁶Al and other short-lived radionuclides, but poor in ⁶⁰Fe. During the integration of the equations of motion of the stars in the cluster, we use the stellar evolution code to determine the amount of ²⁶Al liberated by the Wolf-Rayet stars. The amount of mass in the Wolf-Rayet star is calculated from rotating stellar evolution models (Vuissoz et al. 2004). We fitted the yields for solar composition in their Fig. 2 to

$$\begin{array}{c}\hfill {log}_{10}(m_{\mathrm{Al}}/M_{\odot })=1.70\times {10}^{-8}{(m/M_{\odot })}^{2.29}.\end{array}$$(A.4)

We do not take the time-dependency of the ²⁶Al-yields into account, but adopt a constant mass fraction in the wind assuming that the entire stellar envelope is homogeneously enriched with ²⁶Al at the amount given by Eq. (A.4). The stellar-evolution code provides the appropriate mass-loss rate and wind velocity.

During the integration of the equations of motion, mass lost by any of the Wolf-Rayet stars is accreted by the proto-planetary disks of the other stars. The amount of mass accreted from the wind is calculated using the Bondi-Hoyle-Littleton accretion (Bondi & Hoyle 1944) formalism.

Appendix B: The effect of a supernova explosion on protoplanetary disks

A supernova may have a profound effect on proto-planetary disks in its vicinity. Extensive simulations of this effect are performed by Ouellette et al. (2007, 2010). Similar calculations performed by Portegies Zwart et al. (2018) included the effect of a radiative transfer, by bridging a radiative transfer solver with a hydrodynamics code to measure the effects of the irradiative heating and the subsequent blast wave.

The hydrodynamics is addressed with the smoothed particles hydrodynamics code Fi (Hernquist & Katz 1989; Gerritsen & Icke 1997; Pelupessy et al. 2004), using 10⁵ SPH particles in the disk. The radiative transfer calculations are performed with SPHRay (Ritzerveld & Icke 2006; Altay et al. 2008, 2011), using 10⁶ rays per time step. A self-consistent solution is obtained by coupling the two codes using the interlaced bridge technique (Portegies Zwart & McMillan 2018).

For the disks, we adopt a power-law density profile with exponent −1 (with temperature profile ∝r^−0.5) with Savronov-Toomre Q-parameter q_out = 25 (Safronov 1960; Toomre 1964), for which the disk is everywhere at least marginally stable. The temperature of this disk ranges from 19 K (at the rim) to about 165 K in the central regions. We perform a series of such calculations in a grid in distance between the supernova and the proto-planetary disk of d = 0.05 pc, 0.1, 0.15, 0.2, 0.4 and 0.6 pc and with an incident angle with respect to the disk’s normal of Θ = 15°, 30, 45, 60 and 75°.

B.1. The effect of irradiation by the supernova

We adopted the supernova PS1-aof11 with an energy spectrum of the photons representative for a power law with an index of −3. The peak of luminosity of 1.1  ×  10⁴³ erg s⁻¹ (almost 10^9.8 L_⊙) is reached about 26 days after the supernova explosion, producing a mass in the ejecta of 23.5 M_⊙ (Sanders et al. 2015).

Protoplanetary disks in the vicinity of a supernova will be heated by the radiation. We fitted the result of the calculations discussed above (see also Portegies Zwart et al. 2018). A satisfactory fit, at a 10% accuracy over the entire range in incident angle Θ and distance to the supernova d was obtained:

$$\begin{array}{c}\hfill {log}_{10}(T/K)\simeq 16.47-13.6d^{0.03}cos{(\mathrm{\Theta })}^{-0.02}.\end{array}$$(B.1)

The fit was carried out using the Levenberg-Marquardt algorithm (Levenberg 1944; Marquardt 1963) with the data obtained for the calculated grid as input. In Fig. B.1 we present the mean disk temperature from the simulations overplotted with the fit in Eq. (B.1)

B.2. The effect of the blast-wave impact

After the disk has cooled again, the nuclear blast wave from the supernova hits and ablates the disk, tilts it, and enriches the surviving disk with a small amount of short-lived radionuclides. We calculate the effect of the supernova blast wave on the circumstellar disk by means of hydrodynamical simulations, which we perform using the smoothed particles hydrodynamics code Fi. Both the blast wave as well as the disk are simulated using the same code. We adopted a resolution of 10⁵ SPH particles in the disk as well as for the supernova blast wave.

All hydrodynamical calculations ware performed using supernova PS1-11aof. For the supernova light curve, we adopt the multiple power-law fits to the observed supernova by Sanders et al. (2015). We use fit parameters for supernovae PS1-11aof from their Table 3. The peaks at a luminosity of 1.1 × 10⁴³ erg s⁻¹ (almost 10^9.8 L_⊙) about 26 days after the supernova explosion is initiated. From the fits to the light curves, Sanders et al. (2015) derive masses in the ejecta of 23.5 M_⊙.

This supernova has a rather massive ejecta and we scale its effect on the protoplanetary disk with the amount of mass in the ejecta of the exploding star in the simulations as obtained from the stellar evolution code. We therefore multiply the various derived parameters with the relative ratio of the amount of mass in the ejecta of PS1-a11aof with respect to the amount of mass in the supernova in our simulations m_shell. Instead of simulating the entire impacting shell, we use only the intersecting part of the shell, with a diameter, at the location of the disk, of twice the disk diameter.

Fig. B.1. Temperature of the disk due to the irradiation of supernova PS1-11aof as calculated in Portegies Zwart et al. (2018) at a distance d and with incident angle Θ. The shaded colors in the background give the result of Eq. (B.1), The contours give the results from the grid of radiative-hydro calculations presented in Portegies Zwart et al. (2018).

The effect of the truncation of the disk is calculated in the same grid of parameters as the mean temperature calculations (see Sect. B.1 and also Portegies Zwart et al. 2018). We measured disk sizes 2 years after the first contact between the blast wave and the disk, which gives the blast wave sufficient time to pass the entire disk. The post-supernova disk-size was measured at a distance from the star where the surface density of the disk drops below 2 g cm². We subsequently fit the disk-size as a function of the distance d and incident angle Θ with respect to the supernova, which results in

$$\begin{array}{c}\hfill r_{\mathrm{trunc}}/\mathrm{au}\simeq 66d^{0.63}cos{(\mathrm{\Theta })}^{-0.68}.\end{array}$$(B.2)

The fit, presented in Fig. B.2, is considerably less satisfactorily than the measured disk temperature in Fig. B.1. In particular with a distance to the supernova of between 0.25 and 0.45 pc at an incident angle of Θ ≲ 45° the fit seems to break down. This is a result of the dramatic ablation when the supernova shell hits the disk almost face-on from a short distance. These nearby low-incident-angle impacts tend to lead to a sudden drop in disk density at the measured distance, with a rather low surface-density disk extending to considerable distance from the star. The eventual disk size is calculated with

$$\begin{array}{c}\hfill r_{\mathrm{disk}}=r_{\mathrm{trunc}}\left(\frac{23.5M_{\odot }}{m_{\mathrm{shell}}}\right)·\end{array}$$(B.3)

Figure B.2 presents the disk-size resulting from the calculations using this supernova as a function of incident angle and distance in Fig. B.2.

Fig. B.2. Disk truncation due to the supernova impact (type PS1-11aof) as calculated in Portegies Zwart et al. (2018) in a grid in d between 0.1 pc and 0.6 pc and an incident angle of Θ = 15° to Θ = 75° (contours). The shaded colors in the background give the result of the fit Eq. (B.3).

Fig. B.3. Inclination induced upon the proto-planetary disk due to the supernova blast wave. The shaded colors in the background give the result of fit Eq. (B.4).

As discussed in Wijnen et al. (2017) a disk hit by an external wind will readjust itself to be perpendicular to the wind. This effect is relevant for changing the inclination angle of a disk when interacting with the supernova blast wave. The inclination induced by the supernova blast wave onto the disk is calculated by the median inclination of the truncated disk for each of the simulations in the grid of models (Portegies Zwart et al. 2018). The fit results in

$$\begin{array}{c}\hfill \delta i\simeq 3.8^{°}{d}^{-0.50}cos(2\mathrm{\Theta }-0.5\pi ).\end{array}$$(B.4)

In Fig. B.3 we present the resulting inclination as a function of distance d and incident angle Θ for supernova PS1-11aof.

B.3. The accretion of ⁶⁰Fe and ²⁶Al from the supernova blast wave

The supernova produces ²⁶Al and ⁶⁰Fe, liberated in the supernova blast wave and some may be accreted by the surviving disk. The yield of ⁶⁰Fe in the supernova as a function of the progenitor mass is taken from Fig. 4 of Timmes et al. (1995) and can be fitted by

$$\begin{array}{c}\hfill {log}_{10}(m_{\mathrm{Fe}}/M_{\odot })=1.74{log}_{10}(m/M_{\odot })-6.93.\end{array}$$(B.5)

For a 20 M_⊙ star this results in a yield of ∼2.2  ×  10⁻⁵ M_⊙ in ⁶⁰Fe. For ²⁶Al the yields are

$$\begin{array}{c}\hfill {log}_{10}(m_{\mathrm{Al}}/M_{\odot })=2.43log10(m/M_{\odot })-7.23.\end{array}$$(B.6)

The fraction of mass that is accreted by the disk from a 23.5 M_⊙ supernova shell was calculated by Portegies Zwart et al. (2018), and we fitted their results to

$$\begin{array}{c}\hfill {log}_{10}(f_{\mathrm{acc}})=-4.77d^{0.49}cos{(\mathrm{\Theta })}^{-0.29}.\end{array}$$(B.7)

In Fig. B.4 we present the fraction of mass that is accreted by a protoplanetary disk of 100 au at distance d and incident angle Θ from supernova PS1-11aof. The actual mass of the supernova shell that is accreted by a disk in the simulations is calculated from the supernova shell that intersects with the disk. Here we take the disk (r_disk) into account in order to calculate the actual amount of accreted ⁶⁰Fe.

Fig. B.4. The amount of mass that is accreted by a disk of 100 au, as calculated in Portegies Zwart et al. (2018) in a grid in d between 0.1 pc and 0.6 pc and an incident angle Θ = 15° to θ = 75° (contours). The shaded colors in the background give the result of fit Eq. (B.7).

The fit in Fig. B.4, agrees within a factor of a few with the simulations. It is hard to improve these without performing a new extensive grid of simulations that cover a wider range of parameters, has a higher resolution and takes a wider range of disk morphologies and supernovae into account.

All Tables

All Figures

Fig. 1. Probability density function for acquiring at least one Wolf-Rayet star (solid curves) or at least one star sufficiently massive to experience a supernova (dotted curves) from the adopted broken power-law initial mass function. The red curves give the probabilities when we adopt a single fixed upper limit of 120 M⊙ to the initial mass function. For this mass function, even a very low-mass cluster could still host, by chance, a Wolf-Rayet star or a star sufficiently massive to experience a supernova explosion. The blue curve gives the probabilities of finding a massive star for a mass function with a variable upper mass limit according to Eq. (2). For the latter mass function a cluster should contain at least ∼200 stars to host a star sufficiently massive to experience a supernova, and at least ∼900 stars before it can host a Wolf-Rayet star.

In the text

Fig. 2. Expected number of solar-system equivalents formed per million years in the Galaxy for clusters with a virialized (q = 1) Plummer density distribution. The calculations are performed in a grid of ten in mass (between 900 stars and 104 stars) and radius (between 0.1 pc and 4.0 pc). The sum of the phase-space distances for each cluster is convolved with the star-cluster birth function in mass and radius, the probability of hosting a star of at least 20 M⊙ (see Sect. A.5), and assuming a star-formation rate of 1 M⊙ yr−1 (Robitaille & Whitney 2010). For the cluster-birth mass function we use a Schechter function with α = 1, β = −2.3 and mbreak = 2 × 105 M⊙ (Schechter 1976). For the cluster size distribution we adopt a log-normal distribution with rmean = 5 pc and σ = 3 pc (van den Heuvel & Portegies Zwart 2013).

In the text

	Fig. 3. Probability density function for the final disk size vs. relative inclination, for those systems with a high value of 𝒮sp. The error bar indicates the current parameters of the Sun.
In the text

Fig. 4. Time evolution of the mean disk size (thick orange curve) for all 180 simulations with N = 2500 and rvir = 0.75 pc. The two dotted orange curves give the standard deviation below and above the mean distribution of disk sizes. The green curve (below the orange curve) gives the evolution of mean disk size for stars with a high value of 𝒮sp. The blue curve (axis to the right, also in blue) gives the cumulative distribution of the number of supernovae in these simulations.

In the text

Fig. 5. Final disk radius (at 50 Myr) vs. minimum disk radius induced by a nearby supernova for those systems in which the disk was most severely truncated by the supernovae. These account for the majority (80%) of the cases, in the other ∼20% the supernova was more effective in truncating the disk than a nearby encounter. The iso-density contours are given for equal decadal percentages. The dotted curve gives the equal radius.

In the text

	Fig. 6. Cumulative distributions of the time of the strongest encounter (orange) and the moment of the first preserving supernova (blue) for the systems equivalent to the solar-system. We only plot the distributions for the first 9 supernovae from the 180 simulations with N = 2500 and rvir = 0.75 pc. A similar distribution of supernova distance was also presented in Parker et al. (2014).
In the text

	Fig. 7. Cumulative distributions of the distance between an exploding star and the other stars in the cluster, separated in the order in which the supernovae occur and normalized to the number of runs. We only show the distance distribution for the first 9 supernovae, even though one run experienced as many as 17.
In the text

Fig. 8. Distribution of the preserved abundances in 26Al and 60Fe of proto-planetary disks for which 𝒮sp >  0.01. The equi-composition contours at constant elevation result from the superposition of 180 simulations performed with N = 2500 stars with a virial radius of rvir = 0.75 pc that resulted in 808 solar system equivalents. The error bar indicates the abundance of the solar system, with a size equivalent to the uncertainty in the measured values. The green arrow indicates the 10 Myr change in the abundance in 26Al and 60Fe due to nuclear decay. We note that the distribution is very broad, even the peak of the distribution ranges over more than two orders of magnitude in abundance.

In the text

Fig. 9. Example of the evolution of three (out of seven solar-system equivalent systems) stars in one particular simulation for N = 2500 and rvir = 0.75 pc towards a high value of 𝒮eqiv. The coloured lines (green, orange and blue) each identify the evolution of one star in time. The vertical lines identify the moments when supernova are produced in the cluster; the thickness of these lines is proportional to the amount of mass that was injected in the supernova shell.

In the text

Fig. 10. Expected number of Wolf-Rayet stars (blue) and supernovae (red) in a cluster as function of the number of stars. The red-shaded areas indicate the expectation range (one standard deviation around the mean) of stars sufficiently massive to experience a supernova for the broken power-law initial mass function with an upper mass-limit as provided by Eq. (2). The blue-shaded area gives the same statistics for the number of Wolf-Rayet stars. The dotted curves give the mean of the number of Wolf-Rayet stars (blue) and expected number of supernovae (red) for a mass-function with a fixed upper-limit of 120 M⊙.

In the text

Fig. A.1. Schematic view of the the bar and spiral arms of the Galaxy at the present time. The blue bullet marks the current position of the Sun measured in an inertial system that is fixed at the center of the Galaxy. The axes Xrot and Yrot correspond to a system that corotates with the bar. We note that the spiral arms start at the edges of the bar, and the coordinates (Xrot1, Yrot1) and (Xrot, Yrot) overlap at present.

In the text

	Fig. B.1. Temperature of the disk due to the irradiation of supernova PS1-11aof as calculated in Portegies Zwart et al. (2018) at a distance d and with incident angle Θ. The shaded colors in the background give the result of Eq. (B.1), The contours give the results from the grid of radiative-hydro calculations presented in Portegies Zwart et al. (2018).
In the text

	Fig. B.2. Disk truncation due to the supernova impact (type PS1-11aof) as calculated in Portegies Zwart et al. (2018) in a grid in d between 0.1 pc and 0.6 pc and an incident angle of Θ = 15° to Θ = 75° (contours). The shaded colors in the background give the result of the fit Eq. (B.3).
In the text

	Fig. B.3. Inclination induced upon the proto-planetary disk due to the supernova blast wave. The shaded colors in the background give the result of fit Eq. (B.4).
In the text

	Fig. B.4. The amount of mass that is accreted by a disk of 100 au, as calculated in Portegies Zwart et al. (2018) in a grid in d between 0.1 pc and 0.6 pc and an incident angle Θ = 15° to θ = 75° (contours). The shaded colors in the background give the result of fit Eq. (B.7).
In the text