Radioactive Planet Formation

In the scenario under consideration (Section 2), magnetic fields from the star truncate the circumstellar disk at radius r_X and cosmic rays are generated in the surrounding region. The resulting number flux Φ_CR in cosmic rays in the reconnection annulus is thus given by

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{\mathrm{CR}}\approx \displaystyle \frac{\dot{{ \mathcal N }}}{\pi {r}_{{\rm{X}}}^{2}},\end{eqnarray} \tag{ 8 }$$

where $\dot{{ \mathcal N }}$ is the total number of cosmic rays emitted per unit time.

We first consider a particular spallation reaction operating in the optically thin regime. The rate Γ at which a target nucleus absorbs cosmic rays, and thereby produces a radioactive isotope, can be written in the form

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{k}=\displaystyle \frac{1}{\pi {r}_{{\rm{X}}}^{2}}{\int }_{{E}_{0}}^{\infty }{\dot{{ \mathcal N }}}_{E}{\sigma }_{k}(E){dE}\equiv \displaystyle \frac{{\dot{{ \mathcal N }}}_{k}\langle \sigma {\rangle }_{k}}{\pi {r}_{{\rm{X}}}^{2}},\end{eqnarray} \tag{ 9 }$$

where σ_k (E) is the cross section for the reaction of interest and 〈σ〉_k is the cross section averaged over the distribution of cosmic-ray energies. The index k labels the reaction of interest. Note that the total cosmic-ray luminosity $\dot{{ \mathcal N }}$ (in particle number) includes projectiles of all types, whereas only a fraction ${\dot{{ \mathcal N }}}_{k}$ participate in the reaction.

In this application, we are interested in the amount of ²⁶Al produced, which can be synthesized through a number of spallation reactions. Specifically, reactions involving targets of ²⁷Al, ²⁶Mg, and ²⁸Si are important for proton spallation, and targets of ²⁶Mg and ²⁸Si are relevant for alpha particle radiation. Targets of ²⁴Mg, ²⁵Mg, ²⁷Al, and ²⁸Si are relevant for spallation with helion (³He) particles. These reactions are summarized in Table 1, along with defining parameters of interest (defined below).

Table 1. Nuclear Spallation Reactions for the Production of ²⁶Al

Nuclear Reactions for Spallation
Reaction	Cross Section	Target Abundance	Projectile Fraction
	〈σ〉 (mb)	XT /XH	${\dot{{ \mathcal N }}}_{\chi }/\dot{{ \mathcal N }}$
26Mg(p, n)26Al	417	4.6 × 10−6	0.9
27Al(p, pn)26Al	44	3.5 × 10−6	0.9
28Si(p,2pn)26Al	7.9	3.8 × 10−5	0.9
24Mg(α, pn)26Al	30	3.3 × 10−5	0.1
28Si(α, α pn)26Al	4.4	3.8 × 10−5	0.1
24Mg(3He,p)26Al	216	3.3 × 10−5	0.03
25Mg(3He,pn)26Al	316	4.2 × 10−6	0.03
27Al(3He,α)26Al	33	3.5 × 10−6	0.03
28Si(3He,pα)26Al	60	3.8 × 10−5	0.03

Note. The reactions are grouped according to proton reactions (top), alpha reactions (middle), and helion reactions (bottom). The second column presents the cross section averaged over the normalized cosmic-ray spectrum. The relative abundances of the target nuclei (Lodders 2003) and the projectiles are listed in the last two columns (see text).

Download table as:  ASCIITypeset image

The abundance of ²⁶Al is measured relative to a reference isotope, taken here to be ²⁷Al. Over a given exposure time ${t}_{\exp }$, the abundance ratio is given by the sum

$$\begin{eqnarray}&&{{ \mathcal R }}_{26}=\displaystyle \frac{{X}_{26}}{{X}_{27}}\approx \displaystyle \frac{\dot{{ \mathcal N }}{t}_{\exp }}{\pi {r}_{{\rm{X}}}^{2}}\sum _{k}\langle \sigma {\rangle }_{k}\displaystyle \frac{{X}_{k}}{{X}_{27}}\displaystyle \frac{{\dot{{ \mathcal N }}}_{k}}{\dot{{ \mathcal N }}}.\end{eqnarray} \tag{ 10 }$$

The first factor in the sum accounts for the cross section of the interaction, the second factor is the ratio of abundances of the target nuclei, and the third factor is the ratio of abundances of the projectiles.

Note that the above treatment represents an upper limit on the expected abundance ratio. The estimate does not include the backreaction on the SLR production rate due to the loss of target nuclei or the decay of the radioactive products. If the value of ${{ \mathcal R }}_{26}$ approaches unity, then the ratio will be smaller than assumed here due to the decline in the supply of targets, although this correction is small. In addition, however, if the exposure time ${t}_{\exp }$ becomes comparable to the decay time of the SLR, then the ratio of SLRs will be smaller due to radioactive decay (which occurs at the rate $\lambda =\mathrm{ln}2/{t}_{1/2}$). Taking into account both of these complications, we can generalize the result of Equation (10) to obtain

