RECOVERING THE OBSERVED B/C RATIO IN A DYNAMIC SPIRAL-ARMED COSMIC RAY MODEL

There are three important timescales describing this system.³ The first is the typical time it takes CRs to diffuse out of the galactic halo:

$$\begin{equation} \tau _{{\rm escape}}\approx {Z_h^2 \over 2 D}, \end{equation} \tag{ 1 }$$

where Z_h is the half width of the galactic halo and D is the diffusion coefficient. The second timescale is the typical time it takes CRs to diffuse from the spiral arm to the solar system (s is the distance from the solar system to the nearest spiral arm):

$$\begin{equation} \tau _{s}\approx {s^2 \over 2D}. \end{equation} \tag{ 2 }$$

The two timescales scale with energy through D⁻¹, and therefore, their ratio is energy independent. That is, if some fraction of the CRs are lost from the galactic halo before they reach the solar system, it will be the same fraction at all energies.

Which timescale dominates depends on the question asked. For example, if we care about the average grammage of a steady source originating from spiral arms, then the average grammage will be proportional to max (τ_escape, (τ_escapeτ_s)^1/2) (Shaviv et al. 2009). Namely, the effective diffusion timescale is the larger of the diffusion out of the galactic halo or the geometric mean of the two diffusion timescales.

The third timescale is the time since the last passage of the solar system through a spiral arm⁴:

$$\begin{equation} \tau _{{\rm arm}} \approx {s/ \sin i \over v_{{\rm arm}}}, \end{equation} \tag{ 3 }$$

where v_arm is the pattern speed of the arm relative to the local ISM, i is the pitch angle of the arm,⁵ and s/sin i is approximately the distance that the solar system has traversed from the last arm passage.

At sufficiently low energies, τ_arm is the shortest timescale, and it therefore dominates. Hence, if it takes say 20 Myr for CRs at a given energy to diffuse from the spiral arms to the solar system, but the solar system last passage through the spiral arms was 10 Myr ago, then the CRs in the solar vicinity should be around 10 Myr old, rather than 20 Myr (unless the diffusion out of the galactic halo is significantly longer than the diffusion from the spiral arm).

Given the typical age t of the CRs, which may be determined by diffusion from the arm, out of the galactic halo, or by advection from the spiral arm, we can estimate the amount of secondaries produced as a function of energy. Over a typical time, t, the grammage, g, seen by a propagating CR is:

$$\begin{equation} g \approx x \bar{\rho } \approx t v \bar{\rho }, \end{equation} \tag{ 4 }$$

where x = tv is the physical path length traveled by the CR particle moving at a speed v and $\bar{\rho }$ is the average density along the path.

There are four different regimes for the solution depending on whether the propagation is dominated by advection or by diffusion and whether or not the particle is relativistic.

In the limit dominated by diffusion, one has $t \sim \max (\tau _{{\rm escape}},(\tau _{{\rm escape}} \tau _s)^{1/2}) = \max (Z_h^2,s^2)/2D$. The dependence of the age on energy depends, in turn, on the energy dependence of the diffusivity. We assume that the effective mean free path (mfp) of the CRs depends on the rigidity as ℓ∝R^δ. For Kolmogorov turbulence, we should expect δ ∼ 1/3, but any value up to 0.6 is theoretically feasible. The diffusivity is then of the order of D ∼ ℓv. For relativistic particles, R∝p∝E, and therefore:

$$\begin{equation} g_\mathrm{rel,diff} \propto {\max \left(Z_h^2,s^2 \right) \over E^\delta } \bar{\rho }. \end{equation} \tag{ 5 }$$

For non relativistic particles, $R \propto p \propto \sqrt{E_{{\rm kin}}}$, and:

$$\begin{equation} g_\mathrm{non-rel,diff} \propto {\max \left(Z_h^2,s Z_h \right) \over E^{\delta /2} } \bar{\rho }. \end{equation} \tag{ 6 }$$

At lower energies, in the advection dominated regime, t ∼ τ_arm ≈ s/(v_armsin i). This lifetime is independent of energy. Therefore,

