Modulation of Galactic Cosmic Rays in the Inner Heliosphere, Comparing with PAMELA Measurements

TL;DR

This paper presents a numerical model for the time-dependent modulation of galactic cosmic rays in the inner heliosphere, incorporating a modified magnetic field, a new diffusion coefficient, and observationally-based turbulence dependencies, successfully matching PAMELA data.

Contribution

The study introduces a comprehensive numerical model with a modified Parker magnetic field and a novel diffusion coefficient, improving the understanding of cosmic ray modulation in the heliosphere.

Findings

Model reproduces PAMELA proton spectra during solar minimum.

Analytical turbulence variation matches spacecraft observations.

Radial dependence of turbulence aligns with empirical data.

Abstract

We develop a numerical model to study the time-dependent modulation of galactic cosmic rays (GCRs) in the inner heliosphere. In the model a time-delayed modified Parker heliospheric magnetic field (HMF) and a new diffusion coefficient model, NLGCE-F, from Qin \& Zhang (2014), are adopted. In addition, the latitudinal dependence of magnetic turbulence magnitude is assumed as $\sim (1+\sin^2\theta)/2$ from the observations of Ulysses, and the radial dependence is assumed as $\sim r^S$, where we choose an expression of $S$ as a function of the heliospheric current sheet (HCS) tilt angle. We show that the analytical expression used to describe the spatial variation of HMF turbulence magnitude agrees well with the Ulysses, Voyager 1, and Voyager 2 observations. By numerically calculating the modulation code we get the proton energy spectra as a function of time during the recent solar minimum, it is shown that the modulation results are consistent with the PAMELA measurements.

Equations

$$\frac{\partial f}{\partial t}=-\left({\mathbf{V}}_{sw}+⟨{\mathbf{v}}_d⟩\right)\cdot \nabla f+\nabla \cdot \left({\mathbf{K}}_s\cdot \nabla f\right)+\frac{1}{3}\left(\nabla \cdot {\mathbf{V}}_{sw}\right)\frac{\partial f}{\partial \ln p},$$

$$\delta B^2=\frac{1}{n}\sum _{i=1}^n{\left(B_i-\overline{B}\right)}^2,$$

$$\overline{B}=\frac{1}{n}\sum _{i=1}^n{B}_i.$$

$$\mathbf{B}=\frac{A{B}_0}{r^2}\left({\mathbf{e}}_r-\frac{(r-r_s)\Omega \sin \theta }{V_{sw}}{\mathbf{e}}_{\varphi }\right)\left[1-2H(\theta -{\theta }^{\prime })\right],$$

$$\mathbf{B}=\frac{A{B}_0}{r^2}\left({\mathbf{e}}_r+\frac{r\delta (\theta )}{r_s}{\mathbf{e}}_{\theta }-\frac{(r-r_s)\Omega \sin \theta }{V_{sw}}{\mathbf{e}}_{\varphi }\right)[1-2H(\theta -{\theta }^{\prime })],$$

$$\delta (\theta )=\frac{\delta {}_m}{\sin \theta },$$

$${\mathbf{V}}_{sw}(r,\theta )=V_0\left\{1-exp\left[\frac{40}{3}\left(\frac{r_s-r}{r_0}\right)\right]\right\}\left\{1.475\mp 0.4\tanh \left[6.8(\theta -\frac{\pi }{2}\pm \xi )\right]\right\}{\mathbf{e}}_r,$$

