© E. Jacquet 2019

1 Introduction

Primitive meteorites, or chondrites, are conglomerates of solids that formed and accreted in the solar protoplanetary disk. Among those, the oldest are the refractory inclusions, comprising calcium-aluminium-rich inclusions (CAIs) and amoeboid olivine aggregates (AOA) (e.g., MacPherson 2014; Krot et al. 2004). Thermodynamic calculations predict that CAIs should be the first condensates in a cooling gas of solar composition (Grossman 2010; Davis & Richter 2014). Since the required temperatures would be in the range 1500–2000 K (Wooden et al. 2007), it is widely assumed that refractory inclusions formed close to the Sun, but precisely how close and in what astrophysical setting remains unclear (Wood 2004; Jacquet 2014).

An important clue in this respect may be provided by beryllium-10, a short-lived radionuclide decaying into boron-10 with a half-life of 1.5 Ma, and which, unlike other known extinct radionuclides such as aluminum-26, is not produced by stellar nucleosynthesis but may be formed through spallation by energetic particles (Davis & McKeegan 2014; Lugaro et al. 2018). Indeed, since the original discovery by McKeegan et al. (2000), all CAIs with analytically suitable Be/B have shown ¹⁰B/¹¹B excesses correlated therewith. The slopes of the resulting isochrons translate into initial (i.e., upon the last isotopic equilibration) ¹⁰Be/⁹Be ratios averaging ~6 × 10⁻⁴ in CV chondrite (type A and B) CAIs (Davis & McKeegan 2014), with less than a factor of two spread. Yet systematically lower values around 3 × 10⁻⁴ and 5 × 10⁻⁴ have been found for two CV chondrite FUN CAIs (fractionated and unknown nuclear effects; MacPherson et al. 2003; Wielandt et al. 2012) and platy hibonite crystals (PLAC) in CM chondrites (Liu et al. 2009, 2010), respectively, and, conversely, values up to 10⁻² have been found for Isheyevo CAI 411 (Gounelle et al. 2013) and more recently, fine-grained group II CV chondrite CAIs (Sossi et al. 2017). The overall spread, along with correlated ⁵⁰V excesses (also ascribed to spallation) in the latter objects (Sossi et al. 2017), is difficult to reconcile with simple ¹⁰Be inheritance from the (galactic cosmic ray-irradiated) protosolar cloud (Desch et al. 2004), which should be largely homogeneous, and argues in favor of local production in the disk following flares from the young Sun (Gounelle et al. 2001, 2006, 2013; Sossi et al. 2017), such as those manifested in X-rays by present-day protostars (e.g., Feigelson et al. 2002; Wolk et al. 2005; Preibisch et al. 2005; Telleschi et al. 2007; Güdel et al. 2007; Bustamante et al. 2016) or evidence for enhanced ionization in some protostellar envelopes (Ceccarelli et al. 2014; Favre et al. 2017, 2018).

Since the energetic protons would not penetrate further than ~690 kg m⁻² (Umebayashi & Nakano 1981) into the gas, Sossi et al. (2017) suggested that CAIs spent a few thousand orbits at the inner edge of the protoplanetary disk, similar to earlier studies that had adopted the framework of the X-wind scenario (Lee et al. 1998; Shu et al. 2001; Gounelle et al. 2001). However, among the objections to the latter discussed by Desch et al. (2010), a general issue is that CAI-forming temperatures would be expected significantly further outward (0.1–1 AU), because of local dissipation of turbulence. If nonetheless CAIs all went somehow through this narrow region and escaped accretion to the Sun or ejection from the solar system, they would be essentially “lucky” foreign material in the chondrites that incorporated them. Yet CAIs are present in carbonaceous chondrites at about the abundances predicted by in situ condensation (Jacquet et al. 2012). Also, although (non-CI) carbonaceous chondrites are enriched in refractory elements relative to CI chondrites, this cannot be ascribed to simple CAI addition to CAI-free CI chondritic material for should the CAIs bementally subtracted, these chondrites would have subsolar refractory element abundances (Hezel et al. 2008), not even considering those CAIs that were converted into chondrules (Misawa & Nakamura 1988; Jones & Schilk 2009; Metzler & Pack 2016; Ebert & Bischoff 2016; Jacquet & Marrocchi 2017; Marrocchi et al. 2018). Isotopic systematics for elements of different volatilities also indicate a need for a non-refractory isotopically CAI-like component in chondrites (e.g. Nanne et al. 2019). So there is evidence for a genetic relationship between CAIs and (part of) their host carbonaceous chondrites, which argues against an origin at the disk inner edge, whose contribution to distant chondritic matter would likely be minor.

However, the inner edge of the disk is not the only region where CAIs or their precursors could have been exposed to energetic solar protons. CAIs floating further out in the disk would at times have reached the upper layers of the disk. While such excursions may individually incur modest proton fluences, their cumulative contributions might explain the total ¹⁰Be evidenced in CAIs. Moreover, prior to CAI condensation, the gas exposed at the surface of the disk may also have undergone spallation and passed on the then-produced ¹⁰Be to the later formed condensates. The purpose of this paper is to analytically calculate the amount of beryllium-10 produced by such channels and thence evaluate the viability of having CAIs form at 0.1–1 AU from the young Sun against their irradiation record. Figure 1 provides a sketch of the overall scenario and setting. After several generalities in Sect. 2, I will express the ¹⁰Be/⁹Be ratio in the disk and free-floating solid CAIs in Sect. 3. I then discuss the numerical evaluations thereof and their implications in Sect. 4 before concluding in Sect. 5.

Fig. 1 Sketch of the scenario of 10Be production developed in this paper. Solar cosmic rays cause spallation reactions in the upper layers of the disk, in the gas and/or condensates, within a region (outlined by a dashed separatrix) where magnetic field lines originating from the proto-Sun and its outskirts hit the disk. The field and particle emission geometries are purely illustrative and do not detract from the generality of the formalism.

2 Prolegomena

2.1 Disk model

I consider aprotoplanetary disk in a cylindrical coordinate system with heliocentric distance R. The inner regions can be described under the steady-state approximation so long as the evolution timescales of the system are long compared to their local viscous timescale $$\begin{array}{lll}{t}_{\text{vis}}\left(R\right)\hfill & \equiv \hfill & \frac{R^2}{\nu }\hfill \\ \hfill & =\hfill & 0.02\text{\hspace{0.22em}Ma}{\left(\frac{R}{0.5\text{\hspace{0.22em}UA}}\right)}^{1/2}\left(\frac{1500\text{\hspace{0.22em}K}}{T}\right)\left(\frac{{10}^{-3}}{\alpha }\right),\hfill \end{array}$$(1)

where $\nu =\alpha c_{\text{s}}^2/{\Omega }_{\text{K}}$ is the effective turbulent viscosity, α the dimensionless turbulence parameter (e.g., Jacquet 2013), $c_{\text{s}}=\sqrt{k_{\text{B}}T/m}$ the isothermal sound speed with k_B the Boltzmann constant, T the temperature and m = 3.9 × 10⁻²⁷ kg the mean molecular mass, and Ω_K the Keplerian angular velocity. Since t_vis(R) is also much shorter than the half-life of ¹⁰Be, I will neglect its decay during the CAI formation epoch. This is consistent with the short timescales (no more than a few hundreds of millenia) of CAI formation and (isotopically resetting) re-heating events (e.g., MacPherson 2014), which should have ended after the last production of ¹⁰Be so as to account for the isochron behavior of the Be–B system in each CAI (that is, with ¹⁰ B excesses scaling with Be rather than the target nuclides leading to its parent; see next subsection).

If infall from the parental cloud can be neglected in the inner region, the disk mass accretion rate Ṁ = − 2πRΣu_R (with Σ the disk surface density and u_R the net, turbulence-averaged, gas radial velocity) is uniform and obeys $$\dot{M}\left(1-\sqrt{\frac{R_0}{R}}\right)=3\pi \text{Σ}\nu ,$$(2)

or, equivalently, $$u_{\text{R}}=-\frac{3\nu }{2R\left(1-\sqrt{R_0/R}\right)},$$(3)

where we have assumed that the stress vanishes at the disk inner edge R₀ (e.g., Balbus & Hawley 1998). When calculating ¹⁰Be/⁹Be in later sections, I shall show expressions for a general Ṁ(R) before specializing to the case of uniform mass accretion rate. Appendix A explores the corrections of infall to Eq. (2).

If I inject the temperature due to viscous dissipation of turbulence (e.g., Appendix A of Jacquet et al. 2012) in Eq. (2), I obtain $$\begin{array}{lll}T\hfill & =\hfill & {\left(\frac{3\kappa m{\dot{M}}^2{\left(1-\sqrt{R_0/R}\right)}^2{\Omega }_{\text{K}}{}^3}{128{\pi }^2{\sigma }_{\text{SB}}{k}_{\text{B}}\alpha }\right)}^{1/5}\hfill \\ \hfill & =\hfill & 1200\text{\hspace{0.22em}K}{\left(\frac{\kappa }{{10}^{-2}{\text{m}}^2{\text{kg}}^{-1}}\right)}^{1/5}{\left(\frac{{10}^{-3}}{\alpha }\right)}^{1/5}{\left(1-\sqrt{\frac{R_0}{R}}\right)}^{2/5}\hfill \\ \hfill & \hfill & \times {\left(\frac{0.5\text{\hspace{0.22em}AU}}{R}\right)}^{9/10}{\left(\frac{\dot{M}}{{10}^{-7}{M}_{\odot }/a}\right)}^{2/5},\hfill \end{array}$$(4)

