text:
# Testing two nuclear physics approximations used in the standard leaky box model for the spallogenic production of LiBeB

## 1 Introduction

The production of the isotopes of Lithium, Beryllium and Boron (LiBeB) from the spallation of $`{}_{}{}^{12}\mathrm{C}`$, $`{}_{}{}^{14}\mathrm{N}`$ and $`{}_{}{}^{16}\mathrm{O}`$ in the Inter-Stellar Medium (ISM) by protons and alpha particles in the Cosmic Rays (CRs) outlined in Meneguzzi, Audouze & Reeves (1971) filled the last major gap in the understanding of the origin the elements. This method of manufacture (dubbed Galactic Cosmic Ray Nucleosynthesis (GCRN) or Cosmic Ray Spallation) is thought to be the sole source of $`{}_{}{}^{6}\mathrm{Li}`$, $`{}_{}{}^{9}\mathrm{Be}`$ and $`{}_{}{}^{10}\mathrm{B}`$ and thus observations of their abundances give clues to the contribution of spallogenic $`{}_{}{}^{7}\mathrm{Li}`$ and $`{}_{}{}^{11}\mathrm{B}`$ to the total abundance of these nuclei. The fragility and paucity of these elements renders the derivation of their abundances from stellar spectra challenging but despite this difficulty data now exists over a wide range of metallicity (Molaro *et al.*, 1997; Duncan, Lambert & Lemke, 1992; Duncan *et al.*, 1997; García López et al., 1998; Hobbs, Thorburn & Rebull, 1999; Nissen *et al.*, 1999; Smith, Lambert & Nissen, 1993). After comparing the observations with the predictions of theory many conclusions have been drawn: for example the discrepancy between the observed solar $`{}_{}{}^{11}\mathrm{B}`$/$`{}_{}{}^{10}\mathrm{B}`$ isotopic ratio of 4.05 (Anders & Grevesse, 1989) with the spallogenic prediction of $`2.5`$ (Reeves, 1974) indicates another source of $`{}_{}{}^{11}\mathrm{B}`$ must exist.

The strength of this, and indeed any, conclusion ultimately relies on the validity of the underlying approximations that permeate the GCRN calculation of which there are many. Within the sources of error there are two that emerge from the nuclear physics aspects of the problem. The first is the straight-ahead approximation which assumes that the fragment emerges from a reaction with a velocity equal to that of the projectile, while the second is that the only significant source of LiBeB is the direct, one-step, production from the spalling of the CNO and that production via intermediaries is insignificant. The use of both of these approximations is widespread but remains largely untested. Relaxing either considerably increases the complexity of the calculation yet this is exactly our intention in this paper. In section §2 we present the basics of cosmic ray calculations in the leaky box model, show how the reaction expansion is invoked and how the straight-ahead approximation enters. Then in section §3 we relax the straight-ahead approximation and demonstrate the implications for the production rates and flux of the secondary nuclei. In section §4 we include the two-step processes in the reaction rates and then, motivated by these results, we proceed to section §5 where we use a cascade technique to calculate the production rates exactly.

## 2 CR Basics and the Straight-Ahead Approximation

The theoretical underpinning for the study of cosmic rays is a propagation model that relates a source spectrum to the flux observed at a later time and place. Of the theoretical models used the most popular for the calculation of the rates of production of $`{}_{}{}^{6}\mathrm{Li}`$, $`{}_{}{}^{7}\mathrm{Li}`$, $`{}_{}{}^{9}\mathrm{Be}`$, $`{}_{}{}^{10}\mathrm{B}`$ and $`{}_{}{}^{11}\mathrm{B}`$ is the Leaky-Box Model (see e.g. Cesarsky (1980)) or its variants. The model is simpler and more transparent than the rival diffusion model but it suffers from number of deficiencies, both theoretical and when used to fit some cosmic ray data. Despite this failing the use of the model is widespread because it has empirically proven to be useful in fitting compositional cosmic ray data which is the focus of this paper.

In the leaky box model the flux of a cosmic ray species is given by (Fields, Olive and Schramm, 1994; Ramaty *et al.*, 1997)

$$\varphi (E)=\frac{1}{w(E)\rho }_E^{\mathrm{}}𝑑E^{}Q(E^{})\frac{S(E^{})}{S(E)}$$
(1)
where $`Q`$ represents the source of particles, the mass density of the interstellar medium (ISM) is represented by $`\rho `$ and the integrating factor $`S`$ has a physical interpretation as the probability that a particle produced with energy $`E`$ will survive to form part of the ISM. This survival probability is simply

$$S(E)=\mathrm{exp}\left(_0^E\frac{dE^{}}{w(E^{})\mathrm{\Lambda }(E^{})}\right),$$
(2)
where $`w`$ is the energy per nucleon lost per unit distance traversed and per unit mass density of the medium and $`\mathrm{\Lambda }(E)`$ is the path length or grammage. The quantity $`\varphi w\rho `$ is frequently called the current, the rate of flow of particles from high to low energy, so that in the limit when $`S(E)1`$ all the particles at a given energy will eventually accumulate at $`E=0`$ i.e. they form part of the ISM. This limit is achieved when $`E0`$ so the rate of change of the abundance $`y`$ is therefore

$$\frac{dy}{dt}=\frac{1}{n_H}_0^{\mathrm{}}𝑑E^{}Q(E^{})S(E^{}),$$
(3)
where $`n_H`$ is the number density of hydrogen in the ISM.

