11institutetext: Université Paris-Saclay, CNRS, Institut d’astrophysique
spatiale, 91405, Orsay, France 22institutetext: AIM, CEA, CNRS, Université
Paris-Saclay, Université Paris-Diderot, Sorbonne Paris Cité, F-91191 Gif-sur-
Yvette, France 33institutetext: IRAP, Université de Toulouse, CNRS, UPS, CNES,
Toulouse, France, 14 Avenue Edouard Belin, F-31400 Toulouse, France

# Impact of solar magnetic field amplitude and geometry on cosmic rays
diffusion coefficients in the inner heliosphere

B. Perri 1122 A. S. Brun 22 A. Strugarek 22 V. Réville 33

###### Abstract

Cosmic rays are remarkable tracers of solar events when they are associated
with solar flares, but also galactic events such as supernova remnants when
they come from outside our solar system. Solar Energetic Particles (SEPs) are
correlated with the 11-year solar cycle while Galactic Cosmic Rays (GCRs) are
anti-correlated due to their interaction with the heliospheric magnetic field
and the solar wind. Our aim is to quantify separately the impact of the
amplitude and the geometry of the magnetic field, both evolving during the
solar cycle, on the propagation of cosmic rays of various energies in the
inner heliosphere (within Earth orbit). We focus especially on the diffusion
caused by the magnetic field along and across the field lines. To do so, we
use the results of 3D magnetohydrodynamics (MHD) wind simulations running from
the lower corona up to 1 AU. This gives us the structure of the wind and the
corresponding magnetic field. The wind is modeled using a polytropic
approximation, and fits and power laws are used to account for the turbulence.
Using these results, we compute the parallel and perpendicular diffusion
coefficients of the Parker cosmic ray transport equation, yielding 3D maps of
the diffusion of cosmic rays in the inner heliosphere. By varying the
amplitude of the magnetic field, we change the amplitude of the diffusion by
the same factor, and the radial gradients by changing the spread of the
current sheet. By varying the geometry of the magnetic field, we change the
latitudinal gradients of diffusion by changing the position of the current
sheets. By varying the energy, we show that the distribution of values for
SEPs is more peaked than GCRs. For realistic solar configurations, we show
that diffusion is highly non-axisymmetric due to the configuration of the
current sheets, and that the distribution varies a lot with the distance to
the Sun with a drift of the peak value. This study shows that numerical
simulations, combined with theory, can help quantify better the influence of
the various magnetic field parameters on the propagation of cosmic rays. This
study is a first step towards the resolution of the complete Parker transport
equation to generate synthetic cosmic rays rates from numerical simulations.

###### keywords:

MHD – solar wind – cosmic rays

##  1 Introduction

The Sun possesses a magnetic field that shows a cyclic evolution in time : it
has a cycle of 11 years in amplitude and 22 years in polarity on average, the
shortest cycle observed being of 9 years and the longest one of 14 years
(Hathaway, 2015; Brun and Browning, 2017). During a minimum of activity, the
solar magnetic field has the lowest amplitude of the cycle and its geometry is
mostly dipolar ; during a maximum of activity, its amplitude is at its peak
and its geometry is mostly quadrupolar (DeRosa et al., 2012). This magnetic
field is generated inside the star via a dynamo loop (Moffatt, 1978; Parker,
1993) and fills the whole heliosphere, including the Earth spatial environment
(Owens and Forsyth, 2013).

The heliosphere is filled with the solar wind, a continuous flow of charged
particles ejected from the Sun (Neugebauer and Snyder, 1962). It has a slow
and a fast component at respectively 400 and 800 km/s at Earth orbit, which
corresponds to 1 AU, and hence is transsonic and transalfvénic at this
distance from the Sun (McComas et al., 2003). It was first described using
fluid dynamics (Parker, 1958), then magnetism was taken into account (Weber
and Davis, 1967; Sakurai, 1985). Multiple models of the solar wind have been
designed, from empirical models (Wang and Sheeley, 1990; Arge and Pizzo, 2000)
to MHD numerical simulations in 1D (Lionello et al., 2001; Suzuki et al.,
2013; Pinto and Rouillard, 2017), 2D (Keppens and Goedbloed, 1999; Matt and
Pudritz, 2008; Réville et al., 2015) or 3D (Tóth et al., 2012; Riley et al.,
2015; Réville et al., 2020). The heating of the corona is modeled through a
polytropic approximation (Réville and Brun, 2017) or via Alfvén waves
perturbations (Usmanov et al., 2014). The complete list of phenomena leading
to this heating still eludes our understanding and is a huge current numerical
challenge (Réville et al., 2020).

Cosmic rays (CRs) are highly energetic extra-terrestrial particles with
energies between $10^{2}$ MeV to $10^{11}$ GeV ; they follow a power-law
distribution, except for the low-energy part of the distribution (Reames,
1999; Heber and Potgieter, 2006). They can be emitted by the Sun during sudden
events such as solar flares or coronal mass ejections ; in that case they are
called Solar Energetic Particles (SEPs) and correspond to the low-energy part
of the distribution (up to 1 GeV). They can also be emitted by sudden events
out of our solar system such as gamma-ray bursts or supernova remnants ; in
that case, they are called Galactic Cosmic Rays (GCRs) and correspond to the
high-energy part of the distribution (from 1 GeV) (Shalchi, 2009). As they
progress through the heliosphere, CRs are subject to an adiabatic cooling
while interacting with both the heliospheric magnetic field and the solar
wind, which change significantly their trajectory due to their turbulent
fluctuations (Parker, 1964; Jokipii, 1966). CRs rate is thus influenced by the
cyclic activity of the Sun : SEPs are correlated with solar activity because
sudden solar events are more frequent at maximum of activity ; on the
contrary, GCRs are anti-correlated with solar activity because the magnetic
field at maximum of activity makes it harder for GCRs to penetrate the
heliosphere (Snyder et al., 1963; Heber and Potgieter, 2006). There are also a
certain number of disparities in the CR distribution. The Voyager missions
have suggested the presence of a negative latitudinal gradient in the count
rate of $>70$ MeV protons (Cummings et al., 1987). Ulysses has shown that SEPs
have a North-South asymmetry linked to the one of the magnetic field and the
wind (McKibben, 1998; Perri et al., 2018); it has also shown that GCRs have
larger gradients in their spatial distribution at minimum than at maximum
(Belov et al., 2003). The observed electron to proton ratios, also linked to
the radial and latitudinal gradients, indicate that large particle drifts are
occurring during solar minimum, but diminish significantly towards solar
maximum (Heber and Potgieter, 2006).

To describe the propagation of CRs, the most common approach is the
statistical one using the Parker cosmic ray transport equation from Parker
(1965). One of the biggest challenges in this equation is to determine the
diffusion tensor, especially its dependency in space, time and energy. For the
diffusion parallel to the magnetic field lines, the quasi-linear theory (QLT)
yields good results, especially when extended to take into account time-
dependent and non-linear corrections (Jokipii, 1966; Goldstein, 1976; Bieber
et al., 1994; Dröge, 2003). For the diffusion perpendicular to the magnetic
field lines however, QLT provides only an upper limit using a field line
random walk (FLRW) description (Jokipii, 1966; Forman et al., 1974; Giacalone
and Jokipii, 1999). Various alternate approaches were tested, including the
Taylor-Green-Kubo (TGK) (Taylor, 1922; Green, 1951; Kubo, 1957; Forman, 1977)
or the Bieber and Matthaeus (BAM) (Bieber and Matthaeus, 1997) formulations,
but these methods systematically underestimates the perpendicular diffusion
(Bieber et al., 2004). Finally, the best method to this day is the non-linear
guiding center (NLGC) theory (Matthaeus et al., 2003; Bieber et al., 2004;
Shalchi, 2009) which provides the best agreement with both observations and
simulations. This is due to the assumption that there is a decorrelation
between the diffusive spread of the particle gyrocenters following field lines
and the diffusive spread of those field lines, due to the transverse
complexity of the magnetic field (Matthaeus et al., 2003). There have been
some recent reformulations of this theory such as the Extended NLGC (ENLGC)
theory by Shalchi (2006) to improve the slab contribution, or with the random
ballistic decorrelation (RBD) interpretation (Ruffolo et al., 2012) which
helped further improve the theory, by matching simulations over a wider range
of fluctuation amplitudes ; recently there was also the application of the
Reduced MHD (RMHD) for astrophysics derived by Oughton et al. (2015, 2017) and
the UNLT (Unified Non-Linear Transport) by Shalchi (2017, 2020) as a unified
theory between all the previous approaches.

In most theoretical studies about CR diffusion, prescriptions are used for the
magnetic field and the solar wind ; this usually limits the applications to
certain range of energy or spatial locations. However, thanks to MHD numerical
simulations, it is possible to have global descriptions of these quantities
even in complex configurations. This approach has already been used in the
studies of Luo et al. (2013), Guo and florinski (2014), Wiengarten et al.
(2016) and Kim et al. (2020) to predict the variations of CRs in the complete
heliosphere, coupled with semi-empirical CRs prescriptions. However, we will
focus here only on the very inner heliosphere within Earth orbit, which means
that we do not include the whole dynamics of GCRs coming from outside the
solar system. Similar approach has been used by Chhiber et al. (2017) to study
the diffusion coefficient in the case of a tilted dipole. However in this
study, only the inclination of the dipole was changed, and only meridional
cuts were shown. In this study, we are interested in the correlation with
cyclic activity of CRs, and thus want to study separately the two variations
of the magnetic field over a cycle : variation in amplitude and in geometry.
Thus we will characterize the difference between a reference case being a
dipole of weak amplitude called D1, and a dipole of strong amplitude called
D10 and a quadrupole of weak amplitude called Q1. Finally we will use the 3D
aspect of our simulations by using the same method with realistic
configurations of minimum and maximum of activity using synoptic maps. The
minimum of activity corresponds to October 1995 and the maximum to August 1999
(Carrington rotations of CR 1902 and CR 1954 respectively). Both maps come
from Wilcox Observatory. All results of the simulations and corresponding
post-processing are available at the MEDOC online
facility.111https://idoc.ias.u-psud.fr/wind_predict_cr/

The article is organized as follows. Section 2 presents the wind model used
for the simulations and the post-processing statistical computation of the
cosmic rays diffusion coefficient. Section 3 details the various parametric
studies we performed on the impact of the magnetic field amplitude and
geometry, and the cosmic rays energy. Section 4 introduces our results for
realistic configurations corresponding to a minimum and maximum of activity.
Finally section 5 sums our conclusions on this study and perspectives for
future work.

##  2 Model and equations

In this section we present first the 3D MHD wind model we used to derive the
magnetic field structure and intensity and the solar wind speed in the inner
heliosphere. Then we present the model used in post-processing to compute the
diffusion of cosmic rays along and across magnetic field lines.

###  2.1 Wind model

Our wind model is adapted from Réville et al. (2015); Réville and Brun (2017)
and Perri et al. (2018) using the multi-physics compressible PLUTO code
(Mignone et al., 2007). In these articles the code has shown good agreement
with other wind codes such as the model from Matt and Pudritz (2008), or the
code DIP (Grappin et al., 2010); we have also checked for the conservation of
MHD invariants in Strugarek et al. (2015) and Réville et al. (2015). We solve
the set of the conservative ideal MHD equations composed of the continuity
equation for the density $\rho$, the momentum equation for the velocity field
$\mathbf{v}$ with its momentum written $\mathbf{m}=\rho\mathbf{v}$, the energy
equation which is noted $E$ and the induction equation for the magnetic field
$\mathbf{B}$:

