© R. Vainio et al., Published by EDP Sciences 2014

1. Introduction

Solar energetic particle (SEP) events provide by far most of the >1-MeV proton fluence on interplanetary missions occurring during solar maximum years (e.g., Feynman et al. 1993). In SEP events the highest intensities at MeV energies are often associated with passages of interplanetary shocks driven by fast coronal mass ejections (CMEs). These high-intensity events are, for historical reasons, termed energetic storm particle (ESP) events (Reames 1999). In ESP events, the shock is observed as a continuous source of particles accelerated by multiple crossings of the shock front, facilitated by scattering off magnetic irregularities in the vicinity of the shock. The resulting diffusive shock acceleration (DSA; e.g., Drury 1983) in CME-driven shocks in the solar corona and interplanetary medium is governed by the evolution of the properties of the shock and the ambient turbulence (Bell 1978; Lee 1983, 2005). While magnetohydrodynamic (MHD) simulations are capable of describing the evolution of global shock properties and macro-scale quantities in the solar corona and solar wind during solar eruptions (e.g., Pomoell et al. 2011; Pomoell & Vainio 2012), kinetic scales are not addressed by them. Turbulent fluctuations responsible for the scattering of energetic ions are in gyroresonance with the particles so their scales are of the order of the ion Larmor-radius, beyond the MHD description. Thus, it is necessary to find alternative ways to describe the fluctuations undergoing wave-particle interaction processes in the region ahead the shock wave, called the foreshock.

When modelling the peak fluxes and fluences of the SEP and ESP events, the question of spatial scaling from 1-AU observations to other heliospheric distances is frequently encountered (e.g., Lario & Decker 2011). Especially the new missions to the inner heliosphere (ESA’s Solar Orbiter and NASA’s Solar Probe Plus) have introduced a need for obtaining reliable models of the SEP/ESP radiation environment in the inner heliosphere. Several approaches to tackle the problem have been undertaken. Observations from spacecraft inside 1 AU have been used to obtain scaling laws for fluences and peak fluxes (Lario et al. 2006) as well as peak fluxes observed early in the SEP event before the ESP peak (Reames & Ng 1998). Another way to estimate the spatial scaling is to utilise particle transport modelling (e.g., Aran et al. 2005; Lario et al. 2007; Vainio et al. 2007; Verkhoglyadova et al. 2012). The problem of the modelling approach is that the peak fluxes and fluences at MeV energies are determined to a large extent by the strength of the ESP event. The main complication is that the particle intensity close to the shock is decoupled from the far upstream intensity because of the self-generated turbulence, acting as a trapping region close to the shock. In fact, the strength of the source, as long as it exceeds the threshold of efficient wave generation in the foreshock, has a very weak effect, if any, on the flux that escapes the shock to the far upstream region where the bulk of the observations are made (Ng & Reames 1994; Reames & Ng 1998; Vainio 2003). This effect is known as the streaming limit.

While the remotely observed intensities cannot reveal the true strength of the particle source due to the streaming limit, the maximum proton energy emitted by the shock has a direct connection to the intensity of turbulence upstream of the shock. This is because the rate of diffusive shock acceleration is proportional to the power of the (resonant) magnetic fluctuations in the foreshock (Zank et al. 2000; Rice et al. 2003; Vainio & Laitinen 2007, 2008). The power, in turn, is directly proportional to the trapped particle flux. Thus, by observing the maximum energy of the escaping particles as a function of time allows modelling the intensity of the fluctuations and accelerated particles at the foreshock even when the shock is still in the corona or inner solar wind at a large distance from the observer (Horne et al. 2013), provided that an accurate model of the coupled processes of particle acceleration and wave generation is used.

While analytical solutions to the problem of coupled wave-excitation and particle acceleration have been known since the early development of DSA theory (Bell 1978; Lee 1983), they are typically based on a number of simplifying assumptions: the shock is taken to be a stationary planar discontinuity and the system is assumed to be in a steady state. Neither of these assumptions is true in evolving CME-driven shocks. A more complicated analytical model, taking into account the radial field geometry and evolution of the shock in an adiabatic manner, was constructed by Lee (2005). While this model is able to predict some of the key observational features of large SEP events, the number of analytical approximations in the modelling has to be relatively large to keep the model tractable. Thus, more accurate (numerical) models to assess the energetic particle transport and acceleration through self-amplified turbulence in the foreshock region are necessary. A few numerical simulation codes have been constructed to tackle the problem of coupled particle acceleration/transport and wave generation using a fully kinetic approach (Vainio & Laitinen 2007, 2008; Ng & Reames 2008). These codes describe the physics of particle acceleration at the highest level of accuracy but as a result they are typically rather heavy to run. For that reason they are not yet suitable for operational space-weather use.

A predecessor of the fully kinetic numerical codes was developed at the turn of the century by Zank and co-workers (2000; Li et al. 2003; Rice et al. 2003) and is presently known as the PATH code (Verkhoglyadova et al. 2009, 2010, 2012). It is based on an analytical approximation of the shock-accelerated particle spectrum that can be derived from the MHD properties of the shock and upstream fluid, and only the transport of particles around the shock wave (upstream and downstream) is treated numerically. This approach has its limitations, but a clear advantage is that the particle transport simulation is kept linear (non-linear coupling affects only the analytical part of the solution) and has, therefore, the potential of being applied in real time with relatively reasonable computer resources. In order to obtain realistic estimates for the accelerated particle spectra in PATH-like models, the analytical approximations employed in the methodology would have to be validated using self-consistent numerical computations. Another modelling approach that has great potential in space-weather applications is the semi-empirical source function approach, where the particle source at the shock is parameterised as a function of shock properties and the parameters of the model are then fitted to obtain a good fit to real individual SEP events (Lario et al. 1998; Aran et al. 2005). However, the semi-empirical source function is able to address the ESP population only, when the shock is close to the observer. “Remote sensing” is limited to the escaping ions and, thus, predictions of the ESP events at distances not accessible to observations have to be based on heuristic models of the spatial evolution of the foreshock turbulence.