$$\begin{eqnarray}&&{{ \mathcal R }}_{26}=\displaystyle \frac{\dot{{ \mathcal N }}\langle \sigma {\rangle }_{T}}{\pi {r}_{{\rm{X}}}^{2}}\displaystyle \frac{1}{\lambda -{\rm{\Gamma }}}\left[1-{e}^{-(\lambda -{\rm{\Gamma }}){t}_{\exp }}\right]\unicode{x021DD}\displaystyle \frac{\dot{{ \mathcal N }}\langle \sigma {\rangle }_{T}}{\pi {r}_{{\rm{X}}}^{2}\lambda },\end{eqnarray} \tag{ 11 }$$

where we have defined a total effective cross section

$$\begin{eqnarray}&&\langle \sigma {\rangle }_{T}\equiv \sum _{k}\langle \sigma {\rangle }_{k}\displaystyle \frac{{X}_{k}}{{X}_{27}}\displaystyle \frac{{\dot{{ \mathcal N }}}_{k}}{\dot{{ \mathcal N }}}.\end{eqnarray} \tag{ 12 }$$

Note that we expect λ ≫ Γ. In the limit of short exposure times, specifically ${t}_{\exp }\ll {t}_{1/2}$, we recover Equation (10). In the opposite limit where ${t}_{\exp }\gg {t}_{1/2}$, we find that ${{ \mathcal R }}_{26}\to {\rm{\Gamma }}/(\lambda -{\rm{\Gamma }})\approx {\rm{\Gamma }}/\lambda$. In other words, the effective exposure time becomes ${t}_{\exp }\sim {\lambda }^{-1}\,={t}_{1/2}/\mathrm{ln}(2)$.

This treatment uses the spallation cross sections provided by the TENDL Nuclear Data Library (Koning et al. 2019). The tabulated cross sections listed in Table 1 are averaged over the (normalized) spectrum of cosmic rays, where we use p = 3 as the index. The resulting values are given in Table 1. Note that the experimentally determined cross sections show significant variations in their published values (compare these values with Soppera et al. 2020a, 2020b, 2020c, and with the weighted-average values given in Lee et al. 1998. ⁴ ) The overall uncertainty in the cross sections is less than a factor of 2. Moreover, the values used here (Koning et al. 2019) tend to be smaller than those of the other sources, especially at higher cosmic-ray energy, so that the estimates of this paper represent lower limits. The abundances of the target nuclei are assumed to have solar values (taken from Lodders 2003). Finally, the relative fractions of projectiles (p, α, h) are presented in the last column of the table (see Section 2).

With the specifications given above, the total spallation cross section for ²⁶Al production defined by Equation (12) becomes 〈σ〉_T ≈ 740 mb. It is useful to define the benchmark value of the abundance ratio

$$\begin{eqnarray}&&{{ \mathcal R }}_{26}^{\star }=0.12\left(\displaystyle \frac{\dot{{ \mathcal N }}}{3.75\times {10}^{34}\,{{\rm{s}}}^{-1}}\right){\left(\displaystyle \frac{{r}_{{\rm{X}}}}{0.1\mathrm{au}}\right)}^{-2}\left(\displaystyle \frac{\langle \sigma {\rangle }_{T}}{740\,\mathrm{mb}}\right).\end{eqnarray} \tag{ 13 }$$

At the expected location r_X ∼ 0.1 au, this benchmark value of the nuclear abundance ratio is much larger than the measured value in solar system meteorites. However, this large ratio can only be realized for moderate amounts of rocky material that is exposed to particle radiation under optically thin conditions over long times. In practice, energy losses within the rocky bodies reduce the abundance ratio by a substantial factor. Departures from this idealized limit are considered in the following section.

In general, the target nuclei could reside in a gaseous phase, be locked up in dust grains, or be incorporated in larger rocks. Within the current paradigm of planet formation, we expect most of the targets of interest to condense onto dust grains or larger bodies. The exposure of the targets to the cosmic-ray flux is limited in two stages. First, the cosmic-ray flux can be attenuated on its path to the rocky entities. Second, the cosmic-ray flux is attenuated further as it propagates within the rocky material. In the scenario considered here, the plasma in the reconnection region is optically thin (see the previous section) so that the first type of attenuation occurs through interactions with the sea of rocky bodies. Once the cosmic-ray flux reaches the surface of a given rock, additional attenuation of the cosmic rays takes place within the body itself.

Consider a total mass M_R of rocky material distributed across the reconnection annulus, so that the column density of the region is given by

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{\mathrm{rock}}=\displaystyle \frac{{M}_{R}}{\pi {r}_{{\rm{X}}}^{2}}\approx \,850\,{\rm{g}}\,{\mathrm{cm}}^{-2}\,\left(\displaystyle \frac{{M}_{R}}{{M}_{\oplus }}\right){\left(\displaystyle \frac{{r}_{{\rm{X}}}}{0.1\mathrm{au}}\right)}^{-2}.\end{eqnarray} \tag{ 14 }$$

Here we consider the regime where the rocky bodies with radius R are partially optically thick. As a result, the constituent nuclei are exposed to a fraction ${ \mathcal F }$ of the cosmic-ray flux that strikes the surface, and a corresponding fraction $1-{ \mathcal F }$ of the flux is absorbed. Note that the fraction depends on the rock size, ${ \mathcal F }={ \mathcal F }(R)$, and this function is estimated below. The optical depth of the reconnection annulus is given by