$\frac{\partial}{\partial t}\rho+\nabla\cdot\left(\rho\mathbf{v}\right)=0,$
(1) $\frac{\partial}{\partial
t}\mathbf{m}+\nabla\cdot(\mathbf{mv}-\mathbf{BB}+\mathbf{I}p)=\rho\mathbf{a},$
(2) $\frac{\partial}{\partial
t}E+\nabla\cdot((E+p)\mathbf{v}-\mathbf{B}(\mathbf{v}\cdot\mathbf{B}))=\mathbf{m}\cdot\mathbf{a},$
(3) $\frac{\partial}{\partial
t}\mathbf{B}+\nabla\cdot(\mathbf{vB}-\mathbf{Bv})=0,$ (4)

where $p$ is the total pressure (thermal and magnetic), $\mathbf{I}$ is the
identity matrix and $\mathbf{a}$ is a source term (gravitational acceleration
in our case). We use a polytropic assumption, which yields the following ideal
equation of state: $\rho\varepsilon=p_{th}/(\gamma-1)$, where $p_{th}$ is the
thermal pressure, $\varepsilon$ is the internal energy per mass and $\gamma$
is the adiabatic exponent. This gives for the energy :
$E=\rho\varepsilon+\mathbf{m}^{2}/(2\rho)+\mathbf{B}^{2}/2$.

PLUTO solves normalized equations, using three variables to set all the
others: length, density and speed. If we note with $*$ the parameters related
to the star and with $0$ the parameters related to the normalization, we have
$R_{*}/R_{0}=1$, $\rho_{*}/\rho_{0}=1$ and
$v_{kep}/V_{0}=\sqrt{GM_{*}/R_{*}}/V_{0}=1$, where $v_{kep}$ is the Keplerian
speed at the stellar surface and $G$ the gravitational constant. By choosing
the physical values of $R_{0}$, $\rho_{0}$ and $V_{0}$, one can deduce all of
the other values given by the code in physical units. In our set-up, we choose
$R_{0}=R_{\odot}=6.96\ 10^{10}$ cm, $\rho_{0}=\rho_{\odot}=1.67\ 10^{-16}\
\mathrm{g/cm}^{3}$ and $V_{0}=v_{kep,\odot}=4.37\ 10^{2}$ km/s. We can define
the escape velocity from the Keplerian speed as
$v_{esc}=\sqrt{2}v_{kep}=\sqrt{2GM_{\odot}/R_{\odot}}$.

Control parameter | PLUTO control parameters | Coronal parameters
---|---|---
Density | $\rho_{0}=1.67\times 10^{-16}\ \mathrm{g/cm}^{3}$ | $n=1.0\times 10^{8}\ \rm{cm}^{-3}$
Rotation rate | $v_{rot}/v_{esc}=2.93\times 10^{-3}$ | $\Omega_{0}=2.6\times 10^{-6}\ \rm{s}^{-1}$
Temperature |  $c_{s}/v_{esc}=0.243$ & $\gamma=1.05$ | $1.6\ 10^{6}\ \rm{K}$
Magnetic amplitude | $v_{A}/v_{esc}=0.176$ | $0.5\ \rm{G}$
Table 1: Table of the control parameters of the 3D MHD wind simulation for the
reference case D1. The four control parameters are the density, rotation rate,
temperature and magnetic amplitude. They are expressed in the PLUTO code
normalization system and with the correspondence in physical values.

Our wind simulations are controlled by four parameters : the adiabatic
exponent $\gamma$ for the polytropic wind, the normalized rotation of the star
$v_{rot}/v_{esc}$, the normalized speed of sound $c_{s}/v_{esc}$ and the
normalized Alfvén speed at the equator $v_{A}/v_{esc}$. For the rotation
speed, we take the solar value, which gives $v_{rot}/v_{esc}=2.93\ 10^{-3}$.
We define a reference case for our simulations named D1. For this case, we
choose to fix $c_{s}/v_{esc}=0.243$, which corresponds to a $1.6\ 10^{6}$ K
hot corona for solar parameters and $\gamma=1.05$. This choice of $\gamma$ is
dictated by the need to maintain an almost constant temperature as the wind
expands, which is what is observed in the solar wind. Hence, choosing
$\gamma\neq 5/3$ is a simplified way of taking into account heating in the low
corona, which is not modeled here. However a new model based on heating by
Alfvén waves has been developed and has shown good agreement with the Parker
Solar Probe data, see Réville et al. (2020). The amplitude of the magnetic
field is set by $v_{A}/v_{esc}=0.176$, which corresponds to an amplitude of
0.5 G at the equator for a dipole. All these parameters are summed up in table
1. Some parameters will vary depending on the model discussed, the differences
between the various cases are shown in table 2.

We use the spherical coordinates $(r,\theta,\phi)$. We choose a finite-volume
method using an approximate Riemann Solver (here the HLL solver, cf. Einfeldt
(1988)). PLUTO uses a reconstruct-solve-average approach using a set of
primitive variables $(\rho,\mathbf{v},p,\mathbf{B})$ to solve the Riemann
problem corresponding to the previous set of equations. The time evolution is
then implemented via a second order Runge-Kutta method. To enforce the
divergence-free property of the field, we use a hyperbolic divergence
cleaning, which means that the induction equation is coupled to a generalized
Lagrange multiplier in order to compensate the deviations from a divergence-
free field (Dedner et al., 2002). We use a splitting between the curl-free
background field and the deviation field.

Figure 1: Boundary conditions (on the left) and example of wind simulation (on
the right). For the right panel, we show the relaxed state corresponding to
the reference case D1. The color scale represents the following quantity :
$\mathbf{v}\cdot\mathbf{B}/(c_{s}||\mathbf{B}||)$, which is the solar wind
velocity projected on the magnetic field in units of Mach number. The black
line corresponds to the Alfvén surface where the wind speed equals the Alfvén
speed. White lines correspond to the poloidal magnetic field lines of positive
polarity in solid and negative polarity in dashed lines. We represent only the
10 first solar radii.

The numerical domain dedicated to the wind computation is a 3D sphere with the
radius $r\in[1.001,220]R_{\odot}$, the co-latitude $\theta\in[0,\pi]$ and the
longitude $\phi\in[0,2\pi]$. We use an uniform grid in latitude and longitude
with respectively 256 and 512 points, and a stretched grid in radius with 256
points; the grid spacing is geometrically increasing from $\Delta
r/R_{\odot}=0.001$ at the surface of the star to $\Delta r/R_{\odot}=0.02$ at
the outer boundary. At the latitudinal boundaries ($\theta=0$ and
$\theta=\pi$), we set axisymmetric boundary conditions. At the longitudinal
boundary conditions ($\phi=0$ and $\phi=2\pi$), we set periodic boundary
conditions. At the top radial boundary ($r=220R_{\odot}$), we set an outflow
boundary condition which corresponds to $\partial/\partial r=0$ for all
variables, except for the radial magnetic field where we enforce
$\partial(r^{2}B_{r})/\partial r=0$. Because the wind has opened the field
lines and under the assumption of axisymmetry, this ensures the divergence-
free property of the field. The bottom boundary conditions are shown in the
left panel of figure 1. In the ghost cells (in orange), the density $\rho$ and
pressure $p$ are set to a polytropic profile, the rotation is uniform, the
poloidal speed $V_{\rm{pol}}$ is aligned with the poloidal magnetic field
$B_{\rm{pol}}$; the latter is imposed by a background dipolar field while the
toroidal magnetic field $B_{\phi}$ is linear. In the first point of the
computational domain (in blue), all physical quantities are free to evolve,
except for the poloidal speed $V_{\rm{pol}}$ which is forced to be aligned
with the poloidal magnetic field $B_{\rm{pol}}$ to minimize the generation of
currents at the surface of the star and keep it as close as possible to a
perfect conductor. We initialize the velocity field with a 1D Parker-like
polytropic wind solution and the magnetic field with either a dipole, a
quadrupole or a realistic magnetic field configuration from a synoptic map
depending on the study case.

###  2.2 Cosmic-ray post-processing

####  2.2.1 Diffusion coefficients

We use a statistical approach based of a Fokker-Planck equation, starting from
the Parker cosmic ray transport equation presented in Parker (1965):

$\frac{\partial U}{\partial t}+\frac{\partial}{\partial
x_{i}}\left(Uv_{i}\right)+\frac{\partial}{\partial T}\left(U\frac{\partial
T}{\partial t}\right)-\frac{\partial}{\partial
x_{i}}\left(\kappa_{ij}\frac{\partial U}{\partial x_{j}}\right)=0,$ (5)

where $U(x_{i},T,t)$ is the distribution of cosmic rays depending on their
spatial coordinates $x_{i}$, kinetic energy $T$ and time $t$, $\kappa_{ij}$ is
the diffusion tensor and $v_{i}$ is the wind speed. With the numerical
simulations using the wind model described above, we have a prescription for
the wind speed and the magnetic field structure and amplitude ; the only term
left to compute is the diffusion tensor $\kappa$ which needs to be modeled. We
will focus on this aspect in this article.

The diffusion tensor can be decomposed into three terms (Jokipii and Parker,
1970):

$\kappa_{ij}=\kappa_{\perp}\delta_{ij}+\frac{B_{i}B_{j}}{B^{2}}(\kappa_{||}-\kappa_{\perp})+\epsilon_{ijk}\kappa_{A}\frac{B_{k}}{B},$
(6)

with $\delta_{ij}$ being the Kronecker symbol and $\epsilon_{ijk}$ the Levi-
Civita tensor. $\kappa_{\parallel}$ is the diffusion along the magnetic field
lines and $\kappa_{\perp}$ is the diffusion across the magnetic field lines.
The coefficient $\kappa_{A}$ is the drift coefficient and intervenes mostly
for very energetic particles and strong gradients of the magnetic field
(Jokipii and Levy, 1977). It takes into account the influence of the current
sheet (Jokipii and Thomas, 1981), the solar tilt angle (Lockwood and Webber,
2005) and heliospheric magnetic field polarity (Webber et al., 2005). Drift
effects also contribute to the 22-year cycle observed in the CR intensity and
the CR latitudinal gradients (Heber et al., 1996). Recent studies have shown
that the drift coefficient needs to be reduced to match spacecraft
observations due to turbulence (Manuel et al., 2011), especially for SEPs in
the inner heliosphere (Engelbrecht and Burger, 2015; Engelbrecht et al.,
2017). In this study, we will first focus on the diffusion coefficients which
are better characterized and more prominent for SEPs in the inner heliosphere,
and we will include drift effects in a later study. We can also introduce mean
free path (mfp) $\lambda$, related to the diffusion tensor $\kappa$ by the
following relationship:

$\kappa_{\parallel}=\frac{1}{3}\lambda_{\parallel}\mathbf{v}_{CR},\kappa_{\perp}=\frac{1}{3}\lambda_{\perp}\mathbf{v}_{CR},$
(7)

where $v_{CR}$ is the particle speed.

In this study, we will thus focus on describing the parallel and perpendicular
mean free paths (mfps). To do so, we will not go into details into all the
formulations proposed, but only used the most recent ones which have reached a
general consensus, detailed hereafter; the reader can however find some very
complete reviews in Shalchi (2009) for parallel diffusion and Shalchi (2020)
for perpendicular diffusion. The most efficient geometry to describe such
parameters is the composite geometry described in Bieber et al. (1994) with
80% of 2D geometry (both magnetic fluctuations and wave vectors are
perpendicular to the magnetic field) and 20% of slab geometry (magnetic
fluctuations are perpendicular to the magnetic field, but wave vectors are
parallel to it). It is also supported by wind observations that show a strong
2D component of the turbulence (Matthaeus et al., 1990). From now on,
quantities related to the 2D geometry will be noted with an index $2$ while
quantities related to the slab geometry will be noted with an index $s$.