In this paper, we report on numerical simulations of the foreshock of a coronal/interplanetary shock wave that propagates outward from the outer corona (~5.8 solar radii) to the solar wind (~60 solar radii). We analyse the evolution of the energetic particle intensities and wave power spectra in the foreshock region as a function of distance from the Sun and the encountered seed-particle spectral properties and derive a semi-analytical model that describes the particle scattering properties in the foreshock. The model can be utilised in studies of spatial scaling of the ESP peak intensities as well as in theoretical studies of coronal shock acceleration without the need to make time-expensive computer simulations of the coupled evolution of the turbulence and the accelerated particles.

The structure of the paper is as follows: we briefly describe the simulation code used in the study in Section 2, describe the semi-analytical model to be validated by the simulations in Section 3, present the results of the analysis in Section 4, discuss them in Section 5 and present the conclusions and outlook of the study in Section 6.

2. Simulation model

Coronal Shock Acceleration (CSA; Vainio & Laitinen 2007, 2008; Battarbee et al. 2011; Battarbee 2013) is a numerical simulation code describing the coupled evolution of energetic particles and Alfvénic turbulence in the region upstream of a propagating coronal shock front. The code traces individual ions under the guiding centre approximation in a radial magnetic flux tube that is being traversed by a coronal shock wave with variable speed V _s and normal angle Θ _Bn . The particles are followed under the influence of the large scale magnetic field and superposed turbulent Alfvénic fluctuations. Waves are described by a WKB-based transport equation with additional terms introduced for taking into account wave growth due to accelerated particles and diffusion in frequency that mimics the effect of wave-wave interactions and improves numerical stability.

The wave-particle interactions are described at the quasi-linear level, i.e., particle scattering rates are taken to be proportional to the wave power spectrum at a resonant wavenumber (k _r = Ω/γv = mΩ/p) and wave-growth rates proportional to the particle streaming at the corresponding resonant momentum (p _r = mΩ/k). Here, v, γ, p, Ω =qB/m, m and q are the speed, Lorentz factor, momentum, cyclotron frequency, mass and charge of the ion, respectively, B is the magnetic field and k is the wavenumber of the Alfvén wave. Note that the full quasi-linear resonance condition, k = mΩ/pμ, contains a dependence on the pitch-angle cosine, μ, which is neglected in the CSA code for simplicity. This simplification speeds up the code by at least an order of magnitude.

The output of the code is the SEP intensity, j(x, E, t), as a function of distance from the shock x, energy E, and time t; and the power spectrum of the Alfvén waves, P(x, f, t), as a function of distance, frequency f = Vk/2π, and time. Here, V = u _sw + v _A is the phase speed of the Alfvén waves in the solar frame, u _sw the solar-wind speed and v _A the Alfvén speed. The input file of the code allows the user to specify the initial distance of the shock wave, R _S(0), the temporal evolution of the shock-normal speed (V _S) and obliquity (cosΘ _Bn ), the total simulation time t _max and temporal resolution Δt of the output files, spatial and frequency/speed resolution of the output grids, and the seed-particle distribution as a kappa distribution (Battarbee et al. 2011),$$F_0(R,p)=n\left(R\right)\frac{\mathrm{\Gamma }(\mathrm{\kappa }+1)}{m^3{w}_0^3{\pi }^{3/2}{\mathrm{\kappa }}^{3/2}\mathrm{\Gamma }(\mathrm{\kappa }-\frac{1}{2})}{\left(1+\frac{p^2}{\mathrm{\kappa }{m}^2{w}_0^2}\right)}^{-\mathrm{\kappa }-1},$$(1)with its density scaled to the density of the coronal plasma, and the properties of the solar wind and interplanetary field (density, solar-wind speed and radial magnetic field) at 1 AU. We consider only protons in our simulation. Note that the thermal velocity w ₀ of the kappa distribution is obtained from the kinetic temperature, T, as$$\frac{1}{2}m{w}_0^2=k_{\mathrm{B}}T\frac{\mathrm{\kappa }-\frac{3}{2}}{\mathrm{\kappa }}.$$(2)

The user can specify the ambient distribution of protons as a superposition of two kappa distributions, which adds flexibility to the model. The density profile n(R) (normalised by the user-specified 1-AU density) and the profile of the kinetic temperature T(R) are fixed in the simulation to a semi-empirical model (Cranmer & van Ballegooijen 2005). The solar-wind speed profile u _sw(R) is obtained from the user-specified 1-AU speed using the conservation of mass, i.e.,$$R^2{u}_{\mathrm{s}\mathrm{w}}n=\mathrm{c}\mathrm{onstant}.$$(3)

The modelled magnetic field at the considered distances is radial and set to B(R) = B ₀(R ₀/R)² where B ₀ is the field at R ₀ = 1 AU.

3. Semi-analytical foreshock model