$$\begin{eqnarray}&&{\tau }_{\mathrm{rock}}={{\rm{\Sigma }}}_{\mathrm{rock}}\,\displaystyle \frac{3}{4\rho R}[1-{ \mathcal F }(R)],\end{eqnarray} \tag{ 15 }$$

where ρ ≈ 3.3 g cm⁻³ is the density of the rocks. As a rough approximation, we can determine the fraction of the original cosmic-ray flux that is incident on the rocky bodies by using a simple one-dimensional slab geometry. This fraction is reduced further within the rocks so that the nuclei are exposed to an overall fraction f_τ of the flux given by

$$\begin{eqnarray}&&{f}_{\tau }=\displaystyle \frac{1}{{\tau }_{\mathrm{rock}}}\left[1-{e}^{-{\tau }_{\mathrm{rock}}}\right]{ \mathcal F }(R).\end{eqnarray} \tag{ 16 }$$

The first part of the expression accounts for the attenuation of the cosmic-ray flux by the background sea of rocky material before the flux reaches the surface of a given rock. The final factor ${ \mathcal F }(R)$ accounts for the fraction of the original surface flux that reaches nuclei within the rock.

We thus need to estimate ${ \mathcal F }(R)$. The cosmic rays entering the rocky bodies lose energy in two ways: They can interact with nuclei and drive nuclear reactions such as spallation. The interaction cross sections for such processes are of order σ ∼ 10⁻²⁵ cm². They can also lose energy through Coulomb interactions with electrons in the rock. In this case, the cosmic rays continually lose energy and eventually come to a stop after encountering a sufficient column density of the material. For cosmic rays with incident energy E ∼ 10 MeV, the required stopping column density Σ_s ∼ 0.3 g cm⁻² (Reedy & Marti 1991; Reedy 2015). This stopping criterion corresponds to an effective interaction cross section given by σ_eff ∼ Am_p /Σ_s ∼ 10⁻²² cm², where m_p is the proton mass, and A is the mean atomic mass number of the rocky material. Note that cosmic rays with higher energy have a larger stopping column. Nonetheless, cosmic rays lose their energy within rocks more readily than they interact with nuclei through spallation.

As cosmic rays propagate through the rocky material, the loss of energy is given by the expression

$$\begin{eqnarray}&&\displaystyle \frac{{dE}}{{ds}}=-F{\left(\displaystyle \frac{{E}_{\mathrm{th}}}{E}\right)}^{q},\end{eqnarray} \tag{ 17 }$$

where the power-law index q ≈ 0.7 and the coefficient F depends on the composition of the target and the type of cosmic-ray particle (Reedy & Marti 1991). In addition, the coefficient is scaled such that the benchmark energy E_th is the threshold energy for the reaction of interest. This loss equation can then be integrated to find the energy of the particle remaining as a function of the distance s traveled,

$$\begin{eqnarray}&&E(s)={\left[{E}_{i}^{q+1}-(q+1){{FE}}_{\mathrm{th}}^{q}s\right]}^{1/(q+1)}.\end{eqnarray} \tag{ 18 }$$

Because the cosmic-ray energy spectrum is much steeper than the energy dependence σ(E), we can make the approximation that the cross sections are relatively constant above a threshold energy E_th . In order for a cosmic ray to have sufficient energy to interact after traveling a distance s, E(s) > E_th, and the initial energy E_i must be larger than the value

$$\begin{eqnarray}&&{E}_{i}\gt {E}_{\mathrm{th}}{\left[1+(q+1)(F/{E}_{\mathrm{th}})s\right]}^{1/(q+1)}.\end{eqnarray} \tag{ 19 }$$

Suppose that a flux of particles Φ₀ is incident on the rocky surface. For the assumed cosmic-ray energy spectrum from Equation (3), the flux of particles with energy larger than a given energy E can be written in the form

$$\begin{eqnarray}&&{\rm{\Phi }}(\gt E)={{\rm{\Phi }}}_{0}{\left(\displaystyle \frac{{E}_{0}}{E}\right)}^{p-1}.\end{eqnarray} \tag{ 20 }$$

The flux of particles above the threshold energy, evaluated at a distance s into the rocky material, is thus given by

$$\begin{eqnarray}&&{\rm{\Phi }}(s)={{\rm{\Phi }}}_{0}{\left(\displaystyle \frac{{E}_{0}}{{E}_{\mathrm{th}}}\right)}^{p-1}{\left[1+(q+1)(F/{E}_{\mathrm{th}})s\right]}^{-(p-1)/(q+1)}.\end{eqnarray} \tag{ 21 }$$

The leading factor represents the incident flux of particles with sufficient energy to induce spallation reactions. The factor in square brackets determines how the flux is attenuated with the distance s traveled through the rock.

Next we want to determine the effective cosmic-ray flux impinging on the target nuclei within a rocky body. For simplicity, the rocks of interest can be considered to be uniform density spheres of radius R. Within the reconnection region, the cosmic-ray flux can be considered isotropic. This isotropic flux of particle is thus incident on a spherical rock. For a given position within the rock (r, μ), where r is the spherical radial coordinate and the $\mu =\cos \theta$, the distance coordinate s is given by

$$\begin{eqnarray}&&s=r\mu +{\left[{R}^{2}-{r}^{2}(1-{\mu }^{2})\right]}^{1/2}.\end{eqnarray} \tag{ 22 }$$

