Modeling Transport of Energetic Particles in Corotating Interaction Region -- A case Study

TL;DR

This study models energetic particle transport in Corotating Interaction Regions using a case study from 2008, combining observations and Monte Carlo simulations to understand particle modulation and shock acceleration.

Contribution

It introduces a method to estimate shock locations and particle intensities in CIRs by integrating in-situ data with focused transport simulations assuming no cross-field diffusion.

Findings

Significant modulation of sub-MeV/nucleon particles observed.

Particle acceleration is more efficient at larger heliocentric distances.

Monte Carlo simulations match observations when assuming increased particle acceleration farther from the Sun.

Abstract

We investigate energetic particle transport in Corotating Interaction Regions (CIRs) through a case study. The CIR event we study occurred on $2008$ February $08$ and was observed by both the Advanced Composition Explorer (ACE) and the twin Solar TErrestrial RElations Observatory (STEREO)-B spacecraft. An in-situ reverse shock was observed by STEREO-B ($1.0$ AU) but not ACE ($0.98$ AU). Using STEREO-B observations and assuming the CIR structure does not vary significantly in the corotating frame, we estimate the shock location at later times for both the STEREO-B and ACE observations. Further assuming the accelerated particle spectral shape at the shock does not vary with shock location, we calculate the particle differential intensities as observed by ACE and STEREO-B at two different times by solving the focused transport equation using a Monte-Carlo simulation. We assume that particles move along Parker's field and experience no cross-field diffusion. We find that the modulation of sub-MeV/nucleon particles is significant. To obtain reasonable comparisons between the simulations and the observations by both ACE and STEREO-B, one has to assume that the CIR shock can accelerate more particles at a larger heliocentric distance than at a smaller heliocentric distance.

Tables

Table 1. Table 1: Monte Carlo Results

	Shock Loc	${\left(\delta B/B\right)}^2$	Spacecraft	$\alpha$	$\eta$
Case I	2.55 AU	0.012	STB	6	1.27
Case II	1.39 AU	0.02	ACE	1	0.61

