3. Models and equations

We list below the models present in the USINE package, providing their main parameters, equations, and references. Any propagation model is defined by:

  • the model properties in the propagation volume (geometry, ISM…);
  • the transport coefficients (diffusion \(K\), convection \(V_c\)…);
  • the solutions of the transport equation.

3.1. Diffusion equation

Cosmic-ray propagation in the Galaxy can be described by a diffusion equation (see, e.g., Berezinskii et al. 1990; Schlickeiser 2002; Strong et al. 2002):

\[\begin{split}\frac{\partial n(\vec r, p, t)}{\partial t} &+ \vec\nabla \cdot ( -K\vec\nabla n + \vec{V}_c n ) \\ &+ \frac{\partial}{\partial p} \left[\dot{p} n + \frac{p}{3} \, (\vec\nabla \cdot \vec{V}_c )n\right] - \frac{\partial}{\partial p}\, p^2 K_{pp} \frac{\partial}{\partial p}\, \frac{1}{p^2}\, n\\ &= {\rm Source}(\vec r, p, t) - {\rm Sink}(\vec r, p, t),\end{split}\]

with \(n (\vec r,p,t)\) the CR density per particle momentum \(p\), \(K\) the spatial diffusion coefficient, \(\vec{V}_c\) the convection velocity, \(\dot{p}\equiv dp/dt\) the momentum losses, and \(K_{pp}\) the diffusion coefficient in momentum space for reacceleration. This equation is formally the continuity equation (first line) with the energy current (second line), and source and sink terms (last line).

3.1.1. Source terms

Primary origin:Diffusive shock acceleration (e.g., in supernova remnants) is the favoured mechanism to accelerate all species, with an injection spectrum \(\propto R^{-\alpha}\) with \(\alpha\approx 2\) (\(R=pc/Ze\) is the rigidity). The CRs accelerated in sources are denoted primaries for short (e.g., \(^1{\rm H}\), \(^{16}{\rm O}\), \(^{30}{\rm Si}\)…).
Secondary origin:
 Nuclear interactions of primary CRs on the ISM give rise to secondary particles (\(e^+\), \(e^-\), \(\gamma\), \(\nu\)), anti-nuclei (\(\bar{p}\), \(\bar{d}\)…), and secondary nuclei as fragments of heavier CR nuclei. These CRs are denoted secondaries for short (e.g., Li to B isotopes).
Tertiary origin:
 Non-annihilating inelastic nuclear interactions are reactions in which the particles survive the interactions, but loose some of their energy (e.g., in resonances). This reaction actually both provides a net loss of particles at the energy of interaction, and a gain at lower energies coming from the redistribution of the particles in energy. These redistributed CRs are often denoted tertiaries for short, as they involve re-interactions of secondary CRs.
Radioactive origin:
 Unstable CRs are net sources for their progeny. The CR fraction decaying is set by the competition between decay times (e.g., 1.387 Myr for \(^{10}{\rm Be}\)) and escape time from the Galaxy (tens of Myrs at a few GeV/n): for instance, 15% of \(^{10}{\rm B}\) at 10 GeV/n comes the \(\beta\)-decay of secondary \(^{10}{\rm Be}\). Another decay type is electronic capture (EC) decay, in which CR ions must attach first an electron, the effective half-life being a competition between electron attachment, stripping, and EC-decay.

3.1.2. Sink terms

Destruction:Inelastic interactions on the ISM destroy CRs. The interaction time is constant above GeV/n energies, but gradually becomes a sub-dominant process as the escape time steadily decreases with energy. Also, heavy nuclei are more prone to destruction than lighter ones, as cross sections typically scale as \(A^{2/3}\).
Redistribution:A fraction of anti-nuclei which suffer non-annihilating interactions disappears. This is this very fraction that feeds the tertiary source term described above.
Decay:When unstable species decay, they change their nature. The various decay channels are the ones mirroring those described above in the radioactive source term.

3.2. Generic source terms

3.2.1. Primary source \(q_{\rm prim}^{i}\)