It is useful to make the following definitions:

$$\begin{eqnarray}&&{{\ell }}_{0}=\displaystyle \frac{{E}_{\mathrm{th}}}{(q+1)F},\ \gamma =\displaystyle \frac{p-1}{q+1},\ \mathrm{and}\ \xi =\displaystyle \frac{r}{R},\end{eqnarray} \tag{ 23 }$$

where ℓ₀ ≈ 0.1 cm for E_th = 10 MeV. The target nuclei within the rock receive only a fraction ${ \mathcal F }$ of the original isotropic flux striking the surface, where

$$\begin{eqnarray}\begin{array}{rcl}{ \mathcal F }(R) & = & \displaystyle \frac{3}{2}{\left(\displaystyle \frac{{{\ell }}_{0}}{R}\right)}^{\gamma }{\displaystyle \int }_{0}^{1}{\xi }^{2}d\xi \\ & & \times {\displaystyle \int }_{-1}^{1}d\mu {\left[({{\ell }}_{0}/R)+\xi \mu +{\left[1-{\xi }^{2}(1-{\mu }^{2})\right]}^{1/2}\right]}^{-\gamma }.\end{array}\end{eqnarray} \tag{ 24 }$$

In the limit of small rocky bodies R ≪ ℓ₀, the fraction ${ \mathcal F }\to 1$ as expected. Moreover, the quantity $1-{ \mathcal F }$ $\sim { \mathcal O }(R/{{\ell }}_{0})$ in this regime. In the opposite limit R ≫ ℓ₀, the fraction decreases rapidly with increasing R. In general, we must evaluate the integrals in Equation (24) numerically.

The considerations outlined above can be summarized by writing the nuclear abundance ratio for ²⁶Al in the form

$$\begin{eqnarray}&&{{ \mathcal R }}_{26}=\displaystyle \frac{\dot{{ \mathcal N }}\langle \sigma {\rangle }_{T}}{\pi {r}_{{\rm{X}}}^{2}\lambda }\left[1-{e}^{-\lambda {t}_{\exp }}\right]\,\left[\displaystyle \frac{1-{e}^{-\tau }}{\tau }\right]{ \mathcal F }(R).\end{eqnarray} \tag{ 25 }$$

The leading factor corresponds to the case where the target nuclei experience exposure from an unattenuated cosmic-ray flux in the long-term limit under optically thin conditions. The first dimensionless correction factor in square brackets takes into account departures from saturation due to shorter exposure times. The second factor accounts for the attenuation of the cosmic-ray flux before it reaches the surfaces of the rocky bodies, where the optical depth of the reconnection region is given by τ ≈ τ_rock (see Equation (15)). The final factor results from energy loss of the cosmic flux within the rocky entities, primarily due to Coulomb interactions, as approximated using the fraction ${ \mathcal F }$ from Equation (24).

The resulting estimates for the abundance ratio ${{ \mathcal R }}_{26}$ are shown in Figure 3 as a function of the radius R of the rocky bodies. Results are shown for a range of surface densities Σ_rock, which sets the level of the optical depth and depends on the total mass of rocky material in the reconnection region (see Equation (14)). Curves are shown for total masses corresponding to the planet types of interest, specifically M_R = 0.5M_⊕ (upper red curve), 1M_⊕ (orange), 2M_⊕ (green), 4M_⊕ (blue), and 8M_⊕ (lower magenta curve). These estimates are calculated in the long-time limit so that $\lambda {t}_{\exp }\gg 1$, which is valid for exposure times of order 1 Myr or longer. For comparison, the lower dashed line shows the inferred abundance ratio for our solar system. Although the estimated abundance ratios are substantially larger than those of our solar system, they are also much smaller than the maximum benchmark value given by Equation (13).

In the limit where the rocky bodies are small, the curves in Figure 3 are flat and depend only on the total optical depth of the reconnection region (see Equation (15)). In this regime, the quantity $(1-{ \mathcal F })/R\to$ constant ∼ 1/ℓ₀, so the attenuation of the cosmic rays depends only on the total amount of material (independent of R). Notice also that most of the attenuation results from Coulomb effects, rather than nuclear reactions. As the radius R increases, the optical depth τ_rock of the background rocky material decreases. The reconnection region becomes optically thin for rock sizes R ≳ 200 cm (M_R /M_⊕). At the same time, the optical depth of each individual rock increases, with an increasing fraction of its interior shielded from the incident cosmic-ray flux. The competition between these two effects results in the broad maximum shown for R = 50–100 cm. For even larger rock radii, the abundance decreases sharply as an increasing fraction of the constituent nuclei are shielded. The enrichment levels are less than solar system values for rock sizes R ≳ 10⁴ cm and become negligible for R ≳ 1 km.