So far we have not distinguished between the different nuclei under consideration since equations (1), (2) and (3) are valid for each regardless of how the different species become cosmic rays. There are differences of course, the flux of CNO is almost entirely that accelerated from the ISM whereas the flux of LiBeB is almost entirely produced ‘in-flight’ so to speak. We can thus make a distinction between two types of sources $`Q`$: a primary spectrum $`Q^{(0)}`$ from whatever mechanism is responsible for acceleration of particles from the ISM, and a secondary production term $`Q^{(s)}`$ that is due to particle reactions between primary CRs and the constituents of the ISM. Cosmic rays are classified by whichever is dominant. The production term $`Q^{(s)}`$ is simply the sum of all the reactions that can make the nucleus we are interested in, so we have

$`Q_F(E_F)`$ $`=`$ $`Q_F^{(0)}(E_F)+Q_F^{(s)}(E_F),`$ (4)
$`Q_F^{(s)}(E_F)`$ $`=`$ $`{\displaystyle \underset{P,T}{}}{\displaystyle _0^{\mathrm{}}}𝑑E_P\varphi _P(E_P)n_T\sigma _F^{P,T}(E_P)P_F(E_F|E_P)`$ (5)
where we have used the superscripts/subscript $`P`$, $`T`$ and $`F`$ to denote the quantities associated with the projectile, target or fragment of the reaction $`\mathrm{P}+\mathrm{T}\mathrm{F}+\mathrm{X}`$, $`n_T`$ is the number density of targets, $`\sigma _F^{P,T}`$ is the cross-section for the reaction and $`P_F(E_F|E_P)`$ is the probability distribution for the fragment energy $`E_F`$ at a given projectile energy $`E_P`$. At this point we expanded the projectile flux $`\varphi _P`$ as $`\varphi _P=_{i=0}\varphi _P^{(i)}`$, the reasoning will become clear in a moment, and insert this expression into equation (5) so that the secondary source for the fragment $`Q_F^{(s)}`$ becomes

$`Q_F^{(s)}(E_F)`$ $`=`$ $`{\displaystyle \underset{i=0}{}}{\displaystyle \underset{P,T}{}}{\displaystyle _0^{\mathrm{}}}𝑑E_P\varphi _P^{(i)}(E_P)n_T\sigma _F^{P,T}(E_P)P(E_F|E_P)={\displaystyle \underset{i=1}{}}Q_F^{(i)},`$ (6)
$`Q_F^{(i)}(E_F)`$ $`=`$ $`{\displaystyle \underset{P,T}{}}{\displaystyle _0^{\mathrm{}}}𝑑E_P\varphi _P^{(i1)}(E_P)n_T\sigma _F^{P,T}(E_P)P(E_F|E_P),`$ (7)
which introduces the expansion over $`Q_F^{(s)}`$. Note that this expansion begins at $`i=1`$ not $`i=0`$ because $`Q_F^{(0)}`$ is already in use as the primary source spectrum for nucleus $`F`$. We can now bundle together the primary and secondary sources we introduced in equation (4) into one expression: $`Q_F=_{i=0}Q_F^{(i)}`$. Now we can also expand the fragment’s flux $`\varphi _F`$ in exactly the same fashion as the projectile flux $`\varphi _P`$ so that equation (1), now reads

$$\varphi _F(E_F)=\underset{i=0}{}\varphi _F^{(i)}(E_F)=\underset{i=0}{}\frac{1}{w_F(E_F)\rho }_{E_F}^{\mathrm{}}𝑑E_F^{}Q_F^{(i)}(E_F^{})\frac{S_F(E_F^{})}{S_F(E_F)}.$$
(8)
If we equate terms in $`i`$ we achieve our final result

$$\varphi _F^{(i)}(E_F)=\frac{1}{w_F(E_F)\rho }_{E_F}^{\mathrm{}}𝑑E_F^{}Q_F^{(i)}(E_F^{}).\frac{S_F(E_F^{})}{S_F(E_F)}$$
(9)
The meaning of the expansion now becomes clear, we have expanded the fragment flux $`\varphi _F`$ over the number of reactions that lead from a given projectile $`P`$ to the desired fragment $`F`$ because $`Q_F^{(i)}`$, defined in equation (7), depends upon the fluxes of projectiles $`\varphi _P^{(i1)}`$. Thus $`\varphi _F^{(0)}`$ is the flux that is accelerated from the ISM, $`\varphi _F^{(1)}`$ denotes the flux produced by interactions involving $`\varphi _P^{(0)}`$, $`\varphi _F^{(2)}`$ is that which arises from reactions involving $`\varphi _P^{(1)}`$ and so on.

Inserting the same expansion of $`Q_F`$ into equation (3) leads to a similar expansion of the rate of change of abundance $`dy_F/dt`$ as a sum of terms $`dy_F^{(i)}/dt`$ where each contribution is from the source term $`Q_F^{(i)}`$