A good approximation for parallel mfp is given by Zank et al. (1998):

$\lambda_{\parallel}=6.2742\frac{B^{5/3}}{\langle
b_{s}^{2}\rangle}\left(\frac{P}{c}\right)^{1/3}\lambda_{s}^{2/3}\left[1+\frac{7A_{\parallel}/9}{(1/3+q)(q+7/3)}\right]\
\rm{km},$ (8)

with $B$ the magnetic field norm, $\langle b_{s}^{2}\rangle$ the variance of
the slab geometry fluctuation, $P=\tilde{p}c/Ze$ the particle rigidity
($\tilde{p}$ being the particle momentum and $Ze$ the particle charge), $c$
the speed of light, $\lambda_{s}$ the correlation length for the slab
turbulence, and:

$A_{\parallel}=(1+s^{2})^{5/6}-1,$ (9)
$q=\frac{5s^{2}/3}{1+s^{2}-(1+s^{2})^{1/6}},$ (10)

and :

$s=0.746834\frac{R_{L}}{\lambda_{s}},$ (11)

with $R_{L}=P/Bc$ being the particle Larmor radius. The units are specific in
this formula : as explained in Bieber et al. (1995), the magnetic field $B$ is
in nT, the magnetic fluctuations $\langle b_{s}^{2}\rangle$ in $\rm{nT}^{2}$,
the rigidity $P$ in V, the light speed $c$ in m/s, the correlation length
$\lambda_{s}$ in m, and the final mfp in km. It is valid for rigidities
ranging from 10 MV to 10 GV. We use the formula
$E=\sqrt{(PZ/A)^{2}+E_{0}^{2}}-E_{0}$ to go from rigidity to energy, where $Z$
is the particle charge and $A$ is the mass number. In this study we will
consider only protons with $Z=1$ and $A=1$. Thus the rigidity range between 10
MV and 10 GV is equivalent to energies ranging between 53 keV and 9 GeV.

For the perpendicular diffusion, the formulation that best fits both
observations and data is the nonlinear guiding center (NLGC) theory described
in Bieber et al. (2004). In Shalchi et al. (2004), analytical forms were
derived from NLGC depending on the rigidity of the particle. This formulation
also presents the advantage that the perpendicular diffusion only depends on
the parallel diffusion and the magnetic field properties. Hence, for small
rigidities ($P<10^{2}$ MV):

$\lambda_{\perp}\approx\frac{a^{2}}{2}\frac{\langle
b_{2}^{2}\rangle}{B^{2}}10^{-3}\lambda_{\parallel},$ (12)

with $a^{2}=1/3$ is a numerical factor determined by simulations (Matthaeus et
al., 2003).

For high rigidities ($P>10^{2}$ MV):

$\lambda_{\perp}\approx\left(\frac{2\nu-1}{4\nu}F_{2}(\nu)\lambda_{s}a^{2}\frac{\langle
B^{\prime
2}\rangle}{B^{2}}\frac{2\sqrt{3}}{25}\right)^{2/3}\lambda_{\parallel}^{1/3}\
\rm{m},$ (13)

with the spectral index $\nu=5/6$ and
$F_{2}(\nu)=\sqrt{\pi}\frac{\Gamma(\nu+1)}{\Gamma(\nu+1/2)}$.

In these two formulations, $\langle b_{2}^{2}\rangle$ is in $\rm{nT}^{2}$, $B$
is in nT, $\lambda_{\parallel}$ is in km, $\lambda_{s}$ is in m and the
resulting $\lambda_{\perp}$ is in m.

####  2.2.2 Modeling the turbulence

As the wind model chosen here does not solve yet any equation linked to the
turbulence (see Réville et al. (2020) for a first step in this direction), we
need to model the quantities $\langle b_{s}^{2}\rangle$ and $\lambda_{s}$. We
will use the model described in Chhiber et al. (2017).

To estimate $\langle b_{s}^{2}\rangle$, we combine the two following formulae
linked to the composite geometry used:

$\frac{\langle b_{s}^{2}\rangle}{\langle b_{2}^{2}\rangle}=\frac{1}{4},\langle
b_{2}^{2}\rangle+\langle b_{s}^{2}\rangle=\langle B^{\prime
2}\rangle\Rightarrow\langle b_{s}^{2}\rangle=\frac{1}{5}\langle B^{\prime
2}\rangle,\langle b_{2}^{2}\rangle=\frac{4}{5}\langle B^{\prime 2}\rangle,$
(14)

and use the expression of $\langle B^{\prime 2}\rangle$:

$\langle B^{\prime 2}\rangle=\frac{Z^{\prime 2}}{M_{A}^{\prime}+1}4\pi\rho,$
(15)

with $Z^{\prime 2}=\langle v^{\prime 2}+b^{\prime 2}\rangle$ being the
fluctuation energy and $M_{A}^{\prime}=\langle v^{\prime 2}\rangle/\langle
b^{\prime 2}\rangle$ the Alfvén ratio. From observations (Tu and Marsch,
1995), we can approximate $M_{A}^{\prime}$ by 1 if $r<45R_{\odot}$ and 1/2
beyond $45R_{\odot}$. To model $Z^{\prime 2}$, we can use the Alfvén wave
energy density $\epsilon$ whose expression is $Z^{\prime 2}=2\epsilon/\rho$.
We finally obtain:

$\langle
b_{s}^{2}\rangle=\frac{8}{5}\frac{\epsilon\pi}{M_{A}^{\prime}+1},\langle
b_{2}^{2}\rangle=\frac{32}{5}\frac{\epsilon\pi}{M_{A}^{\prime}+1}.$ (16)

For the Alfvén wave energy density $\epsilon$, we needed an expression that
could adapt to any amplitude and geometry of the magnetic field, as we wanted
to vary these parameters. Instead of solving WKB computations for each case
(similar to Usmanov et al. (2000)), we chose to perform a fit, using an
Alfvén-wave turbulence model (Réville et al., 2020). This model propagates
parallel and anti-parallel Alfvén waves following the WKB theory inside a MHD
wind model similar to what has been described in the previous section. We can
then reconstruct $\epsilon$ using the Elsässer variables $z^{+}$ and $z^{-}$:

$\epsilon=\frac{1}{4}\left((z^{+})^{2}+(z^{-})^{2})\right)\rho.$ (17)

Using the output from this model, we performed a fit of $\epsilon$ depending
on several other physical quantities of the model. According to our tests, the
most relevant quantities are the amplitude of the total magnetic field $B$,
the current $J$ and the poloidal speed $v_{p}$. We obtained the following
formula as best fit:

$\epsilon\approx 0.0104\ B^{1.25}\ J^{0.166}\ v_{p}^{-0.123},$ (18)

with $B$, $J$ and $v_{p}$ in PLUTO code units, so normalized by respectively
$B_{0}=\sqrt{4\pi\rho_{0}}V_{0}$ and $V_{0}=v_{kep}$; the final result
$\epsilon$ is also in code units, and then normalized to fit the amplitude
described in Usmanov et al. (2000) in $\rm{erg}.\rm{cm}^{-3}$. The standard
deviation error associated with each parameter is respectively $1.73\times
10^{-5}$, $2.47\times 10^{-4}$, $2.25\times 10^{-4}$ and $4.01\times 10^{-4}$,
which is reasonable given the parameter values. This fit was obtained for
solar parameters, which means it is valid for magnetic fields between 0.1 and
20 G and for wind speeds between 300 and 800 km/s. We have not tested this fit
for other ranges of input parameters.

Thus we can apply relation 18 to any case and have a turbulence which is self-
consistent with the wind simulation. We show in figure 2 an example of
relation 18 applied to the reference case. The left panel shows the 2D
meridional average of $\epsilon$ in $\rm{erg}.\rm{cm}^{-3}$ from the solar
surface up to 1 AU. The right panel shows the radial profile taken at the
north pole and the latitudinal profile in the northern hemisphere at 1 AU.
Values are very similar to what was found in the case of Usmanov et al. (2000)
for a dipole. The fact that we take into account the current $J$ however
introduces a new effect at the borders of the current sheet, as seen in the
latitudinal profile, with a drop in $\epsilon$ of about 30% near the equator.

Figure 2: 2D (on the left) and 1D (on the right) representation of the Alfvén
wave energy density $\epsilon$ obtained with our fit. For the left panel, we
illustrate the profile of $\epsilon$ for the reference case D1 (see table 2).
For the right panel, we show the radial profile taken at the north pole (top)
and the latitudinal profile taken in the northern hemisphere at 1 AU. Values
are in $\rm{erg}.\rm{cm}^{-3}$.

To estimate $\lambda_{s}$, we use the assumption from Hollweg (1986) which
implies that the correlation length varies as the distance between the
magnetic field lines, which in turn depends on the magnetic field amplitude
(Spruit, 1981), so that $\lambda\propto B^{-1/2}$. Then we use the
relationship $\lambda_{s}=2\lambda_{2}=2\lambda$, observed in the simulations
of Usmanov et al. (2014). Finally :

$\lambda_{s}=4\times 10^{7}B^{-1/2}\ \rm{m}.$ (19)

The proportionality constant is set depending on the boundary conditions to
reach a value of $2\times 10^{7}$ m at the surface of the Sun.

To validate the values provided by these formulae, we have both observational
and computational points of comparison. The observations regarding mfps of
cosmic rays are summarized under the Palmer consensus (Palmer, 1982) : the
values of $\lambda_{\parallel}$ at 1 AU range between 0.08 and 0.3 AU for
rigidities between $0.5$ MV and $5$ GV. For perpendicular diffusion, we can
combine the observations from Chenette et al. (1977) for Jovians electrons and
Burger et al. (2000) for Ulysses measurements of galactic protons : this
yields values for $\lambda_{\perp}$ between 0.003 and 0.01 AU for rigidities
between 1 MV and 10 GV. Numerical simulations have also been performed in
Bieber et al. (2004) for the perpendicular diffusion. Finally another study
described in Chhiber et al. (2017) has used similar methods to derive parallel
and perpendicular mfps for a tilted dipole of amplitude 16 G at the poles. We
discuss in the next section this comparison for various magnetic field
configurations.

##  3 Parametric study

Now that we have a model to compute the parallel and perpendicular mfps from
simulations of a magnetic wind, we will apply this to a parametric study to
understand how various magnetic configurations influence CR diffusion
coefficient. To do so we will focus on 3 cases, described in table 2. Case D1
is the reference case, which parameters were shown in table 1, with a dipole
of amplitude $B_{*}$ of 0.5 G (taken at the surface of the star at the
equator). Case D10 is the same but with a dipole of amplitude $B_{*}$ of 5 G,
so 10 times stronger than D1. Case Q1 is of the same equatorial amplitude as
D1 but with a quadrupolar geometry instead of dipolar. This will allow us to
study separately the effects of amplitude and geometry on the diffusion of
CRs. Finally we will discuss the influence of energy to differentiate SEPs
from GCRs. Because we study in this section only axisymmetric configurations,
we will only focus on the meridional plane corresponding to an azimuthal
average.

