Modeling Ion Acceleration and Transport in Corotating Interaction Regions: The Mass-to-charge Ratio Dependence of the Particle Spectrum

We investigate the role of perpendicular diffusion in shaping the energetic ion spectrum in corotating interaction regions (CIRs), focusing on its mass-to-charge (A/Q) ratio dependence. We simulate a synthetic CIR using the EUropean Heliospheric FORecasting Information Asset and model the subsequent ion acceleration and transport by solving the focused transport equation incorporating both parallel and perpendicular diffusion. Our results reveal distinct differences in ion spectra between scenarios with and without perpendicular diffusion. In the absence of perpendicular diffusion, ion spectra near CIRs show a strong (A/Q) ϵ dependence with ϵ depending on the turbulence spectral index, agreeing with theoretical predictions. In contrast, the incorporation of perpendicular diffusion, characterized by a weak A/Q dependence, leads to similar spectra for different ion species. This qualitatively agrees with observations of energetic particles in CIRs.


Introduction
Corotating interaction regions (CIRs) form when the faster solar wind overtakes slower solar wind.The stream interaction leads to compression waves and eventually the formation of the forward-reverse shock pair.These persistent, large-scale structures are often associated with tens of keV n -1 to several MeV n -1 energetic particles in the inner heliosphere (Richardson 2018).Early work by Fisk & Lee (1980) assumed that CIRs accelerate particles beyond several astronomical units, where shocks are generally formed.Later it was argued that compression regions associated with CIRs can accelerate ions in a similar way to shocks, but without the need of shocks (Giacalone et al. 2002).While Ulysses observations have seen energetic ions at the reverse shock as far as 3 au (Mason et al. 1999), recent observations at 1 au have pointed to a local source for CIR-related energetic ions, specifically those below approximately 1 MeV n -1 (Mason et al. 2008;Ebert et al. 2012;Wijsen et al. 2021).These findings suggest that the energy spectrum observed by L1 spacecraft results from a combination of local acceleration processes and transport effects as particles move inward from distant shocks.Previous work by Mason et al. (2008) has examined the abundance ratios of heavy ions versus their energy/nucleon, revealing a near constancy below 2 MeV n -1 .This constancy is also mirrored in spectral rollover energies (E 0 ).A more recent study by Filwett et al. (2019) found that the rollover energy in stream interaction regions does not show a significant mass-to-charge (A/Q) ratio dependence, a stark contrast to solar energetic particle (SEP) events, where a varying A/Q dependence of ion time profiles and/or spectra are usually seen (Mason et al. 2006;Desai et al. 2016a).
Given that the acceleration and the subsequent transport of energetic ions from CIRs is diffusive in nature, the presence or absence of the A/Q dependence in ion rollover energies can be related to the A/Q dependence of the ion diffusion coefficients.In considering the transport of ions in the solar wind, the ion mean free path (λ) can be expressed as a power law in rigidity (R = A/Qv), λ ∼ R α , implying a diffusion coefficient ) that scales as (A/Q) α v ( α+1) (Ellison & Ramaty 1985;Cohen et al. 2003).Different ions of different energies will have similar time profiles if their diffusion coefficients are the same.This is often referred to as the equal diffusion coefficient condition (e.g., Mason et al. 2006Mason et al. , 2012)).In considering the acceleration process, if κ has a certain A/Q dependence, it will translate to a rollover energy that depends on A/Q (Cohen et al. 2005).Li et al. (2009) suggested that the shock geometry plays a critical role in determining the A/Q dependence in gradual SEP events, with the total diffusion coefficient as a function of shock obliquity.In the work of Li et al. (2009), the rollover energy of shock-accelerated ions is shown to have a (A/Q) ò dependence with ò= − 2 for parallel shocks and ò = −1/5 for perpendicular shocks.Later observations of large SEP events associated with interplanetary shocks (Desai et al. 2016a(Desai et al. , 2016b) ) found the spreading of ò in good agreement with Li et al. (2009).If the observer is not at the acceleration site, transport may complicate the A/Q dependence of the rollover energy.Zhao et al. (2016b) studied how particle transport with pitch-angle scattering might contribute to the A/Q dependence observed in heavy ion spectra during SEP events.These studies have underscored the significance of diffusion coefficient and the shock geometry in understanding the A/Q dependence of the ion acceleration and transport.
Unlike SEP events studied in Desai et al. (2016aDesai et al. ( , 2016b)), where the shock geometry varies as a shock propagates out and where the transport effects also vary with time as the distance between the shock and the observer decreases, CIRs are mostly stationary structures in large scales.Therefore, the roles of ion acceleration and transport in shaping the CIR energetic particle spectrum are less entangled and are easier to understand than SEP events.However, at small scales (several to tens of ion inertial lengths), the shock front can be irregular and nonstationary.Recent observations and kinetic simulations by Kajdič et al. (2019), Guo et al. (2021), and Trotta et al. (2023a, 2023b) have shown that small-scale irregularities are ever present at the shock front, evolving in space and time.These structures can be crucial in understanding the injection of superthermal particles because the injection depends on the shock geometry (Li et al. 2012), and as shown in Trotta et al. (2023a), the local shock obliquity shows high variability along an irregular shock front, leading to an inhomogeneous injection of the superthermal particles.These results suggest that the shock properties in the large-scale magnetohydrodynamic (MHD) simulations, where small-scale irregularities are ignored, can only be regarded as some mean properties of the shock.
Assuming the interplanetary magnetic field takes the form of nominal Parker field, and ignoring perpendicular diffusion and small-scale shock irregularities, then any given observer will connect to a single point at the CIR shock.Under such an assumption, Zhao et al. (2016a) examined how energetic CIR particle spectra depend on the observer's longitude.This picture, however, needs to be revised if we consider the effect of perpendicular diffusion.Particle perpendicular diffusion has two contributions.The first is the meandering of the field line due to the presence of the magnetic fluctuation δB, which contains no A/Q dependence.This has been examined by, e.g., Bian & Li (2021, 2022) and Li & Bian (2023).Because the field line itself is "diffusive," when particles move along these field lines, there will be displacements perpendicular to the unperturbed background field.Therefore, perpendicular diffusion is also referred to cross-field diffusion (with respect to the unperturbed field).Under the zero-gyroradius approximation, particle perpendicular diffusion is the same as that of the meandering field line (Bian & Li 2021, 2022).The second contribution is the finite gyroradius effect.Retaining a nonzero gyroradius, one can evaluate particle's trajectory and compute its perpendicular displacement from the original unperturbed field lines.Matthaeus et al. (2003) developed a nonlinear guiding center (NLGC) theory to compute particle's perpendicular diffusion coefficient where the underlying turbulence is assumed to be of "slab+2D" geometry.This work was later extended to other turbulence geometries (Shalchi 2009;Shalchi et al. 2010).
Previous simulations of CIR events (Zhao et al. 2016a;Wijsen et al. 2021) have not explicitly and systematically considered the effect of perpendicular diffusion.This motivates our current study.In this Letter, we investigate how perpendicular diffusion can influence ion acceleration and transport in a CIR event, and potentially lead to a weak A/Q dependence of ion spectra.In Section 2, we provide a brief discussion of our modeling approach.We first describe the procedure of simulating a synthetic CIR using the EUropean Heliospheric FORecasting Information Asset (EUHFORIA; Pomoell & Poedts 2018) model.We then describe how the particle acceleration and transport is followed from the vicinity of the CIR reverse shock to an observer.This follows closely to recent works by Wijsen et al. (2021Wijsen et al. ( , 2023)).In Section 3, we present a comparative analysis for cases with and without perpendicular diffusion.This comparison focuses on the ion fluxes and spectra at 1 au.The conclusion is given in Section 4.