The semi-analytical foreshock model is based on the 1-D steady-state theory of Bell (1978), which can be given in form (Vainio & Laitinen 2007, 2008; Horne et al. 2013)$$F\left(x,p\right)=F_{\mathrm{S}}\frac{x_0}{x+x_0},$$(4) $$\lambda \left(x,v\right)=\frac{3\left(u_1-v_{\mathrm{A}}\right)}{v}\left(x+x_0\right),$$(5) $$F_{\mathrm{S}}\left(p\right)=\frac{\sigma ϵn_{\mathrm{p}}}{4\pi p_{\mathrm{i}\mathrm{nj}}^3}{\left(\frac{p}{p_{\mathrm{i}\mathrm{nj}}}\right),}^{-\sigma }$$(6) $$x_0\left(p\right)=\frac{2}{\pi ϵ\sigma }\frac{v_{\mathrm{A}}}{{\mathrm{\Omega }}_{\mathrm{p}}}{\left(\frac{p}{p_{\mathrm{i}\mathrm{nj}}}\right)}^{\mathrm{\sigma }-3},$$(7)where F is the particle distribution function, λ is the scattering mean free path, x is the distance from the shock wave towards the upstream region along the mean magnetic field, v and p are the particle speed and momentum, v _A is the Alfvén speed, Ω_p is the proton cyclotron frequency, n _p is the upstream proton density, and u ₁ is the shock speed relative to the upstream medium along the ambient magnetic field. The parameter is the fraction of upstream protons injected into the acceleration process at the injection momentum, p _inj, taken to be p _inj = 2m(u ₁ – u ₂), which results from cold upstream ions making a head-on collision with the downstream scattering centres. The quantity x ₀ is the scale height of the turbulent fluctuation intensity at the shock. Finally, σ = 3r _c/(r _c − 1) is the spectral index of the particle distribution at the shock, and r _c = (u _1n  − v _An )/u _2n  = r _g (1 – M ⁻¹) is the scattering-centre compression ratio assuming that the scattering in the downstream region is governed by frozen-in magnetic fluctuations or strong turbulence with zero cross-helicity. Here, M = u _1n /v _An  = u ₁/v _A is the Alfvénic Mach number of the upstream flow, the subscripts n in the velocities refer to the shock-normal components, and r _g = ρ ₂/ρ ₁ is the gas compression ratio of the shock. Note that this formulation neglects transport perpendicular to the mean magnetic field, so it cannot be applied to quasi-perpendicular shocks.

Since the interplanetary magnetic field is not a homogeneous field but scales with the heliocentric distance as R ⁻² in the inner heliosphere, we have to correct the spatial distribution of particles to take that scaling into account. The radial expansion of the field and the corresponding increase of the flux tube cross-sectional area as is taken into account by considering particle intensity as$$j(x,E)=j_{\mathrm{S}}\left(E\right)\frac{x_0}{x+x_0}{\left(\frac{R_{\mathrm{S}}}{R_{\mathrm{S}}+x}\right)}^2,$$(8)where is the heliocentric distance of the shock and is the intensity as a function of energy at the shock. Instead of the distance of the shock, the intensity can be normalised at another point in the flux tube, e.g., at . Equation (8) is valid as long as the streaming caused by adiabatic focusing is small. This condition can be stated as (Battarbee et al. 2011)$$\mathrm{\lambda }/L\ll 3(u_1-v_{\mathrm{A}})/v,$$(9)where is the focusing length. When, instead,$$\mathrm{\lambda }/L\gtrsim 3(u_1-v_{\mathrm{A}})/v,$$(10)particles can escape from the shock since the adiabatic focusing produces an outward-directed anisotropy that is large enough to overcome isotropisation by scattering. Combining Eqs. (5) and (10) allows us to obtain an estimate on the position of the boundary of the turbulent foreshock (which can be regarded as a free-escape boundary for the ions) as$$x_{\mathrm{f}\mathrm{eb}}\approx R_{\mathrm{S}}-2x_0,$$(11)if the mean free path has acquired a steady-state length. This, however, is not achieved before the time-integrated wave-growth rate, proportional to the time-integrated net flux of particles, well exceeds unity (Vainio 2003). Thus, the estimate (11) applies only to the late phases of the SEP event at energies below the maximum energy in the spectrum.

Our semi-analytical foreshock model is based on the above theoretical formulation, where we regard the parameters and as free, to be fitted to the simulation results. To take account of the finite acceleration time, we use a spectral form that allows an exponential cut-off at high energies. The location of the free-escape boundary, , was analysed by Battarbee et al. (2011) using$$\mathrm{\lambda }(x_{\mathrm{f}\mathrm{eb}},E)/L\left(x_{\mathrm{f}\mathrm{eb}}\right)=3(u_1-v_{\mathrm{A}})/v.$$(12)

They found that the boundary is following the contours of the particle intensity relatively closely, and this holds true also in our simulation modelling. We, therefore, use this as our model assumption.

Instead of momentum, when analysing the simulation results, we use energy as an independent variable and instead of particle distribution function we use particle intensity. Thus, the correspondence with the spectral indices of the power laws and at non-relativistic energies is$$\beta =\frac{\sigma -2}{2}\iff \sigma =2\beta +2.$$(13)

Thus, at non-relativistic energies$$x_0\propto p^{\mathrm{\sigma }-3}\propto E^{\beta -1/2}.$$(14)