$$\begin{eqnarray} g_\mathrm{rel,adv} &\propto & {sc \over {v_{{\rm arm}} \sin i}} \bar{\rho } \propto {\mathrm{c}onst.}, \end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray} g_\mathrm{non-rel,adv} &\propto & {sv\over {v_{{\rm arm}} \sin i}} \bar{\rho } \propto E^{1/2}. \end{eqnarray} \tag{ 8 }$$

In particular, in the advection dominated nonrelativistic case, the grammage is an increasing function of energy. This is one of the cardinal results of this work.

If we combine the results together, we find that there are two net possible energy dependences, depending on whether the diffusion/advection behavior change takes place for nonrelativistic or for relativistic particles. The results are summarized in Figure 1.

Zoom In Zoom Out Reset image size

Under standard diffusion models, the velocity cancels out (e.g., see Strong et al. 2007, p. 295), such that to obtain an increase in the B/C ratio one requires a nonstandard diffusion model, having a diffusion coefficient that depend not only on the rigidity (Di Bernardo et al. 2010). These results should be compared with the standard solutions to the B/C behavior. The B/C behavior at low energies requires a grammage which increases with energy. Therefore, any solution for which the time it takes the CRs to reach the solar system saturates at a finite value that is independent of the energy can recover the B/C behavior. One option is simply to assume that the diffusivity below E_crit/nucleon ∼ 4 GeV is constant. Here, the age saturates at $\tau _{{\rm max}} \sim Z_h^2/D(<E_{{\rm crit}})$. In the solution where a galactic wind with a velocity v_wind is added the age saturates at τ_max ∼ Z_h/v_wind, the time it takes the wind to evacuate the galactic halo from CRs. Another standard solution is that of re-acceleration. Here, low-energy particles are mixed giving them the same average age.

The break energy between the rising B/C and the decreasing B/C can be obtained by comparing the relevant diffusion timescale to the advection timescale, max (τ_escape, τ_s) ≈ τ_arm. Assuming that the break takes place when the particles are relativistic we have:

$$\begin{equation} E_{{\rm break}} \approx E_0 \left(s v_{{\rm arm}} \sin i \over 2 D_0 \right)^{1/\delta }, \end{equation} \tag{ 9 }$$

where we assumed a diffusion coefficient of D = D₀(E/E₀)^δ.

For example, if we reasonably assume that the galaxy has four arms with a 28° pitch angle, rotating with a periodicity (arm passage relative to the solar system) of 145 Myr, and a diffusion coefficient D = D₀(E/E₀)^δ with D₀ = 2 × 10²⁷ cm⁻² s⁻¹, E₀ = 3 GeV and δ = 0.5, we find that E_break ≈ 1.5 GeV if the last spiral arm passage was 5 Myr ago. This break energy will be around 1 GeV and can explain the observed B/C behavior. However, a full numerical calculation is needed to determine the accurate spectrum and the break obtained from a given set of parameters. Moreover, one has to remember that additional factors come into play at this energy range, such as solar modulation and Coulomb cooling.

The source of CRs at the relevant energy range is SNRs. We divide the SN population into degenerate collapse SNe (type Ia) and core collapse SNe (type II and Ibc). Given the observational data of SNe in different spiral armed galaxies, we expect the type Ia contribution to be around 15%–20% (Tammann et al. 1994; van den Bergh & McClure 1994). Here, we take 20%.

Since the progenitors of core collapse SNe live less than 30 Myr, they explode not far from their birth sites, which is primarily the spiral arms of the galaxy (e.g., see Lacey & Duric 2001). Because about 5% of O stars and 15% of H ii regions reside outside of the spiral arms (see Shaviv 2003 and references within), we assume here that 10% of the core collapse SNe are in the "field," i.e., distributed in a homogeneous disk, while the rest are divided equally between two sets of arms (see Section 3.2). We note that this disk is different and generally narrower than the disk-like halo in which CR diffuse.

For the actual radial and vertical distribution of Type Ia and core collapse SNe, we take the respective distributions in Ferrière (2001). For the sets of arms, we take the radial and vertical distribution of the field SNe and add an azimuthal component, as in Shaviv (2003), with spiral arm location and dynamics described below.