For small planets, M_plan ∼ 1M_⊕, Figure 3 indicates that enrichment levels can be more than an order of magnitude larger than those inferred for the solar system (equivalently, a maximum mass of ∼15M_⊕ of material can be enriched to solar system values). In previous studies, which focus on meteoritic enrichment, the exposure times are often smaller than those considered here, with correspondingly smaller yields (see, e.g., Shu et al. 1996 to Sossi et al. 2017). Such studies also require the enriched material to be carried out to larger heliocentric distances. This scenario—which allows for rocky material to be exposed to the largest cosmic-ray fluxes over the longest times—thus maximizes local enrichment. Nonetheless, given the finite number of energetic cosmic rays, there is a limit on the amount of material that can achieve such high abundances of ²⁶Al. Previous studies estimate that only ∼3M_⊕ of material can reach solar system abundances (Duprat & Tatischeff 2007), where this limit can be extended up to ∼10M_⊕ under favorable circumstances (Fitoussi et al. 2008). These results are comparable to the mass of ∼15M_⊕ found here (where we assume long exposure times, so that enrichment is saturated, and use a somewhat different treatment for Coulomb losses). Significantly, all of these estimates (3–15M_⊕) fall below the estimated rocky inventory of the solar system (40–80 M_⊕). This finding suggests that spallation cannot enrich the entire early solar nebula to the high isotopic abundances observed in meteorites.

For typical pebble sizes R ∼ few centimeters, the estimated abundance ratio has the approximate value ${{ \mathcal R }}_{26}\,\approx$ 0.001 ${({M}_{R}/{M}_{\oplus })}^{-1}$. For a planet with mass M_plan ≈ M_R = 1 M_⊕, for example, this enrichment level is ∼15 times larger than that inferred for our solar system. Note that smaller bodies (R ≪ 1 cm) are likely to be swept up into the accretion flow, deposited onto the stellar surface, and be unavailable for planet formation. ⁵ The key issue is thus the length of time ${t}_{\exp }$ that the pebbles remain small (∼1 cm) before being incorporated into larger planetesimals. Next we note that radioactive heating rates, and other physical quantities of interest, depend on the mass fraction X₂₆ of ²⁶Al, where

$$\begin{eqnarray}&&{X}_{26}={X}_{27}{{ \mathcal R }}_{26}\approx 4.5\times {10}^{-6}{\left(\displaystyle \frac{{M}_{R}}{{M}_{\oplus }}\right)}^{-1}.\end{eqnarray} \tag{ 26 }$$

Here we have taken X₂₆ and X₂₇ to be the mass fractions of the rocky material and have used X₂₇ ≈ 4.5 × 10⁻³ for solar system material (estimated from Table 2 of Lodders 2003). Another quantity of interest is the power per unit volume ${ \mathcal H }$ due to radioactivity, given by

$$\begin{eqnarray}\begin{array}{rcl}{ \mathcal H } & = & \lambda \,\epsilon \,{X}_{26}\left(\displaystyle \frac{\rho }{26{m}_{p}}\right){e}^{-\lambda t}\\ & \approx & 0.066\,\mathrm{erg}\,{{\rm{s}}}^{-1}\,{\mathrm{cm}}^{-3}{\left(\displaystyle \frac{{M}_{R}}{1{M}_{\oplus }}\right)}^{-1}\,{e}^{-\lambda t},\end{array}\end{eqnarray} \tag{ 27 }$$

where ≈ 4 MeV is the energy deposited per decay of the ²⁶Al nucleus. The second approximate equality evaluates the volume heating rate for standard values of the parameters.

Given these estimates for radioactive enrichment levels, it is useful to assess how the results depend on the input parameters and other uncertainties. Figure 3 shows how the predicted ²⁶Al abundances vary with the size of the rocky bodies and the total optical depth of the reconnection region (given by the total mass in solids). These curves are calculated in the regime of long exposure times ($\lambda {t}_{\exp }\gg 1$); shorter exposure times can be considered by including the dimensionless correction factor in Equation (25). Additional inputs to the estimates include the reaction cross sections, the target abundances, the relative populations of cosmic-ray species, and the power-law index of the particle radiation spectrum. All of these variables combine to specify the net cross section 〈σ〉_T (see Equation (10)). The cross sections for individual reactions 〈σ〉_k averaged over the cosmic-ray spectrum show little variation with the spectral index because the cross sections are relatively localized in energy (Duprat & Tatischeff 2007; Fitoussi et al. 2008; see also Figure 2). The target abundances scale with the metallicity (where solar values are used here). One uncertain quantity is the fraction of the cosmic-ray flux in helions, taken here to be 0.03. This fraction could vary by factors of ∼3, but the helion contribution to the overall production rate is only 10.7% so that the uncertainty in the ²⁶Al yield is also ∼10%. Perhaps the largest uncertainty is the net luminosity in energetic particles. The value used here, L_p = 10⁻⁴ L_*, provides a good estimate for the expectation value (Feigelson et al. 2005); however, the system to system scatter allows for an order-of-magnitude variation from the mean, with a corresponding variation in the expected radioactive yields.

Zoom In Zoom Out Reset image size

The timescale over which a planetesimal adjusts its internal thermal structure is governed by the Kelvin–Helmholtz time, i.e., the time required for the internal luminosity to radiate away the heat content of the planetesimal. The luminosity is given by L_in = ${ \mathcal H }M/\rho$, where M is the mass of the rocky body. The heat content is given by ${ \mathcal U }$ = CMT, where C is the specific heat capacity of the rock. The characteristic cooling timescale is thus given by