with κ the specific Rosseland mean opacity and σ_SB the Stefan-Boltzmann constant. This confirms that CAI-forming temperatures may be reached at a fraction of an AU for Ṁ ≳ 10⁻⁷ M_⊙∕a.

2.2 ¹⁰Be production rate

The local production rate of ¹⁰Be by spallation normalized to ⁹Be may be written as (e.g., Sossi et al. 2017) $${\left(\frac{\stackrel{•}{{}^{10}\text{Be}}}{{}^9\text{Be}}\right)}_{\text{tg}}={\displaystyle {\int }_0^{+\infty }\text{d}}E{\sigma }_{k,i}(E){\left(\frac{k}{\text{Be}}\right)}_{\text{tg}}{\displaystyle \int \text{d}}\Omega I(i,E,\Omega ),$$(5)

with σ_k,i the ¹⁰Be production cross section for target (tg) nuclide k and cosmic ray species i (with energy per nucleon E) and I(i, E, Ω) the corresponding (orientation-dependent) monoenergetic cosmic ray specific intensity (defined by number and not energy, as in Gounelle et al. 2001). Here and throughout, all isotopic or elemental ratios are atomic and the Einstein convention for summation over the repeated indices k and i is adopted. I assume that the incoming cosmic rays have a solar composition (e.g. ⁴ He/H = 0.1; but see Mewaldt et al. 2007 for details on contemporaneous solar energetic particles), and obey a power law ∝ E^−p for E ≥ E₁₀ ≡ 10 MeV upon arrival on the disk. I will ignore the contribution of secondary neutrons, despite their comparable cross sections (e.g., Leya & Masarik 2009) and larger attenuation columns. This is because a dilute medium such as the surface of the disk will let free decay thwart significant accumulation of neutron flux (Umebayashi & Nakano 1981). However, based on this latter work, this is only marginally true (see Eq. (C.5)), so the production rates given here should be strictly viewed as a lower bounds.

In the following sections, we will be interested in the density (ρ)-weighted vertical average of the production rate at a given heliocentric distance, $${\langle {\left(\frac{\stackrel{•}{{}^{10}\text{Be}}}{{}^9\text{Be}}\right)}_{\text{tg}}\rangle }_{\rho }\equiv \frac{1}{\text{Σ}}{\displaystyle {\int }_0^{\text{Σ}}\text{}}\text{d}{\text{Σ}}^{\prime }{\left(\frac{\stackrel{•}{{}^{10}\text{Be}}}{{}^9\text{Be}}\right)}_{\text{tg}}\approx \frac{2}{\text{Σ}}{\displaystyle {\int }_0^{+\infty }\text{}}\text{d}{\text{Σ}}^{\prime }{\left(\frac{\stackrel{•}{{}^{10}\text{Be}}}{{}^9\text{Be}}\right)}_{\text{tg},+},$$(6)

where Σ′ is the column density integrated vertically from the upper surface and the +subscript in the final equality denotes restriction to cosmic rays from the upper side¹, assuming Σ is much larger than the attenuation column Σ_p,i(E) of the cosmic rays, defined here as $${\text{Σ}}_{p,i}(E)\equiv \frac{1}{I_{+}(i,E,\Omega ;0)}{\displaystyle {\int }_0^{+\infty }\text{d}}{\text{Σ}}^{″}{I}_{+}(i,E,\Omega ;{\text{Σ}}^{″}),$$(7)

where Σ″ is the actually traversed column density. If I envision a cosmic ray beam (with the particles gyrating along a magnetic field line) penetrating the disk with a grazing angle ϕ and ignore secondary particles, Σ″ = Σ′∕(μsinϕ) with μ the cosine of the pitch angle of the cosmic ray with respect to the magnetic field line (Padovani et al. 2018). I thus obtain $${\langle {\left(\frac{\stackrel{•}{{}^{10}\text{Be}}}{{}^9\text{Be}}\right)}_{\text{tg}}\rangle }_{\rho }=\frac{2}{\text{Σ}}{\displaystyle {\int }_0^{+\infty }\text{d}}E{\sigma }_{k,i}(E){\text{Σ}}_{p,i}(E){\left(\frac{k}{\text{Be}}\right)}_{\text{tg}}{\left(\frac{i}{\text{H}}\right)}_{\text{CR}}\frac{\text{d}{F}_H}{\text{d}E}\cdot n,$$(8)

where n is the downward-directed unit vector normal to the disk surface and d F_H ∕dE the incoming differential proton number vector flux. If I scale the latter to the incoming > 10 MeV nucleon⁻¹ energy vector flux F₁₀ and extract O as a proxy for all ¹⁰Be-producing target nuclides, Eq. (8) can be further manipulated to finally yield $${\langle {\left(\frac{\stackrel{•}{{}^{10}\text{Be}}}{{}^9\text{Be}}\right)}_{\text{tg}}\rangle }_{\rho }=\frac{K_{\text{p}}{L}_{10}}{2\pi R\text{Σ}}{\left(\frac{\text{O}}{\text{Be}}\right)}_{\text{tg}}\frac{\partial f(R_0,R)}{\partial R}$$(9)

with $$K_{\text{p}}\equiv \left(p-2\right)E_{10}^{p-2}{\left(\text{}\frac{k}{\text{O}}\text{}\right)}_{\text{tg}}{\left(\text{}\frac{i}{\text{H}}\text{}\right)}_{\text{CR}}{\displaystyle {\int }_0^{+\infty }\text{}}\text{d}E{\sigma }_{k,i}(E){\text{Σ}}_{p,i}(E)E^{-p},$$(10)

which is plotted in Fig. 2 (see Appendix B for further details on its calculation and data sources). Here, f(a, b) is the fraction of the total >10 MeV nucleon⁻¹ energy luminosity L₁₀ that reaches the disk between heliocentric distances a and b. The latter depends on the geometry of energetic particle emission and the magnetic field configuration around the young Sun, which are largely unknown (protosolar corona, star-disk fields, or disk fields; Feigelson et al. 2002; Donati & Landstreet 2009; Padovani et al. 2016). For reference, in the case of a geometrically flared disk, an isotropic point source with straight trajectories would yield f(R₀, R) = H_spall∕R with H_spall the height of the spallation layer above the midplane (see Appendix C for an estimate). For illustration in the plots of the next section, I will use a form of f assuming that the outgoing energy flux density of the energetic particles is uniform over the Sun-centered sphere of radius R₀ and that the magnetic field is approximately dipolar (with a magnetic moment along the rotation axis of the disk) for field lines interior to a maximum irradiation distance R_max. This yields² $$f(R_0,R)=\text{\hspace{0.17em}max}\left(\sqrt{1-\frac{R_0}{R}},\sqrt{1-\frac{R_0}{R_{\text{max}}}}\right).$$

However, in the main text, I will keep general expressions as functions of f. This parameterization of our ignorance will turn out to be quite useful in the calculation.

A means of comparison with previous models invoking irradiation of bare solids (e.g., Lee et al. 1998; Gounelle et al. 2006; Sossi et al. 2017) is to divide the rate given by Eq. (9) by that in a target of the same composition exposed to the same unattenuated proton flux. Assuming isotropic distribution of cosmic rays (in each half-space along the field line), this gives $$\frac{\stackrel{•}{{}^{10}\text{Be}}}{\stackrel{•}{{}^{10}\text{Be}_{\text{bare}}}}=\frac{{\langle {\text{Σ}}_p(E)\rangle }_{\sigma E^{-p}}\text{sin\hspace{0.17em}}\varphi }{\text{Σ}}$$(12)

with $$\begin{array}{lll}{\langle {\text{Σ}}_p(E)\rangle }_{\sigma E^{-p}}\hfill & \equiv \hfill & {\left(\frac{k}{\text{O}}\right)}_{\text{tg}}{\left(\frac{i}{\text{H}}\right)}_{\text{CR}}{\displaystyle {\int }_0^{+\infty }\text{d}}E{\sigma }_{k,i}(E){\text{Σ}}_{p,i}(E)E^{-p}\hfill \\ \hfill & \hfill & /\left[{\left(\frac{k}{\text{O}}\right)}_{\text{tg}}{\left(\frac{i}{\text{H}}\right)}_{\text{CR}}{\displaystyle {\int }_0^{+\infty }\text{d}}E{\sigma }_{k,i}(E)E^{-p}\right],\hfill \end{array}$$(13)

which (for a refractory target) is 90 kg m⁻² for p = 2.5 typical of gradual flares favored by Gounelle et al. (2013) and Sossi et al. (2017), and 5 kg m⁻² for p = 4 commensurate with impulsive flares (with subdominant fluences in the present-day Sun; Desai & Giacalone 2016). Obviously, for Σ ≫ 10⁴ kg m⁻² (the Minimum Mass Solar Nebula at 1 AU; Hayashi 1981) this calls for (intermittent) irradiation timescales 2 ± 1 orders of magnitude above the centuries calculated by Sossi et al. (2017). As alluded to in the Introduction, the transport timescale t_vis may, however, provide the correct order of magnitude. The following section undertakes to calculate more precisely the net ¹⁰Be abundances produced in gas and solids and during radial transport.