$$\frac{dy_F}{dt}=\underset{i=0}{}\frac{dy_F^{(i)}}{dt}=\underset{i=0}{}\frac{1}{n_H}_0^{\mathrm{}}𝑑E_F^{}Q_F^{(i)}(E_F^{})S_F(E_F^{}).$$
(10)
We return our attention to equation (7) and the first term in the reaction source expansion $`Q_F^{(1)}`$. This is the term that corresponds to the formation of $`F`$ via a single interaction of a primary CR and the ISM, i.e. the one-step source. We can break this quantity apart into the sum from each projectile/target pair and likewise expand $`\varphi _F^{(1)}`$ and $`dy_F^{(1)}/dt`$ in the same fashion with each term representing the contribution from the source term $`Q_F^{P,T}`$. Lastly, the primary projectile flux $`\varphi _P^{(0)}`$ may be rewritten as $`\varphi _P^{(0)}=\alpha _P\psi _P`$ so that the integral over the projectile’s energy per nucleon $`E_P`$ yields $`\alpha _P`$, the total number of projectiles per unit time and per unit area. The parameter $`\alpha _P`$ acts as the flux strength and we can therefore remove the trivial dependence upon the chemistry of the CRs. The rate of change of the abundance in the ISM and the flux of nucleus $`F`$ are thus

$`{\displaystyle \frac{dy_F^{(1)}}{dt}}`$ $`=`$ $`{\displaystyle \underset{P,T}{}}{\displaystyle \frac{dy_F^{P,T}}{dt}}={\displaystyle \underset{P,T}{}}\alpha _Py_TR_F^{P,T}`$ (11)
$`\varphi _F^{(1)}(E_F)`$ $`=`$ $`{\displaystyle \underset{P,T}{}}\varphi _F^{P,T}(E_F)={\displaystyle \underset{P,T}{}}{\displaystyle \frac{\alpha _Py_T}{\rho /n_H}}\mathrm{\Psi }_F^{P,T}(E_F)`$ (12)
where

$`R_F^{P,T}`$ $`=`$ $`{\displaystyle _0^{\mathrm{}}}𝑑E_F^{}\left\{{\displaystyle _0^{\mathrm{}}}𝑑E_P\psi _P(E_P)\sigma _F^{P,T}(E_P)P_F(E_F^{}|E_P)\right\}S_F(E_F^{})`$ (13)
$`\mathrm{\Psi }_F^{P,T}(E_F)`$ $`=`$ $`{\displaystyle _{E_F}^{\mathrm{}}}𝑑E_F^{}\left\{{\displaystyle _0^{\mathrm{}}}𝑑E_P\psi _P(E_P)\sigma _F^{P,T}(E_P)P_F(E_F^{}|E_P)\right\}{\displaystyle \frac{S_F(E_F^{})}{w_F(E_F)S_F(E_F)}}.`$ (14)
Equations (13) and (14) serve to define the normalized, one-step reaction rate $`R_F^{P,T}`$ and the normalized, one-step flux $`\mathrm{\Psi }_F^{P,T}`$.

With the basics complete we can fill in some of the mundane details of the calculation. The set of nuclei we will use for both the ISM and CRs are p/$`{}_{}{}^{1}\mathrm{H}`$ and $`\alpha `$/$`{}_{}{}^{4}\mathrm{He}`$ plus all the stable nuclei from $`{}_{}{}^{6}\mathrm{Li}`$ through to $`{}_{}{}^{16}\mathrm{O}`$. In addition we must add $`{}_{}{}^{7}\mathrm{Be}`$ and $`{}_{}{}^{11}\mathrm{C}`$ to our set of nuclei found in the CRs because their decay mode is inner orbital electron capture and so are stable in the ionized state. We will sometimes denote by $`Z`$ the subset of nuclei with mass greater than $`{}_{}{}^{4}\mathrm{He}`$. The stopping power $`w_F(E)`$ is tabulated up to $`12\mathrm{MeV}/\mathrm{nucleon}`$ for each fragment by Northcliffe & Shilling (1970). Above this energy we use the stopping power for protons as tabulated by Janni (1982) and rescale in exactly the same manner as Fields, Olive and Schramm (1994). The important cross-sections needed for our calculation were tabulated by Read & Viola (1984) but there is no cross-sectional data involving $`{}_{}{}^{7}\mathrm{Be}`$ and $`{}_{}{}^{11}\mathrm{C}`$ as a reactant so we have appealed to the fact that each is the mirror of a nucleus for which data does exist and therefore assume that the cross-section for reactions involving $`{}_{}{}^{7}\mathrm{Be}`$ or $`{}_{}{}^{11}\mathrm{C}`$ is the same as for it’s isobaric brother. For the primary CR spectrum we adopt the form presented by Gloeckler & Jokipii (1967) and used by Walker, Mathews & Viola (1985) and Steigman & Walker (1992)

$$\varphi _P^{(0)}(E_P)=\alpha _P\frac{1.6E_0^{1.6}}{(E_P+E_0)^{2.6}}(\mathrm{MeV}/\mathrm{nucleon})^1\mathrm{cm}^2\mathrm{s}^1.$$
(15)
We have set the parameter $`E_0`$ to be the nucleon mass $`m_N`$ though we note Garcia-Munoz *et al.* (1987) preferred a value of $`400\mathrm{MeV}/\mathrm{nucleon}`$.

The last quantity we need to specify is the fragment energy distribution $`P_F(E_F|E_P)`$. This quantity has been historically approximated with a $`\delta `$-function form known as the straight-ahead approximation. The approximation is divided into two, the selection being dependent upon the identity of the projectile $`P`$,

$$P_F(E_F|E_P)=\left\{\begin{array}{cc}\delta (E_F)\hfill & \mathrm{P}\{\mathrm{p},\alpha \}\hfill \\ \delta (E_PE_F)\hfill & \mathrm{P}\{\mathrm{Z}\}\hfill \end{array}\right\},$$
(16)
Note that in the case of $`\mathrm{P}\{\mathrm{p},\alpha \}`$ there is *no* flux of the fragment F.