There are several indications that the Milky Way has two sets of spiral arms (e.g., see Amaral & Lepine 1997; Naoz & Shaviv 2007 and references therein), one with an m = 4 symmetry and one with an m = 2 symmetry. This can also be seen from the work of E. Gavish & N. J. Shaviv (in preparation), in which the two different sets can be from ln (r) − ϕ plots based on the CO maps of Dame et al. (2001). The nominal numbers we take are i = 28° for the four-armed set, and i = 11° for the two-armed set.

We approximate the spiral arms as logarithmic spiral arms, as is expected from the density wave theory (e.g., Binney & Tremaine 1988):

$$\begin{equation} r_j = a_j \exp [ b_j (\phi - \phi _{0,j})], \end{equation} \tag{ 10 }$$

where a_j is the inner radial extent of the arm j, and b_j = tan (i_j) where i_j is the pitch angle of arm j and r, ϕ, and z are cylindrical coordinates centered on the Milky Way galaxy.

Given that there are two sets of spiral arms, it is not surprising that different groups have found different pattern speeds which cluster around two values (see the summary in Shaviv 2003). This is also why an unbiased analysis of open clusters, which allows for multiple patterns, found both pattern speeds (Naoz & Shaviv 2007). Following the latter analysis, we take Ω₂ = 25 (km s⁻¹) kpc⁻¹ for the two-armed set (which is nearly corotating with the solar system; see below) and Ω₄ = 15 (km s⁻¹) kpc⁻¹ for the four-armed set.

We take the rotation curve of Olling & Merrifield (1998) in our model.

We assume that the CRs undergo normal diffusion without additional advective processes relative to the ISM which itself is differentially rotating around the galaxy. This diffusion coefficient is assumed to be location independent up to the galactic halo height Z_h, beyond which CRs can escape without diffusion.

The diffusion coefficient is also assumed to be energy dependent, with the general form $D = D_0 \beta R_{GV}^\delta$, where β = v/c and R_GV is the rigidity in GV. δ is the diffusion coefficient power-law index, which depends on the physical model for the diffusion.

One can show that a Kolmogorov power-law spectrum for the interstellar turbulence having δB ≈ 5 μG gives rise to a diffusion coefficient $\kappa \sim 2 \times 10^{27} \beta R_{GV}^{1/3}\ {\rm cm^{-2}\ s^{-1}}$ (see Strong et al. 2007 and references therein), i.e., δ = 1/3. Other, non-Kolmogorov models for the ISM magnetic turbulence can result in different power laws. For example, Kraichnan-type turbulence produces a steeper spectrum, with δ = 1/2 (Yan & Lazarian 2004).

Therefore, CR diffusion models typically consider a diffusion power law index of δ ∼ 0.3 to 0.6 (Strong et al. 2007). We use a nominal value of δ = 1/2 as it appears to best explain the high energy dependence of B/C.

A major ingredient of the model is the nuclear spallation reactions which the CRs undergo during their propagation. Since we concentrate here on light isotopes, we consider all the spallation and decay reactions of isotopes up to oxygen. Given that the contribution by the heavier nuclei, between Fluorine and Silicon, is only about 11% of the total mass between carbon and silicon, one would at most get a 10% correction if all the above-oxygen elements break up to boron and carbon. We therefore neglect the contribution of the heavier elements in this work, and keep in mind that subsequent analyses will require including these elements if a few percent level accuracy is needed. A summary of the reactions we consider can be found in Garcia-Munoz et al. (1987).

Since the typical timescale for diffusion and for spallation is 10⁷ yr, there is no need to evolve all the short radioactive isotopes separately. That is, when a short lived isotope is formed, it is immediately converted into its daughter products. If there is a branching ratio between different decay processes, then each decay branch is randomly followed with the proper weight. Below oxygen, there is only one long-lived radioactive isotope which is explicitly followed, and that is ¹⁰Be.

For the partial spallation cross-section, we use the yieldx program of Silberberg & Tsao (1990). We use the formulae of Karol (1975) for the total inelastic cross-sections of interactions with helium.

We assume a fixed ratio between oxygen and carbon at the source. As we shall see in Section 4.5, this ratio of 1.45 is obtained by fitting to the observed C/O ratio.

At energies below about 1 GeV, the propagating nucleons lose energy through primarily Coulomb collisions and ionization losses. The expression for Coulomb collisions is approximately (e.g., Strong & Moskalenko 1998)