The canonical power-law spectrum at the shock is obtained only in steady state in an infinite one-dimensional medium. The effects of finite acceleration time and adiabatic focusing lead to a finite maximum energy obtained from the shock. Assuming that scattering of particles in the downstream region is much more efficient than in the upstream region, the acceleration rate at a parallel shock is obtained as (Drury 1983)$$\dot{p}=p\frac{{\tilde{u}}_{1n}({\tilde{u}}_{1n}-{\tilde{u}}_{2n})}{3D_{\mathrm{nn}}}=p\frac{(u_1-v_{\mathrm{A}})(1-r_{\mathrm{c}}^{-1})}{3x_0}=p\frac{u_1(r_{\mathrm{c}}-1)}{3x_0{r}_{\mathrm{g}}},$$(15)where is the diffusion coefficient in the shock-normal direction upstream of the shock (here, the value at the shock is used), is the angle between the upstream magnetic field and the shock-normal, and and denote the shock-frame scattering-centre speed along the shock normal upstream and downstream of the shock, respectively. Here, Eq. (5) has been used to derive the second form from the first. Substituting from Eq. (7) and using , where is the shock-normal speed, gives$$\mathrm{d}{p}^{\mathrm{\sigma }-3}=p_{\mathrm{i}\mathrm{nj}}^{\mathrm{\sigma }-3}\frac{u_1\cos {\mathrm{\Theta }}_{\mathrm{Bn}}}{V_{\mathrm{S}}{r}_{\mathrm{g}}}\frac{\pi ϵ\sigma }{2}\frac{{\mathrm{\Omega }}_{\mathrm{p}}}{v_{\mathrm{A}}}\mathrm{d}R.$$(16)

Integrating this over the distance propagated by the shock gives the cut-off energy. In our model, the shocks propagate in a region where the solar-wind speed is almost constant. Thus, all other parameters on the right hand side are approximately constants except the ion inertial length, , which scales approximately as . Thus,$$p_{\mathrm{c}}^{\sigma -3}=p_{\mathrm{i}\mathrm{nj}}^{\mathrm{\sigma }-3}\left[1+\frac{u_1\cos {\mathrm{\Theta }}_{\mathrm{Bn}}}{V_{\mathrm{S}}{r}_{\mathrm{g}}}\frac{\pi ϵ\sigma }{2}\frac{{\mathrm{\Omega }}_{\mathrm{p}}}{v_{\mathrm{A}}}{R}_{\mathrm{S}}\ln \frac{R_{\mathrm{S}}}{R_{\mathrm{S}0}}\right],$$(17)where the second term in the brackets dominates under normal conditions.

Equation (11) shows, however, that even in a steady state the power law cannot extend to infinity, since turbulent trapping in the foreshock becomes inefficient as . This condition gives another estimate for the cut-off momentum:$$p_{\mathrm{c}}^{\sigma -3}=p_{\mathrm{i}\mathrm{nj}}^{\sigma -3}\frac{\pi ϵ\sigma }{4}\frac{{\mathrm{\Omega }}_{\mathrm{p}}{R}_{\mathrm{S}}}{v_{\mathrm{A}}}.$$(18)

We note that if$$\ln \frac{R_{\mathrm{S}}}{R_{\mathrm{S}0}}=\frac{r_{\mathrm{g}}}{2} \frac{V_{\mathrm{S}}}{u_1\cos {\mathrm{\Theta }}_{\mathrm{Bn}}},$$(19)the two estimates agree. Thus, if the ratio of the final to the initial shock distance is of the order of a few, the shock should have enough time to achieve a steady state. Note, however, that the estimate given for the acceleration rate is actually an overestimate, as it uses the steady-state diffusion coefficient at the shock for the whole upstream region, when in reality the diffusion coefficient increases linearly away from the shock and also evolves (from higher values towards the steady-state value) as a function of time. Thus, numerical factors in the scaling laws are not to be taken as exact and the scaling expected from the theory has to be calibrated using simulations.

In summary, the equations we use in our semi-analytical foreshock model, applicable to non-relativistic energies, are$$j\left(x,E\right)=j_{\mathrm{S}}\frac{x_0}{x+x_0}{\left(\frac{R_{\mathrm{S}}}{R_{\mathrm{S}}+x}\right)}^2,$$(20) $$x_0\left(E\right)=a(E/E_0{)}^{\beta -1/2},$$(21) $$\mathrm{\lambda }(x,E)=\frac{3(u_1-v_{\mathrm{A}})}{v}(x+x_0)e^{x/x_{\mathrm{f}\mathrm{eb}}},$$(22) $$j_{\mathrm{S}}\left(E\right)=C{E}^{-\beta }\exp \left\{-{\left(\frac{E}{E_{\mathrm{c}}}\right)}^{\alpha }\right\},$$(23) $$x_{\mathrm{f}\mathrm{eb}}\left(E\right)=\left[\right(v{t}_{\max }{)}^{-2}+x_1^{-2}{]}^{-1/2},$$(24) $$x_1\left(E\right)=x_0\left(\frac{j_{\mathrm{S}}}{j_{\mathrm{f}\mathrm{eb}}}-1\right),$$(25)where , , , , and are the model parameters to be determined by simulations and is a reference energy chosen from the power-law part of the fitted particle spectrum. The exponential factor in Eq. (22) is introduced to facilitate particle escape at and the first term inside the brackets in Eq. (24) is inserted to make sure that the free-escape boundary is not further from the shock than the fastest propagating escaping particles.

4. Results

4.1. Simulation runs

The simulation runs for the determination of the parameters of the semi-analytical foreshock model were performed for twelve runs that cover the radial distances from 5.8 solar radii up to 60 solar radii. We simulate the shock at four different distances from the Sun so that each run covers approximately a time , where is the initial position of the shock. We fix and to values that correspond to a shock that gradually becomes more oblique mimicking the winding of the interplanetary spiral magnetic field (although the CSA model itself has only a radial field implemented). The solar-wind speed at AU is fixed to , the magnetic field to , and the solar-wind proton density to . In addition, we vary the parameters of the seed-particle population to study the effect of its spectral form and density on the acceleration process. The cases considered are presented in Table 1. Here, we consider a compound seed population with two values of , allowing one to vary the spectral shape of the distribution along with the density. The seed-particle densities and give the densities of the populations that are distributed in velocity space following the kappa distribution with and , respectively. Both are scaled with solar-wind proton density . The table also gives the gas compression ratio , the magnetic compression ratio and the scattering-centre compression ratio , which the simulation code computes from the Rankine-Hugoniot relations.

