Time dependent solar modulation of cosmic rays from solar minimum to solar maximum

TL;DR

This paper models the time-dependent solar modulation of cosmic rays across solar activity cycles, deriving local interstellar spectra and analyzing the effects of heliospheric magnetic field variations.

Contribution

It introduces a two-dimensional modulation model that accounts for temporal changes in the interplanetary environment and derives the LIS during different solar activity periods.

Findings

Rigidity dependence of the diffusion coefficient is key during polarity reversal periods.

A power law relates the diffusion coefficient to the heliospheric magnetic field strength.

The model successfully fits experimental data from Voyager 1, PAMELA, BESS-POLARII, and ACE.

Abstract

We study the time-dependent modulation effect and derive the local interstellar spectra (LIS) for the cosmic ray (CR) proton, helium, boron and carbon. A two-dimensional modulation model including the variation of the interplanetary environment with time is adopted to describe modulation process. The propagation equation of CRs in the heliosphere is numerically solved by the package Solarprop. We derive the LIS by fitting the latest results of several experiments, including Voyager 1, PAMELA, BESS-POLARII and ACE, during low solar activity periods. We further study the modulation in the polarity reversal periods with the PAMELA proton data. We find that the rigidity dependence of the diffusion coefficient is critical to explain the modulation effect during reversal periods. Our results also indicate a power law relation between the diffusion coefficient and the magnitude of the heliospheric magnetic field (HMF) at the Earth.

Figures

Tables

Table 1. Table 1: The parameterization of proton LIS with cubic spline interpolation.

$log\left(E_k/\mathrm{GeV}\right)$	-2.42	-1.41	-0.50	0	0.50	1.00	1.50	2.00
$log\left(I/\left({\mathrm{GeVm}}^2\mathrm{sr}s\right)\right)$	4.3003	4.4676	4.0396	3.4675	2.5702	1.4101	0.0685	-1.3465

Table 2. Table 2: The parameterization of LIS spectrum for helium, boron and carbon with cubic spline interpolation method.

	$log\left(R/\mathrm{GV}\right)$	-1.0	-0.50	0	0.50	1.00	1.50	2.00
He	$log\left(I/\left({\mathrm{GeVm}}^{2}s\mathrm{sr}\right)\right)$	2.4742	2.7412	2.6561	1.7645	0.4564	-0.8996	-2.2882
B	$log\left(I/\left({\mathrm{GeVm}}^{2}s\mathrm{sr}\right)\right)$	-1.5451	-0.6089	-0.3095	-0.7162	-3.5311	-5.0564	-6.6083
C	$log\left(I/\left({\mathrm{GeVm}}^{2}s\mathrm{sr}\right)\right)$	-0.7080	0.2617	0.4494	-0.3155	-1.5787	-2.9170	-4.2698

Equations