We compile the one-step production rates $`R_F^{P,T}`$ as calculated by equation (13) using two, constant values of $`\mathrm{\Lambda }=5\mathrm{g}/\mathrm{cm}^2`$ and $`\mathrm{\Lambda }=10\mathrm{g}/\mathrm{cm}^2`$. In table (1) only those where the reactants are p/$`{}_{}{}^{1}\mathrm{H}`$ and $`{}_{}{}^{16}\mathrm{O}`$. An initial examination of the results shows that loss of the fragments when $`\mathrm{P}\{\mathrm{Z}\}`$ reduces the production rate by a factor $`10`$ relative to the rate when the projectile and target are interchanged. While this may appear to render the inverse reactions unimportant it must be remembered that, at the present time, the abundances of many of the $`Z`$ nuclei in the CR are enhanced relative to their ISM values by approximately equivalent factors. Finally, we also see from the table that increasing the grammage leads to a larger fraction of LiBeB captured from the inverse reactions and so increasing the efficiency of production.

## 3 Relaxing The Straight Ahead Approximation

In section §2 we used the straight-ahead approximation to calculate $`R_F^{P,T}`$ and $`\mathrm{\Psi }_F^{P,T}`$ and we now wish to test its validity. An examination of the approximation was made by Fields, Olive and Schramm (1994) who replaced it with another $`\delta `$-function but with different relation between the fragment and primary energies while Tsao *et al.* (1995) conducted a test for heavy nuclei that included a spread in fragment energy. It is the direction of Tsao *et al.* (1995) that we shall explore.

The momentum distributions of the isotopes produced by fragmentation of various projectile and target nuclei are discussed by Morrissey (1989) and Hufner (1985). The distributions from experiment are found to be normally distributed in the rest frame of the spalled nucleus $`Z`$ with narrow dispersions and a small mean momentum in the direction parallel to the incident projectile. In particular the momentum distributions of the fragments from $`{}_{}{}^{12}\mathrm{C}`$ and $`{}_{}{}^{16}\mathrm{O}`$ projectiles upon various targets ranging in mass from Be to Pb were measured by Greiner *et al.* (1975) who found that their results had no significant correlation with target mass or beam energy. Goldhaber (1974), and more recently Bauer (1989), explain these results in terms of a nuclear model with minimal correlations between the nucleon momenta. They predict the dependence of the momentum per nucleon distribution $`P_F(𝐏_F^{(Z)}|E_P)`$ of a fragment with mass $`A_F`$ to be

$$P_F(𝐏_F^{(Z)}|E_P)=\frac{1}{(2\pi ϵ_F^2)^{3/2}}\mathrm{exp}\left(\frac{(𝐏_F^{(Z)}\overline{𝐏}_F)^2}{2ϵ_F^2}\right)$$
(17)
where $`𝐏_F^{(Z)}`$ is the momentum per nucleon of the fragment in the rest-frame of $`Z`$, $`\overline{𝐏}_F`$ is the mean and $`ϵ_F^2`$ the variance of the distribution. To obtain $`P_F(E_F|E_P)`$ when $`P\{p,\alpha \}`$ we must integrate equation (17) over all momenta $`𝐏_F^{(Z)}`$ that yield $`E_F`$ and make use of the fact that $`𝐏_F^{(Z)}=𝐏_F`$.

$$P_F(E_F|E_P)=d^3𝐏_F^{}P_F(𝐏_F^{}|E_P)\delta (E_FE_F^{})$$
(18)
where $`E_F^{}`$ is the energy associated with $`𝐏_F^{}`$, i.e. $`E_F^{}+m_N=\sqrt{(𝐏_F^{})^2+m_N^2}`$. After a little effort this produces

$$P_F(E_F|E_P)=\frac{1}{\sqrt{2\pi ϵ_F^2}}\left(\frac{E_F+m_N}{\overline{P}_F}\right)\left[\mathrm{exp}\left(\frac{(P_F\overline{P}_F)^2}{2ϵ_F^2}\right)\mathrm{exp}\left(\frac{(P_F+\overline{P}_F)^2}{2ϵ_F^2}\right)\right].$$
(19)
Likewise, when $`\mathrm{P}\{\mathrm{Z}\}`$ we can obtain $`P_F(E_F|E_P)`$ but the relation between $`𝐏_F^{(Z)}`$ and $`E_F`$ is now through a Lorentz transformation and consequently yields the more complicated

$`P_F(E_F|E_P)`$ $`=`$ $`{\displaystyle \frac{1}{P_P\sqrt{2\pi a_F^2}}}\mathrm{exp}\left({\displaystyle \frac{\overline{\gamma _F}^2/\gamma _P^22\gamma _F\overline{\gamma _F}/\gamma _P+1}{2a_F^2}}\right)`$ (20)
$`\left\{\mathrm{exp}(x_{}^2)\mathrm{exp}(x_+^2)+\sqrt{{\displaystyle \frac{\pi }{2}}}{\displaystyle \frac{\overline{\gamma _F}}{a_F}}\left[\mathrm{erf}(x_+)\mathrm{erf}(x_{})\right]\right\}`$
where $`\overline{\gamma _F}=\overline{P_F}/\beta _Pm_N`$, $`a_F=ϵ_F/m_N`$, and $`x_+`$ and $`x_{}`$ are