Table 2. Table 2: Fisk & Lee Fitting

	Shock Loc	${\kappa }_0$	Spacecraft	$J_0$ ( # cm-2 s-2 sr-1 (MeV/n)-1)	$\eta$
Source	1.0 AU	0.03	STB	1227	1.0
Case I	3.63 AU	0.03	STB	5474	3.4
Case II	1.73 AU	0.03	ACE	4167	4.7

Equations

$$v_{sh}=\frac{\beta v_{comp}-v_{sw}}{\beta -1}.$$

$$\frac{CE}{DE}=\frac{v_{sw}}{v_{sw}-v_{sh}}.$$

$$\frac{CD}{DE}=\frac{CE-DE}{DE}=1.92,$$

$$DE=v_{sw}\Delta t,$$

$$\begin{array}{r}\frac{\partial f}{\partial t}+(u_i+v\mu b_i)\frac{\partial f}{\partial x_i}+[\frac{1-3{\mu }^2}{2}{b}_i{b}_j\frac{\partial u_i}{\partial x_j}-\frac{1-{\mu }^2}{2}\frac{\partial u_j}{\partial x_j}-\frac{\mu b_i}{v}(\frac{\partial \mu }{\partial t}+\frac{\partial u_i}{\partial x_j})]v\frac{\partial f}{\partial v}\\ +\frac{1-{\mu }^2}{2}[v\frac{\partial b_i}{\partial x_i}+\mu \frac{\partial u_i}{\partial x_i}-3\mu b_i{b}_j\frac{\partial u_j}{\partial x_i}-\frac{2b_i}{v}(\frac{\partial u}{\partial t}+u_j\frac{\partial u_i}{\partial x_j})]\frac{\partial f}{\partial \mu }=\frac{\partial }{\partial \mu }{D}_{\mu \mu }\frac{\partial f}{\partial \mu }+S-L.\end{array}$$

$$\begin{array}{rl}{B}_r=& B_0(\frac{R_0}{r}{)}^2,\\ B_{\theta }=& 0,\text{mbox}and\\ B_{\varphi }=& B_r\frac{\Omega r\sin \theta }{u}(r>>R_0)\end{array}$$

$$\delta \mu =sign({ϵ}^{\prime }-\frac{1}{2})er{f}^{-1}(ϵ)\sqrt{4D_{\mu \mu }\delta t}+\frac{\partial D_{\mu \mu }}{\partial t}\delta t,$$

$$D_{\mu \mu }=\frac{1}{2}\frac{<(\delta \mu {)}^2>}{\delta t}=\frac{\pi }{4}(1-{\mu }^2){\Omega }_o\frac{kP(k)}{B_o^2},$$

$$P(k)=A_q{\lambda }_c(\delta B{)}^2(1+(k{\lambda }_c{)}^q{)}^{-1}$$

$${\int }_{k_s}^{k_L}P(k)dk=(\delta B{)}^2.$$

$$\delta B^2(r)=\delta B^2(1\mathrm{AU})(\frac{r}{1\mathrm{AU}}{)}^{-3.5}$$

$$\eta =\frac{{\int }_{E_{\min }}^{E_{\max }}dE{J}_{sh}}{{\int }_{E_{\min }}^{E_{\max }}dE{J}_{1\mathrm{AU}}}$$

$$\alpha =\eta \frac{A_{sh}}{A_{1AU}}$$

$$J=J_0(\frac{R}{R_s}{)}^{2/(\beta -1)+V/({\kappa }_{0}v)}{v}^n{e}^{-v/v_0}$$

$$\beta =(\frac{V_{up}^2+{\Omega }^2{R}_s^2}{V_{dn}^2+{\Omega }^2{R}_s^2})\frac{B_{dn}}{B_{up}}$$

$$<\dot{p}>=-\frac{p}{3}\nabla \cdot v_{sw}$$

$$\begin{array}{rl}\dot{p}& =-p{v}_{sw}[\frac{sec\psi }{2L(z)}(1-{\mu }^2)+\cos \psi \frac{d}{dr}sec\psi {\mu }^2]\\ \dot{\mu }& =\frac{v}{2L(z)}[1+\mu \frac{v_{sw}}{v}sec\psi -\mu \frac{v_{sw}v}{c^2}sec\psi ](1-{\mu }^2)-v_{sw}(cos\psi \frac{d}{dr}sec\psi )\mu (1-{\mu }^2)\end{array}$$

$$\begin{array}{rl}{\mathbf{V}}_1=& \Omega R_1{\mathbf{\theta }}_1,\\ {\mathbf{V}}_2=& \Omega R_2{\mathbf{\theta }}_2,\end{array}$$

$${\mathbf{V}}_1-{\mathbf{V}}_2=\Omega R_1{\mathbf{\theta }}_1-\Omega R_2{\mathbf{\theta }}_2.$$

$$\begin{array}{rl}& ({\mathbf{V}}_1-{\mathbf{V}}_2)\cdot (R_1{\mathbf{\theta }}_1-R_2{\mathbf{\theta }}_2)\\ =& \Omega R_1^2{\mathbf{\theta }}_1\cdot {\mathbf{r}}_1+\Omega R_2^2{\mathbf{\theta }}_2\cdot {\mathbf{r}}_2-\Omega R_1{R}_2{\mathbf{\theta }}_1\cdot {\mathbf{r}}_2-\Omega R_2{R}_1{\mathbf{\theta }}_2\cdot {\mathbf{r}}_1.\end{array}$$

$$\begin{array}{rl}{\mathbf{r}}_1\cdot {\mathbf{\theta }}_1& =0,\\ {\mathbf{r}}_2\cdot {\mathbf{\theta }}_2& =0,\\ {\mathbf{r}}_1\cdot {\mathbf{\theta }}_2& =\sin (\Delta \psi ),\\ {\mathbf{r}}_2\cdot {\mathbf{\theta }}_1& =-\sin (\Delta \psi ).\end{array}$$

$$A{D}^2+A{B}^2-2ADAB\cos \alpha =B{D}^2.$$

$$AD=\frac{r\theta \cos \alpha }{1-(\sin \alpha /\cos \psi {)}^2}(1+\tan \alpha tan\psi ).$$

$$\begin{array}{rl}tan\alpha & =\frac{CD}{AC}=\frac{V_{sh}dt}{AC},\\ cot\psi & =\frac{CE}{AC}=\frac{V_{sw}dt}{AC},\end{array}$$

Full text

\setkeys

Gindraft=false \authorrunningheadLULU ZHAO ET AL. \titlerunningheadTRANSPORT OF ENERGETIC PARTICLES IN CIRs \authoraddrCorresponding author: Gang Li, Department of Space Science and CSPAR, University of Alabama in Huntsville, Huntsville, Alabama, USA. (gang.li@uah.edu)

11affiliationtext: Department of Space Science and CSPAR, University of Alabama in Huntsville, Huntsville, Alabama, USA.22affiliationtext: Southwest Research Institute, San Antonio, TX, USA33affiliationtext: Applied Physics Laboratory, Johns Hopkins University, Laurel, MD, USA44affiliationtext: Institute of Space Sciences and School of Space Science and Physics, Shandong University at Weihai, Weihai, China

{article}

1 Introduction

Corotating Interaction Regions (CIRs) are formed when fast solar wind, originating from coronal holes that extend to low latitudes, compresses slow solar wind in front of it. A forward shock can form at the leading edge of the compression region and propagate into the slow solar wind, and a reverse shock can form at the trailing edge of the compression region and propagate into the fast solar wind. Early studies showed that by $3$ to $4$ Astronomical Units (AU), most CIRs are bounded by a forward and reverse shock pair (e.g. Hundhausen and Gosling, 1976). Recently, a statistical study by Jian et al. (2006) found that $31\%$ of CIRs observed at 1 AU are associated with shocks.

CIRs are a major source of energetic particles in the inner heliosphere during solar minimum (e.g. van Hollebeke et al., 1978; McDonald et al., 1976; Richardson et al., 1993). Fisk and Lee (1980) first examined particle acceleration associated with CIRs. They solved the steady state transport equation with a geometry appropriate to CIRs. In their model, particles are accelerated via the first order Fermi acceleration mechanism at either the forward or the reverse shock which are often at a distance of several AU. Energetic particles then propagate along the interplanetary magnetic field back to $1$ AU. Similar to Galactic Cosmic Rays (GCRs) and Anomalous Cosmic Rays (ACRs), the adiabatic cooling can lead to a modulation effect at low energies.

Later, Giacalone et al. (2002) proposed another mechanism which, instead of invoking shocks, considered particle acceleration in gradual and slowly-varying solar wind compression regions. They found that particles can be accelerated up to $10$ MeV/nucleon by a process similar to diffusive shock acceleration, in which particles gain energy by scattering between converging scattering centers. Simulation results similar to observations were obtained for a reasonable set of parameters.

Recent observations by Mason et al. (1997) suggested that the observed energy spectra in CIR events often do not show modulation effects in low energies. Instead, the spectra continue to fall down as a power law from the energy threshold of the instrument, typically $\sim 30~{}$keV/nucleon, and even solar wind energies (Chotoo et al. 2000). More recently, Ebert et al. (2012) studied $73$ CIR-associated suprathermal He intensity enhancements and found that the peak of these sub-MeV He intensities correlate well with the arrival of the compression region trailing edge. These and other recent works (e.g. (Mason et al., 2009; Bučík et al., 2009; Ebert et al., 2012b)) suggest that sub-MeV energetic particles associated with CIRs at 1 AU may be accelerated locally.

While sub-MeV particles may be accelerated locally, higher energy particles suffer less modulation and may well be accelerated beyond 1 AU (e.g. $3$ to $5$ AU) and propagate along Parker field back to $1$ AU. Therefore it is possible that energetic particles observed at 1 AU come from more than one radial distance. This also implies that cross-field diffusion may be important to understand CIR observations. Indeed, there is some evidence that cross-field diffusion may be important in a few CIRs (Dwyer et al., 1997). Because cross field diffusion can drastically reduce the distance a particle travels from several AUs to reach Earth’s orbit, it is therefore more important for low energy particles.

In this work, we study the transport of energetic particles at CIRs using a numerical model and compare our model with an earlier analytical model of Fisk and Lee (1980). To quantify the modulation effect along a single Parker field, we assume particles are tied to single field lines. Since there is no cross-field diffusion, particles on different field lines originate from different locations along the reverse shock. As a case study, we examine one CIR event that occurred between $2008$ February $08$ and $14$. This event was observed by the twin Solar TErrestrial RElations Observatory (STEREO) spacecraft and the Advanced Composition Explorer (ACE) spacecraft. Observations at STEREO-B (hereafter STB) showed an in-situ signature of a reverse shock, which was not seen at ACE and STEREO-A (hereafter STA). This in-situ observation of the shock by STB allows us to estimate the shock location at a later time for ACE and STB observations upstream of the shock.

We further assume the shape of the accelerated particle spectrum at the shock does not vary with heliocentric distance. Note that this is only a working hypothesis. With this hypothesis, we can decouple the acceleration process from the transport process. This hypothesis is bound to introduce some errors in our analysis. However, as we will see, the shock in our event was quite strong at 1 AU, having a compression ratio $\beta=2.17$. Inteplanetary shocks often do not have compression ratios larger than, say, $3$. The shock accelerated particle spectrum for a shock with $\beta=2.17$ is $\sim p^{-5.5}$ while that with a compression ratio $\beta=3.0$ is $\sim p^{-4.5}$. As we will see, the effect of transport (modulation at low energy) is much larger than the uncertainty introduced in this hypothesis. Indeed, considering the ratio of the differential intensities $dJ/dE$(2 MeV)/$dJ/dE$(0.2 MeV) (roughly the energy range for our problem), then the ratio of using $\sim p^{-4.5}$ is larger by a factor of $\sim 3.2$ than that using $\sim p^{-5.5}$. In comparison, the effect of modulation on this ratio is more significant. Therefore, our hypothesis, while crude, is somewhat justified. Also note that the shock location in our event is not too far from 1 AU. For the two periods we consider, it locates at $1.55$ AU and $2.55$ AU, respectively. The plasma parameters at these two shock locations may not differ too much from that at 1 AU and the corresponding compression ratios may be similar to $2.17$. In any case, as we focus on the effects of the transport of energetic particles in this work, we here make the simplest assumption about the shape of the accelerated particle spectrum, that is that it does not vary with heliocentric distance.

To examine the transport process, we use a Monte-Carlo simulation. In particular, we simulate particle spectra at two different times for STB and ACE observations and compare our simulation results with observations. In our simulation, we ignore the gyration degree of freedom of the particles and solve the focused transport equation along a Parker field. Comparing to (Fisk and Lee, 1980), who assumed the motion of particles are diffusive in spatial coordinates and are described by a diffusion coefficient $\kappa$, we retain explicitly the particle pitch angle in our approach. Consequently, we consider explicitly the magnetic focusing effect and the pitch angle diffusion (described by $D_{\mu\mu}$) in our formalism. The Fisk and Lee (1980) approach is justified when the solar wind turbulence is strong. However, when solar wind turbulence is weak, employing $D_{\mu\mu}$ and treating the focusing effect explicitly is more appropriate. To summarize, using a case study where in-situ observation of a CIR-shock is available, we improve upon the Fisk and Lee model by employing a more sophisticated numerical simulation that treats the transport of CIR-associated energetic particles by a focused transport equation.

Our paper is organized as follows. In section 2 we present the multi-spacecraft observations of this event; we discuss the simulation model in section 3 and compare the simulation results with observations in section 4. We conclude in section 5.

2 Event Analysis

The energetic particle measurements we present here are from the Suprathermal Ion Telescope (SIT) (Mason et al., 2008) onboard the STA and STB spacecraft, and the Ultra-low Energy Isotope Spectrometer (ULEIS) (Mason et al., 1998) onboard ACE. The solar wind plasma measurements are obtained from the PLAsma and SupraThermal Ion Composition investigation (PLASTIC) (Galvin et al., 2008) onboard STA & B, and the Solar Wind Electron Proton Alpha Monitor (SWEPAM) (McComas et al., 1998) onboard ACE. We also use the magnetic field data from the Magnetometers on STA & B (Acuña et al., 2008) and ACE (Smith et al., 1998).

The CIR event we study was first observed at STB starting on $2008$ February $08$. It has been reported as Event $25$ by Mason et al. (2009). It has also been studied by Bučík et al. (2011).

Figure 1 shows the in-situ observations (STB, ACE and STA from left to right). The uppermost panel shows the $0.14$-$2.2$ MeV/n He time intensity profiles between $2008$ February $07$ and $2008$ February $17$; the lower panels show, in descending order, the solar wind proton speed, density, temperature and the total magnetic field (magnitude) in the same time period. For the STB observations, the He intensity started to increase at $18:00$ UT, $2008$ February $08$ and reached its maximum at $19:36$ UT, $2008$ February $09$ with a clear peak associated with it. At the time of the peak, discontinuities in the proton speed, magnetic field and temperature, were observed, indicating this time as the passage of the CIR reverse shock. This shock was identified in the STEREO shock list at http://www-ssc.igpp.ucla.edu/\~jlan/STEREO/Level3/STEREO\_Level3_Shock.pdf, maintained by Dr. J. Lan. The in-situ shock passage as observed by STB is shown as the red dashed line in the left panel. The blue and green dashed lines correspond to two observation periods for which simulations were performed. The shock spectrum is obtained by integrating the differential intensity in the period $17:54$ UT to $20:09$ UT, February $09$ and is shown by the red curve in Figure 2. In comparison, the spectrum for a later time period between $13:55$ UT to $15:50$ UT, February $11$ was shown as the blue curve in Figure 2. The modulation at low energies can be clearly seen for this period.

Because the coronal structures from which fast solar wind originated from often evolve slowly, CIRs are considered as a steady state structure in the co-rotating frame to the first order (Mason et al., 2009). One therefore expects to observe the same reverse shock at ACE after $\sim 1.6$ days. Indeed, the He intensity observed by ACE began to increase gradually and reached its maximum at $22:30$ UT, February $10$. However, the proton speed, magnetic field and total pressure showed no discontinuities, so no reverse shock was observed by ACE. This is not surprising since time variations of plasma properties such as density, speed, etc. can cause the location of the shock to vary. The energy spectrum from ACE observations for the period between $09:59$ UT $\sim$ $12:00$ UT, February $11$ is shown as the green curve in Figure 2. During this time, ACE is magnetically connected to the shock at a location not far from 1 AU. Consequently, there is little modulation at low energies.

2.1 Event Geometry

Figure 3 shows the relative locations of three spacecraft at $20:00$ UT, 2008 February $09$, just before the energetic He intensities began to increase on STB. STA, STB and Earth are represented by red (labeled as $A$), blue (labeled as $B$) and green solid circles respectively. The separation angle between STB and Earth (ACE) is $23.656^{\circ}$ and the heliocentric distance of STB is $1.0$ AU and that of Earth is $0.987$ AU.111http://stereo-ssc.nascom.nasa.gov/cgi-bin/make_where_gif Since ACE lies at the L1 Lagrange point between the Sun and the Earth, its heliocentric distance is $0.977$ AU.

The configuration of magnetic field lines as well as the CIR reverse shock with respect to STB & ACE is shown in the cartoon in Figure 4. The Sun is located at $O$. The black half circle is the trajectory of STB and the dashed half circle is the trajectory of ACE in the co-rotating frame. Green, blue and red curves are Parker field lines. The black thick curve extending from $A$ to $D$ represents the CIR reverse shock surface. Point $A$ ($B$) is the intersection between the blue (red) magnetic field line and the STB orbit. Point $H$ ($I$) is the intersection between the green (blue) magnetic field line and the ACE orbit.

When STB is at location $A$, it connects to the local reverse shock, corresponding to observations at $19:36$ UT, $2008$ February $09$. As STB rotates from $A$ to $B$ in the co-rotating frame, its connection point at the shock (i.e. the acceleration site) moves further out. When STB is at location $B$, under the assumption of no cross-field diffusion, it will observe energetic particles that are accelerated by the reverse shock at point $D$.

As shown in Figure 1, the intensity of STB He gradually decreases from $19:36$ UT, February $09$ to $13:55$ UT, February $11$. This decrease is due to two reasons. First, as the acceleration sites moves out, it becomes harder and harder for accelerated particles to propagate back to 1 AU. This is the modulation effect. Second, the density of seed particles, therefore the intensity of the accelerated particles may also vary with heliocentric distance. Note, the composition of the seed particle is still presently under debate. In modeling SEP events, where ions and electrons are accelerated at a CME-driven shock that propagates out from the Sun, it has been argued and assumed that a fraction of the solar wind (e.g. $1\%$-$4\%$ as often used in Li et al. (2003); Li and Zank (2005)) is accelerated at the CME-driven shock. CIR shocks, however, are often quasi-perpendicular shocks and have higher injection energies than quasi-parallel CME-driven shocks. So it is hard for bulk solar wind to be accelerated. Recently, Mason et al. (2012b), through a study of the abundance of 3He and He+, have suggested that the seed particles for CIR shocks is the suprathermal ion pool rather than the bulk solar wind. If the seed particle was the solar wind, then the density of the seed particles will have an $r^{-2}$ dependence. However, if the seed particle was the suprathermal ions, then the radial dependence of the seed particle may be more complicated. For example, the cooling of solar wind may imply faster decay of the seed particles than $r^{-2}$. On the other hand, some continuous particle-wave interaction in the solar wind may suggest a slower decay than $r^{-2}$. Indeed, using Cassini observations, Hill et al. (2009) examined how the intensities of $2$-$60$ keV/nuc suprathermal He*++* vary with $r$. They found that the intensities of suprathermal He*++* decrease slower than $r^{-2}$ (see Figure 2 of (Hill et al., 2009)).

In our event, because STB observed the shock in-situ at 1 AU, if we assume the shock is a steady structure in the co-rotating frame and if shock parameters do not change significantly along the shock surface, then with a particle transport model, we can examine how the intensity of the accelerated particles at the shock varies with radial distance.

From in-situ observation of STB, we find that the average solar wind speed in the fast compression region between 12:00 $\sim$ 18:00 UT, February $09$ is $v_{comp}=620$ km s*-1* and in the fast solar wind right after the shock passage is $v_{sw}=760$ km s*-1* (see Figure 1). Correspondingly, the angle between the upstream magnetic field with the radial direction is $\sim 30^{\circ}$. The reverse shock has a $\theta_{BN}=51^{\circ}$ (see the online list http://www-ssc.igpp.ucla.edu/\~jlan/STEREO/Level3/STEREO_Level3\_Shock.pdf), so the angle between the shock normal and the radial direction, $\alpha$, is $20^{\circ}$. If we approximate the shock normal to be along the radial direction, then we can estimate the shock speed $v_{sh}$ at $1$ AU from

[TABLE]

Here $\beta=2.17$ is the compression ratio of the reverse shock as obtained from in-situ magnetic field data. The calculated shock speed is $500$ km s*-1* in the s/c frame. Note that $v_{sh}<v_{sw}$ since the reverse shock propagates towards the Sun in the fast solar wind frame. Assuming $v_{sh}$ does not vary with $r$, we can obtain the radial distance $OD$ from,

[TABLE]

Using Equation 2, with the above $v_{sw}$ and $v_{sh}$, we find

[TABLE]

and the length of $DE$ can be calculated by

[TABLE]

where $\Delta t=1.84$ days is the time difference between points $A$ and $B$. So the distance between $C$ and $D$ is $1.55$ AU. Therefore, for STB observation at point $B$, the energetic particles are accelerated at point $D$ on the reverse shock, having a heliocentric distance $r_{sh}=2.55$ AU. We remark that the equation 1 is only an approximation and it is only valid when $\alpha$ is small (note that $\cos 20^{\circ}=0.94$). If the shock normal differs from the radial direction substantially (e.g. $\alpha>30^{\circ}$, one can not use equation 1.

Now consider the ACE observations. Given that the sidereal rotation period of the Sun is $24.47$ days and the angle between STB and ACE is $\sim 23.656^{\circ}$, it takes $1.61$ days for ACE to be connected to the exact same portion of the shock structure as STB, implying the enhancement of particle intensity in ACE will be delayed by $1.61$ days from STB. From the observations shown in Figure 1, we find a time delay of the peak intensity between STB and ACE to be $1.14$ days. This difference of $d\tau$ between the observation and by considering solar rotation can be due to a number of reasons. For example, the ACE and the STB have different latitudes, so ACE and STB did not see the same portion of the shock. Finally, the CIR shock may not be completely stationary and its location may vary with time (similar to e.g., the heliospheric termination shock). The effect of this time variation is illustrated in Figure 4. In the cartoon, the dashed thick curve depicts the CIR reverse shock when ACE is at point $H$.

Since ACE did not observe the shock in-situ, it is hard to discern the exact cause of $d\tau$. In the following, we do not consider possible oscillations of the CIR shock and assume that it is given by the thick solid curve as in Figure 4. Following the analysis as in the STB case, we find that for a later time of ACE observation at $10:00$ UT February $11$, the source location at the shock has a $r_{sh}$ of $1.39$ AU.

We perform two numerical simulations for two $2$-hour periods corresponding to STB observations at $13:55$-$15:55$ UT February $11$ and ACE at $10:00$-$12:00$ UT February $11$. The periods are chosen in which the energetic particle intensities have small variations. We note that STA also observed this event about 2 days later than ACE. However, the intensity profile from STA observation showed clear differences from STB and ACE. This can be due to the fact that assumptions of a steady-state CIR reverse shock is only applicable within a short duration. As shown in Mason et al. (2009) (their Figure 6), energetic particle observations of the same CIR event with a few days apart can differ substantially. In our event, the separation between STA and STB is $\sim 45^{\circ}$, $>3$ days apart, and we do not consider STA observations in this study.

As noted before, we assume the shock is in a steady state in the co-rotating frame and that the shock strength does not vary along the shock surface. Assuming the shock parameters do not vary as a function of heliocentric distance allows us to use the 1 AU in-situ shock spectrum from STB as a reference for energetic particle spectra at different times. Therefore, as noted before, we can decouple the transport process from the acceleration process for CIR-associated energetic particles.

Specifically, we will use the in-situ STB 1 AU shock observation of the energetic particle differential current intensity $J(E)$ and scale it by a factor of $\eta$ (defined below) at two different locations at the CIR reverse shock, to calculate the resulting differential intensities as observed at 1 AU by STB and ACE. By comparing simulations with observations we can obtain how $\eta$, which is a measure of the seed particle intensity, varies with $r$. The implication of this radial dependence of $\eta$ is discussed in section 4.

3 Model Description

We use a Monte-Carlo code to study the transport of charged particles. The transport of energetic particles can be described by the focused transport equation (e.g. Skilling, 1971; Isenberg, 1997):

[TABLE]

In the above, the distribution function $f$ is assumed to be gyrotropic; the $b_{i}$ and $u_{i}$ are the components of the unit magnetic field and solar wind speed. $\mu$ is the particle’s pitch angle cosine. $S$ and $L$ are the source and loss terms. The focused transport equation can be derived from the guiding center theory (le Roux and Webb, 2009). The particle’s motion therefore can be regarded as the motion of the guiding center along the background magnetic field line plus a diffusion in pitch angle, governed by $D_{\mu\mu}$.

In the solar wind, the background magnetic field is given by

[TABLE]

where $\theta=90^{\circ}$ corresponds to the magnetic field in the ecliptic plane. In this work, we use $B_{0}=1.53$ Gauss at $R_{0}=1$ $R_{s}$ (solar radius).

Energetic particles are followed in two different frames (see Ruffolo, 1995; Agueda et al., 2005; Kocharov et al., 2003). One is the instantaneous co-rotating frame, the other is the instantaneous solar wind frame. The instantaneous co-rotating frame is a frame which co-rotates with the Sun. The instantaneous solar wind frame is a local inertial frame that co-moves with the solar wind at the particle’s location. In the instantaneous co-rotating frame, the solar wind velocity is $v_{sh}^{co}=v_{sw}/\cos\psi$, where $\psi$ is the angle between the magnetic field direction and the radial direction. The solar wind velocity is then parallel to the local magnetic field, yielding a zero induced electric field. Therefore, the particle’s energy is conserved in this frame. Furthermore, due to the focusing effect, the particle’s pitch angle will change because of the conservation of the particle’s magnetic moment $p_{\perp}^{2}/B$. At each time step, the particle’s location and pitch angle are updated in this instantaneous co-rotating frame. We then transform the particle’s momentum to the instantaneous solar wind frame. This is the frame where the effect of solar wind Magnetohydrodynamic (MHD) turbulence is considered since MHD waves and turbulence are generated locally. Pitch angle diffusion is considered in this frame. After considering pitch angle diffusion, we transform back to the instantaneous co-rotating frame. Finally, we need to transform the particle’s energy and momentum between two instantaneous co-rotating frames at two different times $t$ and $t$ + $\delta t$. As shown in Appendix A, to the order of $(v/c)^{2}$, the particle’s energy is conserved under the Lorentz transformation between these two frames.

In modeling the pitch angle diffusion, we follow Zhao and Li (2014). For each time step $\delta t$, the particle’s pitch angle $\mu$ will change by a small amount $\delta\mu$ given by

[TABLE]

here $\epsilon^{\prime}$ and $\epsilon$ are random numbers uniformly distributed in $[0,1]$. $erf^{-1}$ is the inverse of error function. In Quasi-linear theory (QLT) (Jokipii, 1971; Lee and Lerche, 1974; Wentzel, 1974), $D_{\mu\mu}$ is given by

[TABLE]

where $P(k)$ (given in Equation 9) is the turbulence power spectrum in solar wind. $\Omega_{o}=eB/(\gamma m)$ is particle’s gyrofrequency. $k=\Omega_{o}(v|\mu|)^{-1}$ is particle’s resonant wave number.

[TABLE]

$\lambda_{c}$ is the correlation length and $q$ is the power law index of the turbulence spectrum. $A_{\beta}$ is determined by the normalization condition

[TABLE]

The correlation length $\lambda_{c}$, smallest ($k_{S}$) and largest ($k_{L}$) wave number in $P(k)$ are set to be the typical values at $1$ AU (see Zhao and Li, 2014).

The radial dependence of $B$ is well defined for a Parker field. For $\delta B$, however, the radial dependence is still presently unknown. Some earlier studies by Bruno and Carbone (2005) suggested that a WKB approximation, i.e., $\delta B^{2}\sim r^{-3}$ may be reasonable but slightly underestimated. Following Mason et al. (2012a), we use $\delta B^{2}\sim r^{-3.5}$ in this work, so

[TABLE]

We assume the shape of the source energetic particle spectrum $f(E)\sim J(E)/v$ does not vary with $r$, and is given by a broken power law form as shown in Figure 5. The form of $f(E)$ is chosen to fit the observed differential intensity $J(E)$ at STB immediately downstream of the shock between $17:54\sim 20:09$ UT, February $09$. In Figure 5, the low-energy portion of $f(E)$ (shown in blue) has $f(E)\sim E^{-3.57}$ between $0.094<E/({\rm MeV/nucleon})<0.546$, and the high-energy portion (shown in red) has $f(E)\sim E^{-4.89}$ between $0.546<E/({\rm MeV/nucleon})<2.185$.

For our simulation, we inject a total of $N_{0}=60,000$ protons at location of $r_{sh}$ (which is also the outer boundary) with an initial pitch angle cosine $\mu$ uniformly distributed between $-1$ and [math]. We do not consider the acceleration process in this work. So these particles are injected with a given spectrum. We assume this spectrum, up to a factor, is the same as that from the in-situ observation of the shock at STB. In this way, we attempt to decouple the acceleration process from the transport process. Note that close to the shock, the acceleration process is governed by the diffusion cofficient, which can be substantially smaller than that in the interplanetary medium. So the acceleration process is quite different from the transport process. We release all particles at the same time and follow all of them for a period of $10$ days and obtain the time-integrated spectrum. We then follow the transport of these particles from the shock to 1 AU. Particles leave the simulation domain when they reach either the inner boundary ($r=0.01$ AU) or the outer boundary ($r=r_{sh}$). The particles’ momenta and pitch angles are recorded when they pass $r=1$ AU.

4 Results and Discussion

The simulated proton differential intensities are shown in Figure 6 as the black curves with “diamond” symbols for two cases. The left panel is for case I, corresponding to the STB observation between $13:55$ UT $\sim 15:50$ UT, February $11$. In this case, energetic protons are injected at a heliocentric distance $r=2.55$ AU and the particle differential current is obtained at $1.0$ AU. The observed differential current intensity is shown as the blue curve with the “plus” sign. The right panel is for case II. In this case, energetic protons are injected at a heliocentric distance $r=1.39$ AU and the particle differential intensity is obtained at $0.98$ AU. This case corresponds to the ACE observation between $10:00$ UT $\sim 12:00$ UT, February $11$. The observed differential intensity is shown as the green curve with the “plus” sign. In both cases, we use a two-hour observation window so that enough statistics can be obtained. We assume the shock location does not change during this $2$-hour interval. We do not consider a longer period ($>$ 2 hours) because the acceleration site to which the spacecraft is connected to can change rapidly along the shock surface.

To fit the observed differential intensities at 1 AU in Figure 6, a key parameter is $(\delta B/B)^{2}$ at 1 AU. We calculate $(\delta B/B)^{2}$ using the STB/MAG 1-minute data and the ACE/MAG 4-minute data. In the case of STB observation, the average magnetic field was obtained for a period of 20 hours before the 2-hour observation window and we find $(\delta B/B)^{2}=0.010$. For the ACE observation, the background magnetic field before the 2-hour observation window showed a clear decreasing trend, so we have chosen a 20-hour period after the 2-hour observation window to calculate the average magnetic field and $(\delta B/B)^{2}$. We find $(\delta B/B)^{2}=0.019$ in this case. The magnetic field data and the periods for calculating the average background magnetic field and $(\delta B/B)^{2}$ for both the ACE and STB observations are shown in figure 7. In the figure, the two left panels are for STB and the two right panels are for ACE. In both cases, the upper panels show the total $B$ for an extended period covering the whole event and the lower panels are zoom-in plots of the total magnetic field and $(\delta B/B)^{2}$ for a 20-hour period before (for the case of STB) or after (for the case of ACE) the 2-hour observation window of energetic particles.

Now consider the fitting of STB observation in Figure 6. We vary the ratio of $(\delta B/B)^{2}$ at $1$ AU from the observed in-situ value of $0.01$ to fit the shape of the observed differential intensity at 1 AU. Note that the observed $(\delta B/B)^{2}=0.01$ is an ensemble average of many radially propagating plasma parcels that pass through the spacecraft during a $20$ hour window. These plasma parcels do not consist of a Parker field line on which energetic particles propagate from $2.55$ AU to $1$ AU. So it is only a proxy of the turbulence level along the Parker field line particles propagate on. Nevertheless, the best fit yields a $(\delta B/B)^{2}=0.012$ at $1$ AU, very close to the in-situ observation. The simulated and observed differential current intensities for STB observation agree nicely, as can be seen from the left panel of Figure 6. Next consider the ACE observation. Again, by varying the ratio of $(\delta B/B)^{2}$ at $1$ AU, the best fit yields $(\delta B/B)^{2}=0.02$, close to the in-situ observation of $(\delta B/B)^{2}=0.019$ at ACE. This is shown in the right panel of Figure 6.

After fitting the spectral shape, we next obtain the normalization factor $\eta$ of the injected particles from the fitting. The parameter $\eta$ is defined through,

[TABLE]

where $J_{sh}$ is the differential intensity at the shock ($r=1.39$ AU for the ACE observation and $r=2.55$ AU for the STB observation) and $J_{1{\rm AU}}$ is the differential current when the shock is observed in-situ at 1 AU by STB between $17:54$ $\sim$ $20:09$ UT, February $09$. The parameter $\eta$ reflects how the number of seed particles, i.e., the particles that participated in the shock acceleration process, depends on $r$. In the simulation, we inject particles to a flux tube whose cross section depends on the shock location. Because the footpoint of the flux tube on the shock surface varies with heliocentric distance, we introduce an additional parameter $\alpha$ through,

[TABLE]

where $A_{sh}$ is the cross section of the flux tube that intersects with the shock at $r$, and $A_{\rm 1AU}$ corresponds to the in-situ STB observation of the shock at 1 AU. The factor $A_{sh}/A_{\rm 1AU}$ reflects how the flux tube expands as a function of $r$. Knowing the solar wind speed and the shock speed, the ratio $A_{sh}/A_{1AU}$ can be readily calculated as shown in Appendix B. For our cases, $A_{1.39AU}/A_{1AU}=1.65$, and $A_{2.55AU}/A_{1AU}=4.72$. From the fitting we obtain $\alpha=6$ for case I and $\alpha=1$ for case II. Therefore we obtain $\eta_{\rm caseI}=1.27$ and $\eta_{\rm caseII}=0.61$. Our Simulation results are summarized in Table 1.

These values of $\eta$’s are very important results. Consider first the ACE observation with $r_{sh}=1.39{\rm AU}$. In this case we find that to fit the observation, the number of accelerated particles at $r_{sh}$ is $0.6$ times that at $1$ AU. This decrease yields a radial dependence of $r^{-1.55}$. This is shallower than the radial dependence of solar wind density $r^{-2}$, and is consistent with Hill et al. (2009) where the intensities of suprathermal He*++* was found to decrease slower than $r^{-2}$. Therefore the ACE observation supports the notation that the seed particles for CIR-associated energetic particles are most likely suprathermal ions than the bulk solar wind (Mason et al., 2012b).

Comparing to the ACE observation, the STB observation with $r_{sh}=2.55{\rm AU}$ shows that $J_{sh}(r=2.55{\rm AU})/J_{sh}(r=1{\rm AU})=1.27$. So the number density of the accelerated particles at $2.55$ AU has to be larger than that at 1 AU. This implies that the intensity of the seed particles increases with $r$, instead of decreasing with $r$. It contradicts with the results in (Hill et al., 2009). However, we note that the data points in (Hill et al., 2009) are $0.5$-$1$ AU apart, and most are far beyond the Earth orbit, so they may not provide much constraint about the relatively small shock distances studied here. A larger seed population further out than 1 AU may seem counter-intuitive. However, earlier works (McDonald et al., 1976; van Hollebeke et al., 1978) have shown that CIR intensities increase beyond 1 AU and peak at several AU. Note that in these earlier studies, the enhancement of CIR intensity is largely attributed to the fact that CIR shocks tend to form beyond 1 AU (e.g. Dwyer et al. (1999)). In our case, however, the CIR shock was seen to form in-situ at 1 AU by STB. Note that we have assumed the shock parameters do not vary with the heliocentric distance of the shock. It is possible that the shock strengthens beyond 1 AU in that the injection energy decreases with r, so that more particles can participate in the shock acceleration process. This, of course, still implies that the seed population (i.e. particles participating the shock acceleration process) increases with $r$.

This increase of seed particle intensity at $2.55$ AU is required because the modulation effect for low energy particles is significant. This can be seen from the following arguments. Under the assumption that particles do not diffuse across the field lines, the path length from $2.55$ ($1.39$) AU to $1$ AU is $2.26$ ($0.52$) AU (assuming a solar wind speed of $760$ $km~{}s^{-1}$). If there is no scattering ($\delta B/B_{0}\rightarrow 0$), then as a charged particle moving towards the Sun along the Parker field, the focusing effect (conservation of particle’s magnetic moment) will tend to reverse the particle’s momentum direction. With turbulence included, the pitch angle will undergo both a focusing and a scattering process. To compete with the focusing effect, one may think that a larger $\delta B/B_{0}$ will help. However, if $\delta B/B_{0}$ is too large so that pitch angle scattering dominates focusing, then the motion of particles along the field can be regarded as a diffusion and the time to arrive 1 AU from the shock will be very long. For example, if the mean free path of a $0.5$ MeV proton is $\sim 0.5$ AU (see e.g. equation (27) of Li et al. (2003)), it takes $\sim 10$ mean free paths for an $0.5$ MeV proton to arrive 1 AU from $2.55$ AU, translating to a time period of $21$ hours. Since the adiabatic cooling rate is $\nabla\cdot v_{sw}$, the longer the propagation time, the more deceleration a particle will experience, leading to a stronger modulation. If the mean free path (mfp) $\lambda$ is smaller, since the travel time scales as $\lambda^{-1}$, the adiabatic cooling will lead to a stronger modulation. Finally, for particles of smaller energies, using equation (27) of Li et al. (2003) (with $\alpha=1/3$), one finds that the travel time scales roughly as $E^{-2/3}$, so the modulation will be even stronger.

The above discussion illustrates why a significant number of seed particles has to be present when the shock is further out. Note that one implicit assumption in our scenario is that particles are tied to a single magnetic field line. So we have ignored $\kappa_{\perp}$. The picture can be changed, and perhaps substantially, if we allow particles to diffuse across field lines. This can be seen easily from Figure 4. Within a small longitudinal range we see that the shock goes from $1$ AU (point A) to $2.55$ AU (point D). While it is difficult for low energy particles accelerated at point D ($2.55$ AU) to propagate along the Parker spiral to point B ($1$ AU), if the low energy particles could cross-field diffuse from field lines that connect to the shock at a closer distance than $2.55$ AU (e. g. field lines to the left of $BD$) then their propagation distance would be much smaller . Consequently, if cross-field diffusion plays an important role, then we do not need to require as large a seed population further out. We note that if particles could cross-diffuse in to, e.g. the field line $BD$ as shown in Figure 4, particles on the field line of $BD$ could also cross-diffuse out. However, since low energy particles accelerated at point $D$ will have a hard time to propagate to point $B$, so whether they stay in the same field line or diffuse out to other field lines does not matter much to observations at point $B$.

It is instructive to fit the differential intensities for both case I and case II using the Fisk & Lee model. Following Fisk and Lee (1980), one can show that the upstream differential intensity $J(E)$ is,

[TABLE]

where

[TABLE]

In the above $\beta$ is the compression ratio, $n=-(\beta+2)/(\beta-1)$; $v_{0}=\frac{V(\beta-1)^{2}}{6\kappa_{0}\beta}$; $V$ the solar wind speed; $\Omega$ the angular rotation speed of the Sun; subscripts “up” and “dn” refer to quantities upstream and downstream of the CIR reverse shock; and $J_{0}$ is a fitting constant. At the CIR shock, $J=J_{0}v^{n}e^{-v/v_{0}}$.

There are four free parameters in equation (14). These are $R_{s}$, $J_{0}$, $\beta$ and $v_{0}$. In fitting the observed CIR-associated energetic particle spectrum using the functional form of equation (14), one issue is that for rather different choices of the parameter set ($J_{0}$, $\beta$ and $v_{0}$), one can obtain very comparable fitting results. If, however, we assume that the solar wind speed and the shock compression ratio do not vary as a function of heliocentric distance, then $\beta$ is fixed by the STB in-situ observation. Furthermore, $v_{0}$ is also fixed by the STB in-situ observation. Therefore, the only two free parameters are $J_{0}$ and $R_{s}$.

Figure 8 shows the fitting results for the differential intensities for,

1. STB in-situ observation of the reverse shock during $17:54$ UT $\sim 20:09$ UT, Feburary $09$, 2008 2) ACE observation during $02:30$ UT $\sim 04:30$ UT, Feburary $11$, 2008 and 3) STB observation during $0:00$ UT $\sim 02:00$ UT, Feburary $11$. To obtain the best fitting, hese are slightly different time periods from our simulation.

Fitting the differential intensity at the reverse shock as observed by STB using a compression ratio $\beta=2.13$ from the in-situ measurement, we find $J_{0}$ (in unit of cm*-2* s*-2* sr*-1* (MeV/n)-1) is $1227$ and $v_{0}=0.172$ (MeV/n)1/2 (corresponding to $\kappa_{0}=0.03$). Using these two parameters we then find $J_{0}$(ACE)$=4167$ and $R_{s}=1.73$ AU; $J_{0}$(STB)$=5474$ and $R_{s}=3.63$ AU. The fitting results of $J_{0}$ and $R_{s}$ are summarized in Table 2. These shock locations are to be compared with our simulations where $R_{s}=1.39$ AU for the ACE observation and $R_{s}=2.55$ AU for the STB observation. The arbitrariness of the fitting can be seen from the following: for the ACE fitting, for example, we can obtain an almost equally well fitting with $J_{0}$(ACE)$=2070$ and $R_{s}=1.5$ AU.

In any events, we see that the Fisk & Lee model yields somewhat larger shock locations for both periods. Now considering the differential intensities: the Fisk & Lee model predicts a differential intensity at $R_{s}=1.73$ AU to be $\frac{4167}{1227}=3.4$ times that when $R_{s}=1$ AU (or a differential intensity at $R_{s}=1.5$ AU to be $\frac{2070}{1227}=1.69$ times that when $R_{s}=1$ AU); and a differential intensity at $R_{s}=3.63$ AU to be $\frac{5474}{5160}=4.7$ times that when $R_{s}=1$ AU. Therefore the Fisk and Lee model also suggests that the differential intensity at the shock increases with heliocentric distance of the shock.

5 Conclusion

We have examined energetic particle transport in one CIR event which occurred in $2008$ February. We choose this event because the CIR was observed by both STB and ACE and the CIR-associated energetic particle intensity profiles from these two spacecraft reasonably resemble each other (with a time shift of $\sim 1.4$ days). Furthermore, STB observed the reverse shock in-situ, suggesting that the shock was formed near $1$ AU. Under the assumption of no cross-field diffusion, we develop a Monte-Carlo test particle model to investigate the transport of energetic particles. The model solves the focused transport equation numerically. Comparing to previous analytical work by Fisk and Lee (1980), our model considers explicitly the particle pitch angle evolution. Both the magnetic focusing effect and the pitch angle diffusion process are included. For the cases where the solar wind MHD turbulence is not strong, our approach is more appropriate than that of (Fisk and Lee, 1980).

Assuming the reverse shock can be approximated by a stationary structure in the co-rotating frame, and that the shape of the accelerated particle spectrum at the shock does not vary with heliocentric distance, we calculated the different intensity at $1$ AU for two periods corresponding to an ACE observation and a STB observation. By assuming the accelerated particle spectrum at the shock is given by that observed in-situ at STB and does not vary with heliocentric distance, we avoid considering particle acceleration at the shock explicitly. Presumably, the acceleration at the shock does vary with the heliocentric distance. Under the diffusive shock acceleration framework, the accelerated particle spectrum depends on various shock parameters including the shock geometry and the shock compression ratio. Both can depend on the heliocentric distance. However, since we have only 1 AU observations, no constraints on these parameters can be obtained. We note this as a limitation of the present work.

By assuming a turbulence similar to that given by the WKB approximation, reasonable agreements between the simulation and the observations are obtained for both observations. The best fit of the STB observation yields $(\delta B/B_{0})^{2}=0.012$, similar to the in-situ value of $0.01$. The best fit of the ACE observation yields $(\delta B/B_{0})^{2}=0.02$, also similar to the in-situ value of $0.019$. The ACE fitting suggests that the seed particle density at $r_{sh}=1.39$ AU is $0.61$ times that when $r_{sh}=1$ AU, consistent with the radial dependence of suprathermal He*++* obtained in Hill et al. (2009). However, the STB fitting suggests that the seed particle density at $r_{sh}=2.55$ AU has to be $1.27$ times larger than that when $r_{sh}=1$ AU. This contradicts to Hill et al. (2009). This requirement of a large seed particle density at $r_{sh}=2.55$ AU is due to the fact that the modulation, especially for low energy particles, is strong.

This contradiction may be resolved by including cross-field diffusion. Including cross-field diffusion can effectively negate the modulation effect since with cross-field diffusion included, low energy particles that accelerated at a shock location closer to $1$ AU can diffuse cross-field and contribute to the observed intensity at a location that magnetically connects to the shock at a large distance (e. g. $\sim 2.55$ AU).

In summary, we have developed a numerical Monte-Carlo code to examine energetic particle transport at CIR shocks. We ignored particle cross-field diffusion and consider particles propagating along a single field line. We examined one CIR event which occurred in $2008$ February using our code. We simulated energetic particle spectra for two two-hour windows at both ACE and STB. Very good agreement between the simulation and observation can be obtained if we assume the seed particles at $2.55$ AU is $\sim 1.27$ times that at 1 AU and if the seed particles at $1.39$ AU is $\sim 0.61$ times that at 1 AU. However, we note that this conclusion may be changed if cross field diffusion is included.

Acknowledgements.

We thank the many individuals at the University of Maryland, Johns Hopkins Applied Physics Laboratory, and Goddard Space Flight Center responsible for the construction of the ACE/ULEIS and STEREO/IMPACT/SIT instruments. We acknowledge the use of data from the STEREO and ACE missions. The STEREO/SIT data was obtained from http://www.srl.caltech.edu/STEREO/Level1/SIT\_public.html, while the STEREO/PLASTIC and /MAG data was obtained from http://aten.igpp.ucla.edu/forms/stereo/level2\_plasma\_and\_magnetic\_field.html. The ACE/ULEIS, /SWEPAM, and /MAG data was obtained from the ACE Science Center (http://www.srl.caltech.edu/ACE/ASC/level2/). GL and LZ were supported in part by NSF grant AGS1135432, NASA grants NNX13AE07G, and NNX14AC08G at UAH; RWE, MID, and MAD at SWRI were supported in part by NASA grant NNX13AE07G, NSF grant AGS-1460118 and an internal research grant from Southwest Research Institute. The work at APL was supported by NSF grant AGS-0962653, contract SA4889-26309 and NASA grant NNX13AR20G from the University of California Berkeley. The work at ShangDong University was supported by NNSFC grants 41274175 and 41331068.

Appendix A Instantaneous Co-rotating Frame Transformation

The adiabatic deceleration effect in an expanding solar wind is described by (Parker, 1965; Jokipii and Parker, 1970)

[TABLE]

where $v_{sw}$ is the solar wind speed and $p$ is particle’s momentum in solar wind frame. This expression of average deceleration rate is valid under the assumption that particle’s pitch angle diffusion is nearly isotropic (Ruffolo, 1995). In Ruffolo (1995) the analytical expressions of the particle’s momentum (pitch angle) deceleration rate for an individual particle is given by

[TABLE]

where $\mu$ denotes the particle’s pitch angle in solar wind frame. $\psi$ is the angle between the radial direction and the magnetic field’s tangent direction. $1/L(z)=-B(z)/(\partial B/\partial z)$ is the reciprocal of the scale length of the interplanetary magnetic field. The first equation in Equation 17 describes the adiabatic cooling effect and the second equation describes pitch angle focusing effect.

In this paper, particle’s adiabatic deceleration and focusing effect are not modeled by Equation 17, instead these effects are treated implicitly by a frame transformation approach. Figure 9 shows the local coordinates of two adjacent instantaneous co-rotating frames. In panel ($a$), the black curve is the spiral magnetic field line and $A$, $B$ are two adjacent locations with a time interval $dt$. The radial distance of $A$ is $R_{1}$ and of $B$ is $R_{2}$. $\bm{r_{1}}$ ($\bm{r_{2}}$) is the unit vector in the direction of $R_{1}$ ($R_{2}$). $\bm{\theta_{1}}$ ($\bm{\theta_{2}}$), which is also the direction of instantaneous co-rotating frames at location $A$ ($B$), is unit vector in the direction perpendicular to $\bm{r_{1}}$ ($\bm{r_{2}}$). The four vectors’ relative directions are shown in panel ($b$). $\Delta\psi$ is the angle between $\bm{r_{1}}$ and $\bm{r_{2}}$.

Suppose the Sun rotates with an angular velocity $\Omega$, we get the velocities of the two instantaneous co-rotating frames

[TABLE]

and the relative velocity between two co-rotating frames is

[TABLE]

Take the dot product of ($\bm{V_{1}}-\bm{V_{2}}$) and $\bm{AB}$ and substitute Equation 18, we get

[TABLE]

Clearly from panel ($b$) in Figure 9, $\bm{r_{1}}$ is perpendicular to $\bm{\theta_{1}}$ and $\bm{r_{2}}$ is perpendicular to $\bm{\theta_{2}}$, then

[TABLE]

Combining Equation 20 and 21, the dot product of ($\bm{V_{1}}-\bm{V_{2}}$) and $\bm{AB}$ is zero, which means if the time step is small, the relative velocity between two co-rotating frames is perpendicular to the Parker’s spiral (particle’s trajectory). One can then obtain the particle’s energy and the momentum component parallel and perpendicular to the Parker’s spiral from the Lorentz transformation.

Appendix B Determination of Shock Cross Section

Total number of injected particles ($N_{0}$) at CIR reverse shock plays an important role in the simulation-observation comparison. Assuming the shock strength does vary along shock surface, the total number of injected particles is proportional to the seed particle density (n) and the shock cross section (A) in a flux tube. Figure 10 shows the configuration of the reverse shock. Similar to Figure 4, red and blue curves are Parker spirals, black thick curve extending from $A$ to $D$ is the reverse shock surface. The sun is located at $O$. Red (blue) Parker spiral intersects the radial direction at $A$ and $E$ ($B$ and $D$). The length $OA$ is the same as $OB$ and $OC$, which is $r$. The vectors $\bm{n}$, $\bm{r}$, $\bm{r^{\prime}}$ and $\bm{t}$ are unit vectors in the direction of shock normal, $OA$, $OB$ and tangent direction of the Parker spiral colored in red at point $A$ respectively.

The following statements and calculations are based on the approximation that $\theta$, which is the angle between $OA$ and $OB$, is small. Then, $AB$ is perpendicular to $OA$ and $OB$; $BC$ is perpendicular to $OB$ and $OD$; $AB$ and $BC$ are in the direction perpendicular to the radial direction. Angle between $\bm{r}$ and $\bm{t}$ is $\psi$ and we have $tan\psi=B_{t}/B_{r}$ ($B_{t}$ and $B_{r}$ are given in Equation 6). $\gamma$ is the angle between $\bm{r^{\prime}}$ and $BD$, and $\beta$ is the angle between $BD$ and $CD$. And in the approximation, we have $\gamma=\psi$ and $\beta=\psi$. Angle between shock surface and $AB$ is $\alpha$, and the distance of $CD$ is $dr$.

Follow the law of cosine, in the triangle $ABD$, we get

[TABLE]

The length of $AB$ is given by $r\theta$ and the length of $BD$ is given by $dr/\cos\beta$ (in the right triangle $BCD$), where $dr=AD~{}sin\alpha$ (in the right triangle $ACD$). We, therefore, get

[TABLE]

In the triangle of $ACE$ and $ACD$, the relation between $tan\alpha$ and $tan\psi$ can be expressed as

[TABLE]

where $V_{sh}$ and $V_{sw}$ is the shock speed and fast solar wind speed and $dt$ is time for the sun to rotate by $\theta$ degree. Equation 24 can be simplified into $tan\alpha~{}tan\psi=V_{sh}/V_{sw}$.