Model
A CIR is modeled by using EUHFORIA, which solves the ideal MHD equations from 0.1 to 5 au.Following Wijsen et al. (2019), a steady CIR structure is achieved by setting a uniform solar wind speed of 400 km s −1 at the inner boundary, with the exception of a specific region where the solar wind speed is elevated to 750 km s −1 .This region is defined by the condition , where f and θ are the longitude and latitude, respectively.For simplicity, we assume a monopolar positive magnetic field at the inner boundary to avoid heliospheric current sheet.With this configuration, the interaction between the slow and fast streams leads to the formation of a CIR, characterized by forward and reverse shock waves.
To examine the acceleration and transport of energetic particles, we follow a similar approach to Wijsen et al. (2021) and inject seed ions of 50 keV n −1 with a radial dependence of r −2 to the reverse shock.This, of course, is a simplification, as we remarked earlier that the injection energy and efficiency of seed ions depend on the local shock geometry, and small-scale shock irregularities (Trotta et al. 2023a) can affect the shock geometry.We then follow these particles by solving the threedimensional focused transport equation through the EUHFORIA solar wind domain in the corotating frame (Wijsen et al. 2019), expressed as In this equation, the gyroaveraged particle distribution function is represented by f (x, p, μ, t), which is a function of spatial coordinate x, momentum magnitude p, the pitch-angle cosine μ and time t.Additional elements in Equation (1) include the particle velocity v, the velocity of the solar wind V sw , and the unit vector of the magnetic field b.Both V sw and b are taken from a single steady-state snapshot from the EUHFORIA model, so V sw and b are time independent in the corotating frame.The two terms on the right-hand side of Equation (1) describe particle diffusion.These are the pitch-angle diffusion given by D μμ (related to the parallel diffusion coefficient κ ∥ ) and the perpendicular diffusion, represented by the diffusion tensor κ ⊥ .These terms collectively model the influence of a fluctuating magnetic field, δB, on particle transport.Additionally, we adopt the assumption by Wijsen et al. (2023) for the radial dependence of δB 2 /B 2 .Following Ding et al. (2022) in describing D μμ using the the quasi-linear theory (QLT; Jokipii 1966) and κ ⊥ from κ ∥ using the NLGC theory (Shalchi et al. 2010), we can obtain the Here ò ∥ is decided by the spectral index of the turbulence power spectral index β (I(k) ∼ k β ).Using QLT and NLGC theory, Li et al. (2009) showed that ò ∥ = β + 2 and ò ⊥ = (β + 2)/3.Clearly, κ ⊥ has a weaker A/Q dependence than κ ∥ .A recent statistical study by Park et al. (2023) found that β ranges from −1.4 to −2.0, both downstream and upstream of the fast forward and fast reverse interplanetary shock, with a mean value around −1.7.In this work, we assume β to be −5/3 when considering the transport of energetic ions in the solar wind, away from the reverse shock.This corresponds to a Kolmogorov spectrum.Close to the shock, besides the case of β = − 5/3, we also consider a control case of β = − 2. In this case, one finds both κ ∥ and κ ⊥ are independent of ion species.The turbulence level near a CIR and its reverse/forward shock can be enhanced (Crooker et al. 1999;Richardson 2004).Indeed, wave intensity across interplanetary shocks are often enhanced.A case study of the large 2003 October 29 event (Li et al. 2005a) showed an enhancement of more than a factor of 10.Recent work by Park et al. (2023) suggested that the turbulence level across fast shocks can increase, on average, by 3-5 times.From the rarefaction region, the spectral index β does not differ much from the Kolmogorov value.This is because the CIR shocks are quasi-perpendicular in geometry; therefore there is not much amplification of waves due to streaming ions, as in the case of parallel shocks (Li et al. 2003(Li et al. , 2005b)).
The simulation duration is taken to be 27 days, approximately one solar rotation.At the inner and outer boundary, a reflective and an absorbing boundary condition are used, respectively.The steady-state solution at a desired position is obtained for a continuous time-independent injection by convoluting the Green's function solution at different times.