$`x_+`$ $`=`$ $`{\displaystyle \frac{\gamma _P\gamma _F(1+\beta _P\beta _F)\overline{\gamma _F}}{\sqrt{2}a_F}}`$ (21)
$`x_{}`$ $`=`$ $`{\displaystyle \frac{\gamma _P\gamma _F(1\beta _P\beta _F)\overline{\gamma _F}}{\sqrt{2}a_F}}.`$ (22)
We have used the expressions for $`\overline{𝐏}_F`$ and $`ϵ_F^2`$ given by Morrissey (1989),

$`ϵ_F^2`$ $`=`$ $`ϵ_0^2{\displaystyle \frac{A_ZA_F}{A_F(A_Z1)}},`$ (23)
$`\overline{𝐏}_F`$ $`=`$ $`8{\displaystyle \frac{(A_ZA_F)}{A_Z}}{\displaystyle \frac{\gamma _P+1}{\beta _P\gamma _P}}\widehat{𝐏}_P\mathrm{MeV}/\mathrm{nucleon},`$ (24)
with $`ϵ_0100\mathrm{MeV}`$, $`A_Z`$ the mass number of the nucleus to be spalled (the heavier of $`P`$ and $`T`$) while $`\beta _P=v_P/c`$ and $`\gamma _P=1/\sqrt{1\beta _P^2}`$ are the Lorentz variables of the projectile. We have ignored the decreases in $`ϵ_F`$ that occur for values of E$`{}_{P}{}^{}100`$ MeV/nucleon indicated by Stokstad (1984): at such small energies, and with our choices of $`\mathrm{\Lambda }`$, the value of $`S_F`$ is very close to unity and every fragment isotope is captured by the Galaxy rendering the decrease in $`ϵ_F`$ irrelevant. The expression for $`\overline{𝐏}_F`$ in equation (24) diverges as $`\beta _P`$ becomes small so it must be replaced by another in this limit. The substitute we have used is

$$\overline{𝐏}_F=\gamma _P\beta _P\frac{A_L}{A_P+A_T}m_N\widehat{𝐏}_P,$$
(25)
where $`A_L`$ is the lighter mass of the interacting nuclei $`P`$ and $`T`$. This equation assumes the fragment carries a fraction $`A_F/(A_P+A_T)`$ of the projectile momentum in the frame where the heavier nucleus is at rest and the changeover from equation (24) to (25) again occurs at energies $`100\mathrm{MeV}/\mathrm{nucleon}`$. Cumming, Haustein & Hseuh (1981) state that such a change is expected to occur at roughly this energy.

The two probability distributions are shown in figures (1) and (2) for the interaction of p/$`\mathrm{H}`$ and $`{}_{}{}^{16}\mathrm{O}`$ to produce $`{}_{}{}^{9}\mathrm{Be}`$ when the Lorentz factor of the projectile is $`\gamma _P=2`$. The figures show that the emerging Beryllium energy is peaked below $`5\mathrm{MeV}/\mathrm{nucleon}`$ for the forward case and that the spread in Beryllium energy for the inverse reaction is considerable.

With the probability distributions now under control it is straight forward to calculate the rates of production and the flux by inserting the distributions into equations (13) and (14). We compile the rates again for $`\mathrm{\Lambda }=5\mathrm{g}/\mathrm{cm}^2`$ and $`10\mathrm{g}/\mathrm{cm}^2`$ which we show in table (2) again only for the reactions involving p/$`\mathrm{H}`$ with $`{}_{}{}^{16}\mathrm{O}`$. After comparing with the results in table (1) we can see that straight-ahead approximation accurately predicts the production rate for the forward reactions. This result was to be expected: the distribution of the fragment energy in localized to $`5\mathrm{MeV}/\mathrm{nucleon}`$ as shown in figure (1) and at these energies the range (the distance travelled by a particle before coming to rest per unit mass density of the medium) is much shorter than $`\mathrm{\Lambda }`$ so all fragments are trapped. After examining the inverse reactions we see a slight increase ($`8\%`$) in the production rates relative to the straight-ahead calculations but the differences are small and do not warrant rejection of the approximation. The differences are of a similar magnitude as those found by Tsao *et al.* (1995). The approximation’s success again has a simple explanation: the dispersion of the Beryllium energy $`ϵ_{\mathrm{Be}}`$ is only of order $`10\mathrm{MeV}`$ and so $`𝐏_{\mathrm{Be}}𝐏_\mathrm{O},E_{\mathrm{Be}}E_\mathrm{O}`$.

## 4 TwoStep Reactions

So far we have only considered the first term $`Q_F^{(1)}`$ in the reaction expansion from equation (7) but now we turn our attention to the second term, $`Q_F^{(2)}`$, to determine its contribution to the rate of change the abundance.

If the fragment $`F_1`$ from the reaction $`P+T_1F_1`$ undergoes a subsequent reaction $`F_1+T_2F_2`$ then from equation (7) we have

$$Q_{F_2}^{(2)}(E_{F_2})=\underset{F_1,T_2}{}_0^{\mathrm{}}𝑑E_{F_1}\varphi _{F_1}^{(1)}(E_{F_1})n_{T_2}\sigma _{F_2}^{F_1,T_2}(E_{F_1})P(E_{F_2}|E_{F_1})$$
(26)
Inserting the expansion for $`\varphi _{F_1}^{(1)}`$ and the expression for $`\varphi _F^{P,T}`$ from equation (12) then

