Theoretical Distributions of Short-lived Radionuclides for Star Formation in Molecular Clouds

Giant molecular clouds (GMCs) are complex and varied structures (e.g., Dobbs & Pringle 2013) with evolutionary histories that remain under investigation. Given the current uncertainties concerning the physical characteristics of these environments, we first develop a baseline model that uses standard assumptions about the GMC stellar population and the SLR transport mechanisms. We then consider several variations to our baseline model that explore other viable scenarios in order to gauge the sensitivity of our results to the physical input parameters.

In all cases, we embed a GMC within a well-defined cuboid structure with a 50 pc body diagonal (e.g., Heyer & Dame 2015). The cloud is assumed to contain 10⁶ M_⊙ of gas out of which ${N}_{\mathrm{mc}}^{* }=50,000$ stars form. The three physical dimensions of the cuboid are generated by randomly drawing numbers w_i between e^0.5 and e² (where i = {x, y, z}), and assigning ${{\ell }}_{i}=\mathrm{ln}{w}_{i}$ as relative lengths that are then normalized in accordance to the body diagonal. The resulting GMC volumes for this scheme, which range from 6550 pc³ to 24,100 pc³, have an expectation value 〈V_GMC〉 ≈ 20,000 pc³. We note that a sphere with the same volume would have a radius of r ≈ 17 pc.

In our baseline model, the formation of stars occurs within clusters whose stellar memberships N span from 10² to 10⁴ in accordance to the distribution function (e.g., Lada & Lada 2003):

$$\begin{eqnarray}&&{P}_{c}(N)\propto {N}^{-2}.\end{eqnarray} \tag{ 1 }$$

Clusters are drawn from this distribution function until their collective stellar content reaches N*_mc, with the last drawn cluster's membership reduced as needed to achieve that outcome. Each cluster is assigned a radius through the empirical formula (Adams 2010)

$$\begin{eqnarray}&&{R}_{c}=1\,\mathrm{pc}\,{\left(\displaystyle \frac{N}{300}\right)}^{1/3}.\end{eqnarray} \tag{ 2 }$$

For our set assumptions, the mean cluster membership is 〈N〉 = 465, so that a GMC is expected to have ∼10² clusters with radii R_c = 0.7–3.2 pc. We note, however, that the cumulative probability for finding stars in clusters with membership size N or smaller is given by the expression

$$\begin{eqnarray}&&{P}_{* }(N)=\displaystyle \frac{\mathrm{ln}[N/{N}_{\min }]}{\mathrm{ln}[{N}_{\max }/{N}_{\min }]},\end{eqnarray} \tag{ 3 }$$

indicating that stars populate clusters in a logarithmic sense. In other words, stars can be found with roughly equal probability in each decade of cluster membership size N. For our assumed cluster distribution, half of the GMC stars are expected to belong to clusters with membership N ≤ 10³.

Following the fractal generating scheme of Goodwin & Whitworth (2004), clusters are randomly placed within a cuboid GMC region with spacial fractal dimension d (with d = 2 assumed for our baseline model). The main GMC structure is evenly divided into eight substructures, the centers of which are seeded with a "child" that has a probability P = 2^d−3 of maturing. If a child does not mature, the substructure it was born in becomes inactive, and no clusters are placed within it. If the child matures, the substructure is further divided evenly into eight substructures, and the process is repeated until enough mature children exist to populate the required number of clusters. One complication in the adopted scheme is that the volume of the largest clusters exceeds that of the substructures associated with the final generation of mature children. As such, we stray from the process outlined in Goodwin & Whitworth (2004) by keeping all mature children (as opposed to keeping only the last generation). Larger clusters are placed randomly first at locations with mature children born into substructures that are minimally larger than the cluster volume, and all higher-generation "active" substructures within are removed from the pool of possible cluster sites. The process is repeated until all clusters are placed with no overlapping volume.

In our baseline model, star formation within a given cluster begins at a time chosen randomly between t = 0–10 Myr, and proceeds with equal probability over a time interval of Δt = 2 Myr. The corresponding stellar-birth distribution function for the entire GMC is therefore given by the piecewise function