Note that in our simulations, we consider the densities of the kappa-distributed seed populations to be less than the solar-wind density, the total density of seed protons being either 1.1% (runs a and b) or 2% (runs c) of the total solar-wind proton density. This down-scaling of the distribution is made, because the Monte-Carlo model does not treat the low-energy plasma in a fully self-consistent manner, and this probably leads to an unrealistically high injection probability in quasi-parallel shocks for particles that would belong to the thermal pool (Battarbee et al. 2013).

4.2. Determination of model parameters

4.2.1. Fitting of the particle energy spectrum at the shock

We have fitted the simulated energy spectra at the shock to the curve$$j_{\mathrm{S}}\left(E\right)=C{E}^{-\beta }\exp \{-(E/E_c{)}^{\alpha }\},$$(26)where , , and are fitting parameters. As the cut-off of the spectrum is quite steep (i.e., is larger than 1), we have divided the fitting in two parts. First, we fit the parameters and using the power-law part of the spectrum, for which the points are selected through visual inspection. After that, we fit the remaining parameters using the variable$$z≔\ln \{-\ln (j_{\mathrm{S}}{C}^{-1}{E}^{\beta }\left)\right\}=\alpha y+b,$$(27)where and . Note that the argument of the outer logarithm will be negative for some points in the power-law part of the fit. Therefore, we include to the latter fit only those points, , of the spectrum for which . The selection of points is done, again, by visual inspection. An example of the fit is given in Figure 1. The plot presents differential intensities and the applied unit is abbreviated as .

Fig. 1. An example of the fit to the energy spectrum at the shock. The case considered is Run 2c.

In addition, we have fitted the distribution by holding the power-law spectral index constant, as obtained from the DSA theory, i.e., using the value of at the time , since this fixes one unobservable parameter from the model. The fitted and predicted values are, in general, quite close to each other and the quality of the fits, using the theoretical prediction for , is similar as when treating the power-law index as a free parameter.

4.2.2. Fitting of the spatial particle distribution

The spatial distribution close to the shock was fitted to the curve$$j_E\left(x\right)={\left(\frac{R_0}{R}\right)}^2\frac{j_0{x}_0}{x+x_0},$$(28)where is the heliocentric distance, , is the distance from the shock along the (radial) mean magnetic field and and are fitting parameters. In practice, the result is obtained by fitting the simulated points of to a straight line , where and are fitting parameters. We first plot the simulated points and choose visually that part of the distribution (close to the shock), where the fitting function seems to be obeyed. After the fitting of the points, we also visually check the final result for the part that was not included in the fit. Typically, the distribution function is well represented by the fit close to the shock, but departs from it at larger distances because the regions further out from the shock take a longer time to reach an equilibrium and are also affected by the radially expanding geometry. We fit the spatial distribution using the energy and scale the value of like a power law of energy, as predicted by theory. An example of the fit to the points is given in Figure 2.

Fig. 2. An example of the fit to the spatial distribution ahead of the shock at E = E0. The case considered is Run 2c.

4.2.3. Detemination of injection strength

The simulation model uses a model of injection, in which not all seed particles swept up by the shock are actually injected into the acceleration process. Instead, a considerable part of the low-energy particles are transmitted through the shock and advected downstream to eventually form the downstream thermal population. In order to estimate the injection efficiency of the shock, , we have used the power-law fit of the intensity at the shock, , and compared it with Eq. (6), which gives a relation between and as$$C=\frac{(\beta \mathrm{\prime }+1)ϵn_{\mathrm{p}}}{2\pi m{v}_{\mathrm{i}\mathrm{nj}}}{\left(\frac{1}{2}m{v}_{\mathrm{i}\mathrm{nj}}^2\right)}^{\beta \mathrm{\prime }}$$(29)as long as the injection speed is fixed. As all our shocks are quasi-parallel, we use the assumption that the injection speed (in the upstream rest frame) is , which corresponds to particles entering the downstream region along the field at the fluid speed and making a 180 degree turn in the downstream region while conserving their energy in the downstream rest frame and returning back to the upstream, thus becoming injected into the acceleration process.

4.2.4. Fitting of the foreshock boundary

After fitting the spectral parameters, the foreshock boundary resulting from applying the Eq. (12) to the simulation results is fitted by choosing the contour of the analytically modelled intensity that best matches the criterion (12). This is done by changing the parameter until the curve (25) matches the position of the boundary as well as possible. By considering all the simulated cases, it turns out that the level can be given as:$$j_{\mathrm{f}\mathrm{eb}}=C{E}_{\mathrm{c}}^{-\beta \mathrm{\prime }}\exp (-3.0).$$(30)

An example of the fit to the points is given in Figure 3. The solid curve is the result of Eq. (24) and it follows the behavior at low energies (below  MeV). The dashed curve, giving the foreshock boundary in the simulation, behaves differently at low energies because, in fact, background wave intensities are high enough to trap the particles close to the shock and prevent them from escaping at all.

Fig. 3. An example of the fit to the foreshock boundary. The contours give the particle intensity in pfu. The dashed curve gives the foreshock boundary resulting from the simulation and the solid curve the one resulting from the semi-analytical modelling. The case considered is Run 2c.

4.2.5. Simulated model parameters