$$Q_{F_2}^{(2)}(E_{F_2})=\underset{P,T_1,T_2}{}\frac{\alpha _Py_{T_1}y_{T_2}}{\rho /n_H^2}\underset{F_1}{}_0^{\mathrm{}}𝑑E_{F_1}\mathrm{\Psi }_{F_1}^{P,T_1}(E_{F_1})\sigma _{F_2}^{F_1,T_2}(E_{F_1})P_{F_2}(E_{F_2}|E_{F_1}).$$
(27)
where the flux $`\mathrm{\Psi }_{F_1}^{P,T_1}(E_{F_1})`$ appearing in equation (27) is the normalized one-step flux that was defined in equation (12). Introducing the quantity $`Q_{F_2}^{P,T_1,T_2}`$ as the contribution to equation (27) from each triplet $`P,T_1,T_2`$ and then inserting the result into equation (10) allows us to define a normalized two-step $`R_{F_2}^{P,T_1,T_2}`$ by

$$\frac{\alpha _Py_{T_1}y_{T_2}}{\rho /n_H}R_{F_2}^{P,T_1,T_2}=\underset{F_1}{}_0^{\mathrm{}}𝑑E_{F_2}Q_{F_2}^{P,T_1,T_2}(E_{F_2})S_{F_2}(E_{F_2})$$
(28)
so that

$$\frac{dy_{F_2}^{(2)}}{dt}=\underset{P,T_1,T_2}{}\frac{\alpha _Py_{T_1}y_{T_2}}{\rho /n_H}R_{F_2}^{P,T_1,T_2}.$$
(29)
We have calculated these two-step reaction rates $`R_{F_2}^{P,T_1,T_2}`$ again for $`\mathrm{\Lambda }=5\mathrm{g}/\mathrm{cm}^2`$ and $`10\mathrm{g}/\mathrm{cm}^2`$ using our improvement of the fragment energy distribution from section §3 and in table (3) we show the results for the reactions involving p/$`\mathrm{H}`$ and $`{}_{}{}^{16}\mathrm{O}`$. For the forward cases, i.e. when $`P\{p,\alpha \}`$, we find that the two-step rates vary in magnitude from $`10^{53}\mathrm{mg}\mathrm{cm}^2`$ for $`{}_{}{}^{6}\mathrm{Li}`$ to $`10^{67}\mathrm{mg}\mathrm{cm}^2`$ for $`{}_{}{}^{11}\mathrm{B}`$ whereas the inverse rates, i.e. $`P\{Z\}`$, are all around $`10^{48}\mathrm{mg}\mathrm{cm}^2`$. The forward rates are very sensitive to the threshold energies of the reactions $`F_1+T_2F_2`$ since the fluxes of $`F_1`$ peak at very low energies as indicated by figure (1). In both cases these rates certainly appear much smaller than those in tables (1) and (2) and it is tempting to regard the two-step rates are negligible. However this would be an erroneous conclusion because attention must be given to the pre-factors of equation (28). For an ISM of only $`\mathrm{H}`$ and $`\mathrm{He}`$ with a mass fraction of $`\mathrm{He}`$ of Y the $`\rho /n_H`$ term is

$$\frac{1}{\rho /n_H}=\frac{1Y}{m_H}.$$
(30)
Including this factor we see that the forward rates are at least five orders of magnitude smaller than the one-step rates in table (2). In contrast, when we apply the $`\rho /n_H`$ pre-factor to the inverse rates that we find they are only $`1/10`$ smaller! This same result was found by Ramaty *et al.* (1997) who also saw substantial changes to their calculations of the ratios of production rates when they included the two-step contributions.

This result may seem surprising but it is not entirely unexpected. The inefficient retardation of the CRs means that the flux of intermediaries $`F_1`$ accumulates over a long period of time. The flux of $`F_1`$ is only $`10^210^3`$ smaller than the primary $`\varphi _P^{(0)}`$. After we factor in the smaller average energy of the intermediaries (which will result in a higher fraction of trapped fragments) and the multiplicity of the two-step channels we obtain two-step production rates that are comparable with the one-step.

The significant two-step contribution to the rate of change of abundances means that we must evaluate the higher order terms in the expansion of $`Q_F^{(s)}`$. At the two-step level the number of reaction rates we must compute is considerable, if we find equally large contributions from the three-step reactions and greater then the computational burden becomes excessive. At this point we would have to abandon this approach to spallation calculations, so we must either show that the three-step contribution is negligible or we must find an alternative and more efficient approach to the calculation of production rates.

## 5 A Cascade Calculation

The three-step rates are, in principle, simple to compute since we proceed in the same fashion as the one and two-step equations we derived in sections §2 and §4. In practice there are so many pathways between any given primary projectile $`P_0`$ and the fragment $`F_3`$ under consideration that the procedure becomes increasingly laborious. Instead we adopt the approach of Moskalenko *et al.* (2002) and utilize the simple property that in any spallation reaction the fragment is always lighter than the heavier reactant.