Magnetic parameters | D1 | D10 | Q1
---|---|---|---
Geometry | Dipole | Dipole | Quadrupole
Amplitude (equator) | 0.5 G | 5 G | 0.5 G
Table 2: Magnetic field parameters for the 3 cases D1, D10 and Q1 used in the
parametric study. The case D1 corresponds to a dipole of amplitude 0.5 G, the
case D10 to a dipole of amplitude 5 G and the case Q1 to a quadrupole of
amplitude 0.5 G. The amplitude is specified at the surface of the star at the
equator. For the other physical parameters see table 1.

###  3.1 Influence of amplitude

Over an activity cycle, the surface magnetic field amplitude varies between
minima and maxima of activity, typically by a factor between 4 and 10 (DeRosa
et al., 2012). In order to characterize the impact of a variation of
amplitude, we study how the increase by a factor 10 of the whole amplitude
impacts the diffusion of CRs. To do so, we will focus on cases D1 and D10
described in table 2. We will focus here on protons of rigidity 445 MV, which
corresponds to 100 MeV energy, for comparison with previous studies (Pei et
al., 2010; Chhiber et al., 2017).

Figure 3: Meridional cuts of the parallel mfp $\lambda_{\parallel}$ and the
perpendicular mfp $\lambda_{\perp}$ in AU for a proton of rigidity 445 MV
(which means 100 MeV energy). The two panels on the left correspond to the
case D1, the two panels on the left correspond to the case D10. The same color
scale is used for each mfp for both cases.

Figure 3 shows the 2D meridional cuts of the parallel and perpendicular mfps
in units of AU ; we recall that 1 AU is equal to approximately $1.5\times
10^{8}$ km. On the left we have the case D1 and on the right the case D10. We
used the same color scale for each mfp, independently of the case studied, to
allow better comparison between the figures. We see first that the two mfps
behave differently: the perpendicular mfp is nearly constant in the domain
while the parallel mfp increases with distance. We find that for case D1,
$\lambda_{\parallel}$ varies between 0.05 AU and 0.2 AU depending on the
latitude (which corresponds to $4.5\times 10^{7}$ km), while $\lambda_{\perp}$
varies between $10^{-8}$ AU close to the star in the low corona to $10^{-3}$
AU at Earth orbit (between 1.5 and $1.5\times 10^{5}$ km). They also have
opposite behaviors at the heart of current sheets (see equatorial plane) :
$\lambda_{\parallel}$ decreases while $\lambda_{\perp}$ increases (Bieber et
al., 2004; Chhiber et al., 2017). Finally we can notice the same phenomenon on
the edges of the current sheet : this is where $\lambda_{\parallel}$ actually
reaches its highest value and where $\lambda_{\perp}$ loses one order of
magnitude ; however, it seems that $\lambda_{\parallel}$ is more sensitive to
the edges of the current sheet than $\lambda_{\perp}$, as we see more
structures with better contrast. Since the NLGC formulation implies that
$\lambda_{\perp}$ is proportional to $\lambda_{\parallel}$, it is expected to
find similar structures for the two mfps.

From the 2D figure, we can see clearly that the amplitude of the stellar
magnetic field has an effect near the equator on the spread of the current
sheet. In case D1, the mfps are affected at the equator at 1 AU between
$\theta=87\degree$ and $\theta=93\degree$ with $\lambda_{\parallel}$ being
decreased by 75% of its magnitude from 0.3 AU to 0.08 AU, and
$\lambda_{\perp}$ being increased by 2 orders of magnitude from $10^{-4}$ AU
to $10^{-2}$ AU. We recall that in ideal MHD such as for this model, the
thickness of the current sheet is determined by the numerical resistivity,
which naturally increases with distance in the case of a stretched grid like
we have ; this implies that the expansion of the current sheet with the
distance is most likely a numerical effect. In case D10, the zone of influence
of the current sheet is now from $\theta=75\degree$ to $\theta=105\degree$ at
1AU, which is 5 times bigger. The rest of the numerical domain is isotropic.

Figure 4: Radial cuts of the parallel mfp $\lambda_{\parallel}$ (in blue), the
perpendicular mfp $\lambda_{\perp}$ (in red) and the ratio between them
$\lambda_{\perp}/\lambda_{\parallel}$ (in purple) in AU for a proton of
rigidity 445 MV (which means 100 MeV energy). The solid lines are cuts at the
equator ($\theta=90\degree$), the dashed lines are cuts at $\theta=60\degree$
and dotted lines are cuts near the poles ($\theta=4\degree$). The panel on the
left correspond to the case D1. The panel on the right correspond to the case
D10.

Figure 4 shows radial cuts of $\lambda_{\parallel}$ in blue, $\lambda_{\perp}$
in red and the ratio $\lambda_{\perp}/\lambda_{\parallel}$ in purple at three
different co-latitudes ($\theta=90\degree$ in solid line, $\theta=60\degree$
in dashed line and $\theta=4\degree$ in dotted line). We see more clearly that
the amplitude of the two mfps has changed. With $B_{*}$ 10 times stronger,
$\lambda_{\parallel}$ has increased by 60% at all latitudes while
$\lambda_{\perp}$ loses 1 order of magnitude. However the general radial trend
far from the Sun does not change with the amplitude, especially at
$\theta=60\degree$ and $\theta=4\degree$ : at 1 AU, the parallel mfp evolves
as $\propto r^{0.15}$ and the perpendicular mfp evolves as $\propto r^{1.8}$.
Here also we can clearly see the effect of the current sheet spread. For case
D1, the only latitude at which we see a variation in the trend is the equator
(in solid line). We see a bump at $r=2\times 10^{-2}$ AU ($3\times 10^{6}$ km)
where $\lambda_{\parallel}$ decreases by 30% and where $\lambda_{\perp}$
increases of almost 2 orders of magnitude ; the width of the bump is of
$1.5\times 10^{-2}$ AU ($2.25\times 10^{6}$ km). There is also a slight
increase of $\lambda_{\parallel}$ by 30% right before the bump, as well as a
slight decrease of 10% of $\lambda_{\perp}$. For case D10, the equator is
affected, but also the latitude $\theta=60\degree$ because of the edges of the
current sheet which is now much more extended. Here $\lambda_{\parallel}$ has
no decrease bump, only a slight increase at $r=1.75\times 10^{-2}$ AU, while
we see an increase of 2 orders of magnitude for $\lambda_{\perp}$ from
$10^{-7}$ AU to $10^{-5}$ AU (15 km to $1.5\times 10^{3}$ km). Moreover there
is no sudden bump, the variation starts at $r=2\times 10^{-2}$ AU ($3\times
10^{6}$ km) and is visible until $r=0.9$ AU ($1.35\times 10^{8}$ km). To sum
up, the change in amplitude has an impact on the amplitude of the mfps and on
their radial distribution in the current sheet and on its edges.

The ratio $\lambda_{\perp}/\lambda_{\parallel}$ (purple line) is never greater
than 1 in our simulations : for case D1, it is between $10^{-5}$ and $10^{-2}$
; for case D10, it is between $10^{-6}$ and $10^{-3}$. It is the general
behavior expected, except for some specific regions (Dwyer et al., 1997; Zhang
et al., 2003) which could be local structures involved in the diminishing of
SEPs (Zhang et al., 2009). However this simulation does not seem to have that
kind of structure, contrary to Chhiber et al. (2017) where they had a strong
dipole inclined at $30\degree$ ; this kind of structures were visible when
crossing the current sheet. This may be due to our modeling of the turbulence,
which is not self-consistent with the wind model and may not be realistic
enough to describe properly the current sheet.

###  3.2 Influence of geometry

Now we want to focus on the influence of geometry, as the solar cycle goes
from mainly dipolar at minimum to mostly quadrupolar at maximum (DeRosa et
al., 2012). To do so, we will focus on cases D1 and Q1 described in table 2.
We still focus on protons of rigidity 445 MV, which corresponds to 100 MeV
energy.

Figure 5: Meridional cuts of the parallel mfp $\lambda_{\parallel}$ and the
perpendicular mfp $\lambda_{\perp}$ in AU for a proton of rigidity 445 MV
(which means 100 MeV energy). The two panels on the left correspond to the
case D1, the two panels on the left correspond to the case Q1.

Figure 5 shows the 2D meridional cuts of $\lambda_{\parallel}$ and
$\lambda_{\perp}$ in units of AU for case D1 on the left and case Q1 on the
right. What is immediately striking is that the position of the current sheet
is different : for case D1, there is only one current sheet at the equator,
while for case Q1 there are two current sheets at $\theta=60\degree$ and
$\theta=120\degree$. There is thus a new zone delimited by the two current
sheets where the diffusion is different compared to the high latitude regions
: $\lambda_{\parallel}$ is higher and $\lambda_{\perp}$ is slightly lower. We
will quantify this in more details with the next figure.

Figure 6: Radial cuts of the parallel mfp $\lambda_{\parallel}$ (in blue), the
perpendicular mfp $\lambda_{\perp}$ (in red) and the ratio between them
$\lambda_{\perp}/\lambda_{\parallel}$ (in purple) in AU for a proton of
rigidity 445 MV (which means 100 MeV energy). The solid lines are cuts at the
equator ($\theta=90\degree$), the dashed lines are cuts at $\theta=60\degree$
and dotted lines are cuts near the poles ($\theta=4\degree$). The panel on the
left correspond to the case D1. The panel on the right correspond to the case
Q1.

Figure 6 shows radial cuts of $\lambda_{\parallel}$ in blue, $\lambda_{\perp}$
in red and the ratio $\lambda_{\perp}/\lambda_{\parallel}$ in purple at three
different co-latitudes ($\theta=90\degree$ in solid line, $\theta=60\degree$
in dashed line and $\theta=4\degree$ in dotted line). We can see that the
amplitudes are very similar in both cases : $\lambda_{\parallel}$ varies
between $2\times 10^{-2}$ and $2\times 10^{-1}$ AU while $\lambda_{\perp}$
varies between $10^{-6}$ and $10^{-2}$ AU. However for case Q1
$\lambda_{\perp}$ reaches only $10^{-3}$: we can see with figure 5 that it is
less enhanced in the current sheets than in case D1. The bump is also slightly
closer to the Sun at $r=1.5\times 10^{-2}$ AU. The latitudinal distribution is
very different due to the disposition of the current sheets : the variations
observed at $\theta=90\degree$ are now seen at $\theta=60\degree$. The
equatorial plane in case Q1 is more similar to the polar region due to the
fact that it is now between two current sheets. The radial trend remains the
same as described before. Hence the change in geometry is affecting mostly the
latitudinal distribution of CRs diffusion.

So in the end changing amplitude versus changing geometry are not likely to be
equivalent. It means that there is some hope to combine our simulations to CR
transport models to correlate past CR records on Earth with past evolution of
$B_{\odot}$ in terms of amplitude and geometry.

###  3.3 Influence of energy

We will now vary the rigidity/energy of the protons to see how SEPs and GCRs
might be affected differently.

Figure 7: Rigidity dependency of the parallel mfp $\lambda_{\parallel}$ (in
blue) and the perpendicular mfp $\lambda_{\perp}$ (in red) in the equatorial
plane at 1 AU for the case D1 (on the left) and case D10 (on the right). The
gray area corresponds to the Palmer consensus for reference to observational
data for the parallel mean free path (Palmer, 1982).

