Studies on Cosmic-Ray Nuclei with Voyager, ACE, and AMS-02. I. Local Interstellar Spectra and Solar Modulation

GCRs are modulated by the heliospheric magnetic field carried by solar winds when they enter the heliosphere, resulting in suppression of their fluxes. This solar modulation effect depends on particle energies and is particularly obvious at low energies. In this work, we adopt the force-field approximation of the solar modulation (Gleeson & Axford 1967, 1968), which was actually an approximate solution of Parker's equation (Parker 1965). In this model, the TOA flux is related to the LIS flux as

$$\begin{eqnarray}&&{J}^{\mathrm{TOA}}(E)={J}^{\mathrm{LIS}}(E+{\rm{\Phi }})\times \displaystyle \frac{E(E+2{m}_{p})}{(E+{\rm{\Phi }})(E+{\rm{\Phi }}+2{m}_{p})},\end{eqnarray} \tag{ 1 }$$

where E is the kinetic energy per nucleon, Φ = ϕ · Z/A with ϕ being the solar modulation potential, m_p = 0.938 GeV is the proton mass, and J is the differential flux of GCRs. The only parameter in the force-field model is the modulation potential ϕ. In principle, the force-field model assumes a quasi-steady-state of the solution of Parker's equation. However, the observational GCR fluxes show 11 year variations associated with solar activities. Therefore a time series of ϕ at different epochs is adopted to describe the data.

Usually power-law or broken power-law functions are explored to fit the GCR data. If the observational data cover a wide enough energy range, one can instead use a non-parametric method by means of spline interpolation of GCR fluxes among a few knots (Ghelfi et al. 2016). The spline interpolation is a way to obtain an approximate function smoothly passing through a series of points using piecewise polynomial functions. We use the cubic spline interpolation here, with the highest-order polynomial of three. We work in the $\mathrm{log}(E)-\mathrm{log}(J)$ space of the energy spectrum. The positions of knots of $x=\mathrm{log}(E)$ for helium, boron, carbon, and oxygen are defined as

$$\begin{eqnarray} & & \{{x}_{1},{x}_{2},{x}_{3},{x}_{4},{x}_{5},{x}_{6},{x}_{7},{x}_{8},{x}_{9}\}\\ & = & \{-2.3,-1.6,-0.9,-0.2,0.5,1.2,1.9,2.6,3.3\}.\end{eqnarray} \tag{ 2 }$$

For lithium and beryllium nuclei, the numbers of Voyager-1 data points are very limited, and their number of knots are adopted to be seven, as

$$\begin{eqnarray} & & [{x}_{1},{x}_{2},{x}_{3},{x}_{4},{x}_{5},{x}_{6},{x}_{7}]\\ & = & \{-1.6,-0.83,-0.06,0.71,1.48,2.25,3.3\}.\end{eqnarray} \tag{ 3 }$$

We also check the results through adding or reducing the number of knots, and we find that the results change only slightly in the energy region where no data are available. In the following, y_i parameters at the above fixed x_i knot positions are assumed to be free and are derived through fitting to the data.

The GCR data from AMS-02 (Aguilar et al. 2017, 2018), Voyager-1 (Cummings et al. 2016), and ACE-CRIS⁵ are adopted. For AMS-02 and Voyager-1, the data for helium, lithium, beryllium, boron, carbon, and oxygen nuclei are available, while for ACE-CRIS only the boron, carbon, and oxygen data are available. The AMS-02 data were taken between 2011 May 19 and 2016 May 26. We extract the ACE-CRIS data of the same period from the ACE Science center to derive the LIS spectra. The ACE data of the whole 20 years of operation are then used to study the solar modulation. The uncertainties of the ACE data are the quadratical sum of the statistical ones and the systematic ones, with the latter mainly coming from the geometry factor (2%), the scintillating optical fiber trajectory efficiency (2%), and the spallation correction (1% ∼ 5%) (George et al. 2009). Note that the proton spectra by AMS-02 (Aguilar et al. 2015) and Voyager-1 (Cummings et al. 2016) are not included in this work. This is because the data-taking time for protons of AMS-02 is different from the other nuclei, which may complicate the solar modulation modeling when fitting the LIS spectra. Furthermore, protons are less relevant in the study of GCR propagation compared with the primary and secondary nuclei discussed in this work.

We fit the normalizations of the n spline knots together with the solar modulation potential ϕ. The χ² statistics is defined as