We begin by splitting the flux $`\varphi _A`$ of any nucleus A into the sum of two components, the primary $`\varphi _A^{(0)}`$ accelerated from the ISM and then the secondary flux $`\varphi _A^{(s)}`$. The secondary flux $`\varphi _A^{(s)}`$ is dependent upon the fluxes of all the nuclei B heavier than A which again are sum of primary and secondary components. Thus we have $`\varphi _A=\varphi _A^{(0)}+\varphi _A^{(s)}`$ and $`\varphi _A^{(s)}=_{B>A}f(\varphi _B^{(0)}+\varphi _B^{(s)})`$ where $`f`$ represents the successive application of the linear equations (5) and (1). The secondary flux of the lightest member of $`B`$ is expressible as $`\varphi _B^{(s)}=_{C>B}f(\varphi _C^{(0)}+\varphi _C^{(s)})`$ and so it can be eliminated from the function for $`\varphi _A^{(s)}`$. We continue consecutive elimination of the lightest secondary flux in the expression for $`\varphi _A^{(s)}`$ until we reach the heaviest nucleus we are considering. The flux of this heaviest nucleus does not have a secondary flux by construction so we find that $`\varphi _A^{(s)}`$ is a function only of primary fluxes $`\varphi ^{(0)}`$. We therefore introduce the quantity $`\varphi _A^{P_0}`$ as the contribution to $`\varphi _A^{(s)}`$ from each primary flux which we label by $`P_0`$

$$\varphi _A=\varphi _A^{(0)}+\underset{P_0}{}\varphi _A^{P_0}.$$
(31)
Inserting this expression into equation (5) for $`\varphi _P`$ we obtain

$`Q_F^{(s)}(E_F)`$ $`=`$ $`{\displaystyle \underset{P_0}{}}{\displaystyle \underset{T}{}}{\displaystyle _0^{\mathrm{}}}𝑑E_{P_0}\varphi _{P_0}^{(0)}(E_{P_0})n_T\sigma _F^{P_0,T}(E_{P_0})P_F(E_F|E_{P_0})`$ (32)
$`+{\displaystyle \underset{P_0}{}}{\displaystyle \underset{F^{},T}{}}{\displaystyle _0^{\mathrm{}}}𝑑E_F^{}\varphi _F^{}^{P_0}(E_F^{})n_T\sigma _F^{F^{},T}(E_F^{})P_F(E_F|E_F^{}).`$
where we use the symbol $`F^{}`$ instead of $`P`$ to emphasize that the second term involves reactions between fragments and the ISM. The set of fragments $`F^{}`$ is a subset of the $`Z`$ nuclei so the ISM targets $`T`$ in the second term of equation (32) can only be $`\mathrm{H}`$ or $`\mathrm{He}`$ (we have been ignoring $`Z+Z^{}`$ reactions) and therefore the fragment $`F`$ must be lighter than $`F^{}`$. We introduce the quantities $`Q_F^{P_0}`$ as the contribution to this equation from each term of the expansion over $`P_0`$ and ‘normalize’ the primary flux as $`\varphi _{P_0}^{(0)}=\alpha _{P_0}\psi _{P_0}`$ in exactly the same way as section §2 in order to extract the dependence upon the composition of the CRs. After their introduction we reach the most important expression of this section:

$`Q_F^{P_0}(E_F)`$ $`=`$ $`\alpha _{P_0}{\displaystyle \underset{T}{}}{\displaystyle _0^{\mathrm{}}}𝑑E_{P_0}\psi _{P_0}(E_{P_0})n_T\sigma _F^{P_0,T}(E_{P_0})P_F(E_F|E_{P_0})`$ (33)
$`+{\displaystyle \underset{T}{}}{\displaystyle \underset{F^{}=F+1}{\overset{P_01}{}}}{\displaystyle _0^{\mathrm{}}}𝑑E_F^{}\varphi _F^{}^{P_0}(E_F^{})n_T\sigma _F^{F^{},T}(E_F^{})P_F(E_F|E_F^{}).`$
We have introduced the symbol $`F+1`$ as the lightest fragment $`F^{}`$ from which $`F`$ may be produced in the reaction with $`T`$ while the upper limit to the sum, $`P_01`$, is the heaviest nucleus that can be formed after an interaction of $`P_0`$ with any component of the ISM. We have denoted this nucleus suggestively by $`P_01`$ but we note that in the case of $`P_0\{p,\alpha \}`$ the heaviest intermediary $`F^{}`$ is formed from the heaviest component of the ISM. Equation (33) shows that the source $`Q_F^{P_0}`$ is the sum of the one-step contributions $`P_0+TF`$ plus the reactions involving all secondary nuclei $`F^{}`$ and $`T`$. This equation applies to any fragment $`F`$ and as the mass of $`F`$ increases the number of intermediaries becomes smaller. However a major change occurs when $`\mathrm{F}=\mathrm{P}_01`$ because for this nucleus there are no intermediaries that can produce $`F`$, only the one-step process involving $`P_0`$ and $`T`$ manufacture $`\mathrm{F}=\mathrm{P}_01`$. For this nucleus the source spectrum $`Q_{F=P_01}^{P_0}`$ is calculated with only the first half of equation (33) i.e

$$Q_{F=P_01}^{P_0}(E_{F=P_01})=\alpha _{P_0}\underset{T}{}y_Tn_H_0^{\mathrm{}}𝑑E_{P_0}\psi _{P_0}(E_{P_0})\sigma _{P_01}^{P_0,T}(E_{P_0})P_{P_01}(E_{P_01}|E_{P_0}).$$
(34)
The spectrum $`Q_{F=P_01}^{P_0}`$, and hence the flux $`\varphi _{F=P_01}^{P_0}`$, are proportional to $`\alpha _{P_0}`$ and therefore we introduce the normalized flux $`\mathrm{\Psi }_F^{P_0}`$ in the same fashion as we did in section §2 so that