The primary source term, \(q_{\rm prim}^{i}\), is expressed in unit of \([{\rm \#part}~{\rm (GeV/n)}^{-1}~{\rm m}^{-3}~{\rm Myr}^{-1}]\).

\[q_{\rm prim}^i\equiv \frac{dQ_{\rm prim}^i}{dE_{k/n}} = \frac{A}{Z\beta}\frac{dQ_{\rm prim}^i}{dR}\]

Note

The term \(q_{\rm prim}^{p}\) is generically calculated in TUValsTXYZEVirtual::GetSrcSpectrum().

3.2.2. Secondary source \(q_{\rm sec}^{ij}\)

The secondary source term, \(q_{\rm sec}^{ij}\), is expressed in unit of \([\frac{dn}{dE_{k/n}}~{\rm Myr}^{-1}]=[{\rm \#part}~{\rm (GeV/n)}^{-1}~{\rm m}^{-3}~{\rm Myr}^{-1}]\).

  • Differential production \(\frac{d\sigma^{PF}}{dE_{\rm tot}}~[{\rm mb~GeV}^{-1}]\)

    Integrates on all incoming projectile energies \(E_{k/n}^{\rm in}\), the differential production cross sections of fragment at \(E_{k/n}^{\rm out}\):

    \[\frac{dQ_{\rm sec}^{i\rightarrow j}}{dE_{k/n}^{\rm out}}(E_{k/n}^{\rm out}) \! = \! \int_{0}^{\infty} \!\!\! dE_{k/n}^{\rm in} \;\frac{dn^i}{dE_{k/n}^{\rm in}}(E_{k/n}^{\rm in}) \times \sum_{\rm ISM} \left( n_{\rm ISM} \!\times\! v^{\rm in} \!\times\! A_j \!\times\! \frac{d\sigma^{i+{\rm ISM} \rightarrow j}}{dE_{\rm tot}^{\rm out}}\right).\]

    The extra factor \(A_j\) originates from the conversion from differential kinetic energy per nucleon to differential total energy, which is the format of production cross sections files (see Section 5.6).

  • Straight-ahead approximation for production \(\sigma^{ij}~[{\rm mb}]\)

    For nuclei, it is usually assumed that \(\frac{d\sigma^{i+{\rm ISM} \rightarrow j}}{dE_{k/n}} = \delta(E_{k/n}^{\rm out}-E_{k/n}^{\rm in})\times\sigma^{ij}\), so that

    \[\frac{dQ_{\rm sec}^{i\rightarrow j}}{dE_{k/n}}(E_{k/n}) =\frac{dn^i}{dE_{k/n}}(E_{k/n}) \times \sum_{\rm ISM} \left( n_{\rm ISM} \times v \times \sigma^{i+{\rm ISM} \rightarrow j}\right).\]

Note

The term \(q_{\rm sec}^{ij}\) is generically calculated in TUXSections::SecondaryProduction().

3.2.3. Tertiary source \(q_{\rm ter}^{i}\)

The tertiary source term (relevant for anti-nuclei only), \(q_{\rm ter}^{i}\), is expressed in unit of \([\frac{dn}{dE_{k/n}}~{\rm Myr}^{-1}]=[{\rm \#part}~{\rm (GeV/n)}^{-1}~{\rm m}^{-3}~{\rm Myr}^{-1}]\):

\[\begin{split}\frac{dQ_{\rm ter}^i}{dE_{k/n}^{\rm out}}(E_{k/n}^{\rm out}) &= \int_{E_{k/n}^{\rm out}}^{\infty} \frac{dn^i}{dE_{k/n}^{\rm in}}(E_{k/n}^{\rm in}) \!\times\! \sum_{\rm ISM} \left( n_{\rm ISM} \!\times\!v^{\rm in} \!\times\! A_i \!\times\! \frac{d\sigma^{i+{\rm ISM} \rightarrow i}}{dE_k^{\rm out}}\right) dE_{k/n}^{\rm in}\\ &-\frac{dn^i}{dE_{k/n}^{\rm out}}(E_{k/n}^{\rm out}) \!\times\! \sum_{\rm ISM} \left( n_{\rm ISM} \!\times\! v^{\rm out} \!\times\! \sigma^{\rm ina} (E_{k/n}^{\rm out})\right),\end{split}\]

with \(\sigma^{\rm ina}\) the inelastic non-annihilating cross section. The extra factor \(A_i\) originates from the conversion from total to kinetic energy per nucleon in the differential cross section (see Section 5.6). This equation is solved iteratively (e.g., Donato et al. 2001).

Note

the term \(q_{\rm ter}^{i}\) is generically calculated in TUXSections::TertiaryProduction().

3.2.4. Radioactive source \(q_{\rm rad}^{ij}\)

The radioactive progeny from an unstable nucleus gives rise to a radioactive source term \(q_{\rm rad}^{ij}\), expressed in unit of \([\frac{dn^i}{dE_{k/n}}~{\rm Myr}^{-1}]=[{\rm \#part}~{\rm (GeV/n)}^{-1}~{\rm m}^{-3}~{\rm Myr}^{-1}]\).

  • \(\beta\)-decay

    \[\frac{dQ_{\rm rad}^{i\rightarrow j}}{dE_{k/n}}(E_{k/n}) = \frac{dn^i}{dE_{k/n}}(E_{k/n}) \times \Gamma_{\rm rad}^{ij}\]

with the decay rate \(\Gamma_{\rm rad}^{ij}\) in \([{\rm Myr}^{-1}]\) related to the unstable nucleus half-life:

\[\Gamma_{\rm rad}^{ij} = \frac{\ln(2)}{\gamma \times t_{1/2}}.\]

Note

The term \(\Gamma_{\rm rad}^{ij}\) is generically calculated in TUAxesCrE::GammaRadBETA_perMyr().

  • EC-decay [TODO]

3.3. Solving for E-derivatives

In semi-analytical models considered below, spatial derivatives of the diffusion equation are solved analytically, so that only energy derivatives remain. The corresponding equation on the differential density \(n^{j} \equiv \frac{dn^j}{dE_{k/n}}\), in unit of \([{\rm \#part}~{\rm (GeV/n)}^{-1}~{\rm m}^{-3}~]\), for species \(j\), takes the form, for all models considered in the USINE package,

\[A^{j} n^{j}(E_{k/n}) + \frac{d}{dE} \left( B^{j} n^{j} - C^{j} \frac{dn^{j}}{dE} \right) = Q^{j}(E_{k/n}).\]

Fluxes are mostly power laws and kinetic energy per nucleon is approximately conserved in nuclear fragmentation reactions (straight-ahead approximation). This is the motivation to solve the above second order differential equation on a logarithmic scale in kinetic energy per nucleon. Using \(dE = dEk = \frac{d\ln Ek}{Ek}\) and \(E_{k/n}=E_k/A\), we get \(\frac{d}{dE_{\rm tot}} = \frac{1}{E_k}\frac{d}{d\ln E_{k/n}}\), so that the above equation is rewritten (omitting the indices for the sake of readability)

\[n + \tilde{\cal A} \frac{d}{d\ln E_{k/n}} \left( \tilde{\cal B} n - \tilde{\cal C} \frac{dn}{d\ln E_{k/n}} \right) = \tilde{\cal S}\]

with

\[\tilde{\cal A} = \frac{1}{A \times E_k},\quad \tilde{\cal B} = B,\quad \tilde{\cal C} = \frac{C}{E_k},\quad {\rm and~} \tilde{\cal S} = \frac{1}{A}.\]

This equation is solved using a finite difference scheme with boundary conditions. This amounts to a tridiagonal inversion. A detailed description of the chosen boundary conditions, their associated coefficients in the matrix, and the impact on the solution, as well as the stability of the numerical scheme, is provided in Sect. 3.1, App. C, and App D of Derome et al. (2009) respectively.

Note

Below, we provide for each model \(A^{j}\), \(B^{j}\), \(C^{j}\), and \(Q^{j}\). In the code, these coefficients are used to calculate \(\tilde{\cal A}\), \(\tilde{\cal B}\), \(\tilde{\cal C}\), \(\tilde{\cal S}\) on the energy grid, and the tridiagonal inversion is performed with the function TUNumMethods::SolveEq_1D2ndOrder_Explicit().

3.4. Leaky-box model

3.4.1. Description

The corresponding class is TUModel0DLeakyBox.

Assumptions
  • steady-state
  • no spatial dependence,
  • homogeneous distribution of gas,
  • homogeneous distribution of sources.
References for model and solution

3.4.2. Free parameters

  • \(n_i~[{\rm cm}^{-3}]\): density of the various elements to calculate the mean mass

    \(\displaystyle <n_{\rm ISM}m>~[{\rm g~cm}^{-3}] = \frac{\sum_{i \in {\rm ISM}} (n_i m_i)}{{\cal N}_{\rm Avogadro}}\),

  • \(\lambda_{\rm esc}(E)~[{\rm g~cm}^{-2}]\): escape length, related to the escape rate by

    \(\displaystyle \tau_{\rm esc}~[{\rm Myr}^{-1}] = \frac{\lambda_{\rm esc}}{v <n_{\rm ISM}m>}\)

  • \({\cal V}_A~[{\rm km~s}^{-1}~{\rm kpc}^{-1}]\): pseudo-Alfvénic speed of scatterers for reacceleration, related to a true speed \(V_a\) in diffusion models, see Ptuskin et al. (1994), by (\(h\) and \(L\) are the thin disc and diffusion halo size)

    \(V_a~[{\rm km~s}^{-1}] ={\cal V}_a\times (hL)^{1/2}\).

3.4.3. Equations and solution

The Leaky-box equation can be rewritten as

\[A^{j} n^{j}(E_{k/n}) + \frac{d}{dE} \left( B^{j} n^{j} - C^{j} \frac{dn^{j}}{dE} \right) = Q^{j}(E_{k/n}),\]

with (note that quantities are slightly differently defined in Putze et al. 2009)

\[\begin{split}n^{j}&~[{\rm \#part}~{\rm (GeV/n)}^{-1}~{\rm m}^{-3}~] \equiv \frac{dn^j}{dE_{k/n}}\\ A^{j}&~[{\rm Myr}^{-1}] = \frac{1}{\tau_{esc}} + \sum_{\rm ISM} n_{ISM} v^{j} \sigma_{inel}^{j+\rm ISM} + \frac{1}{\tau_{\beta-{\rm decay}}^{j}},\\ B^{j}&~[{\rm GeV~Myr}^{-1}] = < \frac{dE}{dt}>_{ion,coul.} + (1+\beta^{2}) \frac{{\cal K}_{pp}}{E},\\ C^{j}&~[{\rm GeV}^2~{\rm Myr}^{-1}] = \beta^{2} \times {\cal K}_{pp},\\ Q^{j}&~[{\rm \#part}~{\rm (GeV/n)}^{-1}~{\rm m}^{-3}~{\rm Myr}^{-1}] = q_{\rm prim}^{j} + \sum_{m_i>m_j} q_{\rm sec}^{ij} + q_{\rm tert}^{j} + \sum_{i \in \beta}q_{\rm rad}^{ij}.\end{split}\]

Implementation in USINE (specific to LB)

The default dependence for the escape length and pseudo-diffusion coefficients is:

\[\begin{split}\lambda_{\rm esc}~&[{\rm g~cm}^{-2}] \propto \lambda_0 R^{-\delta},\\ {\cal K}_{pp}~&[{\rm GeV}^2~{\rm Myr}^{-2}]=\frac{4}{3} {\cal V}_a^2\beta^2 E^2 \frac{\tau_{\rm esc}}{\delta(4-\delta^2)(4-\delta)}.\end{split}\]

However, because of implementation constraint in USINE, one must define a diffusion coefficient in space \(K_{00}\) and in momentum \(K_{pp}\) (in USINE parameter file), and we chose to define:

\[\begin{split}K_{00}^{\rm USINE} &\equiv \lambda_{\rm esc}~[{\rm g~cm}^{-2}],\\ K_{pp}^{\rm USINE} &\equiv {\cal K}_{pp} / \tau_{\rm esc},\end{split}\]

Solution

The Leaky-box equation is solved numerically, as explained in Solving for E-derivatives.

3.5. Standard 1D diffusion model

3.5.1. Description

The corresponding class is TUModel1DKisoVc.

Assumptions
  • steady-state
  • 1D geometry: thin plan (half-height \(h\)) + halo (half-height \(L\gg h\))
  • homogeneous distribution of gas and sources in thin plan,
  • isotropic and spatial-independent diffusion in plan and halo,
  • reacceleration in thin plan only,
  • constant wind in halo (implies discontinuity in plan).
References for model and solution

3.5.2. Free parameters

  • \(n_i~[{\rm cm}^{-3}]\): density of the various elements in thin plan
  • \(K(E)~[{\rm kpc}^2~{\rm Myr}^{-1}]\): spatial diffusion coefficient (isotropic and homogeneous)
  • \(K_{pp}(E)~[{\rm GeV}^2~{\rm Myr}^{-1}]\): momentum diffusion coefficient in thin plan
  • \(V_A~[{\rm km~s}^{-1}]\): Alfvénic speed in thin plan for reacceleration
  • \(V_c~[{\rm km~s}^{-1}]\): constant velocity in the halo (perpendicular to plan)
  • \(L~[{\rm kpc}]\): diffusive halo half-height
  • \(r_h~[{\rm kpc}]\): size of local underdense region (flux of radioactive nuclei suppressed)

3.5.3. Equations and solution

The model solution for a species \(j\) is analytical in the halo, and numerical in the thin disc (see its resolution in Solving for E-derivatives):

\[\begin{split}n^j(z) &= n^j(z=0) \times \exp^{(V_{c}z/2K)} \frac{\sinh(S^j(L-z)/2)}{\sinh(S^jL/2)},\\ A^j n^j(z=0) &+ 2h \times \frac{d}{dE} \left( B^j n^j(z=0) - C^j \frac{dn^j}{dE}(z=0) \right) = 2h \times Q^{j}(z=0),\end{split}\]

with (\(z\) is the vertical coordinate so that \(z=0\) in in the thin disc)

\[\begin{split}n^j(z)&~[{\rm \#part}~{\rm (GeV/n)}^{-1}~{\rm m}^{-3}~] \equiv \frac{dn^j}{dE_{k/n}}(z)\\ A^j&~[{\rm kpc~Myr}^{-1}] = V_{c} + 2h\Gamma^{j}_{\rm inel} + K S^{j}\coth\left(\frac{S^{j}L}{2}\right),\\ S^j&~[{\rm Myr}^{-1}] = \sqrt{\frac{V_{c}^{2}}{K^{2}} + 4 \frac{\Gamma_{\rm rad}^{j}}{K}},\\ B^{j}&~[{\rm GeV~Myr}^{-1}] = < \frac{dE}{dt}>_{\rm ion,coul.,adiab.} + (1+\beta^{2}) \frac{K_{pp}}{E},\\ C^j&~[{\rm GeV}^2~{\rm Myr}^{-1}] = \beta^{2} \times {\cal K}_{pp},\\ Q^j(z=0)&~[{\rm \#part}~{\rm (GeV/n)}^{-1}~{\rm m}^{-3}~{\rm Myr}^{-1}] = q_{\rm prim}^{j} + \sum_{m_i>m_j} q_{\rm sec}^{ij} + q_{\rm tert}^{j} + \sum_{i \in \beta}q_{\rm rad}^{ij}.\end{split}\]

Note

Radioactive nuclei also decay in the halo and contribute to a complicated term not described above (see Putze et al. 2010).

Default parameters

The default dependence for the spatial and momentum diffusion coefficients are:

\[\begin{split}K(R)~&[{\rm kpc}^2~{\rm Myr}^{-1}] = \beta K_0 R^\delta\\ K_{pp}~&[{\rm GeV}^2~{\rm Myr}^{-1}]=\frac{4}{3} V_A^2\beta^2 E^2 \frac{1}{\delta(4-\delta^2)(4-\delta)K(R)}.\end{split}\]

3.6. Standard 2D diffusion model

3.6.1. Description

The model is very similar to its 1D counterpart. The corresponding class is TUModel2DKisoVc.

Assumptions
  • steady-state
  • 2D geometry: thin disc (half-height \(h\)) and cylindrical halo (half-height \(L\gg h\)) of radius \(R\),
  • homogeneous distribution of gas in thin disc,
  • radially-dependent distribution of sources in thin disc,
  • isotropic and spatial-independent diffusion in disc and halo,
  • reacceleration in thin disc only,
  • constant wind in halo (implies discontinuity in disc).
References for model and solution

3.6.2. Free parameters

  • \(n_i~[{\rm cm}^{-3}]\): density of the various elements in thin disc
  • \(K(E)~[{\rm kpc}^2~{\rm Myr}^{-1}]\): spatial diffusion coefficient (isotropic and homogeneous)
  • \(K_{pp}(E)~[{\rm GeV}^2~{\rm Myr}^{-1}]\): momentum diffusion coefficient in thin disc
  • \(V_A~[{\rm km~s}^{-1}]\): Alfvénic speed in thin disc for reacceleration
  • \(V_c~[{\rm km~s}^{-1}]\): constant velocity in the halo (perpendicular to disc)
  • \(L~[{\rm kpc}]\): diffusive halo half-height
  • \(r_h~[{\rm kpc}]\): size of local underdense region (flux of radioactive nuclei suppressed)
  • \(R~[{\rm kpc}]\): radius of the cylindrical box

3.6.3. Equations and solution

An expansion along the first order Bessel function is performed

\[n(r,z) = \sum_{i=1}^{\infty} n_i(z) J_0\left(\zeta_i \frac{r}{R}\right),\]

with \(\zeta_i\) the \(i\)-th zero of the zero-th order Bessel function \(J_0\). This automatically ensures the boundary condition \(n(r=R,z)=0\). Each Bessel coefficient \(n_i(z)\) follows an equation very similar to that of the 1D version of the model.

The model solution for Bessel order \(i\) of a species \(j\) is analytical in the halo, and numerical in the thin disc (see its resolution in Solving for E-derivatives)

\[\begin{split}n_i^j(z) &= n_i^j(z=0) \times \exp^{(V_{c}z/2K)} \frac{\sinh(S_i^j(L-z)/2)}{\sinh(S_i^jL/2)},\\ A_i^j n_i^j(z=0) &+ 2h \times \frac{d}{dE} \left( B^j n_i^j(z=0) - C^j \frac{dn_i^j}{dE}(z=0) \right) = 2h \times Q_i^j(z=0),\end{split}\]

with (\(z\) is the vertical coordinate so that \(z=0\) in in the thin disc)

\[\begin{split}n_i^j(z)&~[{\rm \#part}~{\rm (GeV/n)}^{-1}~{\rm m}^{-3}~] \equiv \frac{dn_i^j}{dE_{k/n}}(z)\\ A_i^j&~[{\rm kpc~Myr}^{-1}] = V_{c} + 2h\Gamma^{j}_{\rm inel} + K S_i^j\coth\left(\frac{S^{j}L}{2}\right),\\ S^j&~[{\rm Myr}^{-1}] = \sqrt{\frac{V_{c}^{2}}{K^{2}} + (2\zeta_i/R)^2 + 4 \frac{\Gamma_{\rm rad}^{j}}{K}},\\ B^j&~[{\rm GeV~Myr}^{-1}] = < \frac{dE}{dt}>_{\rm ion,coul.,adiab.} + (1+\beta^{2}) \frac{K_{pp}}{E},\\ C^j&~[{\rm GeV}^2~{\rm Myr}^{-1}] = \beta^{2} \times {\cal K}_{pp},\\ Q_i^j(z=0)&~[{\rm \#part}~{\rm (GeV/n)}^{-1}~{\rm m}^{-3}~{\rm Myr}^{-1}] = q_{i,\,\rm prim}^j + \sum_{m_k>m_j} q_{i,\,\rm sec}^{kj} + q_{i,\,\rm tert}^{j} + \sum_{k \in \beta}q_{i,\,\rm rad}^{kj},\end{split}\]

where \(q_i\) terms are Bessel-Fourier counterparts (or expansions) of the source term.

Note

Radioactive nuclei also decay in the halo and contribute to a complicated term not described above (see Putze et al. 2010).

Default parameters

The default dependence for the spatial and momentum diffusion coefficients are:

\[\begin{split}K(R)~&[{\rm kpc}^2~{\rm Myr}^{-1}] = \beta K_0 R^\delta\\ K_{pp}~&[{\rm GeV}^2~{\rm Myr}^{-1}]=\frac{4}{3} V_A^2\beta^2 E^2 \frac{1}{\delta(4-\delta^2)(4-\delta)K(R)}.\end{split}\]