$$\begin{eqnarray}&&{\chi }^{2}=\displaystyle \sum _{i=1}^{m}\displaystyle \frac{{[J({E}_{i};{\boldsymbol{y}},\phi )-{J}_{i}({E}_{i})]}^{2}}{{{\sigma }_{i}}^{2}},\end{eqnarray} \tag{ 4 }$$

where $J({E}_{i};{\boldsymbol{y}},\phi )$ is the expected flux and J_i(E_i) and σ_i are the measured flux and error for the ith data bin with central energy E_i.

We use the Markov Chain Monte Carlo (MCMC) algorithm to minimize the χ² function, which works in the Bayesian framework. The posterior probability of model parameters ${\boldsymbol{\theta }}$ is given by

$$\begin{eqnarray}&&p({\boldsymbol{\theta }}| \mathrm{data})\propto { \mathcal L }({\boldsymbol{\theta }})p({\boldsymbol{\theta }}),\end{eqnarray} \tag{ 5 }$$

where ${ \mathcal L }({\boldsymbol{\theta }})$ is the likelihood function of parameters ${\boldsymbol{\theta }}$ given the observational data, and $p({\boldsymbol{\theta }})$ is the prior probability of ${\boldsymbol{\theta }}$.

The MCMC driver is adapted from CosmoMC (Lewis & Bridle 2002). We adopt the Metropolis–Hastings algorithm. The basic procedure of this algorithm is as follows. We start with a random initial point in the parameter space and jump to a new one following the covariance of these parameters. The accepted probability of this new point is defined as $\min [p({{\boldsymbol{\theta }}}_{\mathrm{new}}| \mathrm{data})/p({{\boldsymbol{\theta }}}_{\mathrm{old}}| \mathrm{data}),1]$. If the new point is accepted, then repeat this procedure from this new one. Otherwise, go back to the original point. For more details about the MCMC one can refer to Gamerman (1997).

The solar modulation degenerates with the LIS fluxes. To constrain the solar modulation potential as effectively as possible, we jointly fit the boron, carbon, and oxygen data, for which the low energy measurements from both Voyager and ACE data are available. This fit gives ϕ_BCO = 0.696 ± 0.016 GV. Using this value as a prior, we then fit the helium, lithium, and beryllium data and get ϕ_He = 0.657 ± 0.013 GV, ϕ_Li = 0.692 ± 0.016 GV, ϕ_Be = 0.694 ± 0.016 GV. We find that all these fits give ϕ ∼ 0.7 GV for the average solar modulation potentials between 2011 May 19 and 2016 May 26, except for He, which gives a somehow smaller modulation potential. Table 1 gives the best-fitting χ² values and the modulation potentials (with 1σ uncertainties). The probability distribution functions of ϕ_BCO, ϕ_He, ϕ_Li, and ϕ_Be are shown in Figure 1.

Zoom In Zoom Out Reset image size

The best-fit LIS spectra of all of these nuclei are shown by solid lines in Figure 2. We can see that this non-parametric method reproduces reasonably any broad structures of the energy spectrum, such as the breaks at O(1) and O(100) GeV n^₋₁. We use the fitted results of ϕ_He, ϕ_Li, ϕ_Be, and ϕ_BCO to de-modulate the TOA measurements by ACE and AMS-02 to obtain the corresponding LIS fluxes, as shown by the colored data points in Figure 2. The uncertainties associated with the modulation parameter, obtained using the error propagation, are added quadratically to the original (statistical and systematic) uncertainties of the measurements. For kinetic energies smaller than ∼1 GeV n^₋₁, the uncertainties due to the modulation parameter account for ∼10% of the total ones, which become smaller at higher energies. The results of the LIS fluxes are given in Table 4 in the Appendix.

Zoom In Zoom Out Reset image size

In Figure 3, we compare the fitted 2σ results of the LIS fluxes for the primary group (He, C, O) and the secondary group (Li, Be, B), with proper normalizations. For the primary group, the energy spectra of He, C, and O are similar to each other for energies above ∼1 GeV n^₋₁. The low-energy spectrum of He is different from that of C and O, which is possibly due to different energy-loss rates of them in the interstellar medium. Whether there are differences among the injection spectra of these primary nuclei needs detailed studies within specified propagation models. The spectra of secondary nuclei are in agreement with each other within the uncertainties.

Zoom In Zoom Out Reset image size