$`\varphi _{F=P_01}^{P_0}(E_{F=P_01})`$ $`=`$ $`\alpha _{P_0}\mathrm{\Psi }_{F=P_01}^{P_0}(E_{F=P_01})`$ (35)
$`=`$ $`{\displaystyle \frac{1}{w_{P_01}(E_{P_01})\rho }}{\displaystyle _{E_{P_01}}^{\mathrm{}}}𝑑E_{P_01}^{}Q_{P_01}^{P_0}(E_{P_01}^{}){\displaystyle \frac{S_{P_01}(E_{P_01}^{})}{S_{P_01}(E_{P_01})}}.`$
Every quantity on the right-hand side of equation (35) is known so the flux $`\varphi _{F=P_01}^{P_0}`$ is exact. Now that we have $`\varphi _{P_01}^{P_0}`$ we turn our attention to $`Q_{F=P_02}^{P_0}`$, the heaviest fragment that can be produced from $`F^{}=P_01`$. For $`F=P_02`$ the sum over intermediaries $`F^{}`$ in equation (33) involves only the fragment $`F^{}=P_01`$ and we insert our solution for $`\varphi _{F=Z1}^{P_0}`$ from equation (35). Both $`\varphi _{P_0}^{(0)}`$ and $`\varphi _{F^{}=Z1}^{P_0}`$ are already known so the source spectrum $`Q_{F=Z2}^{P_0}`$, and hence the flux $`\varphi _{F=Z2}^{P_0}`$, may also be calculated exactly. As we step down through the nuclei the source spectrum $`Q_F^{P_0}`$ of any fragment $`F`$ is only dependent upon fluxes that are previously derived and so we do not introduce any errors into $`Q_F^{P_0}`$. Even the $`\alpha \alpha `$ fusion reactions that produce $`{}_{}{}^{6}\mathrm{Li}`$ and $`{}_{}{}^{7}\mathrm{Li}`$ can be accommodated in this scheme by first calculating $`\varphi _{{}_{}{}^{7}\mathrm{Li}}^\alpha `$ and then appropriately inserting this solution into the expression for $`Q_{{}_{}{}^{6}\mathrm{Li}}^\alpha `$. For every fragment $`F`$ the source spectrum is proportional to $`\alpha _{P_0}`$ and so we write the rate of change of abundance for $`F`$ for a given primary as

$$\frac{dy_F^{P_0}}{dt}=\alpha _{P_0}R_F^{P_0}=_0^{\mathrm{}}𝑑E_F^{}Q_F^{P_0}(E_F^{})S_F(E_F^{})$$
(36)
which defines the normalized cascade rate $`R_F^{P_0}`$.

Previously we extracted the dependence upon the ISM abundances $`y_T`$ but we have refrained in this cascade calculation because if we do then we run into the same problem as the reaction expansion, namely, vast arrays of rates to calculate for all the possible pathways from the primary $`P_0`$ to $`F`$. In order to calculate the source spectra, the fluxes of the intermediaries and the rate $`R_F^{P_0}`$ we need to specify a composition for the ISM. This introduces a chemical dependence we have been at pains to avoid but the sufficient accuracy of the one-step rates for the forward reactions, when $`P_0\{p,\alpha \}`$, means that we only follow this cascade for $`P_0\{Z\}`$ which is why we called the heaviest intermediary $`P_01`$. In the inverse calculations the heavy-element content of the ISM is irrelevant since all the reactions are with $`\mathrm{H}`$ and $`\mathrm{He}`$ as targets so we only need specify the abundance of Helium in the ISM. The dependence upon the ISM chemistry is therefore very weak since $`y_{\mathrm{He}}`$ is essentially constant. In table (4) we show the cascade rates $`R_F^{P_0}`$ for the case when the primary projectile is Oxygen together with the rate as calculated by using the reaction expansion up to the two-step order. The grammage is set at the constant value of $`\mathrm{\Lambda }=5\mathrm{g}/\mathrm{cm}^2`$, we use the straight-ahead approximation and the helium abundance is set to $`y_{\mathrm{He}}=0.1`$. It is clear from the table that the cascade and reaction-expansion rates agree to $`3\%`$ and thus we conclude that terminating the reaction expansion at the two-step level gives sufficiently accurate rates that the difference between the two methods is negligible.

## 6 Conclusions

The results of any calculation of the spallogenic rates of production for Lithium, Beryllium and Boron often find their application in trying to remove this synthesis mechanism from the observations in order to elucidate the other contributions to the abundances. Similarly, the ratio of the production rates has been frequently used to infer the spallogenic abundance of one of the LiBeB elements after observation of another. Both applications rely on the accuracy of the calculated rates so erroneous results can occur if the estimates for the contribution to the abundance from spallation are wrong. There are many aspects of the calculation where errors can enter but two do not have astrophysical origins: the straight-ahead approximation and the reaction expansion. We find that relaxing the straight-ahead approximation does increase the complexity of the problem because the distribution of the fragment energies must be taken into account but that using a more sophisticated description produces changes that are less than $`5\%`$. This is too small to insist that the straight-ahead approximation be rejected at the present time. Much larger changes were found when we included the two-step terms from the reaction expansion so that our faith in this methodology waned. However, after computing the rates to all order we found changes that the difference between the cascade and two-step rates is only a few percent and therefore sufficiently accurate rates can be calculated by terminating the reaction expansion at the two-step term. We therefore conclude that the errors arising from two assumptions in the nuclear physics aspects of the calculation can be removed leaving only that from the uncertainty in the cross-sections.

The authors would like to thank Gary Steigman for useful discussions. We acknowledge the support of the DOE through grant DE-FG02-91ER40690.