$$\begin{eqnarray}\begin{array}{rcl}{t}_{\mathrm{cool}} & = & \displaystyle \frac{{ \mathcal U }}{{L}_{\mathrm{in}}}=\displaystyle \frac{C\rho T}{{ \mathcal H }}\\ & \approx & 27,000\,\mathrm{yr}{\left(\displaystyle \frac{{ \mathcal H }}{0.066\,\mathrm{erg}\,{{\rm{s}}}^{-1}\,{\mathrm{cm}}^{-3}\,}\right)}^{-1}\left(\displaystyle \frac{T}{2000\,{\rm{K}}}\right).\end{array}\end{eqnarray} \tag{ 28 }$$

The temperature reference scale is chosen because rock melts when T ≳ 2000 K. The resulting cooling time is substantially shorter than the decay timescale for the SLRs (∼1 Myr). Note that this timescale is independent of the mass of the planetesimal.

The internal radioactive luminosity of the rocky body L_in must be transported outwards so that the temperature obeys a diffusion equation of the general form

$$\begin{eqnarray}&&\rho C\displaystyle \frac{\partial T}{\partial t}=\displaystyle \frac{1}{{r}^{2}}\displaystyle \frac{\partial }{\partial r}\left({K}_{{\rm{c}}}{r}^{2}\displaystyle \frac{\partial T}{\partial r}\right)+{ \mathcal H }(r,t),\end{eqnarray} \tag{ 29 }$$

where ${ \mathcal H }(r,t)$ is the heating due to radioactivity and K_c is the thermal conductivity (see Hevey & Sanders 2006 and references therein). We can divide out the density ρ and the specific heat C to obtain

$$\begin{eqnarray}&&\displaystyle \frac{\partial T}{\partial t}=\displaystyle \frac{\kappa }{{r}^{2}}\displaystyle \frac{\partial }{\partial r}\left({r}^{2}\displaystyle \frac{\partial T}{\partial r}\right)+\displaystyle \frac{{{ \mathcal H }}_{0}\kappa }{{K}_{{\rm{c}}}}{e}^{-\lambda t},\end{eqnarray} \tag{ 30 }$$

where ${{ \mathcal H }}_{0}$ is a constant and κ ≡ K_c/C ρ is the thermal diffusivity of the rocky material. If we assume that the initial temperature of the planetesimal is that of the ambient region at T₀, then the initial condition T(r, t = 0) = T₀. Following standard treatments of this problem (Carlslaw & Jäger 1959), the solution has the form

$$\begin{eqnarray}\begin{array}{rcl}T(r,t) & = & {T}_{0}+\displaystyle \frac{{{ \mathcal H }}_{0}\kappa }{\lambda {K}_{{\rm{c}}}}{e}^{-\lambda t}\left[\displaystyle \frac{R\sin {qr}}{r\sin {qR}}-1\right]\\ & & +\displaystyle \frac{{{ \mathcal H }}_{0}{R}^{2}}{{K}_{{\rm{c}}}}\displaystyle \frac{R}{r}\displaystyle \sum _{n=1}^{\infty }{a}_{n}\sin \left(\displaystyle \frac{n\pi r}{R}\right){e}^{-\kappa {n}^{2}{\pi }^{2}t/{R}^{2}},\end{array}\end{eqnarray} \tag{ 31 }$$

where we have defined

$$\begin{eqnarray}&&q\equiv \sqrt{\displaystyle \frac{\lambda }{\kappa }}\ \mathrm{and}\ {a}_{n}\equiv \displaystyle \frac{2}{{\pi }^{3}}\displaystyle \frac{{\left(-1\right)}^{n}}{n\left({n}^{2}-\lambda {R}^{2}/\kappa {\pi }^{2}\right)}.\end{eqnarray} \tag{ 32 }$$

The analytic solution (31) assumes that the thermal conductivity and diffusivity of the rock is constant. However, once the temperature is high enough to melt the rock, T ≳ 2000 K, the molten body can transport heat via convection and the thermal properties change.

Here we use the solution (31) to map out the parameter space for which planetesimals can melt. For a given value of the radioactive heating rate ${{ \mathcal H }}_{0}$, we find the minimum planetesimal size R that is required to achieve a given temperature. For the sake of definiteness, we measure the heating rate in terms of the benchmark value ${{ \mathcal H }}_{{SS}}=6.4\times {10}^{-3}$ erg s⁻¹ cm⁻³ that is often used for planetesimals in our solar system. ⁶ In addition, we require that 90% of the radial extent of the planetesimal reach the given temperature. Figure 4 shows the resulting contours of constant temperature in the plane of size R versus heating rate ${{ \mathcal H }}_{0}$. The curves correspond to temperatures from T = 1000 K (red) up to T = 5000 K (magenta), with the dashed curves showing intermediate temperature values. The expected temperature required for the melting of rock is about 2000 K, as depicted by the green curve in the Figure. As a result, planetesimals are expected to become almost fully molten for the parameter space above the green curve.

For the relatively large radioactive heating rates of interest, t_cool < t_1/2, so that the temperature structure of the planetesimal can reach an equilibrium state. As a result, the steady-state solution of the diffusion Equation (29) is applicable and takes the form

$$\begin{eqnarray}&&T(r)={T}_{S}+{T}_{{\rm{X}}}\left(1-\displaystyle \frac{{r}^{2}}{{R}^{2}}\right),\end{eqnarray} \tag{ 33 }$$