Figure 7 shows the rigidity dependency of the two mfps between 1 and $10^{5}$
MV, which corresponds to between 1 keV and 100 GeV. They are plotted at 1 AU
in the equatorial plane, on the left for case D1 and on the right for case
D10. For reference, we have also plotted the Palmer consensus range (Palmer,
1982) as a grey rectangle, which corresponds to the range observed in data for
the parallel mean free path only. For case D1, $\lambda_{\parallel}$ goes from
$2\times 10^{-2}$ AU at 1 MV to $200$ AU at $10^{5}$ MV, first with a relation
of $\propto P^{0.33}$, then at $2\times 10^{3}$ MV breaks into the relation
$\propto P^{1.31}$. $\lambda_{\perp}$ has a weaker sensitivity to rigidity,
going from $1\times 10^{-3}$ AU at 1 MV to $2\times 10^{-2}$ AU at $10^{5}$
MV, starting with a relation of $\propto P^{0.33}$, then breaking at also
$2\times 10^{3}$ MV to go with the relation $\propto P^{1.36}$. The break at
$2\times 10^{3}$ MV corresponds to the transition between SEPs and GCRs at 1
GeV (Zank et al., 1998). For case D10, the values of $\lambda_{\parallel}$ are
smaller, going from 0.02 Au to 20 AU, and the values of $\lambda_{\perp}$ are
also smaller, going from $3\times 10^{-4}$ AU to $2\times 10^{-3}$ AU. The
break also happens at a different rigidity, here around $2\times 10^{4}$ MV.
The case D1 is closer to the simulations of Bieber et al. (2004) :
$\lambda_{\parallel}$ goes from $2\times 10^{-1}$ AU at 1 MV to $10$ AU at
$9\times 10^{4}$ MV, while $\lambda_{\perp}$ goes from $2\times 10^{-3}$ AU at
1 MV to $1\times 10^{-2}$ AU at $10^{5}$ MV. Case D10 on the other hand is
closer to the results of Chhiber et al. (2017): $\lambda_{\parallel}$ goes
from $0.05$ AU at 1 MV to $20$ AU at $10^{4}$ MV, while $\lambda_{\perp}$ goes
from $2\times 10^{-3}$ AU at 1 MV to $5\times 10^{-3}$ AU at $10^{5}$ MV. We
see more difference with $\lambda_{\perp}$ because we chose a different
modeling from theirs. This shows that the energy dependency is very sensitive
to the intensity of the magnetic field, especially at 1 AU.

Figure 8: Collection of 1D histograms of the distribution of the parallel mfp
$\lambda_{\parallel}$ (on the left) and the perpendicular mfp
$\lambda_{\perp}$ (on the right) depending on rigidity for the case D1. SEPs
and GCRs are separated using a vertical black line at $P=1.69\times 10^{3}$
MV.

Figure 8 shows a collection of 1D histograms of the distribution of the values
of $\lambda_{\parallel}$ (on the left) and $\lambda_{\perp}$ (on the right)
depending on the rigidity of the proton considered. We consider once again
particles with rigidities between 1 and $10^{5}$ MV, which corresponds to
energies between 1 keV and 100 GeV. Thus, on these histograms, SEPs are
located on the left on GCRs on the right, separated by a vertical black line
at $P=1.69\times 10^{3}$ MV (which corresponds to a 1 GeV energy). These
histograms are made considering all the values of the mfps in the azimuthally-
averaged meridional plane between 1 and 220 solar radii and from one pole to
the other. We see for the parallel mfp that the distribution of values for
SEPs is more centered around small values: around 20% of the distribution lies
between 0.004 and 0.05 AU for rigidities under $10^{3}$ MV, and between 10 and
15% of the distribution lies within slightly bigger values between 0.01 and
0.8 AU. For GCRs, the trend is the same, with 20% of the distribution between
0.05 and 0.2 AU, and between 10 and 15% of the distribution between 0.1 and 1
AU. However, the more energetic is the proton, the more spread is the
distribution of values, with a drift of the tail of the distribution. For the
perpendicular mfp, we can see the weaker dependency on rigidity, as SEPs and
GCRs display almost the same behavior : for SEPs, 8% of the distribution is
between $1\times 10^{-7}$ and $3.0\times 10^{-7}$ AU; for GCRs, 8% of the
distribution is between $3.0\times 10^{-7}$ AU and $5.0\times 10^{-7}$ AU.
Compared to $\lambda_{\parallel}$, the distribution is much more spread
towards higher values of mfps. Our models thus allows us to probe the entire
heliosphere to derive statistical information about both SEPs and GCRs, which
will come useful when computing a realistic count of cosmic rays reaching the
Earth.

##  4 Application to real configurations

Figure 9: 3D visualizations of the wind simulations performed with synoptic
maps from Wilcox observatory. On the left it is the minimum of activity of
October 1995, on the right the maximum of activity of August 1999. The color
scale corresponds to the wind speed in the ecliptic plane. The magnetic field
lines are represented as tubes. The purple surface is the contour of the
neutral current sheet.

We will now use synoptic maps from the Wilcox Observatory to see the 3D
configuration of CR diffusion for real magnetic field configurations
(Hoeksema, 2009). We chose a map corresponding to the minimum of activity
reached in October 1995 (with a SSN of 25, 10 erg of energy in the dipole and
0.1 erg of energy in the quadrupole), and a map corresponding to the maximum
of activity reached in August 1999 (with a SSN of 100, 3 erg of energy in the
dipole and 1 erg of energy in the quadrupole) (DeRosa et al., 2012). The
corresponding configuration of the wind up to 1 AU can be seen in figure 9,
especially the position of the current sheet. More information about the maps
can be found in Réville and Brun (2017) where robust simulations have already
been performed and commented. To use the map in the simulation, we project the
magnetic field on spherical harmonics, and then use the decomposition on the
first 15 $\ell$ degrees to reconstruct the magnetic field. It is initialized
with a potential field source surface (PFSS) method (Altschuler and Newkirk,
1969; Schatten et al., 1969; Schrijver and De Rosa, 2003) to fill the whole
corona, and then the relaxation of the wind modifies the magnetic field
according to the MHD equations. These configurations combine the two effects
discussed earlier : at minimum of activity, the magnetic field is mostly
dipolar with a weaker amplitude ; at maximum of activity, it is mostly
quadrupolar with a stronger amplitude. These simulations are also non-
axisymmetric, which means that we expect large longitudinal variations. Please
note that the fit presented in equation 18 was made based on 2.5D simulations,
so without longitudinal variations. To apply this fit to the 3D case, we
actually apply the fit to each longitude to reconstruct the 3D approximation
of $\epsilon$, which may imply some differences to a fit directly in 3D.

Figure 10: Meridional and equatorial cuts at respectively $\phi=0$ and
$\theta=\pi/2$ of the parallel mfp $\lambda_{\parallel}$ (on the left) and the
perpendicular mfp $\lambda_{\perp}$ (on the right) in AU for a proton of
rigidity 445 MV (which means 100 MeV energy). These panels correspond to a
minimum of solar activity reached in October 1995.

Figure 10 shows the 2D meridional and equatorial cuts at respectively $\phi=0$
and $\theta=\pi/2$ of the parallel (on the left) and perpendicular (on the
right) mfps in units of AU for the minimum of activity. The meridional cuts
are indeed pretty similar to case D1 (cf. Figure 3, left panel) : we have one
current sheet located near the equator which decreases $\lambda_{\parallel}$
and increases $\lambda_{\perp}$. The difference is that here the current sheet
is inclined, which is expected with a 30-degree shift to the southern
hemisphere (Mursula and Hiltula, 2003). The equatorial cut shows that the
distribution is rather isotropic, except at $\phi=135\degree$ and
$\phi=270\degree$ where we can see local decrease of the parallel mfp. These
are due to current sheets as well. This shows that even in minimum of
activity, we need to take into account the 3D structure of the magnetic field.

Figure 11: Meridional and equatorial cuts at respectively $\phi=0$ and
$\theta=\pi/2$ of the parallel mfp $\lambda_{\parallel}$ (on the left) and the
perpendicular mfp $\lambda_{\perp}$ (on the right) in AU for a proton of
rigidity 445 MV (which means 100 MeV energy). These panels correspond to a
maximum of solar activity reached in August 1999.

Figure 11 shows the 2D meridional and equatorial cuts at respectively $\phi=0$
and $\theta=\pi/2$ of the parallel (on the left) and perpendicular (on the
right) mfp in AU for the maximum of activity. Here the configuration is very
multipolar, with more than 2 current sheets visible in the meridional cut. The
equatorial plane also shows 4 current sheets crossing the equatorial plane,
but with less angular spread than at minimum of activity. Depending on the
position of the Earth along its orbit at 1 AU, both the parallel and
perpendicular diffusion are thus very different at maximum of activity, with
favored axis of diffusion.

Figure 12: Histograms of the ratio of the parallel and perpendicular mfp in AU
for a proton of rigidity 445 MV (which means 100 MeV energy) in minimum of
activity (on the left) and maximum of activity (on the right) at
$r=50,100,150,200\ R_{\odot}$. The histograms show the distribution in
percentage on a spherical shell at the corresponding radius.

Figure 12 shows 1D histograms of the distribution of values of the ratio of
$\lambda_{\parallel}$ over $\lambda_{\perp}$ in spherical shells at
$r=25,50,100,150,200\ R_{\odot}$ to show the profile of diffusion at a given
radius. The minimum of activity is on the left, the maximum on the right. As
said before, a general trend observed is that the further we go from the Sun,
the more this ratio diminishes, meaning that the perpendicular mfp becomes
more and more important with the distance to the Sun. For the minimum of
activity for example, at 25 $R_{\odot}$, the peak value is at 10, while at 1
AU, the peak value is around 2. At minimum of activity we see a drift for the
peak value of the distribution : at $25R_{\odot}$, the peak is at 10, at
$50R_{\odot}$ it is at 5, at $100R_{\odot}$ at 4, at $150R_{\odot}$ at 3 and
at $200R_{\odot}$ at 2. The spread of the distribution is also different : the
closer we are to the Sun, the more spread is the distribution with the peak
reaching only 15% of the total distribution; the further we go from the Sun,
the more peaked is the distribution with more than 70% of the values between 1
and 2 at 1 AU. This means that at minimum, there is a distinction of behavior
close and far from the Sun : close to the Sun perpendicular mfp is smaller
than parallel mfp but with a more irregular distribution of values. This means
that the diffusion for SEPs and GCRs is different at minimum, because SEPs
originate from the Sun and GCRs from outside the solar system. Such insights
in the low corona can soon be correlated with the results from Parker Solar
Probe to better understand the diffusion of SEPs in the corona. At maximum of
activity, we observe the same behavior, except that the ratio
$\lambda_{\parallel}/\lambda_{\perp}$ is smaller, with the peak of the
distribution reaching only 4 at $25R_{\odot}$. The distribution is also
slightly more spread, with the peak of the distribution corresponding to only
50% of the total distribution at 1 AU. This means that at maximum, because of
the multiple current sheets, it is more probable to have a high perpendicular
diffusion, and thus it is almost as important as the parallel one 25% of the
time at 1 AU.

##  5 Conclusions

In this study, we combine theory and numerical simulations to better
understand the impact of magnetic field on the propagation of CRs between the
Sun and the Earth. We use a 3D MHD polytropic wind simulation based on the
PLUTO code, which can use either analytical formulations or data from synoptic
maps to initialize the bottom boundary condition magnetic field at the surface
of the star. This provides us realistic inputs for the heliospheric magnetic
field and wind speed. We complete this simulation with approximations and fits
for turbulence parameters such as $\lambda_{s}$ and $\epsilon$. We then use
analytical formulations in post-processing to compute the parallel and
perpendicular mfps associated with a proton of a certain rigidity or energy.
For the parallel mfp, we used the formulation of Zank et al. (1998); for the
perpendicular mfp, we used the formulation of Shalchi et al. (2004) derived
from the work of Bieber et al. (2004).