The results of the fits performed on the simulations are given in Table 2. The first two columns give the ion inertial length and de Hoffmann-Teller frame velocity of the upstream plasma in the end of the simulation (location in Table 1). The remaining columns give the fitted parameters of the model.

We note that the simulated spectra are in general slightly softer than the theoretical predictions (). The differences in the spectral indices are, however, not large and, as discussed above, the spectra were also fitted well by holding the power-law spectral index fixed to the theoretical value. The rest of the parameters given in the Table are obtained using the theoretical spectral index . We also note that the results from the runs a and b with are very close to one another. Since the cases a and b have the same seed-particle density with different spectrum, this indicates that shape of the seed-particle spectrum has little influence on the results in quasi-parallel shocks.

4.3. Interdependencies of the parameters

Next, we compare the theoretical values of the key parameters of the foreshock model with simulations. In Figure 4, we present the simulated value of the foreshock scale height (see Eq. (21)) against the theoretical value$$a_{\mathrm{t}\mathrm{h}}=\frac{2}{\pi ϵ\sigma \mathrm{\prime }}\frac{v_{\mathrm{A}}}{{\mathrm{\Omega }}_{\mathrm{p}}}{\left(\frac{E_0}{E_{\mathrm{i}\mathrm{nj}}}\right)}^{\beta \mathrm{\prime }-\frac{1}{2}}.$$(31)

Fig. 4. Comparison of the simulated and the theoretical values of the foreshock parameter a. The solid line gives the power-law fit and the dashed red line the linear approximation between the quantities.

We see that the theoretical and simulated values match with good accuracy. The exponent 1.06 is close to a linear relation between the simulated and theoretical values. A linear law of can be applied without introducing a large error.

In Figure 5 we compare the cut-off momenta obtained from the simulations (i.e., ) against the theoretical prediction of Eq. (17). The value of used in the theoretical prediction is the one from the fit. Here, we see that the simulated cut-off momentum falls short of the theoretical prediction up to about one order of magnitude. However, the points cluster very well around the line$$p_{\mathrm{c},\mathrm{sim}}=0.105 p_{\mathrm{c},\mathrm{th}},$$(32)which allows us to calibrate the theory against simulations. We also fit the two values using a power law (shown as the red curve in Fig. 5). Since the power-law index of the fit is close to unity, we infer that while the simplified theory (with the steady-state at the shock representing the whole upstream region) grossly overestimates the acceleration rate, it still gives a valid representation of the scaling of the cut-off momentum.

Fig. 5. Comparison of the simulated and the theoretical values of the spectral parameter pc.

4.4. Determining the foreshock model parameters

A method to use the observed cut-off energies as an aid to estimate the foreshock parameters remotely now emerges.

1. Determine the empirical cut-off energy, , of the source, e.g., from a fit of the observed particle intensities at 1 AU. Use it to calculate ).

2. Use the inverse relation of Eq. (32), i.e.,$$p_{\mathrm{c},\mathrm{th}}=\frac{p_{\mathrm{c},\mathrm{emp}}}{0.105}=9.52 p_{\mathrm{c},\mathrm{emp}},$$(33)to estimate the corresponding value of the cut-off momentum in the simplified steady-state DSA theory.

3. Invert the Eq. (17) for the theoretical cut-off momentum to get the value of as$$ϵ=\left[{\left(\frac{p_{\mathrm{c},\mathrm{th}}}{p_{\mathrm{i}\mathrm{nj}}}\right)}^{\sigma \mathrm{\prime }-3}-1\right]\frac{2V_{\mathrm{S}}{r}_{\mathrm{g}}{v}_{\mathrm{A}}}{u_1\cos {\mathrm{\Theta }}_{\mathrm{Bn}}\sigma \pi {\mathrm{\Omega }}_{\mathrm{p}}{R}_{\mathrm{S}}\ln (R_{\mathrm{S}}/R_{\mathrm{S}0})},$$(34)where all parameters of the equation can be estimated from (MHD) modelling of the shock. If using MHD simulations of shock propagation, however, note that the correspondence between the simulated and theoretical cut-off energies is calibrated using compression ratios determined from Rankine-Hugoniot equations (with adiabatic index set to 5/3, which might not apply in the MHD simulation model).

4. After has been found, use it to compute$$C=\frac{(\beta \mathrm{\prime }+1)ϵn_{\mathrm{p}}}{2\pi m{v}_{\mathrm{i}\mathrm{nj}}}{E}_{\mathrm{i}\mathrm{nj}}^{\beta \mathrm{\prime }}.$$(35)where and$$a=\frac{1.36}{\pi ϵ(\beta \mathrm{\prime }+1)}\frac{v_{\mathrm{A}}}{{\mathrm{\Omega }}_{\mathrm{p}}}{\left(\frac{E_0}{E_{\mathrm{i}\mathrm{nj}}}\right)}^{\beta \mathrm{\prime }-\frac{1}{2}}.$$(36)

5. Use the values of and in Eqs. (23) and (25) to determine the foreshock boundary.

The values given above fit the simulations well, but when analysing future in-situ observations, they can also be considered as parameters to fit. For example, a lower value of yields a softer cut-off of the spectrum, which may fit the observations better. An energy-dependent value of would allow tuning the modulation effect of the foreshock. For example, the functional form:$$j_{\mathrm{f}\mathrm{eb}}=C{E}^{-\delta }{E}_{\mathrm{c}}^{\delta -\beta \mathrm{\prime }}\exp (-\mathrm{\zeta })$$(37)with would allow modelling situations, where the modulating effect of the foreshock is less severe at low energies than in the CSA simulations. Here, is a fitting parameter as well. This could be considered if ever an empirical model of the foreshock would be constructed based on the same formulation, or if a larger number of simulations with a more diverse set of input parameters was run to have enough data to increase the number of parameters of the semi-analytical model.