$$\left\{\begin{array}[]{l}v_{i}(t-t_{i})\geq r-r_{0}\\ v_{i+1}(t-t_{i+1})<r-r_{0}\\ \end{array}\right.$$

$$\left\{\begin{array}[]{l}v_{0}(t-t_{i})\geq r-r_{0}\\ v_{0}(t-t_{i+1})<r-r_{0}.\\ \end{array}\right.$$

$$⟨{\mathbf{v}}_d⟩=\nabla \times \left({\kappa }_A\frac{\mathbf{B}}{B}\right),$$

$${\kappa }_A=q\frac{P\beta }{3B},$$

$$⟨{\mathbf{v}}_d⟩$$

$$\mathbf{f}$$

$$\eta$$

$$\Gamma$$

$$v_{gc,r}$$

$$v_{gc,\theta }$$

$$v_{gc,\varphi }$$

$$v_{gc,0}=qA\frac{2P\beta r[1-2H(\theta -{\theta }^{\prime })]}{3B_0(1+{\eta }^2+{\Gamma }^2{)}^2}.$$

$${\theta }^{\prime }=\frac{\pi }{2}-\arctan \left[\tan \alpha \sin \left(\varphi +\frac{\left(r-r_s\right)\Omega }{V_{sw}}\right)\right],$$

$$v_{ns,0}=vqA\left[0.457-0.412\frac{d}{R_L}+0.0915{\left(\frac{d}{R_L}\right)}^2\right].$$

$$v_{ns,r}$$

$$v_{ns,\theta }$$

$$v_{ns,\varphi }$$

$${\varphi }^{\prime }=\varphi +\frac{\left(r-r_s\right)\Omega }{V_{sw}},$$

$$\rho =\sqrt{x_1^2+x_2^2+x_3^2}$$

$$x_1$$

$$x_2$$

$$x_3$$

$$\delta B\sim r^S$$

$$\delta B\sim \frac{1+{\sin }^2\theta }{2}.$$

$$\delta B=\delta B_{1\mathrm{AU}}{R}^S\left(\frac{1+{\sin }^2\theta }{2}\right),$$

$$S=-1.56+0.09\ln \frac{\alpha }{{\alpha }_c},$$

$$\ln \frac{{\lambda }_{\sigma }}{{\lambda }_{slab}}=\sum _{i=0}^{n_{\sigma 1}}{a}_i^{\sigma }{\left(\ln \frac{R_L}{{\lambda }_{slab}}\right)}^i$$

Full text

cosmic rays, modulation, turbulence, activity – Sun: heliosphere

1 Introduction

Galactic cosmic rays (GCRs) are modulated by solar wind irregularities while transporting inside the heliosphere. The physical mechanism of cosmic ray transport in the heliosphere is described by a well known Parker transport equation (TPE) (Parker, 1965),

[TABLE]

where $f(\mathbf{r},p,t)$ is the omni-directional cosmic ray distribution function, with $\mathbf{r}$ the position, $p$ the particle momentum and $t$ the time. The distribution function $f(\mathbf{r},p,t)$ is related to the differential intensity $j$ with respect to kinetic energy by $j=p^{2}f$. The terms on the right hand side include all the relevant transport processes, i.e., the outward convection by the solar wind velocity $\mathbf{V}_{sw}$, the pitch-angle averaged drift velocity $\left\langle{\mathbf{v}_{d}}\right\rangle$ caused by irregularity in the global heliospheric magnetic field, $\mathbf{K}_{s}$ is the symmetric part of the diffusion tensor which is diagnal in HMF-aligned coordinate system, and the adiabatic energy loss, which have been successfully illustrated by theoretical and numerical models (see, e.g., Parker, 1965; Zhang, 1999; Pei et al., 2010b; Strauss et al., 2012; Potgieter, 2013; Zhao et al., 2014).

The past solar cycle (Cycle 24) was unusual with a prolonged periodicity, and the solar minimum conditions lasted until early 2010. During the solar minimum of 2006 to 2009, the averaged sunspot number (SSN) was reported to be the lowest since 1914 (Schrijver et al., 2011). The solar polar field strength during the 2003-2004 solar maximum was especially weaker than previous three solar cycles (Svalgaard et al., 2005). The heliospheric magnetic field (HMF) reached $\sim$3 nT at 1 AU in 2009, the lowest value since 1963 (Ahluwalia et al., 2010). The tilt angle ($\alpha$) of the heliospheric current sheet (HCS) had a flatter decline than previous solar minimums (Ahluwalia & Ygbuhay, 2011). It is also reported that the coronal mass ejection (CME) rates and the solar wind dynamic pressure in 2008-2009 were noticeably lower than that in 1997-1998 (Vourlidas et al., 2010; McComas et al., 2008). These extreme solar minimum conditions resulted in a record high-level GCR intensity measured at Earth (Mewaldt et al., 2010; Ahluwalia & Ygbuhay, 2011). Such a prolonged solar minimum provides a good chance to study the modulation of GCRs inside the heliosphere (e.g., Mewaldt, 2013; Adriani et al., 2013b; Potgieter et al., 2014; Zhao et al., 2014).

The Payload for Antimatter-Matter Exploration and Light-nuclei Astrophysics (PAMELA) satellite experiment was designed to study the charged components of cosmic rays (CRs), among which antiparticles are focused. PAMELA has been taking data since it was launched in June 2006, and has brought plentiful scientific results about the heliosphere (Adriani et al., 2014b; Mori et al., 2015; Galper et al., 2017). PAMELA measures precise cosmic ray proton and helium spectra in the rigidity range from $1$ GV to $1.2$ TV (Adriani et al., 2011b). The energy spectra of protons and helium particles show different spectral indexes above 30 GV and both present a spectral hardening at about $230$ GV. The electron and positron energy spectra are measured up to $600$ GeV and $200$ GeV, respectively (Adriani et al., 2011a, 2013a). PAMELA also measures the flux of Boron and Carbon as well as the boron-to-carbon (B/C) ratio which can be utilized to investigate the cosmic ray propagation processes (Adriani et al., 2014a). Adriani et al. (2013b) presented precise galactic proton energy spectra in the range $0.08-50$ GeV for each Carrington rotation from July 2006 to January 2010. The observed proton spectra became progressively softer since July 2006. Later, Adriani et al. (2015) reported precise electron spectra with a six-month interval at the same time period. Such precise CR energy spectra can be used to study physical processes of CR transport in the heliosphere, including convection, diffusion, drifts, and adiabatic energy changes.

The precise cosmic ray spectra measured by PAMELA instrument provides important information for people to understand the origin and propagation of GCRs. Voyager 1 was assumed to have crossed the heliopause (HP) at $\sim 122$ AU on August 2012 (Webber & McDonald, 2013), cosmic ray spectra in the energy between $5$ MeV to $50$ MeV in the very local interstellar medium was then reported by Stone et al. (2013). The cosmic ray spectra data from PAMELA and Voyager 1 can be used to construct the very local interstellar spectrum (LIS, see e.g., Potgieter et al., 2014; Ngobeni & Potgieter, 2014; Potgieter et al., 2015; Vos & Potgieter, 2015; Bisschoff & Potgieter, 2016). The LIS is important as the input spectrum in the solar modulation study. To understand the modulation processes during this unusual solar minimum, three-dimensional (3D) numerical models have been established to solve the Parker transport equation. Zhao et al. (2014) used an empirical diffusion coefficient model according to Zhang (1999) and incorporated a 3D wavy HCS. By reproducing the proton spectra observed by PAMELA and IMP 8 in previous two solar minima, they concluded that increased parallel diffusion and decreased perpendicular diffusion in polar direction caused by low magnetic turbulence might be the possible mechanism for the high GCR intensity in the past solar minimum, which is in contrast to the assumption of enhanced diffusion in the polar regions that was used to explain the observed Ulysses CR gradients (e.g., Potgieter, 2000). Potgieter et al. (2014) used a different diffusion coefficient model, in which modulation parameters vary as a function of time, in their numerical work. They successfully reproduced the PAMELA proton spectra for four selected periods and reported that this solar minimum was “diffusion dominated”, and that the modulation effects of particle drifts were less obvious but still played a significant role. Similar numerical model was also used to study the modulation of galactic electrons (Potgieter et al., 2015) and the modification effects of Parker HMF over the polar region Raath et al. (2016).

Furthermore, due to the precise proton and electron spectra measured by PAMELA from 2006 to 2009 (Adriani et al., 2013b, 2015), one is able to know the different modulation effects between protons and electrons, which is called the charge-sign-dependent modulation. In this past polarity $A<0$ solar minimum, low-energy protons were more sensitive to the changes of heliospheric conditions than low-energy electrons, such phenomenon can be reproduced only by incorporating drifts in the numerical model (Di Felice et al., 2017). The PAMELA proton flux data were also used to study the spatial gradients in the inner heliosphere together with proton flux data from Ulysses COSPIN/KET. de Simone et al. (2011) and Gieseler & Heber (2016) used an empirical approach to calculate the radial and latitudinal gradients of protons during the past solar minimum. It is shown that the radial gradients are always positive while the latitudinal gradients are always negative as expected but with less magnitude than that predicted by earlier works (Potgieter et al., 2001). Following de Simone et al. (2011) and Gieseler & Heber (2016), Vos & Potgieter (2016) also used a numerical model to compute the spatial gradients of protons from 2006 to 2009. They concluded that although the drift effects were weaker than the predictions from those drift-dominated works due to the suppression by the excess diffusion, they still played an important role due to the significant decrease of HMF magnitude until the end of 2009.

Note that numerical models of modulation are usually solved in steady state, interplanetary conditions have to be determined with the solar activity some time before because of the limit of solar wind speed. Ndiitwani et al. (2013) used a time-dependent two-dimensional (2D) numerical model to study CR modulation using PAMELA proton data in this unusual period. Smoothed monthly HMF and HCS, which were embedded in the solar wind plasma, were used to establish a realistic heliospheric conditions. Based on the work of Manuel et al. (2011) and Potgieter et al. (2014), Ndiitwani et al. (2013) established a time-dependent diffusion coefficient model, in which the yearly time-dependent modulation parameters were obtained from the compound model (Manuel et al., 2011) and the empirical model (Potgieter et al., 2014), respectively. However, Ndiitwani et al. (2013) did not consider the variations of solar wind speed. Recently, Boschini et al. (2017) used a 2D heliospheric modulation (HelMod, e.g., Bobik et al., 2012, 2013) model to study the modulation of GCR during solar cycles 23 and 24. In their model, the heliosphere was divided into polar and equatorial regions, the modified Parker spiral HMF (Jokipii & Kota, 1989) and the Parker spatial HMF (Parker, 1958) were used in polar and equatorial regions, respectively. They also re-scaled the heliosphere into 15 radially equally-spaced slices to relate interplanetary conditions with states near Earth. They used a parameter $K_{0}$ to describe the time dependence of diffusion coefficients. Bobik et al. (2012) discussed the relationship between $K_{0}$ and the modulation strength given by the force-field model (FFM, see e.g., Gleeson & Axford, 1968; Gleeson & Urch, 1971), and Boschini et al. (2017) derived $K_{0}$ using modulation strength data from Usoskin et al. (2011). For periods of low solar activity, it was divided into ascending and descending phases for both negative and positive solar magnetic field polarities, and different polynomial equations were used to describe the relationship between $K_{0}$ and the sunspot numbers. Furthermore, they used neutron monitor counting rate to reproduce the variation of $K_{0}$ during periods of high solar activity. Such model was used to study modulation of GCRs with energy approximately larger than 0.5 GeV/nucleon, and the results were consistent with the observations of PAMELA, AMS-02, and Ulysses.

In this paper, we develop a model of GCR modulation in the inner heliosphere to study the GCRs measurements from PAMELA. The paper is organized as follows: In Section 2 we discuss the GCRs modulation model, including the interplanetary conditions input from observations, heliospheric magnetic field and solar wind speed, the particle drifts, the magnetic turbulence throughout the inner heliosphere, the diffusion coefficients, and the heliospheric boundary, from subsections 2.1 to 2.6. In section 3, we describe the numerical methods. In section 4, we show the numerical modulation results and the comparison with the PAMELA observations in recent solar minimum. Conclusions and discussion are shown in section 5.

2 GCRs modulation model

2.1 Interplanetary conditions input from observations

In order to study the time-dependent modulation of GCRs, we need some spacecraft observations near the Earth. Figure 1 illustrates observations of interplanetary conditions as a function of time which is used in our model. Top panel shows the computed tilt angle $\alpha$ until 2015 for the new model from Wilcox Solar Observatory (wso.stanford.edu). Second and third panels show averaged solar wind velocity $V_{sw}$ and HMF magnitude at $1$ AU using the OMNI data (omniweb.gsfc.nasa.gov) for each Carrington rotation. Based on the assumption of isotropic magnetic turbulence, the total variance $\delta B^{2}$ is calculated over Carrington rotation intervals using hourly averages of HMF magnitude from OMNI. We have $n$ ($\sim 655$ here) samples per Carrington rotation, the variance of the total magnetic field magnitude is

[TABLE]

where

[TABLE]

The square root of $\delta B^{2}$ is shown as black line in the bottom panel of Figure 1. Manuel et al. (2011, 2014) (see also, e.g., Strauss & Potgieter, 2010; Strauss et al., 2011b) also calculated the total magnetic field variance over 1 year intervals, and the results are shown as red circles in the bottom panel. It is shown that our calculation is consistent with the results of Manuel et al. (2011, 2014). In the GCRs modulation model, all input parameters are obtained from observations near the Earth shown in Figure 1.

2.2 Heliospheric magnetic field and solar wind speed

The heliospheric magnetic field (HMF), which plays an important role in the modulation of GCRs, is assumed to have an Archimedean spiral due to the solar rotation according to Parker (1958). The Parker spiral HMF can be written as

[TABLE]

where $B_{0}$ is a constant, $A$ is the polarity of HMF whose positive (negative) value represents the magnetic field points outward (inward) in the northern hemisphere of Sun, $\mathbf{e}_{r}$ and $\mathbf{e}_{\phi}$ are unit vectors in the radial and azimuthal directions, respectively, $r$ is heliocentric radius, $\theta$ is the polar angle, $r_{s}$ is the radius of the source surface where the HMF is assumed to be directed radially outwards and we take $r_{s}=r_{\odot}=0.005\textrm{AU}$ with $r_{\odot}$ being the radius of solar surface (Jokipii & Kota, 1989), $\Omega=2.66\times 10^{-6}$ rad s*-1* is the rotation speed of Sun, $V_{sw}$ is the radial solar wind speed, $\theta^{\prime}$ is the heliospheric current sheet (HCS) latitudinal extent, and $H$ is the Heaviside function.

However, the Parker HMF is an oversimplification and gives a low magnetic field intensity at large radial distance in polar heliosphere, which could lead to the too rapid entry of GCRs in the polar regions, so it is necessary to modify the Parker HMF. Jokipii & Kota (1989) suggested superimposing a perturbation field on the Parker spiral HMF since turbulence near the solar surface resulting in a transverse magnetic field at large radial distance in the polar regions. With this modification the Parker spiral HMF becomes

[TABLE]

where $\delta(\theta)$ is the perturbation parameter. In order to have a divergence free magnetic field the perturbation parameter is written as

[TABLE]

here $\delta_{m}$ indicates the perturbation parameter in the equatorial plane. In addition, we use a reflective boundary condition near the poles to avoid singularity, $\theta=2\theta_{0}-\theta$, for $\theta<\theta_{0}$ if $\theta<90^{\circ}$ or $\theta>\theta_{0}$ if $\theta>90^{\circ}$. In this study we set $\delta{}_{m}=2\times 10^{-5}$ (Bobik et al., 2013; Boschini et al., 2017) , $\theta_{0}=2.5^{\circ}$ if $\theta<90^{\circ}$ and $\theta_{0}=177.5^{\circ}$ if $\theta>90^{\circ}$. This modification makes the field decrease as $r^{-1}$ instead of $r^{-2}$ in the polar regions for large $r$ without changing the magnetic field dramatically in the equatorial plane, and it is supported by the observations of Ulysses (e.g., Balogh et al., 1995; Heber & Potgieter, 2006) and tested by numerical models (e.g., Langner, 2004; Bobik et al., 2012, 2013; Raath et al., 2016).

The solar wind speed has a latitudinal dependence during solar minimum, increasing from $\sim 400$ km s*-1* in the equatorial plane to $\sim 800$ km s*-1* in the high latitudes, but during solar maximum, such simple pattern does not exist anymore (McComas et al., 2002; Heber & Potgieter, 2006; Zurbuchen, 2007). Solar activity can be classified in terms of the HCS tilt angle $\alpha$, with $\alpha\leq 30^{\circ}$, $30^{\circ}<\alpha\leq 60^{\circ}$, and $60^{\circ}<\alpha\leq 90^{\circ}$ representing periods of low, moderate, and high solar activity, respectively (Potgieter et al., 2001, 2013). In addition, the solar wind speed accelerates from zero to a constant within 0.3 AU from the Sun according to Sheeley et al. (1997). In this work, we study GCRs in solar minimum, so following Heber & Potgieter (2006) and Potgieter (2013) we express solar wind speed as

[TABLE]

with $V_{0}=400$ km/s, $r_{0}=1$ AU, and $\xi=\alpha+15\pi/180$. The top and bottom sign correspond to the northern and southern hemisphere, respectively. However, if we study GCRs during periods of moderate and high solar activities the solar wind speed can be set as a constant extracted from OMNI data set (Bobik et al., 2012). Note that for simplicity purpose, in solving the TPE Equation (1) numerically in each step we assume the magnitude of solar wind as a constant with the value calculated with Equation (7).

Figure 2 shows particle’s gyro-radius as a function of rigidity (top panel), polar angle (second panel), and radial distance (bottom panel). Interplanetary conditions at $1$ AU are set as $B=5.05$ nT, and $\alpha=15^{\circ}$. The black solid and red dotted lines indicate results from the Parker field and the modified one, respectively. From the figure we can see that generally the modified field agrees with the Parker field, however, for $1$ GV particles, in the polar regions with large solar radial distance, very weak Parker field makes particles’ gyro-radius very large, but the modified model with enhanced field keeps particles’ gyro-radius around several AU.

It is noted that interplanetary conditions at solar radial distance $r$ are related to the states at the source surface $r_{s}$ at some earlier time because of the solar wind flow (Potgieter et al., 2014, 2015), and the heliosphere is dynamic due to the solar activities. In our numerical model, we divide time in days and assume a locally static heliosphere in each day. In the $i^{th}$ Carrington period $t_{i}$, the observation of solar wind velocity at $1$ AU is $v_{i}$. To calculate the interplanetary conditions in solar distance $r$ at time $t$, we can use the input parameters near the Earth (e.g., $V_{sw}$, $B$, $\delta B$, $\alpha$, $A$) at time $t_{i}$ if Equation (8) is satisfied,

[TABLE]

with $r_{0}=1$ AU. For simplification, we use

[TABLE]

with the typical solar wind speed $v_{0}=0.25$ AU/day. It is noted that in the region between two slices of plasma the HMF is not divergence free, but in each step of the numerical solution of the TPE Equation (1), we always keep inside one slice of plasma.

2.3 Particle drifts

Particle drifts play an important role in the solar modulation of GCRs (Jokipii et al., 1977; Jokipii & Kopriva, 1979; Potgieter, 2013), the pitch angle averaged drift velocity caused by irregularity in the HMF is given by

[TABLE]

with $\kappa_{A}$ the drift coefficient. Under the assumption of weak scattering, the drift coefficient is simply written as

[TABLE]

with q the particle charge sign, $P$ the rigidity of particle and $\beta$ the ratio between the speed of particle and that of light. For the modified Parker HMF given in Equation (5), the drift velocity can be written as (Burger & Potgieter, 1989)

[TABLE]

where

[TABLE]

Here, $\delta_{Dirac}$ is the Dirac’s delta function, $\mathbf{v}_{gc}$ is the combination of gradient and curvature drifts, and $\mathbf{v}_{ns}$ is the current sheet drift.

In the following, we show that charge-sign dependent modulation and a 22-year cycle could be caused by gradient and curvature drifts (Potgieter, 2013). During $A<0$ polarity cycles, protons mainly drift inwards along the HCS in the equatorial regions so their intensity can be reduced by the increasing waviness of HCS, therefore, a sharp peak in the temporal profile of GCR intensity is usually observed. However, during the $A>0$ cycles, protons mainly drift inwards from polar regions, therefore, a flatter peak of GCR intensities are usually observed. This effect reverses for negatively charged GCRs. The radial, latitudinal, and azimuthal components of the gradient and curvature drifts are given by

[TABLE]

respectively, where

[TABLE]

The expression for $\theta^{\prime}$ is given by Kota & Jokipii (1983)

[TABLE]

with $\alpha$ the tilt angle. This formula is valid for large tilt angle conditions (Pei et al., 2012; Raath et al., 2015).

In the current sheet, the current sheet drift velocity given by Equation (12) becomes a Dirac function. The singular current sheet drift velocity is not physical and is not easy to deal with in the numerical method (Zhang, 1999). Therefore, we replace the current sheet drift magnitude with a formula shown in Equation (21) by following Burger & Potgieter (1989). With the assumption of Burger & Potgieter (1989), a particle will experience current sheet drift if its distance $d$ to the HCS is less than two gyro radii $2R_{L}$, and the magnitude of $\mathbf{v}_{ns}$ is given by

[TABLE]

It can be derived directly from Equation (12) that the direction of current sheet drift velocity lies in the HCS and is perpendicular to the HMF, and the radial, latitudinal, and azimuthal components of the current sheet drifts can be written as

[TABLE]

respectively, where

[TABLE]

and

[TABLE]

with

[TABLE]

We should note that Equations (22)-(24) are equal to the results of Burger (2012) and Pei et al. (2012) if we use the Parker HMF (i.e., $\eta=0$). Similar expressions of Equations (16)-(18) and (22)-(24) are given by Raath et al. (2016), but using different expressions for the HMF.

Since the Parker HMF is an oversimplification and gives a low magnetic field intensity at large radial distance, especially at high latitudes, drifts become very large over the polar regions of the heliosphere. It is also known that, with the assumption of weak scattering, the particle’s gyro-radius is equivalent to its drift scale, so that Figure 2 also shows that for Parker field in the polar regions with large solar radial distance, $1$ GV particles’ drift speed becomes very large, but for the modified model the drift speed keeps in the similar level as in other regions.

2.4 Magnetic turbulence throughout the inner-heliosphere

The development of magnetic turbulence transport models (TTMs, see, e.g., Zank et al., 1996, 2012, 2017; Breech et al., 2008; Pei et al., 2010a; Oughton et al., 2011; Engelbrecht & Burger, 2013a; Guo & Florinski, 2016), which allows us to have a better scenario of the heliosphere, plays an important role in the ab initio models for cosmic ray diffusion. Especially in the numerical modulation study, the analytical expressions of diffusion is essential, which are directly based on the spatial dependence of magnetic turbulence.

Some theoretical work has been done to study TTMs. For example, Oughton et al. (2011) developed the two-component TTM. Furthermore, Engelbrecht & Burger (2013a, b) solved the TTM of Oughton et al. (2011) for solar minimum interplanetary conditions and it was shown that the results of turbulence quantities agree well with the observations of Ulysses, and consequently they studied the spatial variations of diffusion coefficients by applying the results of TTM in scattering theory, and the results were used in an ab initio model for cosmic ray modulation. Recently, Zank et al. (2017) studied turbulence quantities with the nearly incompressible magnetohydrodynamics (NI MHD) theory. The NI MHD theory can be used to investigate a broader range of solar wind observations, and its solutions given by Adhikari et al. (2017) enhance our understanding of turbulence quantities in the inner heliosphere. However, the TTM from Engelbrecht & Burger (2013a) and Zank et al. (2017) are complicated to some extent, especially for the study of long-term modulation of GCRs. Therefore, some works use analytical expressions to describe the radial and latitudinal dependence of magnetic turbulence (Zank et al., 1996; Burger et al., 2008; Effenberger et al., 2012; Ngobeni & Potgieter, 2014). The magnetic turbulence magnitude in the heliosphere is assumed decreasing as (e.g., Burger et al., 2008; Effenberger et al., 2012; Guo & Florinski, 2014; Ngobeni & Potgieter, 2014; Strauss et al., 2017)

[TABLE]

where $S$ means the radial dependence of the magnetic turbulence magnitude which can be determined later. Based on the observations of Ulysses instruments (Perri & Balogh, 2010), the variance of magnetic field magnitude at high latitude is smaller than that of low latitude, we can give an expression to describe the latitudinal dependence of the magnetic turbulence magnitude as

[TABLE]

Therefore, the analytical expression of magnetic turbulence magnitude can be written as

[TABLE]

where $\delta B_{1\mathrm{AU}}$ represents the observation of turbulence at Earth and $R=r/r_{0}$, $r_{0}=1$ AU.

According to the observations of Ulysses, $S$ should be varying as a function of time (Smith & Balogh, 2008), so it is assumed that $S$ is closely associated with solar activities. Obviously, there exist different expressions to describe such relationship. However, in this work we use the following expression,

[TABLE]

where $\alpha_{c}=1^{\circ}$. Hereafter we denote the turbulence magnitude model from Equations (32) and (33) as TRST (Turbulence magnitude varying as R, S, and Theta) model. Next, it can be shown that the TRST model agrees well with the turbulence magnitude measurements from Ulysses, Voyager 1, and Voyager 2,

Figure 3 represents the comparison between the results of our TRST model and the observations of Ulysses from September 2006 to May 2008. Top panel shows the trajectory of Ulysses, with black and red lines representing the radial distance and heliographic latitude, respectively. During this time period, Ulysses had a fast latitude scan, the radial distance varies from $1.4$ AU to $3.4$ AU and the latitude varies from $-80^{\circ}$ to $80^{\circ}$. In the bottom panel, black circles indicate the square root of magnetic field variance which are computed over 1 day intervals using hourly magnetic field data of Ulysses, red line represents the results of the TRST model, and $\delta B_{1AU}$ is calculated in the same way using magnetic field data of OMNI. Considering the time delayed heliosphere, the Ulysses data has been shifted back to $1$ AU with the typical solar wind speed $v_{0}=0.25$ AU/day. Note that we compute the magnetic field variance over 1 day intervals instead of Carrington rotation intervals due to the large latitude variation during the fast latitude scan. From this figure we can see that the turbulence model TRST provides a good prediction of $\delta B$ observations of Ulysses.

Figure 4 shows magnetic turbulence as a function of radial distance from $2$ AU to $80$ AU. Black circles in the top and bottom panels mean the square root of magnetic field variance which are computed over Carrington rotation intervals using hourly HMF data of Voyager 1 and Voyager 2, respectively. We have also computed the magnetic field variance using methods referred by other works, e.g., Zank et al. (1996); Smith et al. (2006); Isenberg et al. (2010); Adhikari et al. (2015), and the results show little difference. Considering the time delayed heliosphere, the Voyager data has been shifted back to $1$ AU with the typical solar wind speed $v_{0}=0.25$ AU/day. Red lines indicate the results of the TRST model, with $\delta B_{1AU}$ calculated in the same way using hourly HMF data from OMNI. It is shown that the TRST model provides a good prediction of the Voyager 1 and 2 observations. From the top panel we can see that $\delta B$ of Voyager 1 decreases faster during solar minimum to form wave troughs, but from bottom panel the data of Voyager 2 does not show clear wave troughs. It is assumed that the difference between the Voyager 1 and 2 magnetic variance is caused by the latitudinal dependence of turbulence magnitude. Generally speaking, the model TRST provides a good prediction of the turbulence variation properties in the inner heliosphere.

Here, we consider two-component model of turbulence (Matthaeus et al., 1990) in solar wind. It is also necessary for one to know the transport of other properties of turbulence in addition to magnitude. However, for simplicity, we assume only the turbulence magnitude is varying.

2.5 Diffusion coefficients

Turbulent magnetic fields in the solar wind plasma result in diffusion of the cosmic rays parallel and perpendicular to the background HMF, which plays an important role in the modulation processes. In the scattering theory we usually emphasis on the global behavior of diffusion coefficients (or mean free paths). In the field-aligned coordinate, the symmetric and diagonal diffusion tensor $\mathbf{K}_{s}$ is composed of three parts: a parallel diffusion coefficient $\kappa_{\parallel}$ and two perpendicular diffusion coefficients, $\kappa_{\perp r}$ and $\kappa_{\perp\theta}$, the perpendicular diffusion coefficients in the radial and polar directions, respectively. Jokipii (1966) developed the quasi-linear theory (QLT) of the diffusion of cosmic rays, which is considered one of the milestones of the study of cosmic rays. It has been considered that QLT is relatively good to describe the parallel diffusion of cosmic rays, but perpendicular diffusion has long been a puzzle. Therefore, empirical models are used in many studies. For example, in some study of GCR modulations, empirical expressions for the parallel diffusion coefficient are used based on the QLT, and perpendicular ones are set to be proportional to the parallel one (e.g., Potgieter, 2013; Potgieter et al., 2014; Zhao et al., 2014; Potgieter et al., 2015; Vos & Potgieter, 2015; Raath et al., 2016; Guo & Florinski, 2016). The results with this approach are usually consistent with the observations at 1 AU, but some parameters in the diffusion coefficient models have to be decided by comparing numerical results with the observations.

Matthaeus et al. (2003) developed a nonlinear guiding center (NLGC) theory for perpendicular diffusion which agrees well with numerical simulations. Further, Qin (2007) extended the NLGC to describe the parallel diffusion, noted as NLPA. With the combination of the NLGC and NLPA models one gets two implicit integral equations which can be solved simultaneously to obtain the perpendicular and parallel diffusion coefficients. Qin & Zhang (2014) further improved the combination of the NLGC and NLPA with some slight modification, then they obtained a new model, NLGCE-F, by fitting the numerical solution from the improved combination of NLGC and NLPA with polynomials. The model NLGCE-F allows one to calculate diffusion coefficients directly without the iteration solution of integration equations set. It is noted that in order to use this model the properties of HMF and turbulence in solar wind are necessary. In this study, we assume that $\kappa_{\perp r}=\kappa_{\perp\theta}$. The expressions for NLGCE-F are as follows:

[TABLE]

with

[TABLE]

where $\sigma$ indicates $\perp$ or $\parallel$, $\lambda_{\sigma}=\frac{3}{v}\kappa_{\sigma}$, $R_{L}$ means the gyro-radius of the particle, $\lambda_{slab}$ and $\lambda_{2D}$ are the spectral bend-over scales of the slab and 2D components of turbulence, respectively, $E_{total}=\left\langle{\delta B^{2}}\right\rangle$ and $E_{slab}=\left\langle{\delta B^{2}}_{slab}\right\rangle$ are the magnetic turbulence energy from all components and from slab component, respectively, and $\delta B/B$ is the turbulence level. The coefficients $d_{i,j,k,l}^{\sigma}$ and polynomial order $n_{\sigma i}$ are provided by Qin & Zhang (2014), and the computer code with parameters for NLGCE-F can be downloaded in www.qingang.org.cn/code/NLGCE-F. From Qin & Zhang (2014) it is also noted that the model NLGCE-F is valid with the parameters

[TABLE]

In this work we set ${\lambda_{slab}}/{\lambda_{2D}}=10.0$ (Matthaeus et al., 2003), $\lambda_{slab}=0.02r$ with $r$ being the solar distance, and $E_{slab}/E_{total}=0.2$ (Bieber et al., 1994) throughout the heliosphere. The turbulence parameters in solar wind, such as $\lambda_{slab}$ and $\lambda_{2D}$, can only be observed by spacecraft indirectly with complicated theoretical study (e.g., Matthaeus et al., 1990; Adhikari et al., 2017; Zank et al., 2017), so for simplicity purpose we set them in simple forms according to some study for solar wind in 1 AU (e.g., Matthaeus et al., 2003). It is noted that if the particle’s energy is not much more than $10$ GeV and the radial distance is not larger than the distance of termination shock, the values of input parameters in Equation (34) are in the ranges of validation.

Using diffusion model Equation (34) with the turbulence model Equation (32), we are able to establish a time-dependent diffusion coefficients model with all input parameters obtained from the spacecraft observations near Earth. Manuel et al. (2014) also established a time-dependent diffusion model with the time-dependent parameters scaled by HMF magnitude (B) and variance ($\delta B^{2}$) (see also, Ferreira & Potgieter, 2004; Manuel et al., 2011; Potgieter et al., 2014).

Figure 5 shows scenarios of mean free paths as a function of rigidity (top panel), polar angle (second panel), and radial distance (bottom panel). Interplanetary conditions are the same as that in Figure 2, and $\delta B_{1AU}=3.0$ nT. The parallel mean free path $\lambda_{\parallel}$ shows stronger rigidity dependence when the rigidity increases, and $\lambda_{\parallel}$ is larger in the polar region as have been shown in the first and second panels. At lower energy, i.e., when energy is lower than about 3 GV, the parallel mean free path shows the expected $P^{1/3}$ dependence, but at higher energy the model gives a $\sim P^{3/2}$ dependence. The perpendicular mean free path is relatively flat as a function of rigidity and colatitude. As has been shown in the bottom panel, the parallel and perpendicular mean free paths show a gradual increase with the radial distance.

2.6 Heliospheric boundary

Voyager 1 crossed the heliopause at $121.7$ AU in August 2012 and has reached the very local interstellar medium, Zhang et al. (2015) believed that the solar modulation boundary is located a fraction of an AU beyond the heliopause. Therefore, the very local interstellar (LIS, e.g., Potgieter et al., 2014; Vos & Potgieter, 2015) spectrum of GCRs could be used as the input spectrum in our modulation model. However, in our model, we set the outer boundary at a smaller solar distance $r=85$ AU, for simplicity purpose, so that we do not include the termination shock acceleration of GCRs and other complicated phenomenon in outer heliosphere. Furthermore, to keep solar distance $r$ not too large could make sure the diffusion model NLGCE-F always valid in this work. In addition, we assume the GCR source at $r=85$ AU as

[TABLE]

where $J_{0}$ is a constant determined later and $p_{0}=1~{}\mathrm{GeV}/c$ by following Zhang (1999).

3 Numerical Methods

To solve the Parker transport equation, we make use of the time-backward Markov stochastic process method proposed by Zhang (1999). For a pseudo-particle in position $(r,\theta,\phi)$ and momentum $p$, the stochastic differential equations equivalent to Equation (1) have the form (Zhang, 1999; Pei et al., 2010b; Strauss et al., 2011a; Kopp et al., 2012)

[TABLE]

with $i\in{(r,\theta,\phi,p)}$, $x_{i}$ the Ito processes (Zhang, 1999), $s$ the backward time and $\mathrm{d}W_{i}$ satisfy a Wiener process given by the standard normal distribution (Pei et al., 2010b; Strauss et al., 2011a). For the modified Parker HMF used in this work, the matrix components $B_{ij}$ are given by Pei et al. (2010b) (see also Kopp et al., 2012),

[TABLE]

and the components of vector $\mathbf{A}$ are given as follows,

[TABLE]

Therefore the statistical differential equations can be written as

[TABLE]

Note that the diffusion tensor in Equations (55-58) are elements of the symmetric diffusion tensor $\mathbf{K}_{s}$ in spherical coordinates. According to Burger et al. (2008) elements of $\mathbf{K}_{s}$ in spherical coordinates for the modified Parker HMF are written as

[TABLE]

with $\tan\Psi=-B_{\phi}/(B_{r}^{2}+B_{\theta}^{2})^{\frac{1}{2}}$ and $\tan\zeta=B_{\theta}/{B_{r}}$, where $\Psi$ is the HMF winding angle.

In our modulation model we use the time-delayed interplanetary conditions at radius $r$ related to the states at the source surface $r_{s}$ at some earlier time, and the heliosphere is considered dynamic due to the solar activities.

4 Modeling Results

In our GCR modulation model we only need four input parameters which can be obtained from the observations at $1$ AU, i.e., the heliospheric current sheet tilt angle, the solar wind speed, the magnitude of background magnetic field $B$, and the magnetic turbulence magnitude $\delta B_{1AU}$. Figure 6 shows the computed and observed proton spectra for four Carrington rotations in November 2006, December 2007, December 2008, and December 2009 with colors cyan, purple, red, and blue, respectively. Circles means observations of PAMELA, numerical results of modulation modeling are shown as solid lines. The GCR source at $85$ AU is represented by the black line with the constant $J_{0}=1.17\times 10^{4}$ $\mathrm{m}^{-2}\mathrm{s}^{-1}\mathrm{sr}^{-1}\mathrm{(GeV/nuc)}^{-1}$ in Equation (42), and magenta triangles mean Voyager 2 observations at 85 AU reported by Webber et al. (2008). The modulation results show good agreement with the observations of PAMELA.

5 Discussion and Conclusions

In this work, we develop a numerical model to study the time-dependent modulation of cosmic rays in recent solar minimum with PAMELA observations. We use the time-backward Markov stochastic process method (Zhang, 1999) to numerically solve the Parker transport equation. In our GCR modulation model, all the parameters are obtained from the observations of OMNI. We get galactic proton spectra varying as a function of time during the recent solar minimum, which are consistent with the observations of PAMELA.

As the Parker HMF provides a low magnitude in the polar regions at large radial distance, we adopt the modified Parker HMF according to Jokipii & Kota (1989), which can help to avoid the traditional Parker HMF’s problem that the drift speed of GCRs in polar regions are too large. Considering the dynamic phenomena of heliosphere with the solar wind flow from source surface to any solar distance, we use the input parameters observed near the Earth in the earlier time with a typical solar wind speed $v_{0}=0.25$ AU/day. It is noted that this will divide the heliosphere into slices and introduce an additional radial dependence in the HMF magnitude, i.e., $B_{0}=B_{0}(r)$. In the region between two slices of plasma the HMF is not divergence free, but in each step of the numerical solution of the TPE Equation (1), we always keep inside one slice of plasma. In addition, we set the outer boundary at a smaller solar distance $r=85$ AU, for simplicity purpose, so that we do not include the termination shock acceleration of GCRs and other complicated phenomenon in outer heliosphere. The GCR source spectrum we use is consistent with the observations of Voyager 2 at 85 AU reported by Webber et al. (2008). Furthermore, by keeping solar distance $r$ not too large we can make sure the diffusion model NLGCE-F always valid in this work.

The knowledge of transport of magnetic turbulence throughout the heliosphere is very important to determine the diffusion coefficients. According to previous studies (e.g., Zank et al., 1996; Burger et al., 2008; Effenberger et al., 2012; Ngobeni & Potgieter, 2014; Strauss et al., 2017; Perri & Balogh, 2010), we use a model for magnetic turbulence magnitude with Equation (32), i.e., $\delta B\sim r^{S}(1+\sin^{2}\theta)$, in which the latitudinal dependence is assumed from the observations of Ulysses, and the expression of $S$ for the radial dependence is chosen as a function of the heliospheric current sheet tilt angle with Equation (33). We show that the new turbulence magnitude model with Equations (32) and (33), denoted as TRST model, agrees well with the Ulysses, Voyager 1, and Voyager 2 observations. In addition, we assume two-component model of turbulence in solar wind. For simplicity purpose, we only suppose the magnetic turbulence magnitude is varying.

We use the new diffusion model NLGCE-F from Qin & Zhang (2014) which was obtained by fitting the numerical solution from the non-linear parallel and perpendicular diffusion with polynomials. The using of the diffusion model NLGCE-F helps us to get more accurate diffusion coefficients without consuming lots of computing resources. For the drift coefficient, turbulence can provide the suppression (see, e.g., Jokipii, 1993; Fisk & Schwadron, 1995; Giacalone & Jokipii, 1999; Candia & Roulet, 2004; Stawicki, 2005; Minnie et al., 2007; Tautz & Shalchi, 2012). The reduction of drift effects is complicated to be used self-consistently (see, e.g., Bieber & Matthaeus, 1997; Burger & Visser, 2010; Tautz & Shalchi, 2012) or in Ad hoc form (see, e.g., Burger et al., 2000; Potgieter, 2013; Vos & Potgieter, 2016; Nndanganeni & Potgieter, 2016) in modulation works. It is far from complete to understand the effects of turbulence on CR drifts. Therefore, in this work, we use the weak scattering drift coefficient for simplicity purpose.

In the future, we plan to use the modulation model established in this paper to study the 11 and 22 year modulation of GCRs in the inner heliosphere (e.g., McDonald, 1998; Shen & Qin, 2016). If our model works well, we can reproduce the GCR observations by Ulysses, Voyager 1, and Voyager 2 with long period of time. Otherwise, we need to improve our modulation model. Firstly, we could improve the turbulence model by modifying the magnitude model and applying more realistic models for transport of turbulence geometry. Secondly, we could use a more self-consistent dynamic heliosphere model, e.g., a model from MHD simulation. Thirdly, we could include termination shock in the model to study the realistic boundary effects. Fourthly, the GCR source spectrum could be improved. Fifthly, the drift suppression from turbulence could be included.

We are partly supported by grants NNSFC 41374177, NNSFC 41574172, and NNSFC 41125016, and the Specialized Research Fund for State Key Laboratories of China. We used the HCS tilt angle data from the Wilcox Solar Observatory ( wso.stanford.edu ), solar wind speed and magnetic field data from OMNI website ( https://omniweb.gsfc.nasa.gov/ ), and PAMELLA data from Database for Charged Cosmic Ray measurements ( https://tools.asdc.asi.it/CosmicRays/ ). The work was carried out at National Supercomputer Center in Tianjin, and the calculations were performed on TianHe-1 (A).