$$\begin{eqnarray}b(t)=\left\{\begin{array}{ll}t/(20\,{\mathrm{Myr}}^{2}) & 0\leqslant t\leqslant 2\,\mathrm{Myr}\\ 1/(10\,\mathrm{Myr}) & 2\,\mathrm{Myr}\leqslant t\leqslant 10\,\mathrm{Myr}\\ (12\,\mathrm{Myr}-t)/(20\,{\mathrm{Myr}}^{2}) & 10\,\mathrm{Myr}\leqslant t\leqslant 12\,\mathrm{Myr}\\ 0 & t\gt 12\,\mathrm{Myr}.\end{array}\right.\end{eqnarray} \tag{ 4 }$$

The 12 Myr duration of the star formation epoch is comparable to the crossing time of the molecular cloud (for transport speed ∼1 km s⁻¹ and size scale ∼10 pc), where this time is characteristic of molecular-cloud operations (e.g., Elmegreen 2000; Dobbs & Pringle 2013).

Since SLRs are produced exclusively by massive stars (8 M_⊙ ≤ M_* ≤ 120 M_⊙), we adopt an IMF

$$\begin{eqnarray}&&f(m)=0.173\,{m}^{-2.5},\end{eqnarray} \tag{ 5 }$$

(where m = M_*/M_⊙) that mimics the high-mass portion of the observed IMF (e.g., Scalo 1998; Kroupa 2001). The normalization is chosen so that the distribution has a given probability that a star has a mass M_* ≥ 8 M_⊙, i.e.,

$$\begin{eqnarray}&&{\int }_{8}^{120}f(m){dm}=0.005.\end{eqnarray} \tag{ 6 }$$

For completeness, we note that stars near the high end of the allowed mass range can collapse directly to black holes (e.g., see the discussion of Woosley 2017). If black hole formation is efficient, then the SLRs produced during stellar evolution will be incorporated into the black holes instead of being distributed into the molecular cloud. Although the probability of black hole formation is not known, this process would reduce the radioactive yields estimated in this paper, as illustrated below (see Figure 8). ⁵

Massive stars are placed at the center of their respective cluster and evolve into SN events upon reaching an age given by

$$\begin{eqnarray}&&{\mathrm{log}}_{10}\left(\displaystyle \frac{{\tau }_{\mathrm{SN}}}{\mathrm{Myr}}\right)=\displaystyle \frac{1.4}{{\left({\mathrm{log}}_{10}m\right)}^{1.5}},\end{eqnarray} \tag{ 7 }$$

where this form represents an empirical fit (Williams & Gaidos 2007) to results from an ensemble of stellar-evolution simulations (Schaller et al. 1992).

The locations of the remaining "field" stars are then specified within each cluster using a formalism similar to that used to place clusters in the molecular cloud (with fractal dimension d = 2 for our baseline model). However, as stars are effectively point masses, they can all be placed at the active sites of the final generation of mature children, which is done randomly. In addition, the lengths of the cube in which stars are placed for a given cluster is scaled as r_*c = aR_c (with a = 2 for the baseline model), but active sites that are farther than r_*c from the cluster center are eliminated as possible locations (thereby turning a cubic structure into a spherical one). We note that this process for placing field stars does not accurately reflect the dynamics of cluster evolution, but, nevertheless, produces structures that are visually similar to those observed in GMC environments. The structure of a baseline GMC generated using our scheme is shown in the left panel of Figure 1, where blue spheres represent clusters (with radii R_c ), black dots show the location of massive stars, and gray spheres denote the location of each field star. For illustrative purposes, we also show the corresponding column-density map in the right panel of Figure 1, where a Hernquist (Hernquist 1990) density profile

$$\begin{eqnarray}&&\rho (r)=\displaystyle \frac{{M}_{\infty }}{2\pi }\displaystyle \frac{c}{r}\displaystyle \frac{1}{{\left(r+c\right)}^{3}}\end{eqnarray} \tag{ 8 }$$

with c = 0.2 pc centered on each star is used to calculate the density at each point of the region, and M_∞ is set to 20 M_⊙ so that the total mass of the GMC is 10⁶ M_⊙ (ignoring any stellar mass contribution). Noting that the mass contained within a radius r = c for a Henrnquist profile is M_∞/4, we define a characteristic density for star-forming regions as

$$\begin{eqnarray}&&{\rho }_{* }\equiv \displaystyle \frac{{M}_{\infty }/4}{4\pi {c}^{3}/3}=150\,{M}_{\odot }\,{\mathrm{pc}}^{-3}.\end{eqnarray} \tag{ 9 }$$

For ease of comparison to observed values of mass fractions in our solar system, the mass densities of ²⁶Al and ⁶⁰Fe obtained in our work will be presented in terms of ρ_*.

Zoom In Zoom Out Reset image size

Massive stars (M_* ≥ 8 M_⊙) born in GMCs, while quite rare, live only a few million years before exiting in cataclysmic fashion. The resulting SN explosions that follow the death of those massive stars spread the byproducts of the nuclear fusion that powers stellar life, as well as nuclei created through the rapid capture of neutrons during the SN event, into the surrounding GMC environment. The most massive of these stars (M_* ≥ 25 M_⊙) also inject nucleosynthesis products into the surrounding medium via strong winds during a short-lived Wolf–Rayet phase prior to their cataclysmic ends. Our main interest lies in the small fraction of those elements that are unstable and will eventually undergo radioactive decay on timescales ∼1 Myr. Specifically, ²⁶Al, and ⁶⁰Fe are produced with relatively large abundances and have half-lives t_1/2 of 0.72 and 2.6 million years, respectively—long enough for them to reach many of the stellar systems forming within the molecular cloud, but short enough that their ensuing radioactive decay can efficiently heat and ionize the disks within those systems. For completeness, we also consider the production of ²⁶Al via stellar winds during the Wolf–Rayet phase, as well as the production of the lower-yield elements ³⁶Cl (t_1/2 = 0.30 Myr) and ⁴¹Ca (t_1/2 = 0.10 Myr) in SN events.

The integrated abundance of ²⁶Al produced during the Wolf–Rayet phase of massive stars is presented in Gounelle & Meynet (2012) for several different initial masses (see their Figure 3). While the time evolution of this process is complicated, the most important feature for our study is the plateau reached just prior to the ensuing SN event. As a result, we adopt a highly idealized empirical model based on simple fits to the results of Gounelle & Meynet (2012) wherein the mass-dependent integrated abundance of ²⁶Al has the form

$$\begin{eqnarray}&&{M}_{\mathrm{Al};{\rm{w}}}(m)={10}^{-4.7+1.66({\mathrm{log}}_{10}[m]-1.48)}\,{M}_{\odot },\end{eqnarray} \tag{ 10 }$$

and remains constant over a time

$$\begin{eqnarray}&&{\tau }_{w}=\left[0.4+\displaystyle \frac{1.6(m-30)}{90}\right]\,\mathrm{Myr},\end{eqnarray} \tag{ 11 }$$

prior to the SN event for stellar masses M_* ≥ 30 M_⊙. The abundance then diminishes due to radioactive decay following the SN event. Lower-mass stars produce substantially lower yields, and are therefore disregarded. Our highly idealized model is presented in Figure 2, and can be compared directly to Figure 3 of Gounelle & Meynet (2012).

Zoom In Zoom Out Reset image size

The SN yields of the SLRs used in this work are shown in Figure 3 over the range of relevant stellar masses. These yields are taken from the stellar-nucleosynthesis calculations of Chieffi & Limongi (2013); see also Sukhbold et al. (2016) and Limongi & Chieffi (2018). These simulations provide updated yields for the SLRs and include the effects of rotation, compared with the earlier stellar-evolution models of Woosley & Weaver (1995) and Limongi & Chieffi (2006); see also Timmes et al. (1995a, 1995b), Woosley et al. (2002), and Rauscher et al. (2002). The stellar models are computed for discrete mass values, with the yields shown by the solid circles, and with intermediate values determined by interpolation (carried out in log-log space). Since results from the stellar models are listed with a minimum mass of m = 13, yields for smaller progenitor masses are extrapolated assuming constant yields over the range 8 < m < 13. The resulting yields are shown in Figure 3 for the SLRs of interest: ²⁶Al (blue curve), ⁶⁰Fe (red curve), ⁴¹Ca (green curve), and ³⁶Cl (magenta curve). For completeness, Figure 3 also includes the plateau ²⁶Al wind-generated abundances shown in Figure 2, along with the corresponding fit from Equation (10).

Zoom In Zoom Out Reset image size

For a given realization of a GMC environment, the total mass enrichment of a given element as a function of time resulting from all SN events, each occurring at time t_i , is given by

$$\begin{eqnarray}&&{M}_{A}(t)=\sum _{i=1}^{{N}_{\mathrm{SN}}}\,H[t-{t}_{i}]\,{Y}_{{A}_{i}}\,{e}^{-{\lambda }_{A}(t-{t}_{i})},\end{eqnarray} \tag{ 12 }$$

where A denotes the species, λ_A is the corresponding decay constant, and H[x] is the well-known step function. Each realization, in turn, generates a value of M_A (t) that is a member of the parent population which statistically describes mass enrichment for our assumed GMC scenario. The expectation value of the corresponding parent distribution can be calculated by integrating over both the IMF and stellar-birth functions

$$\begin{eqnarray}&&\begin{array}{c} \langle {M}_{A}(t) \rangle =N{\ast }_{{\rm{mc}}}{\displaystyle \int }_{8}^{120}f(m)\,{\displaystyle \int }_{0}^{t}\,b({t}_{\ast })H[t-{t}_{{\rm{SN}}}(m)]\,\\ \times \,{Y}_{A}(m)\,{e}^{-{\lambda }_{A}(t-{t}_{{\rm{SN}}}(m))}\,{{dt}}_{\ast }\,{dm},\end{array}\end{eqnarray} \tag{ 13 }$$

where t_* represents the stellar-birth time, ${t}_{\mathrm{SN}}={\tau }_{\mathrm{SN}}+{t}_{* }$, and b(t_*) is given by Equation (4) for our baseline model. The variance can also be calculated via the expression

$$\begin{eqnarray}&&{\sigma }^{2}(t)=\displaystyle \frac{\langle {M}_{A}{\left(t\right)}^{2}\rangle -\langle {M}_{A}(t){\rangle }^{2}}{\langle {N}_{\mathrm{SN}}(t)\rangle -1},\end{eqnarray} \tag{ 14 }$$

where

$$\begin{eqnarray}&&\begin{array}{c} \langle {M}_{A}{\left(t\right)}^{2} \rangle = \langle {N}_{{\rm{SN}}}(t) \rangle \,N{\ast }_{{\rm{mc}}}{\int }_{8}^{120}f(m)\,{\int }_{0}^{t}\,b({t}_{\ast })H[t-{t}_{{\rm{SN}}}(m)]\,\\ \times \,{\left[{Y}_{A}(m){e}^{-{\lambda }_{A}(t-{t}_{{\rm{SN}}}(m))}\right]}^{2}\,{{dt}}_{\ast }\,{dm},\end{array}\end{eqnarray} \tag{ 15 }$$

and

$$\begin{eqnarray}&&\begin{array}{c} \langle {N}_{{\rm{SN}}}(t) \rangle =N{\ast }_{{\rm{mc}}}{\int }_{8}^{120}f(m){\int }_{0}^{t}\,b({t}_{\ast })\\ H[t-{t}_{{\rm{SN}}}(m)]\,{{dt}}_{\ast }\,{dm}.\end{array}\end{eqnarray} \tag{ 16 }$$

Similar expressions can be obtained for the distribution of ²⁶Al yields produced by winds, for which a given realization of a GMC environment generates a total mass enrichment of

$$\begin{eqnarray}&&\begin{array}{c}{M}_{{\rm{Al}};{\rm{w}}}(t)=\displaystyle \sum _{i=1}^{{N}_{{\rm{SN}}}}\,H[t-{t}_{i}-{\tau }_{w}]H[{t}_{i}-t]\,\\ \,\times \,{M}_{{\rm{Al}};{\rm{w}},{\rm{i}}}+H[t-{t}_{i}]\,{M}_{{\rm{Al}};{\rm{w}},{\rm{i}}}\,{e}^{-{\lambda }_{A}(t-{t}_{i})}.\end{array}\end{eqnarray} \tag{ 17 }$$

The expectation values calculated using the above integral expressions are in complete agreement with mean value obtained from 10⁵ randomly generated realizations of GMCs for our baseline model (with values from each realization calculated using Equation (12)). Likewise, the variances calculated for the parent populations are in good agreement with those obtained from our randomly generated sample, though the latter are slightly wider owing to undersampling of the IMF at high masses. The mean values of the total mass enrichment obtained from our numerical simulations of our baseline GMC model are presented by the solid curves in Figure 4, with the corresponding shading depicting the ±1σ band. For reference, the mean values of the total mass enrichment at 5 Myr intervals of GMC age, along with the corresponding characteristic relative abundance for star-forming regions,

$$\begin{eqnarray}&&{\eta }_{A}\equiv \displaystyle \frac{\langle {M}_{A}\rangle /\langle {V}_{\mathrm{GMC}}\rangle }{{\rho }_{* }},\end{eqnarray} \tag{ 18 }$$

are compiled in Table 1, where the subscripts Al and Al + refer to SN yields only and SN plus wind yields, respectively.

Zoom In Zoom Out Reset image size

Further insight is provided by the skewness and kurtosis of the generated ²⁶Al and ⁶⁰Fe distributions, which are plotted along with the normalized first-central moments 〈M_A 〉/σ as a function of GMC age in Figure 5. Consistent with the results shown in Figure 4, the first-central moments (solid curves) increase with age, indicating that the distributions narrow with age as the death of more stars increases the sample size of SN events which contribute to the yields. The skewness and kurtosis of the distributions decrease with age from fairly large initial values to values of around 0.5 and 3.1, respectively, both of which are slightly larger than their normal counterparts of 0 and 3. These results indicate that the distributions evolve from highly non-Gaussian forms to nearly normal distributions with a median value given by ∼〈M_A 〉 − σ/6. Moreover, the tails are slightly heavier and longer at higher mass-enrichment values than normal distributions. To illustrate this point, we show in Figure 6 the histogram of the ²⁶Al mass-enrichment values generated at a GMC age of 30 Myr that correspond to the results shown in Figures 4 and 5, along with a Gaussian fit to the data. Finally, the probability functions for $\mathrm{log}[{M}_{\mathrm{Al}}]$ (with and without wind contributions) and $\mathrm{log}[{M}_{\mathrm{Fe}}]$ resulting from our analysis are shown at 5 Myr intervals of GMC age in Figure 7. While the means of the ²⁶Al and ⁶⁰Fe distributions reach their highest values at t ≈ 15 Myr when only SN events are considered, the contribution from winds maximizes the mass enrichment of ²⁶Al at t ≈ 10 Myr.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The discussion thus far assumes that stellar progenitors of all masses contribute to the production of SLRs. In some cases, however, SN explosions can stall, so that the stars directly collapse to form black holes. In such cases, the SLRs that are produced via stellar nucleosynthesis are not available for distribution throughout the cloud. Unfortunately, the masses for which black hole production occurs is not well determined. Simulations of SN detonation show complicated behavior, with some progenitor masses producing black holes and comparable masses having successful explosions. Moreover, the trends are not monotonic and different sets of simulations are not in complete agreement (e.g., Sukhbold et al. 2016; Couch et al. 2020; Boccioli et al. 2021). Nonetheless, the likelihood of black hole formation generally increases with progenitor mass. In addition, gravity-wave experiments (Abbott et al. 2016, 2021) are finding black holes in the mass range 20–120 M_⊙. After correcting for observational biases (e.g., Croker et al. 2021), the underlying mass distribution of black holes displays a broad peak at M_bh ≈ 60 M_⊙. The corresponding distribution of progenitor masses (for stars that produce black holes) is expected to peak at somewhat higher mass.

In order to assess the loss of SLRs due to black hole formation, we calculate the total supernova remnant mass enrichment for ²⁸Al and ⁶⁰Fe in our baseline model, but with the contributions from stars with mass M ≥ 100, 80, and 60 M_⊙ removed. The results are shown in Figure 8, along with the curves corresponding to the full mass range and corresponding ±σ boundaries shown in Figure 4. We note that the formation of black holes could have a noticeable effect on our results for GMC ages less than 15 Myr due primarily to the delay in the onset of SLR injection (we note that a 60 M_⊙ star remains on the main sequence for ∼1 Myr longer than a 120 M_⊙ star), but also because the total contribution of SLRs from stars of mass ≥60 M_⊙ is about 14% of the total. Nevertheless, as can be seen by comparing the curves in Figure 8 to the ±σ boundaries of the full mass range contribution, the effect is comparable to (or less than) the statistical fluctuations that arise naturally in the system. However, the cutoff at progenitor mass 60 M_⊙ could be considered conservative (see the aforementioned references). If an even larger fraction of the stellar population produces (large) black holes, then the reduction in SLRs would be correspondingly greater.

Zoom In Zoom Out Reset image size

Determining the distributions of the relative abundances of SLRs in the dense regions where stars form requires an understanding of how radioactive nuclei propagate through the GMC environments. Detailed numerical simulations of SN remnants indicate that the blast wave and ejecta propagate freely for distances of order 1–2 pc and then interact strongly with the background molecular cloud (e.g., Pan et al. 2012). The ejecta are efficiently mixed into the cloud within a relatively thin layer, with thickness ∼0.5 pc, and within a time frame of ∼30,000 yr. Subsequent movement of the radioactive nuclei must take place diffusively. ⁶ The nuclei themselves are tied to the magnetic fields of the cloud, but the gyroradius is small compared to the size scales of interest. As a result, the subsequent propagation of the SLRs occurs primarily through turbulent diffusion.

Since the mixing size and timescales are smaller than those of interest for the propagation of SLRs throughout the cloud, we assume here that each SN immediately delivers its SLRs to a thin shell with radius R = 2 pc, and calculate the ensuing density profiles (including decay) for both ²⁶Al and ⁶⁰Fe relative to the shell center by assuming point diffusion from each part of the shell. Doing so yields an expression

$$\begin{eqnarray}&&\rho (r,t)=\displaystyle \frac{2{M}_{A}{Dt}\,{e}^{-({r}^{2}+{R}^{2})/(4{Dt})}\,{e}^{-{\lambda }_{A}t}\,}{{rR}{[4\pi {Dt}]}^{3/2}}\sinh \left[\displaystyle \frac{2{rR}}{4{Dt}}\right],\end{eqnarray} \tag{ 19 }$$

for the density of the diffusing particles at radius r from the SN event at time t measured relative to that event. The diffusion constant (see Klessen & Lin 2003) is estimated through the expression

$$\begin{eqnarray}&&D=2{v}_{T}{\ell },\end{eqnarray} \tag{ 20 }$$

where v_T is the turbulent transport speed and ℓ is the mean free path. The speed v_T is identified with the observed nonthermal linewidths (Δv) of the molecular clouds, where (Δv) scales with the size L of the observed region according to

$$\begin{eqnarray}&&{v}_{R}\sim ({\rm{\Delta }}v)\sim 1\mathrm{km}{\,{\rm{s}}}^{-1}{\left(\displaystyle \frac{L}{1\,\mathrm{pc}}\right)}^{{\rm{b}}},\end{eqnarray} \tag{ 21 }$$

where the index b = 0.4–0.5 (Larson 1981; Jijina et al. 1999). Equation (21) holds over a range of size scales L = 1–30 pc, with a corresponding range of velocity scales ∼1–5 km s⁻¹. These values imply that the effective diffusion constant lies in the range D ∼ 1–150 (km s⁻¹)pc. Given the uncertainty within this possible range, we adopt D = 10 (km s⁻¹)pc for our baseline model (roughly, the geometric mean of the range), but will also explore how our results vary for other values.

In contrast to the SN ejecta, the ²⁶Al produced during the Wolf–Rayet phase of a massive star is carried by stellar winds out to the boundary where shocks stall the winds. Given that wind speeds are typically ∼100 km s⁻¹, we assume that this process instantaneously distributes the produced ²⁶Al uniformly within a sphere of radius R_w = 5 pc centered on the massive star (Gounelle & Meynet 2012; Deharveng et al. 2010).

For a given realization of our baseline model for a GMC, the densities ρ_Al and ρ_Fe at every field star location are determined at a given GMC age by summing over the contributions of all SN events. The process is then repeated for a total of 10⁴ GMC realizations (thereby building up a total sample of ≈5 × 10⁸ field point locations). Corresponding values for the relative abundances of SLRs are obtained by dividing the calculated densities by the characteristic density of a star-forming region ρ_* = 150 M_⊙pc⁻³ (see Equation (9)). The resulting probability functions for $\mathrm{log}[{\rho }_{\mathrm{Al}}/{\rho }_{* }]$ (with and without wind contributions) and $\mathrm{log}[{\rho }_{\mathrm{Fe}}/{\rho }_{* }]$ at 5 Myr intervals of GMC age are shown in Figure 9. The mean values for each distribution are indicated by the corresponding long vertical lines. For comparison, the characteristic relative abundances presented in Table 1 are shown by the short vertical lines at the top of the figure.

Zoom In Zoom Out Reset image size

The distributions depicted in Figure 9 show a number of trends. Early on, SLRs injected into the GMC have not had much time to propagate, resulting in large spatial density fluctuations characterized by broad distributions with mean values larger than the expected characteristic values. Over time, diffusion smooths out the fluctuations, and the relative abundance distributions become somewhat narrower. Taken as a whole, the distributions of SLRs span several orders of magnitude, with the bulk of the results falling in the range 10⁻¹¹–10⁻⁸ for both radioactive species of interest (where this range is consistent with previous estimates; see Vasileiadis et al. 2013; Kuffmeier et al. 2016). Significantly, the high end of the range includes solar system benchmark values (X₂₆ ∼ 4 × 10⁻⁹ and X₆₀ ∼ 10⁻⁹), although the typical mass fractions are somewhat lower.

In addition, we note that diffusion competes with radioactive decay to remove SLRs from the GMC environment with a characteristic timescale that depends on the diffusion constant,

$$\begin{eqnarray}&&{\tau }_{D}=\displaystyle \frac{{R}_{\mathrm{GMC}}^{2}}{4D}=7\,\mathrm{Myr}\,{\left(\displaystyle \frac{{R}_{\mathrm{GMC}}}{17\,\mathrm{pc}}\right)}^{2}\,{\left(\displaystyle \frac{D}{10\,(\mathrm{km}\,{{\rm{s}}}^{-1})\mathrm{pc}}\right)}^{-1},\end{eqnarray} \tag{ 22 }$$

and which is comparable to the GMC age for our baseline model. Since τ_D is larger than the half-lives of our SLRs, diffusion is not expected to have a significant effect on lowering the SLR content for the GMC for our standard scenario (notice, however, that we consider more-efficient SLR propagation out of the cloud in the following section). We also note that the difference between the distribution means and expected characteristic values become more pronounced at later ages for the longer-living ⁶⁰Fe.

In addition to the absolute values of the radioactive abundances, it is also interesting to consider the ratio of different SLRs. We complete our analysis by plotting the probability distributions for the abundance ratios $\mathrm{log}[{\rho }_{\mathrm{Fe}}/{\rho }_{\mathrm{Al}}]$ (with and without winds) calculated in the star-forming regions at our six representative ages in Figure 10. The corresponding ratios for $\mathrm{log}[{\rho }_{\mathrm{Cl}}/{\rho }_{\mathrm{Al}}]$ and $\mathrm{log}[{\rho }_{\mathrm{Ca}}/{\rho }_{\mathrm{Al}}]$ are shown in Figure 11.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The distributions of SLR abundance ratios for ⁶⁰Fe/²⁶Al shown in Figure 10 display a peak at ∼2.5 and extend down to ∼0.3. This lower value is roughly comparable to the benchmark value X₆₀/X₂₆ ≈ 0.25 that we would expect for the early solar system, if ⁶⁰Fe enrichment occurs at levels corresponding to the upper end of the measured range. For comparison, the abundance ratio inferred from gamma-ray observations (Diehl et al. 2006; Diehl 2013) is ∼0.6, which corresponds to a galactic scale and steady-state average. Although this value falls within the expected range, it is below both the peak and the mean of the distribution. We also note that these distributions are generally above the upper limits used in many considerations of the nascent solar system Gounelle (2015). One possible interpretation of this finding is that another mechanism (e.g., spallation) could have contributed to the production of ²⁶Al when our solar system formed, thereby lowering the ratio.

One limiting case for the star formation history of a GMC is to have all of the stars form at once. In this scenario, all stars form at t = 0, so that the stellar-birth distribution function is given by b(t) = δ(t). All other model parameters are the same as for our baseline model. As can be seen from Figure 12, the immediate onset of star formation injects SLRs into the GMC at earlier times, leading to the yields of both ⁶⁰Fe and ²⁶Al peaking at around 4–5 Myr (compare with Figures 2 and 3 from Voss et al. 2009).

Zoom In Zoom Out Reset image size

The resulting probability functions for $\mathrm{log}[{\rho }_{\mathrm{Al}}/{\rho }_{* }]$ (with and without wind contributions) and $\mathrm{log}[{\rho }_{\mathrm{Fe}}/{\rho }_{* }]$ at 5 Myr intervals of GMC age are shown in Figure 13. The mean values for each distribution are indicated by the corresponding lower vertical lines, and can be compared directly to their baseline model counterparts, which are shown by the upper vertical lines. Not surprisingly, instantaneous star formation produces higher SLR densities for the first 10 Myr, and lower SLR densities at later times.

Zoom In Zoom Out Reset image size

In this scenario, we consider the possibility that star formation sweeps across the GMC (e.g., Elmegreen & Lada 1977). For the sake of definiteness, we implement this scenario by correlating the time that star formation begins within a given cluster to the relative center-of-mass position of that cluster along the longest axis of the GMC, and set the sweeping time to 10 Myr. All other model parameters are the same as for our baseline model. Once star formation begins within a cluster, it proceeds with equal probability over a time interval Δt = 2 Myr, resulting in the same stellar-birth distribution function for the entire GMC as our baseline model. As such, the total SLR mass enrichment as a function of time is identical to that shown in Figure 4. However, as shown in Figure 14, the correlation between stellar position and stellar birth generally leads to slightly lower values for the probability functions of $\mathrm{log}[{\rho }_{\mathrm{Al}}/{\rho }_{* }]$ and $\mathrm{log}[{\rho }_{\mathrm{Fe}}/{\rho }_{* }]$ at star-forming sites (the exception being at 5 Myr).

Zoom In Zoom Out Reset image size

In the standard model, stars are placed within clusters, which in turn are placed within the GMC using a scheme that creates a structure with fractal dimension d = 2. We now consider what effect cloud structure has on our results by considering opposite extremes. Specifically, we first consider a model where the same scheme is used to place clusters within the GMC and stars within clusters as for the baseline model, but with a fractal dimension of d = 1.6 for placing both the clusters in the cloud and the stars within the cluster. This scenario leads to an increase in the substructure within the GMC. At the other extreme, we consider a model where all stars are placed randomly within the GMC in absence of any clusters.

For the lower-fractal-dimension scenario, star formation evolves in the same manner as our baseline model. Specifically, star formation within a given cluster begins at a time chosen randomly between t = 0–10 Myr, and proceeds with equal probability over a time interval of Δt = 2 Myr. For consistency, stellar birth in the random placement model is assigned through the same stellar-birth distribution function that is used to describe the scenario with lower fractal dimension (which in turn is the same as for our baseline model). All other model parameters are the same as for our baseline model.

As illustrated in Figure 15, lowering the fractal dimension from d = 2 to d = 1.6 has a minimal effect on how SLRs are distributed throughout the GMC. The resulting probability distributions are nearly the same. In contrast, the scenario with randomly placed stars (shown in Figure 16) leads to lower SLR abundances at stellar sites. In this random model, the stellar locations are spread out more evenly over the entire GMC and the SNe themselves are less clustered. As a result, the random distribution of stellar locations leads to less concentration of the SLRs and a corresponding deficit of high abundances.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

While the cluster distribution function chosen for our baseline model is based on our current understanding of star formation, local clusters (within 1 − 2 kpc) have been observed to have memberships that range from N = 30 to 2000 (Lada & Lada 2003; Porras et al. 2003). Motivated by this local sample, we also consider the scenario in which the cluster membership N spans a smaller range from 10² to 10³. All other model parameters are the same as for our baseline model.

Within the smaller range of stellar membership size N, clusters are drawn from the (usual) distribution function given by Equation (1). As before, we sample the distribution until the collective stellar content reaches N*_mc, with the last drawn cluster's membership reduced as necessary to achieve that outcome. For this set of assumptions, the mean cluster membership size is only 〈N〉 = 256, so that a GMC is expected to have ∼200 clusters with radii R_c = 0.7–1.5 pc. We note that for our assumed cluster distribution, half of the GMC stars are expected to belong to clusters with membership N ≤ 316.

As shown in Figure 17, the resulting distributions are shifted toward lower values compared to the standard case. The smaller clusters, which are more populous, act to distribute the stars more evenly across the cloud (but not as evenly as in the case of the random-distribution model of Section 3.3). As a result, fewer locations are near multiple SNe, and the SNe themselves are more spread out, so that fewer locations have high abundances of SLRs.

Zoom In Zoom Out Reset image size

As noted in Section 2.3, particle transport in GMC environments is poorly constrained both observationally and theoretically. We can address this issue by considering varying values for the diffusion constant, which determines the manner in which the SLRs are transported across the molecular cloud. Given the timescales of interest (t ∼ 10 Myr) and the standard value of the diffusion constant (D = 10 (km s⁻¹)pc), the typical transport distances are comparable to the typical sizes of molecular clouds ($r\sim \sqrt{{Dt}}\sim 10$ pc). As a result, larger values of the diffusion constant allow for a substantial fraction of the SLRs to be carried out of the GMC and hence become unavailable for enriching forming stars at later epochs. In contrast, smaller values lead to much less transport and allow for enhanced localized enrichment.

Here, we first consider the case where the diffusion constant has the smaller value D = 1 (km s⁻¹)pc, as shown in Figure 18. The resulting distributions of SLRs show enhanced abundances due to the smaller diffusion constant. In this scenario, the SLRs stay relatively close to their birth clusters and are effective at enriching stars forming in the vicinity. This enhancement effect is larger for ⁶⁰Fe compared to ²⁶Al due to its longer half-life.

Zoom In Zoom Out Reset image size

For larger values of the diffusion constant, D = 100 (km s⁻¹)pc, Figure 19 shows that the SLR abundances are significantly lower than in the standard baseline model. This result is expected, as the larger diffusion constant allows for a substantial fraction of the SLRs to be swept out of the cloud.

Zoom In Zoom Out Reset image size

Note that of all the parameters varied in this section, the value of the diffusion coefficient is the least constrained. With its possible range of two orders of magnitude, this quantity has the largest effect on the resulting distributions of SLRs abundances.

The use of a single diffusion coefficient to describe particle transport in a GCM is highly idealized. Indeed, recent insights about the evolution of molecular clouds indicate that massive-star feedback erodes clouds within a few Myr, and that wind-blown bubbles are important sites where radioisotopes spread. Far from homogeneous, GMCs are highly nonuniform, and exhibit thin filamentary structures extending over ∼100 pc. The spherical injection and propagation models adopted in our work are therefore too simplistic to fully describe the transport of matter in such a hot and stirred massive-star environment.

While creating a particle-transport model that takes into account the structural complexities of GMCs is beyond the scope of our work, we model the effects on nonuniformity in the GMC environment by selecting the diffusion coefficient used to determine the density at a given stellar site (via Equation (19)) from a distribution function F(D). In order to generate this distribution function, we first calculate the column density along a 4 pc line of sight extending in six directions from the position of each field star for one realization of our GMC. Doing so yields ≈300,000 values of column densities whose distribution characterizes the nonuniform structure of our GMC environments. Since the diffusion constant is expected to be inversely proportional to the mean free path, we then generate a corresponding distribution of diffusion constants through the relation D_i = 10〈N〉/N_i pc, where 〈N〉 is the mean of the column-density distribution, and N_i is the ith member of the column-density distribution. The resulting distribution of diffusion constants F(D) is shown in Figure 20.

Zoom In Zoom Out Reset image size

The spread in F(D) indicates that particle diffusion will not proceed uniformly in all directions, but rather will vary in accordance with the diffusion-constant distribution. Note that the width of the distribution is approximately a factor of 2 on either side of the mean. With this result, we perform a calculation using our baseline-model parameters, but with the diffusion constant D used to calculate the SLR densities at each star formation site (via Equation (19)) randomly selected from our generated distribution function F(D). The results of this procedure are shown in Figure 21 (solid curves) along with the corresponding distributions from the baseline model (dotted curves), for which D = 10 km s⁻¹ pc (as presented in Figure 9). The inclusion of nonuniform propagations slightly broadens the density-distribution functions (as expected) with slightly lower mean values. However, as can be seen by comparing the solid and dotted curves, the overall effect for our model is quite modest. We also note that filamentary clouds allow for more complicated transport behavior, wherein gas streaming away from one filament can reach other filaments. As a result, future work should explore this issue further, including a dynamical cloud medium shaped by feedback of winds and explosions.

Zoom In Zoom Out Reset image size

Our main conclusions can be summarized as follows:

• 1.  
  Distributed enrichment of SLRs from SN sources produces a wide distribution of abundances of radioactive nuclei across molecular clouds. For the isotope ⁶⁰Fe, we find abundance levels in the approximate range ρ_SLR/ρ_* ∼ 10⁻¹¹–10⁻⁸, where such enrichment is expected during the first ∼10 Myr of star formation. The corresponding enrichment levels for ²⁶Al are somewhat smaller but roughly comparable (see Figure 9). The distributions of SLR abundances are time dependent, and decrease steadily after the epoch of star formation has finished (here, for times later than 12 Myr).
• 2.  
  The ratio of ⁶⁰Fe/²⁶Al also displays a wide distribution (Figure 10), with typical values somewhat greater than unity (∼1–2). Significantly, SNe lead to greater enrichment levels for ⁶⁰Fe, compared to ²⁶Al, which is opposite to the trend observed/inferred for our solar system (and for the Galaxy as a whole). This trend arises because the two isotopes are produced with roughly comparable abundances in SNe, whereas ⁶⁰Fe has a significantly longer half-life (2.6 Myr, compared to 0.72 Myr for ²⁶Al). These results correspond to the typical values sampled from the entire distributions found by this paper. Any particular realization still could have more ²⁶Al than ⁶⁰Fe. For example, the scenario with efficient SLR propagation (large diffusion constant) and evaluated at late evolution times allows for larger values of ²⁶Al/⁶⁰Fe (see Figure 19).
• 3.  
  The distributions for the absolute enrichment levels indicate that this GMC-wide enrichment scenario could explain the absolute values of the abundances of SLRs inferred for the early solar system. The probability of attaining the observed solar system mass fractions depends on time (as determined by the star formation history of the cloud). For our standard scenario, forming stars have a ∼10% chance of reaching ²⁶Al enrichment levels for cloud times t = 10–15 Myr, with considerably lower probabilities outside that span of time. On the other hand, the distributions of abundance ratios ⁶⁰Fe/²⁶Al are generally not consistent with current estimates.
• 4.  
  In addition to our baseline-enrichment scenario (Figure 9), we have varied the input parameters of the model. Specifically, this paper provides SLR distributions for instantaneous star formation (Figure 13), sequential star formation (Figure 14), different fractal dimensions for cloud structure (Figure 15), uniform–random distribution of star formation sites (Figure 16), smaller clusters (Figure 17), varying SLR propagation efficiencies (Figures 18 and 19), and for the effects of nonspherical propagation (Figures 20 and 21). Over most of this parameter space, the distributions of SLRs display only modest variations. The most important variable—in terms of changing the distributions of SLRs abundances—is the diffusion coefficient for the propagation of SLRs through the molecular cloud. If the diffusion constant becomes sufficiently large (D ≳ 30 pc² Myr⁻¹), then the enrichment levels decrease substantially. The key issue required for maintaining significant levels of nuclear enrichment is the retention of the SLRs within the molecular cloud.

In terms of the absolute SLR abundances that affect forming stars and planets, we find intermediate results. The expected nuclear enrichment is large enough to provide significant sources of heating and ionization for circumstellar disks and planets, but is not large enough to completely dominate the picture (see the discussion of Lichtenberg et al. 2019, Reiter 2020, and references therein). Moreover, the distributions of SLRs abundances are wide, roughly comparable to their mean values. As a result, enrichment levels must be described in terms of probability distributions, i.e., we cannot assign a single, typical value to the expected degree of radioactive enrichment.

Although the focus of this paper is to provide the distributions of SLRs for the entire population of stars forming within a cloud, it is useful to put our solar system in context. Here, our results present a somewhat complicated picture for nuclear-enrichment scenarios. On one hand, the probability of reaching the abundance levels for ⁶⁰Fe and ²⁶Al inferred for our solar system are high enough that the observed values are not problematic. In this sense, our solar system is not out of the ordinary (consistent with previous work, from Jura et al. 2013 to Young 2020). On the other hand, the abundance ratio of iron to aluminum, ⁶⁰Fe/²⁶Al, is generally larger than unity, rather than smaller than unity as observed for the solar system (where ⁶⁰Fe/²⁶Al ≲ 0.25). This ratio is also smaller than unity on galactic scales, ⁶⁰Fe/²⁶Al ∼0.4–0.9, as measured by gamma-ray lines (Wang et al. 2020).

This study considers the enrichment provided by the transport of SLRs across molecular clouds, where the sources include SNe and stellar winds. Additional sources of radioactive material could also play a role, including spallation in the interstellar medium (Desch et al. 2010) and spallation from local sources of cosmic rays (Shu et al. 1997). Significantly, spallation can provide an additional source for ²⁶Al, whereas only stellar nucleosynthesis is known to produce ⁶⁰Fe. As a result, these additional sources of radioactive nuclei act to decrease the abundance ratio ⁶⁰Fe/²⁶Al.

Taken together, the results of this paper suggest that the expected abundances of both ⁶⁰Fe and ²⁶Al in molecular clouds are large enough to affect star and planet formation. Specifically, enrichment levels comparable to those inferred for the early solar nebula can be attained with reasonable probability, although such values fall toward the high end of the distribution. In any case, the expected enrichment levels provide significant contributions to ionization and heating in forming stars and planets. On the other hand, enrichment through SNe and stellar winds alone does not provide a full explanation for all the solar system SLRs, including their exact abundance ratios. In particular, the mass ratios ⁶⁰Fe/²⁶Al predicted from this study are generally larger than estimates for both the early solar system and the Galaxy as a whole (although the tail of the distribution includes smaller values). Additional work is thus necessary to obtain a full understanding of the SLR abundances for the early solar system, the Galaxy, and the current population of forming stars and planets.