Results
All analyses are done in the equatorial plane, as in the recent work by Wijsen et al. (2023).Figure 1(a) displays the radial velocity of the solar wind, V r .Slow and fast winds can be seen clearly in the plot.Figure 1(b) plots ∇ • V sw .Since the solar wind propagates radially out, ∇ • V sw is positive if V sw is a constant, but can become negative near stream interface.In Figure 1(b), these negative divergence regions include both the forward and reverse shocks.The gray shaded regions indicate expanding solar wind areas.In this study, we limit the injection only at the reverse shock as this allows for a clean analysis of the particle spectra.We note that the reverse shock is more efficient in accelerating particles (Lee & Fisk 1982).Figure 1(c) shows the parallel mean free path λ ∥ for 50 keV protons.Away from the reverse shock, λ ∥ increases with the heliocentric distance, reflecting the gyrofrequency of particles and their wave-particle resonance conditions with the background magnetic field (Wijsen et al. 2023).Note that in the rarefaction region particles exhibit a larger parallel mean free path.Figure 1(d) shows the resulting perpendicular mean free path λ ⊥ for 50 keV protons, computed from λ ∥ (see details in Ding et al. 2022).The similar azimuthal variation of λ ⊥ is due to its dependence on  l 1 3 .Near the reverse shock, the mean free path decreases.This is because near the compression regions the wave activity increases due to stream interactions.We note that EUHFORIA simulations do not capture how the upstream wave intensity evolves.To model the transport of energetic particles, a prescription of the wave intensity, which amounts to a description of particle diffusion coefficient, is needed (Wijsen et al. 2021(Wijsen et al. , 2023)).Recently, in examining two energetic storm particle events observed by Solar Orbiter, Ding et al. (2024) assumed that the wave intensity upstream of the shock decays exponentially from the shock.In this work, we follow a similar approach to describe the wave intensity upstream of the CIR shocks.By way of example, we approximate the wave intensity I(k, x) upstream of the reverse shock by where k is the wavenumber; x is the position upstream of the shock and x 0 is the closest point to x on the shock; I 0 is the ambient wave intensity and a is the amplitude of enhanced wave intensity at the shock, taken to be 5 in our simulation; L, the diffusion length scale, is assumed to be 0.05 au.Note that L could be energy dependent and radial dependent, similar to the case considered at coronal mass ejection-driven shocks (Li et al. 2005b(Li et al. , 2021;;Ding et al. 2022Ding et al. , 2024)).Upstream of the shock, I(k) reduces exponentially to the ambient level in a distance of ∼L.Downstream of the shock, where ∇ • V sw < 0, we assume ), which is a constant.Since we follow the evolution of the particle distribution function with the focused transport equation, acceleration is implicitly taken care of when particles cross the shock.In fact, as shown in Wijsen et al. (2019Wijsen et al. ( , 2021Wijsen et al. ( , 2023)), shock is not necessary for accelerating particles, as only a negative divergence is needed (Giacalone et al. 2002).Our simulation follows closely to Wijsen et al. (2019Wijsen et al. ( , 2021Wijsen et al. ( , 2023)), but with the consideration of the additional turbulence near the compression region, as shown in Equation (2).We remark that the ion acceleration efficiency and their maximum energies positively correlate with the turbulence level near the CIR.It could be important to understand some large CIR events with maximum energy up to ∼20 MeV n -1 (Richardson 2004;Mason et al. 2008).
Figure 2 shows the divergence of the solar wind velocity, ∇ • V sw , and the shock obliquity, θ BN , as functions of radial distance along the reverse shock.Notably, the magnitude of ∇ • V sw initially increases with radial distance, reaching a maximum around 3.5 au, before gradually decreasing.This implies a radial variation in the strength of the reverse shock, suggesting a more efficient particle acceleration occurring around 3.5 au.This phenomenon could potentially correlate with the observed radial dependence of ion intensity during CIR events, where a positive correlation between ion intensities and radial distances has been recently noted (Allen et al. 2021).A more comprehensive exploration of this relationship will be pursed in a future study.The value of θ BN shows a gradual increase from 78°to 88°as one moves from 1 to 5 au, indicative of the shock geometry becoming more perpendicular further out.The obliquity of the CIR shock is crucial for understanding ion acceleration.As discussed by Li et al. (2009), the total diffusion coefficient κ can be generally expressed as where κ ∥,0 and κ ⊥,0 represent the parallel and perpendicular diffusion coefficients for protons at speed v 0 , γ ∥ and γ ⊥ are the indices denoting the dependence on proton speed and ò ∥ and ò ⊥ refer to the A/Q dependence.In the QLT and NLGC theory, one can obtain γ ∥ = (β + 3) and γ ⊥ = (β + 5)/3 (Ding et al. 2022).
To understand the impact of perpendicular diffusion on the A/Q dependence during CIR events, we compare simulation results at 1 au with and without the inclusion of perpendicular diffusion.Figure 3 presents the time profiles of solar wind speed and number density, and spectrograms of normalized ion flux intensity as functions of energy and time, including proton (H), helium (He), oxygen (O), and iron (Fe).The values of A/Q for H, He, O, and Fe are chosen to be 1/1, 4/2, 16/6, and 56/ 12, respectively (Desai et al. 2016a).To better compare the spectrogram and energy spectra, the injection of different ions is assumed to be the same.The observer is fixed in the inertial frame.The left panels of Figure 3 display simulation results for the case without perpendicular diffusion.In this case, a clear peak in ion intensity is observed near the compression region, followed by a gradual decay.The decay phase reflects the fact that the observer connects to further and further parts of the shock (the CIR structure is rotating counterclockwise in the inertial frame) where the injection drops.Note that ions accelerated beyond 1 au need to propagate back to 1 au.Because no perpendicular diffusion is included, a fixed observer in the corotating frame only connects to a single point at the shock.A fixed observer in the inertial frame, of course, connects to different parts of the shock at different times, yet the connection point moves along the shock with a constant rotation.This is vastly different from the case when we include the perpendicular diffusion, which is shown in the right panels of Figure 3.We see now the peaks near the compression regions are smeared out.Energetic ions appear ∼40 hr prior to the compression, and over a total of 120 hr (30 hr before and 90 hr after the compression), the intensities maintain relatively plateau-like shape.This suggests that the inclusion of the perpendicular diffusion allows an observer to "probe" a large portion of the shock surface at one time instead of only "probe" a single point on the shock surface at a time.
Figure 4 compares the time-integrated ion spectra near the compression region: from 45 to 60 hr.This period is chosen such that the observed particles are predominantly accelerated near 1 au, reducing the transport effect in the solar wind.In the left panels, where perpendicular diffusion is not included, the ion spectra exhibit good power-law-like behavior with an exponential decay.The spectra toward the higher-energy end show clear A/Q dependence, which is also reflected in the fitted rollover energy.We fit each ion spectra using a form of a single power law with an exponential tail, - ) , where γ is the spectral index and E 0 is the rollover energy.Dashed lines represent the fitting results.In this case, the rollover energies for different ions are MeV n −1 , which shows E 0 ∼ (A/Q) −0.51±0.04 .Since no perpendicular diffusion is included, this A/Q dependence is governed solely by κ ∥ and it is directly linked to the turbulence spectral index β, taken to be −5/3 here.The theoretical prediction of rollover energy with only κ ∥ is E 0 ∼ (A/Q) −2(β+2)/(β+3) indicated by Equation (3).Our fitting result is consistent with the prediction, E 0 ∼ (A/Q) −1/2 when β = −5/3.In the case of β = −2, the theory predicts no A/ Q dependence.This is shown in the inset of left panel.We do not show the corresponding spectrogram of ion intensity in Figure 3 for β = −2 since these are less informative.For β > − 5/3, the A/Q dependence becomes more significant (not shown here).The right panels in Figure 4 show the spectra with perpendicular diffusion included.In the main panel, where β = − 5/3, using the NLGC theory, we have κ ⊥ ∼ (A/Q) 1/9 .Following Equation (3), this leads to a E 0 ∼ (A/Q) −1/5 if only perpendicular diffusion dominates the acceleration process.Different ion spectra in the main panel show only minor differences above 1 MeV n −1 comparing to the left main panel.The fitted rollover energies shows a E 0 ∼ (A/Q) −0.27±0.07 .This is slightly stronger than (A/Q) −1/5 , since κ || is also playing a role.We also examine the case where κ ⊥ is A/Q independent (i.e., ò ⊥ = 0).Such a choice corresponds to considering only the field line meandering contribution to κ ⊥ .The resulting spectra are shown in the inset.As expected, it leads to an even less pronounced A/Q dependence, characterized by (A/Q) −0.18±0.02 .Again, this A/Q dependence is attributed to κ ∥ .As β → −2, then A/Q dependence diminishes, as indicated in the left subpanel.
Figure 4 is the most important result of this work.It demonstrates clearly the crucial role of perpendicular diffusion in regulating the accelerated ion spectrum near CIR shocks.It is interesting to note that the maximum energy of ions when perpendicular diffusion is included is lower than that when perpendicular diffusion is ignored.This is because perpendicular diffusion can effectively move particles away from the shock region and reduce the acceleration efficiency.This is similar to the behavior of the seed population when perpendicular diffusion is included (Wijsen et al. 2023).

Conclusion
In this study, we have explored the influence of perpendicular diffusion on ion acceleration and transport in CIRs, particularly focusing on the implications for the A/Q dependence in ion rollover energies.In the absence of perpendicular diffusion, ion spectra near the compression region exhibit a pronounced A/Q dependence in rollover energies, which is related to the turbulence power spectral index.However, the inclusion of perpendicular diffusion largely removes such a A/Q dependence, leading to a markedly different behavior.The spectra of different ion species do not differ much and the rollover energies of ions are similar.Our findings highlight the crucial role of perpendicular diffusion in the ion acceleration near CIR shocks.It implies that CIR environment could be the best place to examine perpendicular diffusion.Observations of CIR energetic particle spectrum (Mason et al. 2008;Filwett et al. 2019) showed no strong A/Q dependence.We suggest that this is because of a weak A/Q dependence of κ ⊥ , and the fact that CIR shocks are largely quasi-perpendicular.Our results provide a model basis for interpreting the A/Q dependence of ion spectra in CIR observations.

Figure 1 .
Figure 1.The solar wind radial speed and the assumed parallel and perpendicular mean free paths for 50 keV protons in the solar equatorial plane.Panel (a) shows the radial solar wind speed.White dashed lines are magnetic field lines from the simulation.Panel (b) shows the divergence of the solar wind velocity.The compression waves can be recognized as the colorful spiral-shaped structures.Panels (c) and (d) show, respectively, the parallel and perpendicular mean free paths for 50 keV protons.The simulation domain is from r = 0.1 au to r = 5 au.

Figure 2 .
Figure 2. The divergence of the solar wind velocity ∇ • V sw and the shock obliquity θ BN as a function of radial distance along the CIR reverse shock.

Figure 3 .
Figure 3.The left and right panels show from top to bottom: solar wind speed, solar wind number density, the spectrogram of normalized energetic ion flux intensities vs. energy and time (proton, helium, oxygen, iron) for an observer at 1 au.Left (right) panels are for the case without (with) perpendicular diffusion.

Figure 4 .
Figure 4.The normalized time-integrated ion spectra for an observer at 1 au.The left and right panels represent the case without and with perpendicular diffusion , respectively.In the left main panel β = − 5/3, and in the inset β = −2; in the right main panel ò ⊥ = 1/9, and in the inset ò ⊥ = 0.The dash lines show the fitted spectra.See the text for details.