As an example of this procedure, we present the modeled foreshock boundary for all 12 simulated cases in Figure 6. The values of the foreshock model are obtained using the fitted value of the for each simulation but deriving the rest of the model parameters following the steps 1–5 above. As can be seen, the model reproduces the simulations quite accurately for all cases considered, demonstrating the ability of the simple model to give a valid prescription of the foreshock under different conditions, based on a single remotely observable parameter, . In case we allow tuning of the parameters and , we even receive a slightly better fit using and .

Fig. 6. The location of the foreshock boundary as obtained from the semi-analytical model (solid curves) compared with the location obtained from each simulation (dashed curve). The three adjacent panels present the comparison for the runs 1a, 1b and 1c (top row) to runs 4a, 4b and 4c (bottom row). The colour scale is the same as in Figure 2.

4.5. Modelled mean free path

Figure 7 presents a comparison between the simulation result and the semi-analytical model of the proton mean free path. It is seen that there is an excellent agreement close to the shock. At larger distances, the simulated values are larger than the modelled ones, as expected, since the ambient mean free path is assumed to be infinite in the semi-analytical model. Inside the foreshock boundary, the curves differ by less than 20%.

Fig. 7. The modelled mean free path. The dashed curve gives the CSA-simulated mean free path of 2.84-MeV protons and the solid curve the result from the semi-analytical model. The vertical dot-dashed line gives the position of the foreshock boundary. The considered case is Run 2c.

4.6. Modelled peak intensity at the shock

Figure 8 presents the peak intensity of an ESP event at 1 MeV, using the values of shock and solar-wind parameters of our simulation model in Eq. (29). We have included only the simulations where the total seed-particle density (combining the two -distributions) corresponds to 1.1% (i.e., 0.1% + 1.0%) of the proton density (i.e., excluding simulation runs 1c, 2c, 3c and 4c, where the density is 2%). We see that the scaling of the intensity towards the Sun is rather steep, on average over the simulated distance range. These results should not be taken as a definite prediction of our simulation model, since the case for shock propagation considered here does not correspond to a generic or even an average case in any sense. However, Eq. (29) can be easily used to predict the ESP peak intensity at the shock, once the parameter has been determined using observations. This is because the stationary state, as described by the DSA theory, is well achieved in our simulations close to the shock at energies well below the spectral cut-off.

Fig. 8. The modelled scaling of an ESP event. The points give the intensity of the ESP event at 1-MeV proton energy as a function of radial distance resulting from the simulations.

5. Discussion

We have developed a semi-analytical model for the foreshock of a travelling CME-driven shock that stays quasi-parallel during its propagation from the outer corona to the interplanetary medium. The model is based on the theory of diffusive shock acceleration and on earlier simulation results, which allow us to estimate the mean free path of protons in the foreshock and the boundary of the foreshock region. The model parameters can be estimated from shock properties and from the observable value of the cut-off energy in the particle source function at the shock.

The simulations show that the foreshock scale height is well represented by the DSA theory but that the cut-off energy of the accelerated particle spectrum is considerably lower than the simple theoretical prediction. A well-ordered linear relation between the theoretical (over-)estimate and the simulation results was, however, found. This relation allows us to renormalise the observed cut-off energies to the theoretical ones and thus extract parameters of the injection using the observed cut-off momentum produced by the shock during its evolution. This is an encouraging result for the remote sensing of the foreshock region.

The good fit obtained between the simulated and theoretical cut-off momentum and the fact that the cut-off momentum is significantly lower than values from Eq. (18), imply that none of the simulation runs were yet limited in maximum energy by escape facilitated by focusing. Note, however, that if a shock propagates in a medium with greatly varying properties, like the solar-wind acceleration region, the acceleration to the highest energies might occur at low altitudes. This would result in the later phases of shock evolution being controlled by the steady-state escape condition with a decreasing value of the cut-off energy. But even in that case, the cut-off energy can be directly related to the value of and, thus, an estimate of the evolution of the foreshock strength as a function of radial distance can be obtained.

The cases presented in this study are based on a set of snapshot simulations of shock propagating over a limited distance between 5.8 and 60 solar radii. One of the features of the semi-analytical model is that the cut-off energy of the energy spectrum keeps increasing almost linearly as the shock propagates outwards. This is a result of the constant injection efficiency of the shock. It would be of interest to study cases where the shock properties change over the acceleration time as this would lead to a variable value of . This requires re-integration of the equation describing the increase of the cut-off momentum as a function of radial distance. A power-law dependence of , for example, gives a prediction for the cut-off momentum of$$\mathrm{d}{p}_{\mathrm{c}}^{\sigma -3}=p_{\mathrm{i}\mathrm{nj}}^{\sigma -3}\frac{u_1\cos {\mathrm{\Theta }}_{\mathrm{Bn}}}{V_{\mathrm{S}}{r}_{\mathrm{g}}}\frac{\pi {ϵ}_0\sigma }{2R_0^{\xi -1}}\frac{{\mathrm{\Omega }}_{\mathrm{p}0}}{v_{\mathrm{A}0}}{R}^{\xi -1}\mathrm{d}R,$$(38)where the subscript 0 refers to values at 1 AU. Thus,$$p_{\mathrm{c}}^{\sigma -3}=p_{\mathrm{i}\mathrm{nj}}^{\sigma -3}\frac{u_1\cos {\mathrm{\Theta }}_{\mathrm{Bn}}}{V_{\mathrm{S}}{r}_{\mathrm{g}}}\frac{\pi {ϵ}_0\sigma }{2\xi R_0^{\xi -1}}\frac{{\mathrm{\Omega }}_{\mathrm{p}0}}{v_{\mathrm{A}0}}(R_{\mathrm{S}}^{\xi }-R_{\mathrm{S}0}^{\xi }),$$where is the starting radial distance of the shock. Clearly, a negative value of would result in a cut-off momentum that first increases as a function of distance and then tends to a constant, given by$$p_{\mathrm{c}}^{\sigma -3}=p_{\mathrm{i}\mathrm{nj}}^{\sigma -3}\frac{u_1\cos {\mathrm{\Theta }}_{\mathrm{Bn}}}{V_{\mathrm{S}}{r}_{\mathrm{g}}}\frac{\pi {ϵ}_0\sigma }{2(-\xi )}\frac{{\mathrm{\Omega }}_{\mathrm{p}0}{R}_0}{v_{\mathrm{A}0}}{\left(\frac{R_0}{R_{S0}}\right)}^{-\xi }.$$