$$\begin{equation} \left. {\it dE} \over dt \right|_{\mathrm{Coul}} \approx -{4 \pi r_e^2 c^3 m_e Z^2 n_e \ln \Lambda \over \beta }, \end{equation} \tag{ 11 }$$

where ln Λ ∼ 40–50 is the Coulomb logarithm and r_e is the classical electron radius. This expression assumes that the typical velocity βc is much higher than the thermal speed of electrons in the ISM.

Equation (11) can be written as an energy loss with grammage, once we consider that dg ≈ (n_pm_pβc/X)dt. Here, X is the hydrogen mass fraction.

$$\begin{eqnarray} &&\left. d(E/A) \over dg \right|_{\mathrm{Coul}} \; \approx \; -{4 \pi X \left(n_e \over n_p\right) \left(m_e \over m_p\right) \left(Z^2 \over A \right) {r_e^2 c^2 \over \beta ^2} \ln \Lambda }. \nonumber\\ \end{eqnarray} \tag{ 12 }$$

When written in this form, it is clear that for each particle, the Coulomb "cooling" is proportional only to β⁻². This expression can be easily incorporated into the simulation, as described below.

The ionization losses have a similar expression to the Coulomb losses. The difference is that instead of ln Λ, there is a more complicated expression that we use (see Strong & Moskalenko 1998). For typical ISM conditions, the ionization losses below 1 GeV are roughly one-fifth of the Coulomb collisional losses.

The CR spectrum observed at Earth is different from the spectrum outside the solar system because of interaction with the solar wind. To a good approximation, this interaction can be quantified through the solar modulation parameter ϕ, which is the energy per proton lost by the CR particle (e.g., Usoskin et al. 2011). The energy per nucleon lost is

$$\begin{equation} (E/A)_{{\rm earth}} = (E/A)_{{\rm ISM}} - (Z/A) \phi. \end{equation} \tag{ 13 }$$

Typically, ϕ ranges between 300 MV and 500 MV in the periods when the measurements are taken (e.g., Usoskin et al. 2011), which is the range we take here.

The CRs spatial distribution from the axisymmetric components does not depend on the rotation of the underlying ISM. However, when considering the diffusion from the spiral arms, the rotation of the arms and the rotation of the ISM should be considered.

We solve the diffusion from each spiral set (as well as from the galactic disk) separately. For each set, we solve in the frame of reference of the moving spiral arm set. On one hand, it implies that we need to add the differential rotation of the galactic disk as a radially dependent drift velocity, V_arm(r) − V_disc(r), which may appear as a complication. On the other hand, it turns the time-dependent problem into a time-independent one. This is because the boundary conditions in the frame of reference which is rotating with the spiral arm are fixed in time.

Thus, at a given time, the CR density distribution is obtained by rotating the distribution from each set by the appropriate rotational phase and adding the distribution from the axisymmetric sources.

For the initial relative abundances for the primary elements, we take 40% carbon, 2% nitrogen, and 58% oxygen. The relative abundances of each isotope for a given element are assumed to be solar. As we shall see in Section 4.5, these abundances recover the observed C/O and N/O ratios.

For the ISM components we take the distributions given in Ferrière (2001). The total ISM density distribution is required for the determination of various CR/ISM interactions.

Our code is a Monte Carlo simulation. This implies that we solve the full diffusive propagation of a large sample of particles, having different initial energies, different origins (arms or disk), and representing different compositions. The particles are followed from their formation until they either leave the galaxy, break into lighter particles than beryllium or cool below 100 MeV. For each particle, we also follow its various interactions, including reactions which change their identity, that is, spallation on the ISM and radioactive decay. When this occurs, we follow the lighter particle that formed. By simulating many particles with different initial energies, we can construct the CR spectrum at different locations in the galaxy and at different times.

The methodology we chose to follow, that of a Monte Carlo simulation, is different from most present day simulations (such as galprop, Strong & Moskalenko 1998, and dragon, Di Bernardo et al. 2010), which solve the diffusion partial differential equations (PDEs). The main advantage of the latter type is that the accuracy increases linearly with the number of grid points. The disadvantage is that it is much harder to include new physical ingredients, such as time dependence or different geometries.