Fig. 2 Kp as a function of the power-law exponent of the incoming cosmic ray energy distribution, for a solar (without possible implantation contributions) and a rocky target (see Sects. 3.1 and 3.2, respectively).

3 ¹⁰Be abundances

3.1 Bulk ¹⁰Be/⁹Be in the disk

In this subsection, I calculate the average ¹⁰Be/⁹Be of the disk as a function of heliocentric distance, irrespective of its (temperature-dependent) physical state (condensation fraction) to which spallation reactions are insensitive. I assume solar abundances throughout the inner disk owing to tight coupling of the condensates with the gas (e.g., Jacquet et al. 2012).

A steady-state gradient of ¹⁰Be/⁹Be should arise in the disk because production of ¹⁰Be is balanced by loss to the Sun by accretion. The transport equation reads $$\begin{array}{l}\frac{1}{R}\frac{\partial }{\partial R}[R(\text{Σ}{u}_{\text{R}}{c}_{\text{Be}}{\left(\frac{{}^{10}\text{Be}}{{}^9\text{Be}}\right)}_{\text{disk}}-D_{\text{R}}\text{Σ}\frac{\partial }{\partial R}\left(c_{\text{Be}}{\left(\frac{{}^{10}\text{Be}}{{}^9\text{Be}}\right)}_{\text{disk}}\right))\hfill \\ =\text{Σ}{c}_{\text{Be}}{\langle {\left(\frac{\stackrel{•}{{}^{10}\text{Be}}}{{}^9\text{Be}}\right)}_{\text{disk}}\rangle }_{\rho }=\frac{K_{p,\text{disk}}{L}_{10}{c}_{\text{O}}}{2\pi R}\frac{\partial f(R_0,R)}{\partial R},\hfill \end{array}$$(14)

with c_{Be, O} the Be,O concentration in number per unit mass of solar gas and D_R the turbulent diffusion coefficient. The K_p,disk factor here includes contributions from k = ¹²C, ¹⁶ O, and ¹⁴ N (even though the last one is unimportant, with N/O = 0.14, compared to C/O = 0.5; Lodders 2003) and i = ¹H, ⁴ He. In addition, it includes contributions from indirect reactions (Desch et al. 2004), i.e. those where the heavy nuclei originate from the cosmic rays instead of the target, since, although the hereby produced ¹⁰ Be will retain the momentum of the incoming particles, it should be stopped further downstream in the disk. Since Galilean invariance dictates σ_i,k = σ_k,i, this amounts to multiplying each “direct” reaction contribution in the right-hand side of Eq. (10) by $1+\left[{(k/i)}_{\text{CR}}/{(k/i)}_{\text{disk}}\right]{\text{Σ}}_{p,k}(E)/{\text{Σ}}_{p,i}(E)=1+A_{\text{k}}/Z_{\text{k}}^2$, where the second expression uses the assumption that the disk has a (solar) composition, identical to that of the cosmic rays (see also Appendix B). The resulting K_p,disk is plotted asthe “solar target” curve in Fig. 2.

Another “indirect” contribution would be solar wind implantation of ¹⁰ Be produced onthe proto-Sun (Bricker & Caffee 2010). Although Bricker & Caffee (2010) originally envisioned implantation on bare solids, one can equally envision implantation on the gas disk surface (followed by vertical mixing); this would amount to adding $(L\text{\_{^}10\text{Be}}/2\pi R)\partial f_{\text{SW}}/\partial R$ to the right-hand side of Eq. (14), with $L\text{\_{^}10\text{Be}}$ the young Sun’s ¹⁰Be (number) production rate and f_SW the counterpart of f. This would amount to an effective increase of K_p,disk of $(\text{d}{f}_{\text{SW}}/\text{d}f)L\text{\_{^}10\text{Be}}/(L_{10}{c}_{\text{O}})$. If I set $\text{d}{f}_{\text{SW}}/\text{d}f(L\text{\_{^}10\text{Be}}/L_{10})$ equal to the present-day ratio of long-term average 1 AU ¹⁰Be (Nishiizumi & Caffee 2001)³ and >10 MeV proton (Reedy 1996) energy omnidirectional fluxes, this evaluates to 6 × 10⁻²⁰ s² m⁻². From Fig. 2, at face value, solar wind implantation thus appears negligible compared to local (in-disk) spallation (except for steep p ≳ 4), but since we do not really know how toextrapolate ¹⁰Be production on the Sun to its very active early times, the nominal one order-of-magnitude deficit may not warrant definitive conclusions yet.

Returning to Eq. (14), the requirement that ¹⁰Be/⁹Be vanishes at infinity leads to the first integration $$\left(\frac{\partial }{\partial R}-\frac{u_{\text{R}}}{D_{\text{R}}}\right)\text{left(frac{^}10\text{Be}{^}9\text{Be}right)\_disk}=\frac{K_{p,\text{disk}}{L}_{10}}{2\pi R\text{Σ}{D}_{\text{R}}}{\left(\frac{\text{O}}{\text{Be}}\right)}_{\text{disk}}f(R,+\infty ).$$(15)

Provided the radial Schmidt number Sc_R ≡ ν∕D_R has a finite lower bound, the requirement that ¹⁰Be/⁹Be does not diverge at the inner edge of the disk (since u_R∕D_R would not be integrable there from Eq. (3)) leads to the unique solution $$\text{begin{align}left(frac{^}10\text{Be}{^}9\text{Be}right)\_disk\&}=\text{\&}{K}_{p,\text{disk}}{L}_{10}{\left(\frac{\text{O}}{\text{Be}}\right)}_{\text{disk}}[\frac{f(R,+\infty )}{\dot{M}}\text{nn\&\&}-{\displaystyle {\int }_{R_0}^R\text{exp}}\left({\displaystyle {\int }_{R^{\prime }}^R\frac{u_{\text{R}}}{D_{\text{R}}}}\text{d}{R}^{″}\right)\frac{\partial }{\partial R}\left(\frac{f(R^{\prime },+\infty )}{\dot{M}}\right)\text{d}{R}^{\prime }]\text{nn\&}=\text{\&}\frac{K_{p,\text{disk}}{L}_{10}}{\dot{M}}{\left(\frac{\text{O}}{\text{Be}}\right)}_{\text{disk}}{f}_{\text{eff}},\text{nnend{align}}$$(16)

where the second equality assumes uniform Ṁ and Sc_R, with $$f_{\text{eff}}\equiv f(R,+\infty )+{\displaystyle {\int }_{R_0}^R{\left(\frac{{R^{\prime }}^{1/2}-R_0^{1/2}}{R^{1/2}-R_0^{1/2}}\right)}^{3{\text{Sc}}_{\text{R}}}}\frac{\partial f(R_0,R^{\prime })}{\partial R}\text{d}{R}^{\prime },$$(17)

which is plotted in Fig. 3.

Interestingly, the result is weakly sensitive to the details of the protoplanetary disk (e.g., turbulence level) with the reference to heliocentric distance being only implicit in f_eff. This is because the approximate ∝ R⁻² decrease of the flux density is essentially compensated by the R² factor in the radial transport timescale t_vis. The above expression can be interpreted as a sum of contributions advected from outer regions (see next subsection, in particular Eq. (19)) and contributions diffused back outward. ¹⁰Be/⁹Be decreases monotonically outward, falling off roughly as R^−3Sc_R∕2 outside the irradiated region, conforming to the equilibrium gradient of a passive scalar (Clarke & Pringle 1988; Jacquet & Robert 2013). In the limit Sc_R → 0, it converges pointwise toward the inner edge value f(R₀, +∞) and a flat profile, but Sc_R is probably of order unity in magnetohydrodynamical turbulence (Johansen et al. 2006), not to mention the possibility of laminar wind-driven accretion (and thus even higher Sc_R) further out (e.g., Bai 2016). So contrary to earlier statements (Gounelle et al. 2013; Kööp et al. 2018a), ¹⁰ Be production in the gas would entail no spatial uniformity of the ¹⁰Be/⁹Be ratio.

Fig. 3 feff (Eq. (17)) as a function of heliocentric distances for two values of Rmax for my toy parameterization of the particle flux distribution (see end of Sect. 2.2) and two values of the radial Schmidt number, with the lower one (higher relative diffusivity) in dotted lines. The parameter feff is proportional to the 10Be/9Be ratio (whose absolute magnitude is discussed at Eq. (22)) in the bulk (gas+solids) disk for given mass accretion rate and proton output and energy distributions.

3.2 In situ ¹⁰Be production in free-floating solids

I now consider a (sub-) millimeter rocky solid, say a CAI, formed (or more precisely, whose beryllium last equilibrated with the gas) at an heliocentric distance R₁. It will certainly inherit the ¹⁰Be/⁹Be ratio of the local reservoir calculated above, but it should also acquire additional ¹⁰ Be so long as itwanders through the irradiated region of the disk. The purpose of this subsection is to estimate this contribution.

Since the inner disk (where irradiation may occur) is dense, I assume the solid to be tightly coupled to the gas (that is, a gas-grain decoupling parameter S ≪ 1; Jacquet et al. 2012; this is consistent with the lack of bulk chemical fractionation assumed previously). It is also small enough for the thin target approximation (e.g., Sossi et al. 2017) to apply, that is, for Eq. (5) to apply at the scale of the whole particle as a function of the local (outside) monoenergetic intensities. Since the radial transport timescale t_vis is much longer than the vertical transport timescale H²∕ν with H = c_s∕Ω_K the pressure scale height, I can use theaverage rate of ¹⁰Be production given by Eq. (9). That is, I view the random vertical motion of the solids due to turbulence as an ergodic process (see, e.g., Fig. 8 in Ciesla 2010) and identify the time average of this rate for a given solid to the instantaneous average over the long-term probability distribution of same-size solids, which here follows the distribution of the gas (but see Appendix C for the case of finite settling)⁴. For this refractory target, the sum on the right-hand side of Eq. (10) is essentially restricted to k = O and i = ¹H, ⁴ He (with no “indirect reaction” contribution). This corresponds to the “rocky target” curve in Fig. 2.

As a result of the turbulent diffusion superimposed on advection by the gas mean velocity, the radial motion of an individual solid has a stochastic character. I will be content to set the “typical timescale” spent crossing a radial bin of width ΔR ≪ R to |ΔR∕u_R|. Indeed, even when the flow is opposite to the transport envisioned (as in this subsection), the fraction of a population of solids at R at t = 0, which has diffused upstream beyond that distance at time t ($~\text{erfc}\left((\Delta R-u_{\text{R}}t)/\sqrt{4D_{\text{R}}t}\right)/2$, with erfc the complementary error function), is maximum for t = − ΔR∕u_R. The total ¹⁰Be produced in the solid during transport between an heliocentric distance R₁ and its exit from the irradiated region is then given by $$\text{begin{align}left(frac{^}10\text{Be}{^}9\text{Be}right)\_CAI \&}=\text{\&}\int \text{\_}{R}_1\text{^}+\infty \text{leftlangleleft(frac{overset}•\text{{^}10\text{Be}}{^}9\text{Be}right)\_tgrightrangle\_}\rho \frac{\text{d}R}{\left|u_{\text{R}}\right|}\text{nn \&}=\text{\&}{K}_{p,\text{CAI}}{L}_{10}{\left(\frac{\text{O}}{\text{Be}}\right)}_{\text{CAI}}{\displaystyle {\int }_{R_1}^{+\infty }\frac{\partial f(R_0,R)}{\partial R}}\frac{\text{d}R}{\dot{M}}\text{nn \&}=\text{\&}\frac{K_{p,\text{CAI}}{L}_{10}}{\dot{M}}{\left(\frac{\text{O}}{\text{Be}}\right)}_{\text{CAI}}f(R_1,+\infty ),\text{end{align}}$$(19)

where the last equality assumes a uniform Ṁ(R). Figure 4 plots f(R, +∞). Unlike the inherited disk value (which benefited from outward diffusion), nonzero values do not extend to source radii outside the irradiated region.

Fig. 4 f(R, +∞) as a functionof heliocentric distances for two values of Rmax for my toy parameterization of the particle flux distribution (see end of Sect. 2.2). This is proportional to the 10Be/9Be ratio produced by in situ irradiation of a refractory solid as a function of its heliocentric distance of origin.

4 Discussion

4.1 Magnitude of ¹⁰Be production

Provided a given CAI underwent isotopic equilibration after the last production of ¹⁰ Be, the ¹⁰Be/⁹Be indicated by the slope of an internal isochron is the sum of that inherited from the disk gas and that produced in situ in the solid (Eqs. (16) and (19), respectively), which may be written as $$\text{left(frac{^}10\text{Be}{^}9\text{Be}right)\_final}=\frac{K_{\text{p},\text{disk}}{L}_{10}{f}_{\text{sum}}}{\dot{M}}{\left(\frac{\text{O}}{\text{Be}}\right)}_{\text{disk}}$$(20)

with $$\begin{array}{lll}{f}_{\text{sum}}\hfill & \equiv \hfill & f(R_1,+\infty )\left(1+\frac{K_{\text{p},\text{CAI}}{\left(\text{O}/\text{Be}\right)}_{\text{CAI}}}{K_{\text{p},\text{disk}}{\left(\text{O}/\text{Be}\right)}_{\text{disk}}}\right)\hfill \\ \hfill & \hfill & +{\displaystyle {\int }_{R_0}^{R_1}{\left(\frac{{R^{\prime }}^{1/2}-R_0^{1/2}}{R_1^{1/2}-R_0^{1/2}}\right)}^{3{\text{Sc}}_R}}\frac{\partial f(R_0,R)}{\partial R}\text{d}{R}^{\prime }.\hfill \end{array}$$(21)

The contribution of in situ production (the term with the “CAI” subscripts in the above equation) is subdominant with respect to spallation in the solar gas (that is, f_sum ≈ f_eff(R₁)) since, not to mention the ignored incipient settling (see Appendix C), (i) K_p,disk outweighs K_p,CAI by a factor ~3 (Fig. 2) as it includes additional production channels (e.g., target ¹² C, and indirect reactions) and (ii) Be, as a refractory element (half-condensation temperature of 1452 K according to Lodders 2003), should be concentrated relative to O in a CAI with respect to the overall disk. Indeed Be concentrations typically are of order 0.1–1 ppm in the data of McKeegan et al. (2000) and Gounelle et al. (2013), compared to the CI chondrite value of 25.2 ppb (the latter amounting to O∕Be = 10⁷ ≈ half the solar value; Palme et al. 2014). Nevertheless, the fine-grained CAIs of Sossi et al. (2017) exhibit concentrations down to 6 ppb⁵. So, at least for some Be-poor CAIs having spent a random walk in the irradiated region of the disk longer than the “typical timescale” used in Sect. 3.2, in situ production may not be entirely negligible. The deviation of the initial boron isotopic composition from solar in some CAIs, too large for an irradiated solar gas, may be a collateral effect of such a contribution (Liu et al. 2010; Gounelle et al. 2013); the same applies for cosmogenic ³ He and ²¹ Ne excesses found by Kööp et al. (2018a) in the same type of inclusions. This indicates that, at least for some objects, heliocentric distance R₁ was inside the proton-irradiated region of the disk.

In order to numerically evaluate Eq. (20), I will, as previous authors (Lee et al. 1998; Gounelle et al. 2001; Sossi et al. 2017), use the X-ray luminosity L_X as a proxy for L₁₀, since only the former can be measured from distant young stellar objects (although Ceccarelli et al. (2014) estimated a proton output equivalent to L₁₀ ≳ 10²⁷ W, for L_X < 7 × 10²⁴ W, from the ionization level seen toward the single source OMC-2 FIR 4). Since their luminosities are higher than the most powerful flares of the contemporaneous Sun (e.g., Feigelson et al. 2002), solar flares are the least improper analogs in that calibration. Lee et al. (1998) derived a scaling L₁₀ = 0.09L_X by ratioing the total proton and X-ray fluences of the impulsive flares of solar cycle 21 (peaking around 1980). However, while impulsive flares correlate best with X-ray outputs (Lee et al. 1998), gradual flares actually dominate the proton fluences at 1 AU (Desai & Giacalone 2016), whatever the ratio at the X point considered by Lee et al. (1998) may be. This warrants re-examination of the scaling between X-rays and solar energetic particles (SEP) without restricting to impulsive flares. Emslie et al. (2012) evaluated the energy outputs of 38 solar eruptive events between 2002 and 2006. Those 21 with reported SEP outputs totaled a SEP output of 1.5 × 10²⁵ J (dominated by >10 MeV protons since the energy spectra is shallower below this; Mewaldt 2006) and 2.6 × 10²⁴ J in X-rays, hence a ratio of average L₁₀∕L_X ≈ 6 which I adopt as a normalizing value⁶. Another relevant – if more indirectly for our purpose – dataset is the catalog of 314 SEP events between 1984 and 2013 of Papaioannou et al. (2016). They represent a total >10 MeV proton 1 AU energy fluence of 5.1 × 10³ J m⁻² and total X-ray fluence of 54 J m⁻² hence (L₁₀∕L_X)∂f∕∂ ln R∕sin ϕ ≈ 23 (which averages different geometries of the flare sources vis-à-vis the Earth) assuming isotropic distribution of SEP at 1 AU, suggestive of the same order of magnitude for L₁₀∕L_X. Equation (20) then becomes $$\text{begin{align}left(frac{^}10\text{Be}{^}9\text{Be}right)\_final\&}=\text{\&}6\times {10}^{-4}\left(\frac{K_{p,\text{disk}}}{9\times {10}^{-19}{\text{s}}^2{\text{m}}^{-2}}\right)\left(\frac{{\left(\text{O/Be}\right)}_{\text{disk}}}{2.2\times {10}^7}\right)\left(\frac{f_{\text{sum}}}{0.1}\right)\text{nn\&\&}\times \left(\text{}\frac{L_{10}/L_X}{6}\text{}\right)\left(\frac{L_X}{3\times {10}^{23}\text{\hspace{0.22em}W}}\right)\left(\frac{{10}^{-7}{M}_{\odot }/a}{\dot{M}}\right).\text{nnend{align}}$$(22)

The nominal value obligingly coincides with the average one measured for regular CAIs in CV chondrites (Davis & McKeegan 2014). This, however, should not hide the considerable uncertainties in several factors, in addition to the L₁₀ ∕L_X discussed above. In particular, K_p,disk, normalized to the value for p = 2.5 typical of gradual flares favored by Sossi et al. (2017), may lose one order of magnitude or so (depending on the possible implantation contributions) for a steeper energy distribution; on the other hand, the normalization value for f_sum, inspired from the case of ballistic emission, could be an underestimate if the cosmic rays are focused toward the disk (see Fig. 3), but it depends on the essentially unknown magnetic field configuration and turbulent diffusion efficiency. While the predictive value of the model should not thus be overrated, this result does show that, in the current state of our knowledge, the evidence of extinct ¹⁰ Be in refractory inclusions does not mandate an origin at the very inner edge of the disk, and formation over a wider range of heliocentric distance, say 0.1–1 AU, can be envisioned. This is the main point of this work.

4.2 Spatio-temporal variations of the ¹⁰Be/ ⁹Be ratio

In this model, ¹⁰Be/⁹Be is proportional to the L_X∕Ṁ that is an observable. Since L_X has only a shallow dependence on time (e.g., ∝t^−0.36 for one solar mass according to Telleschi et al. 2007), compared to the decrease of the accretion rate (∝ t^−1.4 according to Hartmann et al. 1998), L_X∕Ṁ should increase over time (as t^1.04 if I combine these examples although derived from different data sources; see Fig. 5). It may indeed be verified in Fig. 6 that L_X∕Ṁ anticorrelates with Ṁ, if with a fair amount of scatter, which reminds us of the elusive determinants of L_X. Since with decreasing Ṁ isotherms should recede toward the Sun (see Eq. (4)), f_sum may also be expected to increase for a given T(R₁) (e.g., a condensation front ) for a fixed magnetic field configuration. So, barring a strong decrease of L₁₀ ∕L_X, ¹⁰Be/⁹Be should increase with time for a given formation and/or equilibration temperature. That is, mass loss from the disk would increase the relative importance of the surficial layers prone to spallation, overcoming the decline in proton luminosity.

This expectation would provide an explanation for the lower ¹⁰Be/⁹Be of FUN CAIs in CV chondrites (MacPherson et al. 2003; Wielandt et al. 2012) and PLACs in CM chondrites (Liu et al. 2010) compared to “regular” CAIs. Indeed these CAIs are characterized by nucleosynthetic stable isotopic anomalies (e.g., in Ti, Ca; MacPherson 2014; Kööp et al. 2018b) much in excess of those seen in bulk chondrules (Gerber et al. 2017) or meteorites (Trinquier et al. 2009). This suggests they formed early, perhaps during infall of an isotopically heterogeneous protosolar cloud, before turbulent mixing essentially suppressed the heterogeneities (e.g., Boss 2012). Perhaps the L_X ∕Ṁ was so low that their ¹⁰Be/⁹Be was in fact dominated by the background inherited from the protosolar cloud (with perhaps some irradiation from the protostar directly on the envelope before arrival on the disk; Ceccarelli et al. 2014) rather than later dominating local spallation (Liu et al. 2010; Wielandt et al. 2012).

At the other end of the spectrum, a relatively late formation might explain the high ¹⁰Be/⁹Be (7 × 10⁻³) of fine-grained CV chondrite CAIs analyzed by Sossi et al. (2017) or that of Isheyevo CAI 411 (10⁻²; Gounelle et al. 2013). Indeed the former CAIs have group II rare earth element patterns that are commonly ascribed to condensation in a reservoir previously depleted in an ultrarefractory condensate (e.g., Boynton 1989). One may then speculate that this preliminary fractionation took some time to complete (given the tight coupling of millimeter-sized solids and gas) so group II CAIs had to be of relatively late formation. However, while Isheyevo CAI 411 shows ²⁶ Al/²⁷Al < 10⁻⁶ – an ²⁶Al depletion common to many CH chondrite CAIs –, the fine-grained CV chondrite CAIs tend to exhibit ²⁶ Al/²⁷Al ratios no lower than their melted counterparts (MacPherson et al. 2012; Kawasaki et al. 2019), unlike what a later formation would suggest in view of the ²⁶Al half-life of 0.72 Ma (e.g., MacPherson et al. 2012). This conclusion is, however, predicated on an homogeneous ²⁶ Al distribution, which is far from a foregone assumption in early times. In fact, the nucleosynthetic anomaly-bearing CAIs discussed above show reduced ²⁶ Al/²⁷Al (~ 10⁻⁵; MacPherson et al. 2014; Park et al. 2017) indicating an increase of that ratio in the CAI-forming region between their formation epoch and that of melted CAIs, and perhaps somewhat beyond. Interestingly also, some ²⁶ Al should be produced by spallation along with ¹⁰Be; for their fine-grained CAIs, Sossi et al. (2017) estimated an increase of ²⁶ Al/²⁷Al of (0.3 − 1.2) × 10⁻⁵, comparable to the spread resolved by MacPherson et al. (2012) and Kawasaki et al. (2019). At any rate, this calls into question the usefulness of ²⁶ Al as a chronometer during CAI formation, although it may have become more homogenized afterward. Alternatively, the coarse-grained CAIs may have acquired ¹⁰Be/⁹Be lower than their fine-grained counterparts because of equilibration with their surroundings at higher heliocentric distances, during the localized heating events that melted them.

This very dependence of ¹⁰Be/⁹Be on heliocentric distance also implies that the comparatively younger chondrules (or other meteoritic material), which likely formed several AUs away from the Sun (so probably outside the range for in situ irradiation), need not have high ¹⁰Be/⁹Be, not to mention effects of finite time of transport or radioactive decay, as timescales for the chondrule formation epoch (0–3 Ma after CAIs; Nagashima et al. 2018; Connelly & Bizzarro 2018) or their transport are comparable to the half-life of ¹⁰ Be. For a Schmidt number of order unity, one might expect a 2 ± 1 orders of magnitude depletion with respect to the CAI-forming region (one order of magnitude closer to the Sun), but it is difficult to narrow down. Still, no chondrule ¹⁰Be/⁹Be values are known yet since silicate crystallization there hardly fractionates Be from B (Davis & McKeegan 2014).

Fig. 5 LX∕Ṁ for Taurusas a function of age (data from Güdel et al. 2007). Although a least-square fit (dashed line) shows a systematic increase with age, scatter is evident, in part owing to the diversity of stellar masses. Overplotted is the curve (red, continuous) corresponding to one solar mass (Telleschi et al. 2007; see beginning of Sect. 4.2). For reference, the normalizations used in Eq. (22) amount to LX ∕Ṁ = 107.7 J kg−1.

Fig. 6 LX∕Ṁ for young stellar objects as a function of Ṁ. Taurus and Orion data are taken from Güdel et al. (2007) and Bustamante et al. (2016), respectively (where LX ∕Ṁ also shows little systematic trend with stellar mass).

4.3 Implications for vanadium isotopic systematics

While it is beyond the scope of the current work to calculate the effects of proton irradiation on other isotopic systems, some comments on the reinterpretation of the vanadium-50 excesses (up to 4.4o) reported by Sossi et al. (2017) are in order. In this framework, spallogenic ⁵⁰V would also be dominated by inheritance from the nebula, rather than in situ production, but to a lesser extent as the target nuclei/element ratio (essentially Ti/V playing the role of O/Be; Sossi et al. 2017), between two refractory elements, would not markedly differ between a solar gas and a condensate (see second paragraph of Sect. 4.1). So the contention by Kööp et al. (2018a) that lack of evidence of in situ production of noble gases by solar energetic production in fine-grained CAIs analyzed by Vogel et al. (2004)⁷ prevents a spallogenic origin for their ⁵⁰V excesses andfavors mass-dependent fractionation upon condensation is not necessarily valid. Since differential attenuation of the energetic particles modifies their energy distribution at finite penetration columns, the production ratio (independent of the fluence, and only dependent on p) of two nuclides would be affected. In particular, for a given p, since the reactions producing ⁵⁰V are efficient at lower energies than for ¹⁰Be (Lee et al. 1998), where attenuation is strongest, ⁵⁰V production would be comparatively lower. Thus, for given (measured) ¹⁰Be/⁹Be and δ⁵¹V, a steeper slope than in the formalism of Sossi et al. (2017) would be indicated (for a given target composition, which however would have to be revised to solar). Since, if we ratio Eq. (8) with its ⁵⁰ V counterpart, the only difference with the Sossi et al. (2017) formalism is the presence of Σ_p (E) in the integrals in the denominator and the numerator, one may surmise that a shift of p comparable to the power-law exponent of Σ_p(E) in the Continuous Slowing Down Approximation (about 1.82; see Appendix B) would largely cancel out its effect and restore the observed ratio. Although a steeper p would diminish ¹⁰ Be production for a given L₁₀, this would hardly affect our conclusions given the other uncertainties discussed at the end of Sect. 4.1. This also does not take into account the neutron contribution alluded to in Sect. 2.2, which may alter inferences on the energy distribution. Clearly, a dedicated study on the simultaneous effects expected for different isotopic systems in this scenario, whose proportions are independent of the uncertainties on the fluences, would be worthwhile.

4.4 X-rays and D/H fractionation

As mentioned above, X-ray emission should accompany cosmic ray flares from the Sun. These could also leave isotopic fingerprints (though unrelated to spallation) in meteorites. Indeed Gavilan et al. (2017) linked deuterium enrichment of chondritic organic matter (previously ascribed by Remusat et al. 2006, 2009 to ionizing radiation) to X-ray irradiation near the surface of the disk. The energy fluence F expected for organic matter wandering in the outer disk can be calculated with the same formalism as above if I replace ${\sigma }_{k,i}(E){(i/\text{H})}_{\text{CR}}{(k/\text{Be})}_{\text{tg}}$ by E, so I end this discussion with this short aside. For an attenuation column ${\text{Σ}}_{\text{X}}(E)=0.15\text{\hspace{0.22em}kg\hspace{0.17em}}{\text{m}}^{-2}{(\text{E/1\hspace{0.22em}keV})}^n$ with n = 2.485 and an energy distribution $E{I}_{\text{X}}\propto \text{exp}\left(-E/(k_{\text{B}}{T}_{\text{X}}\right))$ with T_X the emission temperature (Igea & Glassgold 1999), the energy-weighted average of the former (the equivalent of K_p ) is Σ_X (k_BT_X)Γ(n + 1). For distant part of the (flared) disk, I can use f(R₀, R) = H_X∕R (with H_X the X-ray absorption height, likely a few times H by analogy with Appendix C) so that the net fluence has the fairly well-constrained expression (modified after Eq. (19)) $$\begin{array}{lll}F\hfill & =\hfill & \Gamma (n+1){\text{Σ}}_{\text{X}}\left(k_{\text{B}}{T}_{\text{X}}\right)\frac{L_{\text{X}}}{\dot{M}}\left[\frac{H_{\text{X}}}{R}\right]\hfill \\ \hfill & =\hfill & 2\times {10}^7\text{\hspace{0.22em}J\hspace{0.17em}}{\text{m}}^{\text{−2}}\left(\frac{\Gamma (n+1)}{3}\right)\left(\frac{{\text{Σ}}_X(1\text{\hspace{0.22em}keV})}{0.15\text{\hspace{0.22em}kg\hspace{0.17em}}{\text{m}}^{-2}}\right){\left(\frac{k_{\text{B}}{T}_{\text{X}}}{1\text{\hspace{0.22em}keV}}\right)}^n\hfill \\ \hfill & \hfill & \times \left(\frac{L_{\text{X}}}{3\times {10}^{23}\text{\hspace{0.22em}W}}\right)\left(\frac{{10}^{-8}{\text{M}}_{\odot }/a}{\dot{M}}\right)\left(\frac{\left[H_{\text{X}}/R\right]}{0.1}\right),\hfill \end{array}$$(23)

with $\left[H_{\text{X}}/R\right]$ the difference of H_X∕R between the locus of formation of the organic matter (or the location inward of which settling ceases to be important, i.e. the S = 1 line of Jacquet et al. 2012; see Appendix C) and that of accretion in the chondrite parent body. This is five orders of magnitude short of the critical fluence of 5 × 10²⁷ eV cm⁻² = 8 × 10¹² J m⁻² indicated by the experiments of Gavilan et al. (2017) for 0.5–1.3 keV photons (and is also much shorter than the astrophysical estimate of Gavilan et al. 2017, which erroneously used X-ray attenuation at surface layers as representative of the bulk). While the steady-state approximation may overestimate the surface density in the outer regions of interest (which nevertheless should not allow efficient settling of the grains), alleviating it would hardly bridge the gap. Nevertheless, the quantitative assessment of isotopic effects of irradiation, whether X-rays (in particular at higher energies) or other parts of the spectrum, is still in its infancy and this additional calculation is essentially intended for future applications in similar scenarios.

5 Conclusion

I have analytically investigated the production of short-lived radionuclide beryllium-10 in surface layers of the disk irradiated by protosolar flares. I found that ¹⁰Be production in the gas outweighs ¹⁰Be production in solids after condensation because the gas contains a greater breadth of suitable target nuclides (e.g., ¹² C, ¹ H) and incurs less dilution by stable refractory Be. Taking into account incipient settling and possible implantation of solar wind-borne ¹⁰ Be would further widen the difference. Although many uncertainties remain on the magnetic field configuration, the scaling of cosmic rays with X-ray luminosities and other factors, it does appear that this model can reproduce ¹⁰Be/⁹Be ratios measured in CAIs. Therefore, the past presence of ¹⁰ Be does not require (at least at present) that CAIs formed at the inner edge of the disk and allows formation at a fraction of an AU, as thermal models would suggest, and more in line with evidence for a genetic link with their host carbonaceous chondrites (abundance,fraction of the refractory budget; Jacquet et al. 2012). If this model holds true, an interesting corollary (barring strong variations of the energetic protons/X-ray ratio) is that the oldest CAIs should have the lowest ¹⁰Be/⁹Be ratios, which would explain those of nucleosynthetic anomalies-bearing CAIs ((F)UN, PLAC). This would also suggest that the fine-grained group II CAIs in CV chondrites measured by Sossi et al. (2017) were a relatively late generation of refractory inclusions, a possibility that remains to be explored. This does not mean that chondrules, which formed at comparatively much larger heliocentric distances, should have high ¹⁰Be/⁹Be, since this ratio decreases outward, following passive diffusion outside the irradiated inner disk. I finally note that the same formalism allows an estimate of fluences of X-rays (produced in the same protosolar flares) or other types of radiations, on aggregates freely floating in the disk which can be, for example, compared to experimental evidence of D/H fractionation in meteoritic organic matter by irradiation.

Acknowledgments

I am grateful to the organizers and participants of the workshop “Core to Disk”, a program of the ψ2 initiative that took place at the Institut d’Astrophysique Spatiale in Orsay from May 14 to June 22, 2018, in particular Patrick Hennebelle and Matthieu Gounelle, who elicited my interest in the subject. I also thank Dr Ming-Chang Liu, Manuel Güdel, and Thomas Preibisch who kindly answered my requests for information. Comments by the anonymous reviewer greatly improved the clarity of the result sections as well as the discussion of other isotopic systems. This work is supported by ANR-15-CE31-004-1 (ANR CRADLE).

Appendix A Steady-state disk with infall

In this appendix, I study the effects of infall on the steady-state solution for the disk. I adopt the Cassen & Moosman (1981) expression for the infall source term per unit area ${\dot{\text{Σ}}}_{\text{in}}$ for a cloud in solid-body rotation, assuming bipolar jets do not suppress it in the inner disk. The mass accretion rate is no longer constant but obeys $$-\frac{\partial \dot{M}}{\partial R}=2\pi R{\dot{\text{Σ}}}_{\text{in}}=\frac{\stackrel{.}{M_{\text{in}}}}{2R_{\text{C}}\sqrt{1-R/R_{\text{C}}}}\theta (R_{\text{C}}-R),$$(A.1)

with R_C the instantaneous centrifugal radius, Ṁ_in the total infall rate on the star+disk system, and θ the Heaviside function. Treating all matter arriving inside R = R₀ as directly accreted by the star, this can be integrated as $$\dot{M}={\dot{M}}_{\star }-{\dot{M}}_{\text{in}}\left(1-\sqrt{1-\frac{R}{R_{\text{C}}}}\right),$$(A.2)

with Ṁ_⋆ the mass accretion rate of the star. If the whole disk can be treated as in steady state, we should have Ṁ_⋆ = Ṁ_in but I refrain from making this identification immediately as the steady-state approximation may only hold locally in the inner regions and/or in the absence of infall (Ṁ_in = 0) altogether (as assumed in the main text).

Angular momentum conservation reads (Cassen & Moosman 1981) $$\dot{M}=6\pi R^{1/2}\frac{\partial }{\partial R}\left(R^{1/2}\text{Σ}\nu \right)+2\pi R^2\text{Σ}\frac{{\dot{M}}_{\star }}{M_{\star }}+4\pi R^2{\dot{\text{Σ}}}_{\text{in}}\left(1-\sqrt{\frac{R}{R_{\text{C}}}}\right),$$(A.3)

with M_⋆ the stellar mass. If the latter is much bigger than the mass of the disk, the middle term on the right-hand side can be neglected upon integration (after division by R^1∕2), which yields $$\left({\dot{M}}_{\star }-{\dot{M}}_{\text{in}}\right)\left(1-\sqrt{\frac{R_0}{R}}\right)+{\dot{M}}_{\text{in}}{\left(\frac{R_{\text{C}}}{R}\right)}^{1/2}{\left[\sqrt{1-x}\left(x^{1/2}-\frac{x+2}{3}\right)\right]}_{R_0/R_{\text{C}}}^{R/R_{\text{C}}}=3\pi \text{Σ}\nu .$$(A.4)

For R ≪ R_C this approximates to $$3\pi \text{Σ}\nu \approx \left({\dot{M}}_{\star }-{\dot{M}}_{\text{in}}\frac{R+\sqrt{R{R}_0}+R_0}{2R_{\text{C}}}\right)\left(1-\sqrt{\frac{R_0}{R}}\right)\approx \dot{M}\left(1-\sqrt{\frac{R_0}{R}}\right)\left(1-\frac{{\dot{M}}_{\text{in}}}{{\dot{M}}_{\star }}\frac{\sqrt{R{R}_0}+R_0}{2R_{\text{C}}}\right),$$(A.5)

which isthus a modest correction to Eq. (2).

For R ≫ R₀ and Ṁ_⋆ = Ṁ_in, I have $${\dot{M}}_{\text{in}}[\sqrt{1-\frac{R}{R_{\text{C}}}}(1-\frac{{\left(R/R_{\text{C}}\right)}^{1/2}}{3})+\frac{2}{3}{\left(\frac{R_{\text{C}}}{R}\right)}^{1/2}\left(1-\sqrt{1-\frac{R}{R_{\text{C}}}}\right)]=3\pi \text{Σ}\nu .$$(A.6)

The factor between brackets decreases from 1 to 2/3 when R increases from 0 to R_C. Beyond the centrifugal radius, Ṁ = 0 and, from Eq. (A.3), R^1∕2Σν becomes constant so that $$3\pi \text{Σ}\nu =\frac{2}{3}{\dot{M}}_{\text{in}}{\left(\frac{R_{\text{C}}}{R}\right)}^{1/2}.$$(A.7)

For bounded αT, 2πRΣ is not integrable, but this steady-state solution, which merely sets an upper bound to the disk viscous expansion, can only hold for distances of t_vis (R) shorter than the timescale of evolution of the system⁸.

Appendix B Calculation of K_p

I take the attenuation column for protons as the minimum between the continuous slowing-down approximation (CSDA) and diffusion regimes of Padovani et al. (2018), after integrating their Eqs. (E.3) and (21)⁹, respectively: $$\begin{array}{lll}\hfill & \hfill & {\text{Σ}}_{p,H}(E)=m\text{min}\left[\frac{E^{1+s}}{A\left(p-1\right)},\frac{1}{2(p-1)}\sqrt{\frac{E\left({\alpha }_p-1\right)}{3{\sigma }_{MT}{L}_p(E)}}\frac{\Gamma \left(\frac{p-1}{{\alpha }_p-1}+1\right)}{\Gamma \left(\frac{p-1}{{\alpha }_p-1}+\frac{1}{2}\right)}\right]=\frac{1}{p-1}\text{min}[\frac{\Gamma \left(\frac{p-1}{{\alpha }_p-1}+1\right)}{\Gamma \left(\frac{p-1}{{\alpha }_p-1}+\frac{1}{2}\right)}1.2\text{\hspace{0.22em}kg\hspace{0.17em}}{\text{m}}^{-2}{\left(\frac{E}{10\text{\hspace{0.22em}MeV}}\right)}^{1.82},\hfill \\ \hfill & \hfill & 6\times {10}^2\text{\hspace{0.22em}kg\hspace{0.17em}}{\text{m}}^{-2}{\left(\frac{1\text{\hspace{0.22em}GeV}}{E}\right)}^{0.14}\frac{\Gamma \left(\frac{p-1}{{\alpha }_p-1}+1\right)}{\Gamma \left(\frac{p-1}{{\alpha }_p-1}+\frac{1}{2}\right)}],\hfill \end{array}$$(B.2)

where A = 1.77 × 10⁻⁶ eV^1+sm², s = 0.82, α_p = 1.28 (noted α in Padovani et al. 2018), σ_MT = 10⁻³⁰ m², L_p (E) = 10⁻²¹ m² eV ${(E/\text{GeV})}^{{\alpha }_p}$ are defined in Padovani et al. (2018) and Γ is Euler’s gamma function. For the other species, I adopt ${\text{Σ}}_{p,i}={\text{Σ}}_{p,H}{A}_i/Z_i^2$, with Z_i, A_i its atomic and mass numbers of species i (e.g., Gounelle et al. 2001). The transition between the CSDA and the diffusion regime thus occurs at $$E_{\text{crit}}=0.33\text{\hspace{0.22em}GeV}{\left(\frac{\Gamma \left(\frac{p-1}{{\alpha }_p-1}+1\right)}{\Gamma \left(\frac{p-1}{{\alpha }_p-1}+\frac{1}{2}\right)}\right)}^{0.51}.$$(B.3)

As to the cross sections, similar to Desch et al. (2004), I adopt the Yanasak et al. (2001) formulation $${\sigma }_{k,i}(E)={\sigma }_0\text{min}\left({\left(\frac{E}{E_0}\right)}^x,1\right)\theta (E-E\text{th}),$$(B.4)

with σ₀, Eth, E₀, x constants (for a given (k, i)) and θ is again the Heaviside function. I fit those for i = ⁴He using the Lange et al. (1995) data (unknown to Yanasak et al. 2001), with power-law fitting below 30 MeV nucleon⁻¹ and the flat regime being assigned the maximum cross section measured (rather than a fit beyond the aforementioned threshold, as this regime is only incipient in the energy range studied). The resulting parameters are presented in Table B.1.

Then I have for Eth ≤ E₀ ≤ E_crit $${\displaystyle {\int }_0^{+\infty }{\sigma }_{k,i}}(E){\text{Σ}}_{p,i}{E}^{-p}\text{d}E={\sigma }_0\frac{{\text{Σ}}_{p,i}^{\text{CSDA}}(E_0)}{E_0^{p-1}}(\frac{1-{(E_{\text{th}}/E_0)}^{2+s+x-p}}{2+s+x-p}+\frac{{(E_{\text{crit}}/E_0)}^{2+s-p}-1}{2+s-p}+\frac{{(E_{\text{crit}}/E_0)}^{2+s-p}}{p+({\alpha }_p-3)/2}),\text{}$$(B.5)

and for E_th ≤ E_crit ≤ E₀ $$\begin{array}{lll}{\displaystyle {\int }_0^{+\infty }{\sigma }_{k,i}}(E){\text{Σ}}_{p,i}{E}^{-p}\text{d}E\hfill & =\hfill & {\sigma }_0\frac{{\text{Σ}}_{p,i}^{\text{CSDA}}(E_0)}{E_0^{p-1}}(\frac{{(E_{\text{crit}}/E_0)}^{2+s+x-p}-{(E\text{th}/E_0)}^{2+s+x-p}}{2+s+x-p}+\frac{{(E_{\text{crit}}/E_0)}^{s+({\alpha }_p+1)/2}-{(E_{\text{crit}}/E_0)}^{2+s+x-p}}{x-p-({\alpha }_p-3)/2}\hfill \\ \hfill & \hfill & +\frac{{(E_{\text{crit}}/E_0)}^{s+({\alpha }_p+1)/2}}{p+({\alpha }_p-3)/2}),\hfill \end{array}$$(B.6)

whose contributions can be summed to yield K_p (Eq. (10)).

Appendix C Height ofspallation layer and effect of settling on ¹⁰Be production in solids

The spallation layer column Σ_spall may be defined as $${\text{Σ}}_{\text{spall}}\equiv \text{ frac}1\text{{left(overset}•\text{{^}10\text{Be}}{/}^9\text{Beright)\_tg},\text{+}(0)\text{}}\int \text{\_}0\text{^}+\infty \text{left(frac{overset}•\text{{^}10\text{Be}}{^}9\text{Be}right)\_tg},\text{+d}{\text{Σ}}^{\prime }={\langle {\text{Σ}}_p(E)\rangle }_{\sigma E^{-p}}\frac{\text{sin\hspace{0.17em}}\varphi }{2}.$$(C.1)

Assuming a vertically isothermal disk (and thus a gaussian gas density $\rho (z)=\text{Σ}\text{exp}(-{(z/H)}^2/2)/(\sqrt{2\pi }H)$), the height H_spall of this layer obeys $${\text{Σ}}_{\text{spall}}=\frac{\text{Σ}}{2}\text{erfc}\left(\frac{H_{\text{spall}}}{H\sqrt{2}}\right).$$(C.2)

For Σ ~ 10^5±1 kg m⁻² and Σ_spall = 10^1±1 kg m⁻², I have χ ≡ H_spall∕H ≈ 3 ± 1, so that $$\frac{H_{\text{spall}}}{R}=0.16\left(\frac{\chi }{3}\right){\left(\frac{R}{0.5\text{\hspace{0.22em}AU}}\right)}^{1/2}{\left(\frac{T}{1500\text{\hspace{0.22em}K}}\right)}^{1/2}.$$(C.3)

Since, for x ≳ 2, I have¹⁰ $$\text{erfc}(x)\approx \frac{\text{exp}(-x^2)}{\sqrt{\pi }x},$$(C.4)

the density at z = H_spall may be approximated as $$\rho (H_{\text{spall}})=\frac{\chi {\text{Σ}}_{\text{spall}}}{H}=8\times {10}^{-9}\text{\hspace{0.22em}kg\hspace{0.17em}}{\text{m}}^{\text{−3}}\left(\frac{\chi }{3}\right)\left(\frac{{\text{Σ}}_{\text{spall}}}{10\text{\hspace{0.22em}kg\hspace{0.17em}}{\text{m}}^{-2}}\right){\left(\frac{1500\text{\hspace{0.22em}K}}{T}\right)}^{1/2}{\left(\frac{0.5\text{\hspace{0.22em}AU}}{R}\right)}^{3/2}.$$(C.5)

The nominal value (corresponding to a number density of 2 × 10¹² cm⁻³) is halfway in the range studied by Umebayashi & Nakano (1981) so the neutron contribution (ignored in this paper) can be seen in their Figs. 2 and 5 to be limited. Then, the settling parameter (Jacquet et al. 2012) for a solid of density ρ_s and radius a there is $$\begin{array}{lll}{S}_z(H_{\text{spall}})\hfill & =\hfill & \sqrt{\frac{\pi }{8}}\frac{{\rho }_{\text{s}}a}{\rho (H_{\text{spall}})c_{\text{s}}{\delta }_z}\approx \sqrt{\frac{\pi }{8}}\frac{{\rho }_{\text{s}}a}{{\text{Σ}}_{\text{spall}}{\delta }_z\chi }=0.2\left(\frac{{\rho }_{\text{s}}a}{1\text{\hspace{0.22em}kg\hspace{0.17em}}{\text{m}}^{\text{2}}}\right)\left(\frac{10\text{\hspace{0.22em}kg\hspace{0.17em}}{\text{m}}^{-2}}{{\text{Σ}}_{\text{spall}}}\right)\left(\frac{{10}^{-1}}{{\delta }_z}\right)\left(\frac{3}{\chi }\right),\hfill \end{array}$$(C.6)

with δ_z the vertical diffusion coefficient normalized to $c_{\text{s}}^2/{\Omega }_K$ (which should be of order α).

Simulations of ideal MHD turbulence indicate δ_z ~ α ∝ ρ⁻¹ over a large extent of the disk thickness, although it should stall in the corona (e.g., Jacquet 2013, and references therein). Since S_z may be expected to be ≲1 there from the above evaluation, it may not be a dramatic underestimate to treat it as constant throughout the vertical thickness of the disk. Under that assumption, a population of identical solids has a density (Jacquet 2013) $${\rho }_{\text{c}}(z)=\frac{{\text{Σ}}_{\text{c}}\sqrt{1+S_z}}{\sqrt{2\pi }H}\text{exp}\left(-\frac{S_z+1}{2}{\left(\frac{z}{H}\right)}^2\right),$$(C.7)

with Σ_c the surface density of the population. The ¹⁰Be production rate averaged over this distribution is then $$\text{nn leftlangleleft(frac{overset}•\text{{^}10\text{Be}}{^}9\text{Be}right)\_tgrightrangle\_}{\rho }_{\text{c}}\equiv \frac{2}{{\text{Σ}}_{\text{c}}}{\displaystyle {\int }_0^{+\infty }\text{d}}{\text{Σ}}_{{\text{c}}^{\prime }}\text{left(frac{overset}•\text{{^}10\text{Be}}{^}9\text{Be}right)\_tg},\text{+}\left({\text{Σ}}^{\prime }({\text{Σ}}_{p^{\prime }})\right),$$(C.8)

with ${\text{Σ}}_{{\text{c}}^{\prime }}$ and Σ′ the columns of the solids and gas, respectively, integrated vertically from the surface. If I let $$\text{nn left(frac{overset}•\text{{^}10\text{Be}}{^}9\text{Be}right)\_tg},\text{+}\equiv \text{ left(frac{overset}•\text{{^}10\text{Be}}{^}9\text{Be}right)\_tg},\text{+}(0)g\left(\frac{{\text{Σ}}^{\prime }}{{\text{Σ}}_{\text{spall}}}\right),$$(C.9)

and since $${\text{Σ}}^{\prime }=\frac{\text{Σ}}{2}\text{erfc}\left(\frac{{\text{erfc}}^{-1}\left(2{\text{Σ}}_{{\text{c}}^{\prime }}/{\text{Σ}}_{\text{c}}\right)}{\sqrt{1+S_z}}\right),$$(C.10)

this becomes $$\text{begin{align}nnleftlangleleft(frac{overset}•\text{{^}10\text{Be}}{^}9\text{Be}right)\_tgrightrangle\_}{\rho }_{\text{c}}\text{}=\text{}\frac{\text{Σ}}{2{\text{Σ}}_{\text{spall}}}\text{}{\displaystyle {\int }_0^{+\infty }\text{}}g\left(\text{}\frac{\text{Σ}}{2{\text{Σ}}_{\text{spall}}}\text{erfc}\left(\frac{{\text{erfc}}^{-1}(x)}{\sqrt{1+S_z}}\right)\right)\text{d}x\text{ leftlangleleft(frac{overset}•\text{{^}10\text{Be}}{^}9\text{Be}right)\_tgrightrangle\_}\rho \text{}\approx \text{}\frac{\text{Σ}}{2{\text{Σ}}_{\text{spall}}}\text{}{\displaystyle {\int }_0^{+\infty }\text{}}\text{}g\left(\frac{\text{Σ}}{2{\text{Σ}}_{\text{spall}}}\frac{x^{1/(S_z+1)}\sqrt{S_z+1}}{{\left(\sqrt{\pi }{\text{erfc}}^{-1}(x)\right)}^{S_z/(S_z+1)}}\right)\text{d}x\text{ leftlangleleft(frac{overset}•\text{{^}10\text{Be}}{^}9\text{Be}right)\_tgrightrangle\_}\rho \text{},\text{nn end{align}}$$(C.11)

where I have used approximation (C.4), with the corollary ${\text{erfc}}^{-1}(z)\approx \sqrt{-\text{ln}(\sqrt{\pi }z{\text{erfc}}^{-1}(z))}$ for z ≪ 1. Roughly speaking, the integrand only contributes (~1) for a g argument ≲ 1, that is $x\lesssim {(\sqrt{\pi }{\text{erfc}}^{-1}(x))}^{S_z}{(2{\text{Σ}}_{\text{spall}}/\text{Σ}\sqrt{S_z+1})}^{S_z+1}$ so that $$\text{nnleftlangleleft(frac{overset}•\text{{^}10\text{Be}}{^}9\text{Be}right)\_tgrightrangle\_}{\rho }_{\text{c}}~\frac{1}{\sqrt{S_z+1}}{\left(\sqrt{\pi \text{ln}(\text{Σ}/{\text{Σ}}_{\text{spall}})}\frac{2{\text{Σ}}_{\text{spall}}}{\text{Σ}}\right)}^{S_z}\text{leftlangleleft(frac{overset}•\text{{^}10\text{Be}}{^}9\text{Be}right)\_tgrightrangle\_}\rho .$$(C.12)

So ¹⁰Be production should quickly dropas S_z approaches unity and become negligible for S_z ≳ 1.

References

All Tables

All Figures

Fig. 1 Sketch of the scenario of 10Be production developed in this paper. Solar cosmic rays cause spallation reactions in the upper layers of the disk, in the gas and/or condensates, within a region (outlined by a dashed separatrix) where magnetic field lines originating from the proto-Sun and its outskirts hit the disk. The field and particle emission geometries are purely illustrative and do not detract from the generality of the formalism.

In the text

	Fig. 2 Kp as a function of the power-law exponent of the incoming cosmic ray energy distribution, for a solar (without possible implantation contributions) and a rocky target (see Sects. 3.1 and 3.2, respectively).
In the text

Fig. 3 feff (Eq. (17)) as a function of heliocentric distances for two values of Rmax for my toy parameterization of the particle flux distribution (see end of Sect. 2.2) and two values of the radial Schmidt number, with the lower one (higher relative diffusivity) in dotted lines. The parameter feff is proportional to the 10Be/9Be ratio (whose absolute magnitude is discussed at Eq. (22)) in the bulk (gas+solids) disk for given mass accretion rate and proton output and energy distributions.

In the text

	Fig. 4 f(R, +∞) as a functionof heliocentric distances for two values of Rmax for my toy parameterization of the particle flux distribution (see end of Sect. 2.2). This is proportional to the 10Be/9Be ratio produced by in situ irradiation of a refractory solid as a function of its heliocentric distance of origin.
In the text

Fig. 5 LX∕Ṁ for Taurusas a function of age (data from Güdel et al. 2007). Although a least-square fit (dashed line) shows a systematic increase with age, scatter is evident, in part owing to the diversity of stellar masses. Overplotted is the curve (red, continuous) corresponding to one solar mass (Telleschi et al. 2007; see beginning of Sect. 4.2). For reference, the normalizations used in Eq. (22) amount to LX ∕Ṁ = 107.7 J kg−1.

In the text

	Fig. 6 LX∕Ṁ for young stellar objects as a function of Ṁ. Taurus and Orion data are taken from Güdel et al. (2007) and Bustamante et al. (2016), respectively (where LX ∕Ṁ also shows little systematic trend with stellar mass).
In the text