$$V_{sw}=\{\begin{array}{ll}{V}_{max},& \theta \le 30°or\theta \ge 150°.\\ V_{min}(1+|\cos \theta |),& 30°<\theta <150°.\end{array}$$

$$\vec{B}=A\frac{B_0}{r^2}({\vec{e}}_r-\tan \psi {\vec{e}}_{\varphi })(1-2H(\theta -{\theta }_{ns})),$$

$$\frac{\partial f}{\partial t}=-({\vec{V}}_{sw}+{\vec{V}}_{drift})\cdot \nabla f+\nabla \cdot [K^s\cdot \nabla f]+\frac{\nabla \cdot {\vec{V}}_{sw}}{3}\frac{\partial f}{\partial \ln p},$$

$$K_{\parallel }=$$

$${\vec{V}}_{ns}^w=$$

Full text

1 Introduction

After accelerated in sources, CRs are injected and propagated in the Galactic interstellar space. When entering the heliosphere, the intensities of CRs at low energies are significantly affected by several local effects, such as the interactions with the outward solar wind with an embedded magnetic field [1]. Therefore, the observed CR spectra are modulated with the solar activity cycle, and are different from those outside the heliosphere, namely the local interstellar spectra. The study of the solar modulation is essential, as it is indispensable for reproducing the low energy ($<30\,\mathrm{GeV}$) LIS and can also help us to understand the physical process in the CR-heliosphere interaction. The LIS are critical to determine the injection information and the propagation model of Galactic CRs [2, 3, 4, 5, 6, 7, 8], which are also closely related to some new physics studies, such as the indirect detection of dark matter [9, 10, 11, 12, 13, 14, 15].

The current data of CR experiments provide unprecedentedly good opportunities for the research of solar modulation. The Voyager 1, which has crossed the boundary of heliosphere (heliopause) since August 2012, can give direct measurements of CR LIS from a few to hundreds MeV/nucleon [16, 17, 18]. Among the experiments, the PAMELA experiment [19, 20] is particularly compelling for the study of solar modulation 111Investigating the solar modulation with AMS-02 data is in preparation.. The PAMELA collaboration has published eight years of CR data (2006/07-2014/02), continuously recording the variation of the CR proton spectrum from the late declining phase of solar cycle 23 to the maximum phase of solar cycle 24 [21, 22]. Moreover, PAMELA performs precise measurements of the proton spectrum in a wide energy range of $80\,\mathrm{MeV}$ $-$ $50\,\mathrm{GeV}$, partly overlapping with the energy range of Voyager 1. Then combining the results of PAMELA and Voyager 1, we can give good constraints to the model of solar modulation and obtain a reasonable CR proton LIS.

The most widely used model of solar modulation is the force field approximation [23]. It is oversimplified to deal with all the current precise data. The more detailed models have been developed to interpret the observations [24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37]. In this work, we consider a time-dependent 2D modulation model including the diffusion, convection, drift and adiabatic energy loss processes, to study the modulation effect over different solar activity periods and derive the CR LIS. Based on the public code Solarprop 222http://www.th.physik.uni-bonn.de/nilles/people/kappl/ [38], we numerically solve the propagation equation of CRs in the heliosphere with stochastic differential equation approach. The typical parameters related to the interplanetary medium environment, such as the magnitude of the heliospheric magnetic field (HMF), the solar wind speed, and the tilt angle of the heliospheric current sheet (HCS), are taken from the observations. The scale factor of the diffusion coefficient is set to be a time-dependent free parameter to accommodate the observations.

In order to obtain the proton LIS in the form of cubic spline interpolation, we simultaneously fit the Voyager 1 data and the PAMELA data [18, 21, 39] in some low solar activity periods. Then we attempt to explain the PAMELA results in the polarity reversal periods, which is more challenging compared with the case of less active periods. Besides, ACE and BESS experiments provide time-dependent spectra for CR nuclei [40, 41, 42, 43], and Voyager 1 also measures the low energy LIS of CR nuclei [18]. We then derive the LIS of helium, boron and carbon in the same framework with the study of proton.

The paper is organized as follows. In Section  2, we present the important ingredients describing the heliosphere environment. In Section  3, we briefly describe the dominant physical mechanisms in the solar modulation effect. In Section  4, we drive the proton LIS using the PAMELA data and compare the calculated spectra with the PAMELA results over different solar activity periods. We also give an empirical relation between diffusion coefficient and heliospheric magnetic field strength at the Earth. In Section  5 we present the LIS of helium, boron and carbon. Finally, we give the summary in Section  6.

2 The global characteristics of the heliospheric environment

The CR propagation in the heliosphere is affected by the solar activity. The main factors that affect the CR transportation, namely solar wind speed, the magnitude and orientation of heliospheric magnetic field and the inclination between the heliospheric current sheet and the equatorial plane named tilt angle, are all correlated with the solar activity.

The solar wind speed is variable in both latitude and time. During the solar minimum, the typical solar wind speed at the low latitude is about $400\,\mathrm{km}/\mathrm{s}$, while it increases by almost a factor of two at the high latitude [44]. With the increase of solar activity, the boundary of slow and fast solar wind rises rapidly [45]. The solar wind speed is described as [46]

[TABLE]

where $V_{max}$ is taken to be $760\,\mathrm{km}/\mathrm{s}$ and $V_{min}$ is the observation value near the Earth. For simplicity, we take the latitude average value as an approximation in this work.

The solar wind carries the Sun’s magnetic field into the interplanetary space and forms the heliospheric magetic field (HMF) with an Archimedean spiral structure given by [47]:

[TABLE]

where $A=+$ ($A=-$) indicates the magnetic field polarity for the solar magnetic field lines pointing outward (inward) in the northern hemisphere and inward (outward) in the southern hemisphere, $r$ is the distance from the Sun, $B_{0}\sqrt{(1+\tan^{2}\psi(r=1\mathrm{AU},\theta=\pi/2))}$ is the magnetic filed magnitude at the Earth, $H$ is Heaviside function, $\theta$ is the polar angle and $V_{sw}$ is the solar wind speed. The spiral angle $\psi$ is defined as $\tan{\psi}=\frac{\Omega r\sin\theta}{V_{sw}}$, where $\Omega=2.866\cdot 10^{-6}~{}\rm{rad}/\rm{s}$ is the rotation speed of the Sun. $\theta_{ns}$ determines the position of heliospheric current sheet (HCS), which divides the heliosphere into two regions with opposite polarities. The structure of the HCS is parameterized by its tilt angle and is well related to the solar activity. Some modifications of the Parker HMF have been proposed in [48, 49, 50], but remain inconclusive. We adopt the standard Parker HMF model in this work.

In Figure  1 we show the time profile of the parameters characterizing the global heliospheric environment in our model. The first and second panels show the averaged solar wind speed ($V_{sw}$) and HMF magnitude at $1\,\mathrm{AU}$ ($B_{E}$) using the data from OMNI website interface333omniweb.gsfc.nasa.gov for each Carrington rotation (about 27.28 days), respectively. The third panel illustrates the variation of the tilt angle ($\alpha$) taken from the WSO website 444wso.stanford.edu with the “new” model. The estimated periods of changeover of the solar magnetic polarity [51, 52] are represented by the red bands. During the time period considered in this work, the solar wind speed and the magnetic field magnitude near the Earth vary in the ranges of 328–610 $\mathrm{km}/\mathrm{s}$ and 3.1–9.1 $\mathrm{nT}$, respectively. The tilt angle have a obvious variation from $4.5^{\circ}$ to $74.5^{\circ}$. Both the HMF strength and tilt angle show a clear nearly 11-year cycle.

As shown in Figure  1, the description for the heliosphere environment with fixed parameters is not realistic. In the original Solarprop, not only the magnetic field strength but also the tilt angle change with time. We further extend the code to allow the variation of solar wind speed. In this work, we study the solar modulation effect including a smooth time correlation to the solar activity. As a valid approximation, the input parameters are averaged in 7–12 months, which denote the time scale of the solar wind propagation from the Sun to the modulation boundary.

3 The CR propagation in the heliosphere

The CR propagation within heliosphere is dominantly affected by four effects, including the diffusion resulted from scattering magnetic irregularities, the convection caused by outward solar wind, the drift induced by the irregularity of the global heliosphere magnetic field, and the adiabatic energy loss [53]. The modulation effect in the heliosheath is neglected here. For the discussion of this possible effect, we refer the reader to [54, 55, 56, 57]. The review of the solar modulation can be found in [58, 53].

The CR propagation within heliosphere can be described by the Parker transport equation [1]:

[TABLE]

where $f(\vec{r},p,t)$ is the omni-directional distribution function, $\vec{r}$ is the position in a heliocentric spherical coordinate system, $p$ is the particle momentum, $\vec{V}_{sw}$ is the solar wind speed, $\vec{V}_{drift}$ is the drift speed, and $K^{s}$ is the symmetric part of the diffusion tensor. The differential intensity related to the distribution function is given by $I=p^{2}f$.

It is far from completely understanding the parameterization of the diffusion tensor [59, 60, 30, 61]. In this work, a simple spatial and rigidity dependence of the parallel diffusion coefficient $K_{\parallel}$ is adopted as the following form [38]

[TABLE]

where $k\equiv k_{0}\cdot 3.6\times 10^{22}\mathrm{cm}^{2}/\mathrm{s}$ is a scale factor, which describes the time dependence of the diffusion coefficient and reflects the variability of interplanetary medium properties, $B_{E}$ is the strength of HMF at the Earth, and $B$ is the strength of HMF at the particle position. This rigidity dependence is suggested for modulation during the solar minimum in [62] based on the quasilinear theory [63]. For particle rigidity above a threshold value the linear rigidity dependence of diffusion coefficient is commonly adopted in many works [23, 64, 65, 66, 67, 68, 69]. The diffusion coefficient perpendicular to the large scale HMF $K_{\perp}$ is taken to be $K_{\perp}=0.02K_{\parallel}$ according to the test particle simulation [70]. This form of $K_{\perp}$ is widely adopted in literatures [71, 31].

The drift effect leads to a charge-sign dependence and a 22-year cycle in the solar modulation effect [72, 73, 74]. The drift speed from the gradient and curvature of the HMF is written as $V_{gc}=q\frac{\beta R}{3}\,\nabla\times\frac{\vec{B}}{B^{2}}$ [73]. The HCS drift is caused by the change of the field direction at the crossing of the HCS. In this study, we describe the HCS drift following [75], where a thick, symmetric transition region determined by the tilt angle is used to simulate a wavy neutral sheet. Combining the $V_{gc}$ and effective wavy neutral sheet drift speed $V_{ns}^{w}$ together show that the drift speed is divergence-free in the region of $\pi/2-\alpha-\theta_{\triangle}<\theta<\pi/2+\alpha+\theta_{\triangle}$, where $\delta_{\triangle}\approx\frac{2RV_{sw}}{AB_{0}\Omega\cos{\alpha}}$. The $V_{ns}^{w}$ is given by

[TABLE]

where $q$ is the charge sign and $v$ is the particle speed. The product of $qA$ determines the drift direction. During the $A<0$ cycle the positive charge particles drift inwards mainly through the HCS near the equatorial regions. Otherwise, during the $A>0$ cycle, they mainly drift inwards from polar regions.

In this work, we use the public monte carlo code Solarprop [38] to numerically solve the transport equation. The computation is based on the equivalence of a set of stochastic differential equations to the Parker equation. For the details of the numerical method, we refer the reader to [76, 77, 78, 38]. We take the termination shock as the modulation boundary and assume that it is $100\,\mathrm{AU}$ from the Sun. In our default calculation, the only free parameter is the time dependent scale factor of the diffusion coefficient $k_{0}$. Other input parameters are obtained from observations, such as the solar wind speed, the magnitude of magnetic field and the HCS tilt angle.

4 The solar modulation for proton

The PAMELA experiment performed a systematic measurement of the CR proton spectrum in the period 2006-2014 from the late declining phase of solar cycle 23 to the maximum of cycle 24. The detailed comparison between the calculated energy spectra and the observations for different solar activity levels can improve our understanding about the modulation process.

4.1 The local interstellar spectrum for proton

The Voyager 1 crossed the heliopause on the August 2012 and provided the CR LIS at very low energies. Additionally, the CR spectra measured by PAMELA and AMS-02 with rigidity above a few tens of $\mathrm{GV}$ are not affected by the modulation. However, until now, there have been no experiments to measure the LIS in the gap. To derive the LIS we parameterize the LIS with the cubic spline interpolation method [79, 80]. We obtain the proton LIS by fitting the calculated spectra to the Voyager 1 observation and a series of PAMELA data. The GNU Scientific Library (GSL)555https://www.gnu.org/software/gsl/ is used to perform the least-sqaures fitting. In order to avoid the influence of the polarity reversal occurred in the late of 2012, the PAMELA sample data are chosen between Jul. 2006 to Feb. 2012. For considerable saving in computing time we use 12 sets of data to construct the LIS. The corresponding Carrington rotation number of the data used in the fit are 2045, 2052, 2058, 2064, 2070, 2076, 2082, 2088, 2093, 2107, 2114, 2121.

The knots and the corresponding proton intensities are listed in Table  1. Figure  2 shows the obtained proton LIS and the experimental results, including the Voyager 1, PAMELA, and AMS-02 data [18, 81, 82]. The LIS agrees with the proton flux measured by Voyager 1 outside the heliosphere at low energies below $300\,\mathrm{MeV}$, and is consistent with the data measured by PAMELA and AMS-02 at high energies above a few $10\,\mathrm{GeV}$. The LIS shows that it is less affected by the solar modulation for proton above $10\,\mathrm{GeV}$.

4.2 Comparing the calculated proton spectra with the observations before the polarity reversal

After deriving the proton LIS, we can calculate the modulated proton spectra in different periods, and compare them to the spectra observed by PAMELA. We firstly focus on the modulation before the polarity reversal. The northern and southern polar fields reversed in 2012 November and 2014 March, respectively [52]. In every period, the diffusion coefficient is adjusted to reproduce the observed PAMELA spectrum in the range of 0.08–40 $\mathrm{GeV}$.

In Figure  3, we show that our calculated proton spectra are consistent with the corresponding PAMELA data over a wide energy range in 8 different periods. The time evolution of proton flux is close related to the solar activity. The proton flux gradually increased from 2006 to 2009 until reaching the maximum value during 2009 December, and then decreased from 2010 to 2012. We also show the time profile of reduced-$\chi^{2}$ ($\chi^{2}/{d.o.f.}$) in the fit to the PAMELA data on a solar rotation period basis in Figure  4. We find that the reduced-$\chi^{2}$ in all the fit is much smaller than 1. It indicates that our calculated spectra agree remarkably well with observations. Figure  4 also indicates that $\chi^{2}$ after 2011 has larger values than those before 2011.

4.3 Modulation in the polarity reversal period

Although the calculated proton spectra show a good agreement with the PAMELA observations in the period 2006/07–2012/10, there are large discrepancies in the subsequent polarity reversal period. The polarity reversal often occurs near the solar maximum.

It is a challenge to model the modulation effect in polarity reversal period. The gradual reversal process and the frequent solar events disturb the interplanetary medium, therefore the magnetic filed structure becomes more complex in this period. The diffusion and drift coefficients related to the Parker magnetic field model might not be appropriate during the polarity reversal period. To account for these fact, we modify the diffusion coefficient by introducing a power law rigidity dependence. The diffusion coefficient is described as $k_{\perp}\propto k_{\parallel}\propto R^{\delta}$ (DC $\propto R^{\delta}$) for the particle rigidity above $0.1\,\mathrm{GV}$. We take $\delta$ as a free parameter rather than a constant 1 in our default case. Some studies argue that the drift effect vanishes during the solar maximum [83, 84]. We attempt to turn off the drift effect in our model.

We attempt some assumptions for the modulation effect in the polarity reversal periods. There are two assumptions for the rigidity dependence of the diffusion coefficient: $DC\propto R$ and $DC\propto R^{\delta}$. For the drift effect, we consider three cases: the polarity is positive ($A>0$), the polarity is negative ($A<0$) and no drift effect. Thus there are total six assumptions for the diffusion coefficient and drift effect. The force-field approximation model is also included in our comparison. In Figure 5 we show the resulting $\chi^{2}$ from the fits to the PAMELA proton data for different scenarios. We find that the assumption with a $DC\propto R^{\delta}$ and without drift provides the best fit. It is evident that adopting a variable power law rigidity dependence of the diffusion coefficient can significantly reduce the $\chi^{2}$.

The time profile of $\delta$ and reduced-$\chi^{2}$ for the best fit are shown in Figure 6. As the diffusion coefficient is related to the magnetic field power spectrum, the property of the HMF turbulence during polarity reversal should be different in the quiet epoches.

4.4 An empirical relation between the diffusion coefficient and the magnitude of HMF at the Earth

It is well known that the diffusion coefficient is anti-correlated with HMF strength. In this work, the time dependence of the diffusion coefficients is described by a scale factor $k_{0}(t)$. The time variations in magnetic field strength are introduced by averaging $B_{E}$ over the time taken the solar wind to reach the modulation boundary from the Sun. Figure 7 shows the time profile of $k_{0}(t)$ and backward time average of HMF at the Earth $\langle B_{E}\rangle$.

We follow the previous works [85, 86] and assume $k_{0}=\left(\frac{B_{c}}{\langle B_{E}\rangle}\right)^{n}$, where $B_{c}$ and $n$ are free parameters. Figure 8 shows the relation between $k_{0}$ and ${\langle B_{E}\rangle}$ during 2006/07 to 2013/02. For the period of 2013/06–2014/02, the correlation between $k_{0}$ and ${\langle B_{E}\rangle}$ is weak (see Figure 7), so we do not take into account this period. Obviously there is a discrepancy between $B_{c}$ in the declining phase (2006/07–2010/03) and increasing phase (2010/10–2013/02) of observed cosmic ray intensity level, while we find that the power $n$ approximates 2 and slightly varies with time. This result is consistent with the conclusion in [87]. From the empirical relation found in the fit, we can use the estimated magnetic field strength to obtain the diffusion coefficient and get predictions for the modulated spectra.

5 LIS for helium, carbon and boron

In order to derive the LIS for other nuclei, we investigate the modulation effect for the CR helium, boron, and carbon. Here, we also use the cubic spline interpolation method to construct the LIS, and derive the helium LIS by minimizing the weighted differences between the calculated spectra and the observations of Voyager 1 and BESS-POLARII. The same method is used but the combination of data from the Voyager 1, PAMELA and ACE boron (carbon) observation in 2009 to obtain the boron (carbon) LIS. The parameters of the LIS are summarized in Table  2.

In Figure  9 we show the calculated CR helium spectra and compare them with the AMS-01, BESS98, BESS-POLARI, and BESS-POLARII results [88, 42, 43]. The LIS of helium and the Voyager 1 data are also shown. In order to reproduce the spectra observed by AMS-01, a factor 1.25 is adopted to scale down the LIS. Figure  10 and Figure  11 show the LIS and calculated spectra for boron and carbon, respectively. For comparison, the ACE and PAMELA results [40, 41, 89] are also shown. It is shown that our calculated spectra can well reproduce the observations. Note that in our calculation, the CR spectra for different nuclei in the same period share modulation parameters. This is an important improvement compared with the force field approximation model, where the potential parameter should be specified for each CR specie [41].

6 conclusion

In this work, we study the solar modulation of CRs and derive the new CR LIS with a time-dependent modulation model and the latest CR experiments. The parameter describing the global characteristics of heliospheric environment, such as the solar wind speed,magnetic field magnitude and tilt angle, are all obtained from observations. In out default calculation, the only free parameter is the scale factor of the diffusion coefficient.

We adopt the long-term PAMELA observation to derive the LIS and investigate the modulation effect of the CR proton. We utilize some data samples of PAMELA during the low solar activity periods and the result of Voyager 1 to obtain the proton LIS, then all the PAMELA data before the polarity reversal can be well reproduced. Modeling the modulation effect during the polarity reversal period is challenging, since the theory of CR propagation in the heliosphere during this period is poorly understood. The complex magnetic field configuration increases the uncertainties of diffusion and drift effects. In order to reproduce the observations, we change the linear relation between diffusion coefficient and rigidity to $DC\propto R^{\delta}$ and assume there is no significant drift effect in the polarity reversal period. We find that the diffusion coefficient is anti-correlated to the HMF strength at the Earth. An empirical relation can be described as $DC\propto{\langle B_{E}\rangle}^{-n}$, where $n$ is $\sim 2$ and slightly varies with time.

We also study the modulation effect and derive the LIS for the CR helium, boron and carbon. In the calculation, the parameters in the modulation model are taken to be same for different CR species in the same period. Since the calculated spectra can well explain several experimental results, our approach is a good description for dealing with the modulation effect. Using the LIS derived here, uncertainties in the study of the CR propagation in the Galaxy can be reduced.

Acknowledgement

This work is supported by the National Key R&D Program of China (No. 2016YFA0400200), the National Natural Science Foundation of China (Nos. U1738209 and 11851303).

Appendix A Comparison with the Force-Field approximation

In order to compare our results with the Force-Field approximation, we show the difference of $\chi^{2}$ from the fits to the long-term PAMELA proton observations in two scenarios in Figure 12. We can see that in most periods the Force-Field approximation give a larger $\chi^{2}$ than our results.