Given the LIS fluxes of CRs, we can then obtain the time evolution of the solar modulation potentials using the long-term measurements of ACE. The ACE data in each Bartels rotation period (27 days) from 1997 to 2016 are extracted. Using the LIS spectra of boron, carbon, and oxygen nuclei, we can derive monthly values of the solar modulation potential. A Bayesian approach is adopted to properly take into account the uncertainties of the LIS spectra. The posterior probability of ϕ is given by

$$\begin{eqnarray}&&p(\phi )\propto \int {{ \mathcal L }}_{{ACE}}(\phi ,{\boldsymbol{y}})p({\boldsymbol{y}})d{\boldsymbol{y}},\end{eqnarray} \tag{ 6 }$$

where ${{ \mathcal L }}_{{ACE}}\propto \exp (-{\chi }_{{ACE}}^{2}/2)$ is the likelihood of model parameters ($\phi ,{\boldsymbol{y}}$) and $p({\boldsymbol{y}})$ is the prior probability distribution of ${\boldsymbol{y}}$ that is obtained through the fit in Section 3.1. The above integration is simply calculated through adding the parameter sets of the last 50% of the Markov chains together, weighted by their stopping numbers at each point.

The posterior mean values (solid lines) and the associated 1σ and 2σ bands (thick and thin shaded regions) of ϕ for the fittings to boron, carbon, and oxygen nuclei individually and simultaneously are shown in Figure 4. We find that the carbon and oxygen data give very close results to the modulation potential, while the boron data give slightly larger results. Since the fluxes of boron are lower than that of carbon and oxygen, the corresponding uncertainties of ϕ derived from the boron data are also larger. Within the uncertainties, these results are consistent with each other. We tabulate the 27 day time series of ϕ_BCO and the associated lower and upper limits in Table 3 in the Appendix. We show the change of TOA fluxes of boron, carbon, and oxygen with time, accompanied by the time series of ϕ_BCO in Figure 7 in the Appendix.

Zoom In Zoom Out Reset image size

Figure 5 compares our results (red curve and associated 68% and 95% bands) of the modulation potential for the joint fit with previous results. The gray line and shaded band show the monthly results⁶ from neutron monitors given in Ghelfi et al. (2017), and the yellow line represents the results⁷ also derived from the neutron monitor data given in Usoskin et al. (2011). Other data points are derived from the studies of various GCR data (Corti et al. 2016; Ghelfi et al. 2016, 2017). The results from different analyses show rough consistency with each other. Quantitatively, they may differ by as large as 50%, in particular for the periods of solar maximum around 2001 or minimum around 2010. The difference may come from different energy ranges of relevant data sets and/or assumptions of the LIS spectra of GCRs adopted in different works. One improvement in our work is the use of the Voyager-1 data taken outside of the solar system to constrain the LIS spectra of GCR nuclei, which makes our LIS spectra less uncertain compared with most of the previous studies.

Zoom In Zoom Out Reset image size

Another indicator of solar activities is the sunspot number. Observational evidence shows a strong correlation between the sunspot numbers and solar activities. Figure 6 shows the relationship between the solar modulation potential ϕ_BCO(t) obtained in this work and the sunspot numbers⁸ at time t − t₀, where t₀ represents a time delay from the solar activity to the modulation of GCRs. We assume a linear correlation between them, as

$$\begin{eqnarray}&&\phi (t)={\phi }_{1}+{\phi }_{2}\times \displaystyle \frac{N(t-{t}_{0})}{100},\end{eqnarray} \tag{ 7 }$$

which is shown by the solid line in Figure 6. The fitting parameters are given in Table 2. Note that the sunspot numbers fluctuate significantly, and thus the uncertainties of the parameters are statistically meaningless. The fit gives a time delay of ∼0.9 years, which can be understood as the time for solar winds traveling across the solar system (∼100 au) with a typical speed of ∼500 km s⁻¹ (Yuan et al. 2017). The results without time delay are also shown in Figure 6 for comparison. We can see that the scattering of data points are clearly larger in the case of no time delay. A similar time delay was also found in previous works (e.g., Kuznetsov et al. 2017; Tomassetti et al. 2017). Tomassetti et al. (2017) found a time delay of 0.68 ± 0.10 years, which is consistent with ours within a 2σ level. Different time delays in the even (∼0.5 years) and odd (∼1.3 years) solar cycles were suggested in Kuznetsov et al. (2017), whose average is fairly consistent with our result.

Zoom In Zoom Out Reset image size