As the solar cycle evolves in time, the magnetic field is going from weak and
mostly dipolar at minimum of activity, to strong and mostly quadrupolar at
maximum of activity. We decoupled these two effects to study separately the
influence of the amplitude and the geometry of the magnetic field on CR
diffusion. We have thus run 3 cases : case D1 with a weak dipole, case D10
with a strong dipole and case Q1 with a weak quadrupole. The comparison
between cases D1 and D10 shows that the magnetic field amplitude affects the
amplitude of the diffusion coefficients : increasing $B_{*}$ by a factor 10
enhances $\lambda_{\parallel}$ of a factor 10 and decreases $\lambda_{\perp}$
by a factor 100. It also affects the spread of the current sheet, which means
that the mfps present variations on a wider zone for a stronger field; hence
the amplitude also affects the radial distribution of the mfps. However
numerical effects can also enhance this tendency, so this result must be
treated carefully. The comparison between cases D1 and Q1 shows that the
magnetic field geometry does not affect the amplitude of the mfps but changes
the location of the current sheet, which changes the latitudinal distribution
of the mfps with more variations at the equator for the dipole, and more at
$\theta=60\degree$ for the quadrupole. Finally we separated SEPs and GCRs by
changing the energy of the particle. Our tests on cases D1 and D10 shows that
the rigidity/energy dependency of the diffusion is actually very sensitive to
the magnetic field amplitude. GCRs have higher values of mfps than SEPs with a
more spread distribution of values. We have thus demonstrated that amplitude
and geometry have a different impact on SEPs and GCRs, which could be a first
step for retrieving information about the past of the solar-terrestrial
interactions using CRs readings over the last centuries (Finley et al., 2019).

We then applied our post-processing to configurations computed using a
synoptic map. We simulated the state of the inner heliosphere within Earth
orbit corresponding to the minimum of activity of October 1995 and the maximum
of activity of August 1999. This allows us to have for the first time 3D maps
of the CR diffusion between the Sun and the Earth at specific dates. This
shows that the diffusion is highly non-axisymmetric for real configurations.
We also show that for real configurations, there are different behaviors close
and far from the Sun (below and above $150R_{\odot}$), especially at minimum
of activity. We also show that perpendicular diffusion is not to be neglected
at 1 AU, especially at maximum of activity when the current sheet
configuration can be very complex. Thanks to this study, we thus provide
useful and easy tools to derive CR diffusion from any wind simulation.

This study is a first step towards bridging theory and simulations of wind and
CRs. We focused here on the diffusion coefficient, because it is one of the
most difficult term to model in the Parker CR transport equation (Parker,
1965) ; it still needs to be completed with the drift coefficient to take into
account the full influence of the magnetic field. We plan however to go beyond
by modeling the other terms of the equation, which is rather straight-forward
using the wind and magnetic field parameters provided by our 3D MHD
simulation, and finally solving the Fokker-Planck equation. Thus we would be
able to have the spatial distribution of CRs in the heliosphere to compare it
with data from the Earth neutron monitors (Heber and Potgieter, 2006), from
the Moon data (Poluianov et al., 2018) or even Venus (Lorenz and Lawrence,
2015) or Mars data (Lee et al., 2006). In particular it would be interesting
to add the effects of adiabatic cooling as it has a major influence on the
modulation of CRs (Jokipii and Kopriva, 1979). We can also change the
population of particles with our formulations, seeing the effect of positive
or negative charge, in particular in regard to the polarity of the magnetic
field, in order to see the influence of the 22-year cycle on the modulations
of CRs (Heber and Potgieter, 2006). We also plan to apply such formulations to
a more realistic wind model, because as we have explained it, the modeling of
the turbulence presented here has limits; it is suited to have the turbulence
evolving with the MHD model. Such work for a turbulence-based wind model with
Alfvén wave heating is currently being undertaken, see Réville et al. (2020).
CRs propagation through the heliosphere, and especially in the lower corona,
is bound to be an important subject in the years to come thanks to the
combined efforts of Parker Solar Probe and Solar Orbiter; here we focused on
the magnetic field, but the wind structures can also locally affect the
propagation of CRs (Cohen et al., 2020; McComas et al., 2019). Finally, to
fully understand the complete dynamic of GCRs, we would need to extend our
model to go beyond 1 AU, which could be done using more ressourceful numerical
methods such as AMR (Adaptative Mesh Refinement). For space weather purposes,
the model would need to be adapted to be time-dependent like in Kim et al.
(2020).

###### Acknowledgements.

We thank Steve Tobias for the original motivation to perform this study. We
thank Karl-Ludwig Klein and Sophie Masson for useful discussions. This work
was supported by a CEA ’Thèse Phare’ grant, by CNRS and INSU/PNST program, by
CNES SHM funds and by the ERC Synergy grant WholeSun. Computations were
carried out using CEA CCRT and CNRS IDRIS facilities within the GENCI 20410133
allocation, and a local meso-computer founded by DIM ACAV+. We thank Eric
Buchlin and the MEDOC facility for hosting the simulation outputs.