The price one pays when simulating using the Monte Carlo method is that the accuracy increases only as the square root of the number of particles used. However, there are several major advantages. First, it is almost straightforward to include any physical effect. In our case, we require spatially and temporally dependent CR sources associated with the spiral arms. The code is also easily parallelizable because each simulated CR particle is independent of other CRs.

Another very important advantage is that the different CR variables, including artificial ones, can be easily recorded. The best example is the CR age and the amount of grammage that the CR has traversed. This allows obtaining the path length and age distributions. In fact, if one artificially switches off cooling, it is possible to see at low energies multiple passages of the spiral arm in the path-length distribution (PLD). This is impossible to obtain when solving the diffusion PDEs.

We have used 10¹⁰ CR particles in the simulations. For each CR, we randomly choose an initial isotopic identity, as described in Section 3.8. We then continue to randomly choose the initial location in the galaxy, where each particle was accelerated, with a distribution described in Section 3.1.

Instead of simulating the small steps corresponding to the physical mfp of the CRs, around 1 pc, our MC simulation uses a larger effective mfp, which lumps together many small effective physical steps. This corresponds to a larger effective time step chosen to be τ_esc/100. For our nominal model, this corresponds to an effective mfp of 25 pc. To avoid numerical error, we keep the effective mfp smaller than the typical length scale over which the gas density varies. To maintain accuracy, we decrease the time step 1000 fold when the CR particle is in the vicinity of the solar system (within a sphere of radius 0.2 kpc). This roughly corresponds to an effective mfp of 1 pc, which is of the order of the real mfp. At each time step, we check whether the particle had left the galaxy, and if it did not, we calculate, using the cross-sections described in Section 3.4, the grammage that the particle traversed in this time step. We then calculate the probability that the particle had undergone a nuclear reaction with ISM nuclei. If the particle did break into a secondary particle, we then follow each secondary particle with the same methodology until it either breaks into a particle lighter than beryllium until it escapes the galaxy or until it cools below 100 MeV.

Although the energy and location of each simulated particle is continuous, both are recorded with a set of four-dimensional arrays (location + energy), for each of the stable isotopes between ⁹Be to ¹⁶O, and the long lived ¹⁰Be.

The CR particles are not restricted to a grid, but the densities are recorded on a coarse grid covering the galaxy, and a finer grid within a sphere of radius 0.1 kpc around the solar system. The coarse grid spans a volume of 30 kpc × 30 kpc × 2Z_h, where Z_h is the height of the galactic halo, which depends on the chosen model parameters. It is composed of (0.1 kpc)³ sized cells. Within the finger grid includes (1 pc)³ sized cells. The energies are recorded in 10 bins per decade, covering the range of 0.1–1000 GeV.

We begin with "model A," by comparing our results with the simplest homogeneous model, to verify the validity of the numerical methods.

We first simulate the B/C ratio obtained in a model for which the spallation cross-sections are artificially fixed to their values at 1 GeV/nucleon. This allows us to verify that the energy dependences are the power-laws we expect (see Section 2.1). Figure 2 demonstrates that the Monte Carlo simulation agrees well with the analytical approximation. Specifically, we obtain the expected broken power-law driven by the diffusion law, with a break at 1 GeV/nucleon associated with the change from the nonrelativistic to relativistic regime. At nonrelativistic velocities, the power law is still less than δ/2 because of the saturation effect of the spallation.

Zoom In Zoom Out Reset image size

The next model we simulate, Model B, includes the energy-dependent nuclear cross-sections and cooling, such that we recover the results of galprop's simplest model (specifically, model 01000 of Strong & Moskalenko 1998, having no break in the diffusion-law, no wind, and no re-acceleration). The model has the following parameters: D₀ = 7 × 10²⁷ cm⁻² s⁻¹, Z_h = 1 kpc, and δ = 0.6. D₀ is normalized at a rigidity of 3 GeV/nucleon. The results are plotted in Figure 3, which should be compared with Figure 3 of Strong & Moskalenko (1998). Although there are several small differences, the general agreement is apparent. The inconsistency at low energies can be explained given that the two models still have small differences. These include different cross-section tables as well as different spatial distributions for the SNe and the ISM density. Unlike the galprop (Strong & Moskalenko 1998) and dragon (Di Bernardo et al. 2010) simulations, our present simulations can be used to obtain the path length distributions. These are plotted in Figure 4. As expected, the distributions are exponentials. This also explains why the leaky box model is generally a good approximation for CR diffusion in the homogeneous disk.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