Under such conditions, focusing-driven particle escape (see Eq. (18)) will eventually limit the cut-off momentum, as that decreases strongly as a function of distance from the Sun. For the variable- model Eq. (18) becomes:$$p_{\mathrm{c}}^{\sigma -3}=p_{\mathrm{i}\mathrm{nj}}^{\sigma -3}\frac{\pi {\epsilon }_0\sigma }{4}\frac{{\mathrm{\Omega }}_{\mathrm{p}0}{R}_0}{v_{\mathrm{A}0}}{\left(\frac{R_0}{R_{\mathrm{S}}}\right)}^{-\xi }.$$(39)

Assuming a value of and that the acceleration-time-limited cut-off momentum estimate is one order of magnitude too large (as indicated by our simulations), the focusing-driven particle escape will start to limit the cut-off momentum when$$\frac{R_{\mathrm{S}}}{R_{\mathrm{S}0}}\gtrsim \frac{V_{\mathrm{S}}{r}_{\mathrm{g}}}{0.2 u_1\cos {\mathrm{\Theta }}_{\mathrm{Bn}}}\succsim 20,$$but this estimate would have to be confirmed by further simulations with variable injection efficiency.

Finally, future improvements of the model should also consider a more realistic form of the background magnetic field. By using an Archimedean spiral instead of the radial field we could also consider distances further out in the heliosphere. This would allow comparison of the model with in-situ observations.

6. Conclusions and outlook

We have studied a simple semi-analytical model of a foreshock region of a propagating interplanetary shock driven by a fast CME. Our modelling is based on the steady-state diffusive shock acceleration theory laid out by Bell (1978) but modified to take into account time dependence (of the particle cut-off momentum) and finite extent of the foreshock region. The modifications are parametrised in an ad hoc manner, but calibrated using a simulation model computing the coupled evolution of the energetic protons and Alfvén waves in the medium ahead of the shock.

Our simulations indicate that the simple theoretical estimates of the cut-off momentum, based on the self-generated turbulence levels at the shock, predict one order of magnitude too large cut-off momentum. The most plausible reason for the gross overestimate is that the simulated diffusion coefficient is neither in a steady state nor spatially homogeneous. Thus, modifications to the SEP models depending on similar theoretical estimates (e.g., Zank et al. 2000; Li et al. 2003; Rice et al. 2003; Verkhoglyadova et al. 2009, 2010, 2012) are needed.

By calibrating the relation between the theoretical and simulated cut-off momentum, our model allows for remote sensing of the injection efficiency of the shock by observing the cut-off momentum of particles accelerated at the shock at different heliospheric distances. This allows us to construct analytical models for the ESP peak intensities and foreshock turbulence, which greatly facilitate future modelling efforts of large SEP events in the corona and inner heliosphere.

Acknowledgments

The research leading to these results has received funding from the European Union Seventh Framework Programme (FP7/2007-2013) under Grant Agreement No. 262468. CSC – IT Center for Science Ltd. is acknowledged for providing the computation resources for the project.

References

All Tables

All Figures

	Fig. 1. An example of the fit to the energy spectrum at the shock. The case considered is Run 2c.
In the text

	Fig. 2. An example of the fit to the spatial distribution ahead of the shock at E = E0. The case considered is Run 2c.
In the text

	Fig. 3. An example of the fit to the foreshock boundary. The contours give the particle intensity in pfu. The dashed curve gives the foreshock boundary resulting from the simulation and the solid curve the one resulting from the semi-analytical modelling. The case considered is Run 2c.
In the text

	Fig. 4. Comparison of the simulated and the theoretical values of the foreshock parameter a. The solid line gives the power-law fit and the dashed red line the linear approximation between the quantities.
In the text

	Fig. 5. Comparison of the simulated and the theoretical values of the spectral parameter pc.
In the text

	Fig. 6. The location of the foreshock boundary as obtained from the semi-analytical model (solid curves) compared with the location obtained from each simulation (dashed curve). The three adjacent panels present the comparison for the runs 1a, 1b and 1c (top row) to runs 4a, 4b and 4c (bottom row). The colour scale is the same as in Figure 2.
In the text

	Fig. 7. The modelled mean free path. The dashed curve gives the CSA-simulated mean free path of 2.84-MeV protons and the solid curve the result from the semi-analytical model. The vertical dot-dashed line gives the position of the foreshock boundary. The considered case is Run 2c.
In the text

	Fig. 8. The modelled scaling of an ESP event. The points give the intensity of the ESP event at 1-MeV proton energy as a function of radial distance resulting from the simulations.
In the text