where T_S is the surface temperature and where the temperature scale T_x is determined by the rock properties

$$\begin{eqnarray}&&{T}_{{\rm{X}}}\equiv \displaystyle \frac{{ \mathcal H }{R}^{2}}{6{K}_{{\rm{c}}}}\approx 51\,{\rm{K}}\left(\displaystyle \frac{{ \mathcal H }}{{{ \mathcal H }}_{{SS}}}\right){\left(\displaystyle \frac{R}{1\,\mathrm{km}}\right)}^{2}.\end{eqnarray} \tag{ 34 }$$

In order for most of the planetesimal volume to be molten, the temperature scale T_x must be significantly above the melting temperature of rock. For example, if we require T > 2000 K at r/R = 0.90, the scale T_x ∼ 5000 K. For solar system heating rates, one thus requires R ≳ 10 km. As the SLR abundances increase, the critical radius for melting scales (approximately) as ${R}_{\mathrm{crit}}\sim {{ \mathcal H }}^{-1/2}$ (consistent with the slope of the curves on the right side of Figure 4).

Once a planet has formed out of the planetesimals, the internal luminosity will be enhanced due to the radioactive decay of ²⁶Al (and other SLRs). This luminosity contribution can be written in the form

$$\begin{eqnarray}\begin{array}{rcl}{L}_{\mathrm{in}} & = & \lambda \epsilon {X}_{26}\left(\displaystyle \frac{{M}_{\mathrm{plan}}}{26{m}_{p}}\right)\,{e}^{-\lambda t}\\ & \approx & 1.2\times {10}^{26}\,\mathrm{erg}\,{{\rm{s}}}^{-1}\left(\displaystyle \frac{{M}_{\mathrm{plan}}}{{M}_{R}}\right){{\rm{e}}}^{-\lambda t}.\end{array}\end{eqnarray} \tag{ 35 }$$

Here the time variable is measured from the production of the planetesimals, as radioactive enrichment becomes highly inefficient when rocky bodies reach that size. Because the planet, with mass M_plan, is made from the rocky material in the reconnection annulus, with mass M_R , we expect the two mass scales to be comparable. As a result, this scenario defines a characteristic planetary luminosity L ∼ 3 × 10⁻⁸ L_⊙, which is independent of the planetary mass (if planet formation is efficient so that M_plan ∼ M_R ).

For comparison, the power intercepted by the planet from its host star takes the form

$$\begin{eqnarray}\begin{array}{rcl}{P}_{* } & = & {L}_{* }{\left(\displaystyle \frac{{R}_{\mathrm{plan}}}{2r}\right)}^{2}\\ & \approx & 1.8\times {10}^{26}\,\mathrm{erg}\,{{\rm{s}}}^{-1}{\left(\displaystyle \frac{{R}_{\mathrm{plan}}}{{R}_{\oplus }}\right)}^{2}{\left(\displaystyle \frac{r}{0.1\,\mathrm{au}}\right)}^{-2},\end{array}\end{eqnarray} \tag{ 36 }$$

where we have used L_* = L_⊙. For Earth-sized planets, the external power is roughly comparable to the internal power for the expected abundances of SLRs. For larger planets, the external power is relatively larger. This ordering only holds for newly formed planets, however, as the internal luminosity decreases exponentially with time. Finally, note that the residual heat left over from planet formation provides another power source for young planets.

The magnetic truncation radius of T Tauri star/disk systems falls in a range centered on r_X ∼ 0.1 au. This region is coincident with a substantial population of observed exoplanets found in multiplanet systems (Figure 1).

The planet-forming material in the reconnection region is exposed to intense cosmic-ray radiation that will synthesize radioactive nuclei through spallation reactions. The resulting nuclear isotope ratio for ²⁶Al has a fiducial value ${{ \mathcal R }}_{26}\approx {10}^{-3}$, a factor of ∼20 larger than the isotope ratio inferred for our solar system (Figure 3).

Zoom In Zoom Out Reset image size

The predicted isotope ratio ${{ \mathcal R }}_{26}$ scales inversely with the total amount of rocky material M_R in the reconnection region. If we assume that the planets form from the same reservoir of material, then ${{ \mathcal R }}_{26}\propto {M}_{\mathrm{plan}}^{-1}$. The cosmic radiation is attenuated by propagating through rocky material, where Coulomb losses provide the dominant source of energy loss. This treatment does not include spallation produced via secondary particles (e.g., Pavlov et al. 2014), so the abundances presented here represent lower limits.

With the enhanced abundances of SLRs, the minimum radius necessary for planetesimals to become completely molten decreases (Figure 4). For the expected isotope ratio corresponding to M_R = 1M_⊕, planetesimals are predicted to melt for radii greater than R ≳ 3 km. For comparison, the radius needed for planetesimals to melt is about ∼20 km for the standard abundance of radioactive material.

Zoom In Zoom Out Reset image size

Because the planetesimals are expected to be entirely molten, the composition of planets forming near r ∼ 0.1 au can be different from those forming farther out. The volatile component of the planet will be determined by those materials that are soluble in magma (see the review of Chao et al. 2021). In particular, water is efficiently outgassed in molten planetesimals (e.g., Lichtenberg et al. 2019) so that planets forming in the reconnection region are predicted to be severely dehydrated. These effects should be elucidated in future work.