After having validated our numerical simulations by recovering previous results for models with a homogeneous galactic source distribution, we can proceed to present the results of our model, where a significant fraction of the CR is formed near the galactic spiral arms.

When choosing model parameters, one has to remember that previously considered "canonical" values are not binding. The reason is that the different path length distribution expected in the spiral arm model, which lacks short path lengths, requires a smaller typical diffusivity and halo size than the standard values obtained in the homogeneous models if the model is to fit observations (Shaviv et al. 2009). Our goal in the present work is not to carry out an extensive parameter study but instead to show that model parameters exist for which the different isotopic spectra are recovered (see the discussion in Section 5).

The parameters of Model C that we use are D = 5.5 × 10²⁶ cm² s⁻¹, Z_h = 250 pc, and δ = 0.5. Following actual spiral arm passage about 20 Myr ago, we assume that the peak CR acceleration took place 6 Myr ago. The four arms have a pitch angle of i = 28°, while the two-arm set has 11°. Other parameters are as described in Section 3.

Figure 5 depicts the CR density at the galactic plane for two energies. Clearly, the inhomogeneous source distribution gives rise to an inhomogeneous CR distribution in the galactic disk. Moreover, the effects of the different timescales is evident as well. At low energies, the spiral arm passage and escape timescales are comparable. This implies that the advection relative to the spiral arms becomes important, and one can see the "smearing" associated with it. At high energies, the escape is much faster and one can consider the spiral arms to be effectively at rest (cf. Section 2.1).

Zoom In Zoom Out Reset image size

The path length distributions that we obtain at several energies are depicted in Figure 6. One of the unique features of the present model is that the distributions are not exponential. Instead, they center around the typical grammage required for diffusion from the nearby arm. As we shall elaborate upon in the discussion, this nonexponential behavior has various interesting ramifications. The difference between the PLD in Model B (Figure 4) and the PLD in Model C (Figure 6) is the underlying reason why we do not need a large halo in our model.

Zoom In Zoom Out Reset image size

Another interesting result is the different spectral indices of CR inside and outside the spiral arms. A nondynamic model can give no spectral difference. Here, however, Figure 7 reveals that the ratio between the density at different energies becomes location-dependent, with the ratio becoming "harder" in the spiral arms. For example, the ratio between 100 GeV CRs and 1 GeV CRs in the wake of the arms can be as large as 10^0.5 lower. In other words, the spectral index can soften by about 0.25 when leaving the arms. This is consistent with observations of γ-ray emission, which show that the spectrum in the direction of the orion arm has a spectral index harder by 0.4 ± 0.2 than in a high galactic latitude direction (Rogers et al. 1988; Bloemen 1989).

Zoom In Zoom Out Reset image size

As predicted in Section 2.1, adding the dynamic spiral arms saturates the age at low energies, and gives a B/C ratio that increases with energy. Figure 8 depicts the B/C for our nominal model.

Zoom In Zoom Out Reset image size

The ¹⁰Be/⁹Be provides information on the typical age of the CRs at different energies. The actual relation between the average age and the ¹⁰Be/⁹Be ratio depends on the path length distribution. Therefore, we compare the model predictions with the ¹⁰Be/⁹Be ratio directly. Figure 9 depicts the ¹⁰Be to ⁹Be ratio we obtain and its consistency with observations.

Zoom In Zoom Out Reset image size

To verify that the initial ratio at the sources of C to O is consistent, we plot the carbon to oxygen ratio we obtain in Figure 10 and compare it with the observed ration.

Zoom In Zoom Out Reset image size

Because our model considers all isotopes up to ¹⁶O, we can also predict the nitrogen to oxygen ratio (note that nitrogen in CRs is mostly secondary). The behavior is as expected similar to the B/C ratio (see Figure 11). However, it should be noted that the slope of the predicted C/O and N/O ratios is larger than the observed slope. Namely, there is some unexplained inconsistency between the power law index δ required to best fit the B/C ratio and the δ which best fits the C/O and N/O ratios. This inconsistency exists in the standard model and it remains here as well.

Zoom In Zoom Out Reset image size