## References

  * Altschuler and Newkirk (1969) Altschuler, M. D., and G. Newkirk, 1969.  Magnetic Fields and the Structure of the Solar Corona. I: Methods of Calculating Coronal Fields.  _Solar Physics_ , 9(1), 131.  10.1007/BF00145734.
  * Arge and Pizzo (2000) Arge, C. N., and V. J. Pizzo, 2000.  Improvement in the prediction of solar wind conditions using near-real time solar magnetic field updates.  _Journal of Geophysical Research_ , 105(A5), 10,465.  10.1029/1999JA000262.
  * Belov et al. (2003) Belov, A. V., E. A. Eroshenko, B. Heber, V. G. Yanke, A. Raviart, R. Müller-Mellin, and H. Kunow, 2003.  Latitudinal and radial variation of >2 GeV/n protons and alpha-particles at solar maximum: Ulysses COSPIN/KET and neutron monitor network observations.  _Annales Geophysicae_ , 21(6), 1295.  10.5194/angeo-21-1295-2003.
  * Bieber et al. (1995) Bieber, J. W., R. A. Burger, and W. H. Matthaeus, 1995.  The Diffusion Tensor throughout the Heliosphere.  In International Cosmic Ray Conference, vol. 4 of _International Cosmic Ray Conference_ , 694.
  * Bieber and Matthaeus (1997) Bieber, J. W., and W. H. Matthaeus, 1997.  Perpendicular Diffusion and Drift at Intermediate Cosmic-Ray Energies.  _The Astrophysical Journal_ , 485(2), 655.  10.1086/304464.
  * Bieber et al. (2004) Bieber, J. W., W. H. Matthaeus, A. Shalchi, and G. Qin, 2004.  Nonlinear guiding center theory of perpendicular diffusion: General properties and comparison with observation.  _Geophysical Research Letters_ , 31(10), L10,805.  10.1029/2004GL020007.
  * Bieber et al. (1994) Bieber, J. W., W. H. Matthaeus, C. W. Smith, W. Wanner, M.-B. Kallenrode, and G. Wibberenz, 1994.  Proton and Electron Mean Free Paths: The Palmer Consensus Revisited.  _The Astrophysical Journal_ , 420, 294.  10.1086/173559.
  * Brun and Browning (2017) Brun, A. S., and M. K. Browning, 2017.  Magnetism, dynamo action and the solar-stellar connection.  _Living Reviews in Solar Physics_ , 14(1), 4.  10.1007/s41116-017-0007-8.
  * Burger et al. (2000) Burger, R. A., M. S. Potgieter, and B. Heber, 2000.  Rigidity dependence of cosmic ray proton latitudinal gradients measured by the Ulysses spacecraft: Implications for the diffusion tensor.  _Journal of Geophysical Research_ , 105(A12), 27,447.  10.1029/2000JA000153.
  * Chenette et al. (1977) Chenette, D. L., T. F. Conlon, K. R. Pyle, and J. A. Simpson, 1977.  Observations of Jovian electrons at 1 AU throughout the 13 month Jovian synodic year.  _The Astrophysical Journal Letters_ , 215, L95–L99.  10.1086/182487.
  * Chhiber et al. (2017) Chhiber, R., P. Subedi, A. V. Usmanov, W. H. Matthaeus, D. Ruffolo, M. L. Goldstein, and T. N. Parashar, 2017.  Cosmic-Ray Diffusion Coefficients throughout the Inner Heliosphere from a Global Solar Wind Simulation.  _The Astrophysical Journal Supplement Series_ , 230(2), 21\.  10.3847/1538-4365/aa74d2, 1703.10322.
  * Cohen et al. (2020) Cohen, C. M. S., E. R. Christian, A. C. Cummings, A. J. Davis, M. I. Desai, et al., 2020.  Energetic Particle Increases Associated with Stream Interaction Regions.  _The Astrophysical Journal Supplement Series_ , 246(2), 20\.  10.3847/1538-4365/ab4c38, 1912.08244.
  * Cummings et al. (1987) Cummings, A. C., E. C. Stone, and W. R. Webber, 1987.  Latitudinal and radial gradients of anomalous and galactic cosmic rays in the outer heliosphere.  _Geophysical Research Letters_ , 14(3), 174.  10.1029/GL014i003p00174.
  * Dedner et al. (2002) Dedner, A., F. Kemm, D. Kröner, C.-D. Munz, T. Schnitzer, and M. Wesenberg, 2002\.  Hyperbolic Divergence Cleaning for the MHD Equations.  _Journal of Computational Physics_ , 175(2), 645.  10.1006/jcph.2001.6961.
  * DeRosa et al. (2012) DeRosa, M. L., A. S. Brun, and J. T. Hoeksema, 2012.  Solar Magnetic Field Reversals and the Role of Dynamo Families.  _The Astrophysical Journal_ , 757(1), 96.  ArXiv: 1208.1768, 10.1088/0004-637X/757/1/96.
  * Dröge (2003) Dröge, W., 2003.  Solar Particle Transport in a Dynamical Quasi-linear Theory.  _The Astrophysical Journal_ , 589, 1027–1039.  10.1086/374812.
  * Dwyer et al. (1997) Dwyer, J. R., G. M. Mason, J. E. Mazur, J. R. Jokipii, T. T. von Rosenvinge, and R. P. Lepping, 1997.  Perpendicular Transport of Low-Energy Corotating Interaction Region–associated Nuclei.  _The Astrophysical Journal_ , 490(1), L115.  10.1086/311003.
  * Einfeldt (1988) Einfeldt, B., 1988.  On Godunov-Type Methods for Gas Dynamics.  _SIAM Journal on Numerical Analysis_ , 25(2), 294.  10.1137/0725021.
  * Engelbrecht and Burger (2015) Engelbrecht, N. E., and R. A. Burger, 2015.  A comparison of turbulence-reduced drift coefficients of importance for the modulation of galactic cosmic-ray protons in the supersonic solar wind.  _Advances in Space Research_ , 55, 390–400.  10.1016/j.asr.2014.09.019, URL http://adsabs.harvard.edu/abs/2015AdSpR..55..390E.
  * Engelbrecht et al. (2017) Engelbrecht, N. E., R. D. Strauss, J. A. le Roux, and R. A. Burger, 2017.  Toward a Greater Understanding of the Reduction of Drift Coefficients in the Presence of Turbulence.  _The Astrophysical Journal_ , 841, 107.  10.3847/1538-4357/aa7058, URL http://adsabs.harvard.edu/abs/2017ApJ...841..107E.
  * Finley et al. (2019) Finley, A. J., S. Deshmukh, S. P. Matt, M. Owens, and C.-J. Wu, 2019.  Solar Angular Momentum Loss over the Past Several Millennia.  _The Astrophysical Journal_ , 883, 67.  10.3847/1538-4357/ab3729.
  * Forman (1977) Forman, M. A., 1977.  The velocity correlation function in cosmic-ray diffusion theory.  _Astrophysics and Space Science_ , 49, 83–97.  10.1007/BF00647077.
  * Forman et al. (1974) Forman, M. A., J. R. Jokipii, and A. J. Owens, 1974.  Cosmic-Ray Streaming Perpendicular to the Mean Magnetic Field.  _The Astrophysical Journal_ , 192, 535–540.  10.1086/153087.
  * Giacalone and Jokipii (1999) Giacalone, J., and J. R. Jokipii, 1999.  The Transport of Cosmic Rays across a Turbulent Magnetic Field.  _The Astrophysical Journal_ , 520, 204–214.  10.1086/307452.
  * Goldstein (1976) Goldstein, M. L., 1976.  A nonlinear theory of cosmic-ray pitch-angle diffusion in homogeneous magnetostatic turbulence.  _The Astrophysical Journal_ , 204, 900–919.  10.1086/154239.
  * Grappin et al. (2010) Grappin, R., J. Léorat, S. Leygnac, and R. Pinto, 2010.  Search for a self-consistent solar wind model.  In M. Maksimovic, K. Issautier, N. Meyer-Vernet, M. Moncuquet, and F. Pantellini, eds., Twelfth International Solar Wind Conference, vol. 1216 of _American Institute of Physics Conference Series_ , 24–27.  10.1063/1.3395848.
  * Green (1951) Green, M. S., 1951.  Brownian Motion in a Gas of Noninteracting Molecules.  _Journal of Chemical Physics_ , 19, 1036–1046.  10.1063/1.1748449.
  * Guo and florinski (2014) Guo, X., and V. florinski, 2014.  Galactic Cosmic-Ray Modulation near the Heliopause.  _The Astrophysical Journal_ , 793(1), 18.  10.1088/0004-637X/793/1/18, URL https://iopscience.iop.org/article/10.1088/0004-637X/793/1/18.
  * Hathaway (2015) Hathaway, D. H., 2015.  The Solar Cycle.  _Living Reviews in Solar Physics_ , 12(1), 4.  10.1007/lrsp-2015-4.
  * Heber et al. (1996) Heber, B., W. Droege, P. Ferrando, L. J. Haasbroek, H. Kunow, R. Mueller-Mellin, C. Paizis, M. S. Potgieter, A. Raviart, and G. Wibberenz, 1996.  Spatial variation of $>40$MeV/n nuclei fluxes observed during the ULYSSES rapid latitude scan.  _Astronomy and Astrophysics_ , 316, 538–546.
  * Heber and Potgieter (2006) Heber, B., and M. S. Potgieter, 2006.  Cosmic Rays at High Heliolatitudes.  _Space Science Reviews_ , 127(1-4), 117.  10.1007/s11214-006-9085-y.
  * Hoeksema (2009) Hoeksema, J. T., 2009.  Evolution of the large-scale magnetic field over three solar cycles.  _Proc. IAU_ , 5(S264), 222–228.  10.1017/S1743921309992675.
  * Hollweg (1986) Hollweg, J. V., 1986.  Transition region, corona, and solar wind in coronal holes.  _Journal of Geophysical Research_ , 91(A4), 4111.  10.1029/JA091iA04p04111.
  * Jokipii (1966) Jokipii, J. R., 1966.  Cosmic-Ray Propagation. I. Charged Particles in a Random Magnetic Field.  _The Astrophysical Journal_ , 146, 480.  10.1086/148912.
  * Jokipii and Kopriva (1979) Jokipii, J. R., and D. A. Kopriva, 1979.  Effects of particle drift on the transport of cosmic rays. III - Numerical models of galactic cosmic-ray modulation.  _The Astrophysical Journal_ , 234, 384–392.  10.1086/157506, URL http://adsabs.harvard.edu/abs/1979ApJ...234..384J.
  * Jokipii and Levy (1977) Jokipii, J. R., and E. H. Levy, 1977.  Effects of particle drifts on the solar modulation of galactic cosmic rays.  _The Astrophysical Journal_ , 213, L85.  10.1086/182415, URL http://adsabs.harvard.edu/doi/10.1086/182415.
  * Jokipii and Parker (1970) Jokipii, J. R., and E. N. Parker, 1970.  on the Convection, Diffusion, and Adiabatic Deceleration of Cosmic Rays in the Solar Wind.  _The Astrophysical Journal_ , 160, 735.  10.1086/150465.
  * Jokipii and Thomas (1981) Jokipii, J. R., and B. Thomas, 1981.  Effects of drift on the transport of cosmic rays. IV - Modulation by a wavy interplanetary current sheet.  _The Astrophysical Journal_ , 243, 1115.  10.1086/158675, URL http://adsabs.harvard.edu/doi/10.1086/158675.
  * Keppens and Goedbloed (1999) Keppens, R., and J. P. Goedbloed, 1999.  Numerical simulations of stellar winds: polytropic models.  _Astronomy and Astrophysics_ , 343, 251.
  * Kim et al. (2020) Kim, T. K., N. V. Pogorelov, C. N. Arge, C. J. Henney, S. I. Jones-Mecholsky, et al., 2020.  Predicting the Solar Wind at Parker Solar Probe Using an Empirically Driven MHD Model.  _The Astrophysical Journal Supplement Series_ , 246(2), 40\.  ArXiv: 1912.02397, 10.3847/1538-4365/ab58c9, URL http://arxiv.org/abs/1912.02397.
  * Kubo (1957) Kubo, R., 1957.  Statistical-Mechanical Theory of Irreversible Processes. I.  _Journal of the Physical Society of Japan_ , 12, 570–586.  10.1143/JPSJ.12.570.
  * Lee et al. (2006) Lee, K., L. Pinsky, V. Andersen, C. Zeitlin, T. Cleghorn, F. Cucinotta, P. Saganti, W. Atwell, and R. Turner, 2006.  Helium cosmic ray flux measurements at Mars.  _Radiation Measurements_ , 41, 1123–1125.  10.1016/j.radmeas.2006.07.011.
  * Lionello et al. (2001) Lionello, R., J. A. Linker, and Z. Mikić, 2001.  Including the Transition Region in Models of the Large-Scale Solar Corona.  _The Astrophysical Journal_ , 546(1), 542.  10.1086/318254.
  * Lockwood and Webber (2005) Lockwood, J. A., and W. R. Webber, 2005.  Intensities of galactic cosmic rays of 1.5 GV rigidity at Earth versus the heliospheric current sheet tilt.  _Journal of Geophysical Research: Space Physics_ , 110(A4), A04102.  10.1029/2004JA010880, URL https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2004JA010880.
  * Lorenz and Lawrence (2015) Lorenz, R. D., and D. J. Lawrence, 2015.  Gamma rays and cosmic rays at Venus: The Pioneer Venus gamma ray detector and considerations for future measurements.  _Planetary and Space Science_ , 109-110, 129–134.  10.1016/j.pss.2015.02.009.
  * Luo et al. (2013) Luo, X., M. Zhang, H. K. Rassoul, N. V. Pogorelov, and J. Heerikhuisen, 2013.  Galactic Cosmic-Ray Modulation in a Realistic Global Magnetohydrodynamic Heliosphere.  _The Astrophysical Journal_ , 764(1), 85.  10.1088/0004-637X/764/1/85, URL https://iopscience.iop.org/article/10.1088/0004-637X/764/1/85.
  * Manuel et al. (2011) Manuel, R., S. E. S. Ferreira, M. S. Potgieter, R. D. Strauss, and N. E. Engelbrecht, 2011.  Time-dependent cosmic ray modulation.  _Advances in Space Research_ , 47(9), 1529–1537.  10.1016/j.asr.2010.12.007, URL http://www.sciencedirect.com/science/article/pii/S0273117710007945.
  * Matt and Pudritz (2008) Matt, S., and R. E. Pudritz, 2008.  Accretion-powered Stellar Winds. II. Numerical Solutions for Stellar Wind Torques.  _The Astrophysical Journal_ , 678(2), 1109.  10.1086/533428.
  * Matthaeus et al. (1990) Matthaeus, W. H., M. L. Goldstein, and D. A. Roberts, 1990.  Evidence for the presence of quasi-two-dimensional nearly incompressible fluctuations in the solar wind.  _Journal of Geophysical Research_ , 95, 20,673.  10.1029/JA095iA12p20673.
  * Matthaeus et al. (2003) Matthaeus, W. H., G. Qin, J. W. Bieber, and G. P. Zank, 2003.  Nonlinear Collisionless Perpendicular Diffusion of Charged Particles.  _The Astrophysical Journal_ , 590(1), L53.  10.1086/376613.
  * McComas et al. (2019) McComas, D. J., E. R. Christian, C. M. S. Cohen, A. C. Cummings, A. J. Davis, et al., 2019.  Probing the energetic particle environment near the Sun.  _Nature_ , 576, 223–227.  10.1038/s41586-019-1811-1.
  * McComas et al. (2003) McComas, D. J., H. A. Elliott, N. A. Schwadron, J. T. Gosling, R. M. Skoug, and B. E. Goldstein, 2003.  The three-dimensional solar wind around solar maximum.  _Geophysical Research Letters_ , 30(10), 1517.  10.1029/2003GL017136.
  * McKibben (1998) McKibben, R. B., 1998.  Three-dimensional Solar Modulation of Cosmic Ray and Anomalous Components in the Inner Heliosphere.  _Space Science Reviews_ , 83, 21.  10.1023/A:1005044707585.
  * Mignone et al. (2007) Mignone, A., G. Bodo, S. Massaglia, T. Matsakos, O. Tesileanu, C. Zanni, and A. Ferrari, 2007.  PLUTO: A Numerical Code for Computational Astrophysics.  _The Astrophysical Journal Supplement Series_ , 170(1), 228\.  10.1086/513316.
  * Moffatt (1978) Moffatt, H. K., 1978.  Magnetic field generation in electrically conducting fluids.  Cambridge University Press.
  * Mursula and Hiltula (2003) Mursula, K., and T. Hiltula, 2003.  Bashful ballerina: Southward shifted heliospheric current sheet.  _Geophysical Research Letters_ , 30(22), 2135.  10.1029/2003GL018201.
  * Neugebauer and Snyder (1962) Neugebauer, M., and C. W. Snyder, 1962.  Solar Plasma Experiment.  _Science_ , 138(3545), 1095.  10.1126/science.138.3545.1095-a.
  * Oughton et al. (2017) Oughton, S., W. H. Matthaeus, and P. Dmitruk, 2017.  Reduced MHD in Astrophysical Applications: Two-dimensional or Three-dimensional?  _The Astrophysical Journal_ , 839(1), 2.  10.3847/1538-4357/aa67e2, URL https://iopscience.iop.org/article/10.3847/1538-4357/aa67e2.
  * Oughton et al. (2015) Oughton, S., W. H. Matthaeus, M. Wan, and K. T. Osman, 2015.  Anisotropy in solar wind plasma turbulence.  _Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences_ , 373(2041), 20140,152.  Publisher: Royal Society, 10.1098/rsta.2014.0152, URL http://royalsocietypublishing.org/doi/10.1098/rsta.2014.0152.
  * Owens and Forsyth (2013) Owens, M. J., and R. J. Forsyth, 2013.  The Heliospheric Magnetic Field.  _Living Reviews in Solar Physics_ , 10(1), 5.  10.12942/lrsp-2013-5.
  * Palmer (1982) Palmer, I. D., 1982.  Transport coefficients of low-energy cosmic rays in interplanetary space.  _Reviews of Geophysics and Space Physics_ , 20, 335.  10.1029/RG020i002p00335.
  * Parker (1958) Parker, E. N., 1958.  Dynamics of the Interplanetary Gas and Magnetic Fields.  _The Astrophysical Journal_ , 128, 664.  10.1086/146579.
  * Parker (1964) Parker, E. N., 1964.  The Scattering of Charged Particles by Magnetic Irregularities.  _Journal of Geophysical Research_ , 69, 1755–1758.  10.1029/JZ069i009p01755.
  * Parker (1965) Parker, E. N., 1965.  The passage of energetic charged particles through interplanetary space.  _Planetary and Space Science_ , 13(1), 9.  10.1016/0032-0633(65)90131-5.
  * Parker (1993) Parker, E. N., 1993.  A Solar Dynamo Surface Wave at the Interface between Convection and Nonuniform Rotation.  _The Astrophysical Journal_ , 408, 707.  10.1086/172631.
  * Pei et al. (2010) Pei, C., J. W. Bieber, B. Breech, R. A. Burger, J. Clem, and W. H. Matthaeus, 2010.  Cosmic ray diffusion tensor throughout the heliosphere.  _Journal of Geophysical Research (Space Physics)_ , 115(A3), A03103.  10.1029/2009JA014705.
  * Perri et al. (2018) Perri, B., A. S. Brun, V. Réville, and A. Strugarek, 2018.  Simulations of solar wind variations during an 11-year cycle and the influence of north-south asymmetry.  _Journal of Plasma Physics_ , 84(5), 765840,501.  10.1017/S0022377818000880.
  * Pinto and Rouillard (2017) Pinto, R. F., and A. P. Rouillard, 2017.  A Multiple flux-tube Solar Wind Model.  _The Astrophysical Journal_ , 838(2), 89.  10.3847/1538-4357/aa6398.
  * Poluianov et al. (2018) Poluianov, S., G. A. Kovaltsov, and I. G. Usoskin, 2018.  Solar energetic particles and galactic cosmic rays over millions of years as inferred from data on cosmogenic 26Al in lunar samples.  _Astronomy and Astrophysics_ , 618, A96.  10.1051/0004-6361/201833561.
  * Reames (1999) Reames, D. V., 1999.  Particle acceleration at the Sun and in the heliosphere.  _Space Science Reviews_ , 90, 413.  10.1023/A:1005105831781.
  * Riley et al. (2015) Riley, P., R. Lionello, J. A. Linker, E. Cliver, A. Balogh, et al., 2015.  Inferring the Structure of the Solar Corona and Inner Heliosphere During the Maunder Minimum Using Global Thermodynamic Magnetohydrodynamic Simulations.  _The Astrophysical Journal_ , 802(2), 105.  10.1088/0004-637X/802/2/105.
  * Ruffolo et al. (2012) Ruffolo, D., T. Pianpanit, W. H. Matthaeus, and P. Chuychai, 2012.  Random Ballistic Interpretation of Nonlinear Guiding Center Theory.  _The Astrophysical Journal Letters_ , 747(2), L34.  10.1088/2041-8205/747/2/L34.
  * Réville and Brun (2017) Réville, V., and A. S. Brun, 2017.  Global Solar Magnetic Field Organization in the Outer Corona: Influence on the Solar Wind Speed and Mass flux Over the Cycle.  _The Astrophysical Journal_ , 850(1), 45.  10.3847/1538-4357/aa9218.
  * Réville et al. (2015) Réville, V., A. S. Brun, S. P. Matt, A. Strugarek, and R. F. Pinto, 2015.  The Effect of Magnetic Topology on Thermally Driven Wind: Toward a General Formulation of the Braking Law.  _The Astrophysical Journal_ , 798(2), 116.  10.1088/0004-637X/798/2/116.
  * Réville et al. (2020) Réville, V., M. Velli, O. Panasenco, A. Tenerani, C. Shi, et al., 2020.  The Role of Alfvén Wave Dynamics on the Large-scale Properties of the Solar Wind: Comparing an MHD Simulation with Parker Solar Probe E1 Data.  _The Astrophysical Journal Supplement Series_ , 246(2), 24\.  10.3847/1538-4365/ab4fef.
  * Sakurai (1985) Sakurai, T., 1985.  Magnetic stellar winds: a 2-D generalization of the Weber-Davis model.  _Astronomy and Astrophysics_ , 152, 121.
  * Schatten et al. (1969) Schatten, K. H., J. M. Wilcox, and N. F. Ness, 1969.  A model of interplanetary and coronal magnetic fields.  _Solar Physics_ , 6(3), 442.  10.1007/BF00146478.
  * Schrijver and De Rosa (2003) Schrijver, C. J., and M. L. De Rosa, 2003.  Photospheric and heliospheric magnetic fields.  _Solar Physics_ , 212(1), 165.  10.1023/A:1022908504100.
  * Shalchi (2006) Shalchi, A., 2006.  Extended nonlinear guiding center theory of perpendicular diffusion.  _Astronomy and Astrophysics_ , 453(3), L43–L46.  10.1051/0004-6361:20065465, URL http://www.aanda.org/10.1051/0004-6361:20065465.
  * Shalchi (2009) Shalchi, A., 2009.  Nonlinear Cosmic Ray Diffusion Theories, vol. 362.  Springer.  10.1007/978-3-642-00309-7.
  * Shalchi (2017) Shalchi, A., 2017.  Time-dependent perpendicular transport of energetic particles in magnetic turbulence with transverse complexity.  _Physics of Plasmas_ , 24, 050,702.  10.1063/1.4982805, URL http://adsabs.harvard.edu/abs/2017PhPl...24e0702S.
  * Shalchi (2020) Shalchi, A., 2020.  Perpendicular Transport of Energetic Particles in Magnetic Turbulence.  _Space Science Reviews_ , 216(2), 23.  10.1007/s11214-020-0644-4, URL https://doi.org/10.1007/s11214-020-0644-4.
  * Shalchi et al. (2004) Shalchi, A., J. W. Bieber, and W. H. Matthaeus, 2004.  Analytic Forms of the Perpendicular Diffusion Coefficient in Magnetostatic Turbulence.  _The Astrophysical Journal_ , 604(2), 675–686.  10.1086/382128.
  * Snyder et al. (1963) Snyder, C. W., M. Neugebauer, and U. R. Rao, 1963.  The Solar Wind Velocity and Its Correlation with Cosmic-Ray Variations and with Solar and Geomagnetic Activity.  _Journal of Geophysical Research_ , 68, 6361.  10.1029/JZ068i024p06361.
  * Spruit (1981) Spruit, H. C., 1981.  Magnetic flux tubes.  In S. Jordan, ed., The Sun as a Star, Monograph Series on Nonthermal Phenomena in Stellar Atmospheres, vol. 450, 385–413. NASA Special Publications.
  * Strugarek et al. (2015) Strugarek, A., A. S. Brun, S. P. Matt, and V. Reville, 2015.  Numerical Aspects of 3D Stellar Winds.  In 18th Cambridge Workshop on Cool Stars, Stellar Systems, and the Sun, vol. 18 of _Cambridge Workshop on Cool Stars, Stellar Systems, and the Sun_ , 589–602.  1410.3537.
  * Suzuki et al. (2013) Suzuki, T. K., S. Imada, R. Kataoka, Y. Kato, T. Matsumoto, H. Miyahara, and S. Tsuneta, 2013.  Saturation of Stellar Winds from Young Suns.  _Publications of the Astronomical Society of Japan_ , 65, 98.  10.1093/pasj/65.5.98.
  * Taylor (1922) Taylor, G. I., 1922.  The Motion of a Sphere in a Rotating Liquid.  _Proceedings of the Royal Society of London Series A_ , 102, 180–189.  10.1098/rspa.1922.0079.
  * Tu and Marsch (1995) Tu, C.-Y., and E. Marsch, 1995.  Magnetohydrodynamic Structures Waves and Turbulence in the Solar Wind - Observations and Theories.  _Space Science Reviews_ , 73(1-2), 1.  10.1007/BF00748891.
  * Tóth et al. (2012) Tóth, G., B. van der Holst, I. V. Sokolov, D. L. De Zeeuw, T. I. Gombosi, et al., 2012.  Adaptive numerical algorithms in space weather modeling.  _Journal of Computational Physics_ , 231(3), 870.  10.1016/j.jcp.2011.02.006.
  * Usmanov et al. (2000) Usmanov, A. V., M. L. Goldstein, B. P. Besser, and J. M. Fritzer, 2000.  A global MHD solar wind model with WKB Alfvén waves: Comparison with Ulysses data.  _Journal of Geophysical Research_ , 105(A6), 12,675.  10.1029/1999JA000233.
  * Usmanov et al. (2014) Usmanov, A. V., M. L. Goldstein, and W. H. Matthaeus, 2014.  Three-fluid, Three-dimensional Magnetohydrodynamic Solar Wind Model with Eddy Viscosity and Turbulent Resistivity.  _The Astrophysical Journal_ , 788(1), 43.  10.1088/0004-637X/788/1/43.
  * Wang and Sheeley (1990) Wang, Y.-M., and N. R. Sheeley, 1990.  Magnetic flux Transport and the Sunspot-Cycle Evolution of Coronal Holes and Their Wind Streams.  _The Astrophysical Journal_ , 365, 372.  10.1086/169492.
  * Webber et al. (2005) Webber, W. R., B. Heber, and J. A. Lockwood, 2005.  Time variations of cosmic ray electrons and nuclei between 1978 and 2004: Evidence for charge-dependent modulation organized by changes in solar magnetic polarity and current sheet tilt.  _Journal of Geophysical Research: Space Physics_ , 110(A12).  10.1029/2005JA011291, URL https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2005JA011291.
  * Weber and Davis (1967) Weber, E. J., and L. Davis, 1967.  The Angular Momentum of the Solar Wind.  _The Astrophysical Journal_ , 148, 217.  10.1086/149138.
  * Wiengarten et al. (2016) Wiengarten, T., S. Oughton, N. E. Engelbrecht, H. Fichtner, J. Kleimann, and K. Scherer, 2016.  A GEeneralized Two-Component Model of Solar Wind Turbulence and ab initio Diffusion Mean-Free Paths and Drift Lengthscales of Cosmic Rays.  _The Astrophysical Journal_ , 833(1), 17.  10.3847/0004-637X/833/1/17, URL https://iopscience.iop.org/article/10.3847/0004-637X/833/1/17.
  * Zank et al. (1998) Zank, G. P., W. H. Matthaeus, J. W. Bieber, and H. Moraal, 1998.  The radial and latitudinal dependence of the cosmic ray diffusion tensor in the heliosphere.  _Journal of Geophysical Research_ , 103(A2), 2085.  10.1029/97JA03013.
  * Zhang et al. (2003) Zhang, M., J. R. Jokipii, and R. B. McKibben, 2003.  Perpendicular Transport of Solar Energetic Particles in Heliospheric Magnetic Fields.  _The Astrophysical Journal_ , 595(1), 493.  10.1086/377301.
  * Zhang et al. (2009) Zhang, M., G. Qin, and H. Rassoul, 2009.  Propagation of Solar Energetic Particles in Three-Dimensional Interplanetary Magnetic Fields.  _The Astrophysical Journal_ , 692(1), 109–132.  10.1088/0004-637X/692/1/109.