With the expected enhancement of SLRs, planets with masses M_plan ∼ 1M_⊕ generate significant internal luminosities due to radioactive decay. If the planets form rapidly, this internal luminosity can be larger than the power intercepted from the host stars (at r = 0.1 au and at early times). As a result, assessments of atmospheric evaporation must include internal energy sources due to radioactivity.

For the most common stars in the galaxy, with masses M_* ∼ 0.25 M_⊙, the reconnection region for cosmic-ray exposure (r_X ∼ 0.1 au) corresponds to the habitable zone (after the stars reach the main sequence). This paper shows that the raw materials that make up this class of potentially habitable planets are naturally irradiated by local cosmic rays (Figure 1), leading to spallation and radioactive enrichment (Figure 3). One potential consequence of this chain of events is that habitable-zone planets orbiting red dwarfs are depleted in water and other volatiles.

For the scenario developed in this paper, the predictions for enhanced levels of radioactive enrichment depend on the timing of planet-forming events. The rocky bodies—in the form of pebbles—must be large enough (R ≳ 1 cm) that they can reside in the reconnection region and avoid being swept into the star (keep in mind that this region is depleted in gas). At the same time, the rocks must remain small enough (R ≲ 300 cm) so that they can be efficiently irradiated and enriched. The subsequent production of planetesimals must then occur rapidly in order for the radioactive species to remain active and melt the larger rocky bodies. The exposure time ${t}_{\exp }$ of the pebbles must thus be somewhat longer than the half-life t_1/2 ∼ 0.72 Myr, whereas the time t_pls required to produce planetesimals must be shorter. The required ordering of timescales (t_pls < t_1/2 < t_exp) is typically satisfied in planet formation models that operate via the streaming instability (Drazkowska et al. 2016; Simon et al. 2016; Johansen et al. 2021). However, if departures occur, for example, shorter exposure times while the rocky bodies are small, the level of radioactive enrichment would be correspondingly smaller. A related issue is that the gas in the reconnection must be optically thin (Σ ≲ 0.01 g cm⁻²) for maximum exposure to cosmic rays, whereas the streaming instability is most efficient for comparable densities of gas and solids (Youdin & Goodman 2005). Because the surface density of solids is larger in this scenario, the growth of planetesimals by streaming instability is correspondingly slower. On the other hand, the solids are already concentrated and can collapse to form planetesimals via gravitational instability (Goldreich & Ward 1973; Youdin & Shu 2002).

Another possible way to reduce enrichment is for planets to form at more distant radial locations in the disk and then migrate inwards. Because the degree of radioactive heating can be much larger for the case of in situ formation, the composition of planets, and the corresponding the mass–radius relation, could be different for the two scenarios. Future work should thus outline the detailed chemical compositions resulting from melted planetesimals with varying levels of radioactive enrichment. Such studies could then predict how planetary properties are different for planets that migrate versus those that form locally.

In the present treatment, the solution to the diffusion equation for the planetesimal temperature assumes that the body is made of solid rock (so that the thermal diffusivity is constant). The resulting temperature profiles provide a good estimate for when and where the rocky planetesimals become hot enough to melt. After the body is sufficiently molten, however, convection carries the heat through the planetary layers where temperatures are high. Future work should thus construct evolutionary models (building on the work of Hevey & Sanders 2006) for planetesimals where the internal heat is transported convectively. The time-dependent temperature structure can then be used to determine the degree of volatile losses and hence the chemical composition of the objects. In addition, the largely molten nature of the planetesimals could affect their subsequent assembly into larger planetary bodies. Collisions between the resulting liquid entities should thus be studied. Once planets are formed, the internal luminosity from SLRs will affect the surface temperatures and the rate of atmospheric mass loss. This process should also be modeled in greater detail.

This paper has focused on the isotope ²⁶Al, which is expected to have the largest impact on planetesimals and forming planets. Other SLRs can be made through spallation and should be included in future work. Consider, for example, the case of ⁵³Mg, with half-life 3.7 Myr. The solar system abundance of this species is a factor of 10 smaller than that of ²⁶Al. If the two isotopes are produced with initial abundance ratios similar to those of our solar system, then the relative power contribution will vary with time according to L_Mn53/L_Al26 = (λ₅₃/λ₂₆) (X₂₆/X₅₃) $\exp [-({\lambda }_{26}-{\lambda }_{53})t]$ ≈ 0.02 $\exp [t/1.3\mathrm{Myr}]$. The manganese luminosity should thus dominate after ages of ∼5 Myr, although by that time the overall power will have decreased by a factor of ∼100. In addition to SLRs, long-lived radioactive nuclei can also be produced by spallation. For example, protons of sufficiently high energy (∼100 MeV) can produce ⁴⁰K through spallation reactions with iron (mostly ⁵⁶Fe, which is relatively abundant).

This paper shows that planets forming in close proximity to their host stars are exposed to enormous fluxes of particle radiation and can become substantially enriched in SLRs. The degree of radioactive enrichment will vary with the planet's location, the timing of planet-forming events, and the details of the local magnetic field structure. The resulting variations in planet properties—due to heating, differentiation, outgassing, evaporation, and other processes—will thus add to the diversity found in the observed population of extra-solar planets.