Research Article
A semianalytical foreshock model for energetic storm particle events inside 1 AU
^{1}
Department of Physics, University of Helsinki, PO Box 64, FI 00014, Finland
^{2}
Deparment of Physics and Astronomy, University of Turku, FI 20014, Finland
^{3}
Finnish Meteorological Institute, PO Box 503, 00101
Helsinki, Finland
^{4}
Jeremiah Horrocks Institute, University of Central Lancashire, PR1 2HE
Preston, Lancashire, UK
^{*} Corresponding author: rami.vainio@utu.fi
Received:
16
September
2013
Accepted:
3
February
2014
We have constructed a semianalytical model of the energeticion foreshock of a CMEdriven coronal/interplanetary shock wave responsible for the acceleration of large solar energetic particle (SEP) events. The model is based on the analytical model of diffusive shock acceleration of Bell (1978), appended with a temporal dependence of the cutoff momentum of the energetic particles accelerated at the shock, derived from the theory. Parameters of the model are recalibrated using a fully timedependent selfconsistent simulation model of the coupled particle acceleration and Alfvénwave generation upstream of the shock. Our results show that analytical estimates of the cutoff energy resulting from the simplified theory and frequently used in SEP modelling are overestimating the cutoff momentum at the shock by one order magnitude. We show also that the cutoff momentum observed remotely far upstream of the shock (e.g., at 1 AU) can be used to infer the properties of the foreshock and the resulting energetic storm particle (ESP) event, when the shock is still at small distances from the Sun, unaccessible to the insitu observations. Our results can be used in ESP event modelling for future missions to the inner heliosphere, like the Solar Orbiter and Solar Probe Plus as well as in developing acceleration models for SEP events in the solar corona.
Key words: Energetic particle / SEP / Shocks / Heliosphere / Interplanetary medium
© R. Vainio et al., Published by EDP Sciences 2014
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
1. Introduction
Solar energetic particle (SEP) events provide by far most of the >1MeV 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 highintensity 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 CMEdriven 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 macroscale 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 Larmorradius, beyond the MHD description. Thus, it is necessary to find alternative ways to describe the fluctuations undergoing waveparticle 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 1AU 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 selfgenerated 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 waveexcitation 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 CMEdriven 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 selfamplified 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 spaceweather use.
A predecessor of the fully kinetic numerical codes was developed at the turn of the century by Zank and coworkers (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 shockaccelerated 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 (nonlinear 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 PATHlike models, the analytical approximations employed in the methodology would have to be validated using selfconsistent numerical computations. Another modelling approach that has great potential in spaceweather applications is the semiempirical 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 semiempirical 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 seedparticle spectral properties and derive a semianalytical 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 timeexpensive 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 semianalytical 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 WKBbased 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 wavewave interactions and improves numerical stability.
The waveparticle interactions are described at the quasilinear 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 wavegrowth 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 quasilinear resonance condition, k = mΩ/pμ, contains a dependence on the pitchangle 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 solarwind 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 shocknormal 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 seedparticle distribution as a kappa distribution (Battarbee et al. 2011),(1)with its density scaled to the density of the coronal plasma, and the properties of the solar wind and interplanetary field (density, solarwind speed and radial magnetic field) at 1 AU. We consider only protons in our simulation. Note that the thermal velocity w _{0} of the kappa distribution is obtained from the kinetic temperature, T, as(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 userspecified 1AU density) and the profile of the kinetic temperature T(R) are fixed in the simulation to a semiempirical model (Cranmer & van Ballegooijen 2005). The solarwind speed profile u _{sw}(R) is obtained from the userspecified 1AU speed using the conservation of mass, i.e.,(3)
The modelled magnetic field at the considered distances is radial and set to B(R) = B _{0}(R _{0}/R)^{2} where B _{0} is the field at R _{0} = 1 AU.
3. Semianalytical foreshock model
The semianalytical foreshock model is based on the 1D steadystate theory of Bell (1978), which can be given in form (Vainio & Laitinen 2007, 2008; Horne et al. 2013)(4) (5) (6) (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 _{1} 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 _{1} – u _{2}), which results from cold upstream ions making a headon collision with the downstream scattering centres. The quantity x _{0} 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 ^{ −1}) is the scatteringcentre compression ratio assuming that the scattering in the downstream region is governed by frozenin magnetic fluctuations or strong turbulence with zero crosshelicity. Here, M = u _{1n }/v _{An } = u _{1}/v _{A} is the Alfvénic Mach number of the upstream flow, the subscripts n in the velocities refer to the shocknormal components, and r _{g} = ρ _{2}/ρ _{1} 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 quasiperpendicular shocks.
Since the interplanetary magnetic field is not a homogeneous field but scales with the heliocentric distance as R ^{−2} 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 crosssectional area as is taken into account by considering particle intensity as(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)(9)where is the focusing length. When, instead,(10)particles can escape from the shock since the adiabatic focusing produces an outwarddirected 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 freeescape boundary for the ions) as(11)if the mean free path has acquired a steadystate length. This, however, is not achieved before the timeintegrated wavegrowth rate, proportional to the timeintegrated 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 semianalytical 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 cutoff at high energies. The location of the freeescape boundary, , was analysed by Battarbee et al. (2011) using(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 nonrelativistic energies is(13)
Thus, at nonrelativistic energies(14)
The canonical powerlaw spectrum at the shock is obtained only in steady state in an infinite onedimensional 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)(15)where is the diffusion coefficient in the shocknormal direction upstream of the shock (here, the value at the shock is used), is the angle between the upstream magnetic field and the shocknormal, and and denote the shockframe scatteringcentre 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 shocknormal speed, gives(16)
Integrating this over the distance propagated by the shock gives the cutoff energy. In our model, the shocks propagate in a region where the solarwind 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,(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 cutoff momentum:(18)
We note that if(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 steadystate 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 steadystate 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 semianalytical foreshock model, applicable to nonrelativistic energies, are(20) (21) (22) (23) (24) (25)where , , , , and are the model parameters to be determined by simulations and is a reference energy chosen from the powerlaw 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 freeescape 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 semianalytical 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 solarwind speed at AU is fixed to , the magnetic field to , and the solarwind proton density to . In addition, we vary the parameters of the seedparticle 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 seedparticle 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 solarwind proton density . The table also gives the gas compression ratio , the magnetic compression ratio and the scatteringcentre compression ratio , which the simulation code computes from the RankineHugoniot relations.
Simulation parameters. The first column gives the label of the run, the second column the shocknormal speed, the third the shocknormal angle, the fourth the initial position of the shock, the fifth the simulation time, the sixth the final distance of the shock, the seventh and eight columns the parametrisation of the seed population, the ninth and tenth columns give the shock’s gas compression ratio and magnetic compression ratio and the eleventh column its scatteringcentre compression ratio.
Note that in our simulations, we consider the densities of the kappadistributed seed populations to be less than the solarwind density, the total density of seed protons being either 1.1% (runs a and b) or 2% (runs c) of the total solarwind proton density. This downscaling of the distribution is made, because the MonteCarlo model does not treat the lowenergy plasma in a fully selfconsistent manner, and this probably leads to an unrealistically high injection probability in quasiparallel 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(26)where , , and are fitting parameters. As the cutoff 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 powerlaw part of the spectrum, for which the points are selected through visual inspection. After that, we fit the remaining parameters using the variable(27)where and . Note that the argument of the outer logarithm will be negative for some points in the powerlaw 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 powerlaw 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 powerlaw 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(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 = E_{0}. 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 lowenergy 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 powerlaw fit of the intensity at the shock, , and compared it with Eq. (6), which gives a relation between and as(29)as long as the injection speed is fixed. As all our shocks are quasiparallel, 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:(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 semianalytical 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 HoffmannTeller 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.
Simulation results. The unit of intensity is .
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 powerlaw 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 seedparticle density with different spectrum, this indicates that shape of the seedparticle spectrum has little influence on the results in quasiparallel 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(31)
Fig. 4.
Comparison of the simulated and the theoretical values of the foreshock parameter a. The solid line gives the powerlaw 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 cutoff 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 cutoff momentum falls short of the theoretical prediction up to about one order of magnitude. However, the points cluster very well around the line(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 powerlaw index of the fit is close to unity, we infer that while the simplified theory (with the steadystate at the shock representing the whole upstream region) grossly overestimates the acceleration rate, it still gives a valid representation of the scaling of the cutoff momentum.
Fig. 5.
Comparison of the simulated and the theoretical values of the spectral parameter p_{c}. 
4.4. Determining the foreshock model parameters
A method to use the observed cutoff energies as an aid to estimate the foreshock parameters remotely now emerges.

Determine the empirical cutoff 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.,(33)to estimate the corresponding value of the cutoff momentum in the simplified steadystate DSA theory.
3. Invert the Eq. (17) for the theoretical cutoff momentum to get the value of as(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 cutoff energies is calibrated using compression ratios determined from RankineHugoniot 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(35)where and(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 insitu observations, they can also be considered as parameters to fit. For example, a lower value of yields a softer cutoff of the spectrum, which may fit the observations better. An energydependent value of would allow tuning the modulation effect of the foreshock. For example, the functional form:(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 semianalytical 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 semianalytical 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 semianalytical 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 semianalytical model. Inside the foreshock boundary, the curves differ by less than 20%.
Fig. 7.
The modelled mean free path. The dashed curve gives the CSAsimulated mean free path of 2.84MeV protons and the solid curve the result from the semianalytical model. The vertical dotdashed 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 solarwind parameters of our simulation model in Eq. (29). We have included only the simulations where the total seedparticle 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 cutoff.
Fig. 8.
The modelled scaling of an ESP event. The points give the intensity of the ESP event at 1MeV proton energy as a function of radial distance resulting from the simulations. 
5. Discussion
We have developed a semianalytical model for the foreshock of a travelling CMEdriven shock that stays quasiparallel 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 cutoff 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 cutoff energy of the accelerated particle spectrum is considerably lower than the simple theoretical prediction. A wellordered linear relation between the theoretical (over)estimate and the simulation results was, however, found. This relation allows us to renormalise the observed cutoff energies to the theoretical ones and thus extract parameters of the injection using the observed cutoff 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 cutoff momentum and the fact that the cutoff 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 solarwind 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 steadystate escape condition with a decreasing value of the cutoff energy. But even in that case, the cutoff 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 semianalytical model is that the cutoff 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 reintegration of the equation describing the increase of the cutoff momentum as a function of radial distance. A powerlaw dependence of , for example, gives a prediction for the cutoff momentum of(38)where the subscript 0 refers to values at 1 AU. Thus,where is the starting radial distance of the shock. Clearly, a negative value of would result in a cutoff momentum that first increases as a function of distance and then tends to a constant, given by
Under such conditions, focusingdriven particle escape (see Eq. (18)) will eventually limit the cutoff momentum, as that decreases strongly as a function of distance from the Sun. For the variable model Eq. (18) becomes:(39)
Assuming a value of and that the accelerationtimelimited cutoff momentum estimate is one order of magnitude too large (as indicated by our simulations), the focusingdriven particle escape will start to limit the cutoff momentum whenbut 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 insitu observations.
6. Conclusions and outlook
We have studied a simple semianalytical model of a foreshock region of a propagating interplanetary shock driven by a fast CME. Our modelling is based on the steadystate diffusive shock acceleration theory laid out by Bell (1978) but modified to take into account time dependence (of the particle cutoff 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 cutoff momentum, based on the selfgenerated turbulence levels at the shock, predict one order of magnitude too large cutoff 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 cutoff momentum, our model allows for remote sensing of the injection efficiency of the shock by observing the cutoff 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/20072013) under Grant Agreement No. 262468. CSC – IT Center for Science Ltd. is acknowledged for providing the computation resources for the project.
References
 Aran, A., B. Sanahuja, and D. Lario, Fluxes and fluences of SEP events derived from SOLPENCO, Ann. Geophys., 23, 3047–3053, 2005. [CrossRef] (In the text)
 Battarbee, M., Acceleration of Solar Energetic Particles in coronal shocks through selfgenerated turbulence, PhD Thesis, University of Turku, Annales Universitatis Turkuensis AI 475, 126 pp, 2013 (In the text)
 Battarbee, M., T. Laitinen, and R. Vainio, Heavyion acceleration and selfgenerated waves in coronal shocks, A&A, 535, A34, 2011. [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
 Battarbee, M., R. Vainio, T. Laitinen, and H. Hietala, Injection of thermal and suprathermal seed particles in coronal shocks of varying obliquity, A&A, 558, A110, 2013. [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
 Bell, A.R., The acceleration of cosmic rays in shock fronts. I, Mon. Not. R. Astron. Soc., 182, 147–156, 1978. (In the text)
 Cranmer, S.R., and A.A. van Ballegooijen, On the generation, propagation, and reflection of Alfvén Waves from the solar photosphere to the distant heliosphere, Astrophys. J. Suppl., 156, 265–293, 2005. [NASA ADS] [CrossRef] (In the text)
 Drury, L. O’C., An introduction to the theory of diffusive shock acceleration of energetic particles in tenuous plasmas, Rep. Prog. Phys., 46, 973–1027, 1983. [NASA ADS] [CrossRef] (In the text)
 Feynman, J., G. Spitale, J. Wang, and S. Gabriel, Interplanetary proton fluence model–JPL 1991, J. Geophys. Res., 98, 13281–13294, 1993. [CrossRef] (In the text)
 Horne, R., S. Glauert, N. Meredith, H. Koskinen, R. Vainio, et al., Forecasting the Earth’s radiation belts and modelling solar energetic particle events: Recent results from SPACECAST, J. Space Weather Space Clim., 3, A20, 2013. [CrossRef] [EDP Sciences] (In the text)
 Lario, D., and R.B. Decker, Estimation of solar energetic proton missionintegrated fluences and peak intensities for missions traveling close to the Sun, Space Weather, 9, S11003, 2011. [CrossRef] (In the text)
 Lario, D., B. Sanahuja, and A.M. Heras, Energetic particle events: Efficiency of interplanetary shocks as 50 keV < E < 100 MeV proton accelerators, Astrophys. J., 509, 415–434, 1998. [CrossRef] (In the text)
 Lario, D., M.B. Kallenrode, R.B. Decker, E.C. Roelof, S.M. Krimigis, A. Aran, and B. Sanahuja, Radial and longitudinal dependence of solar 4–13 MeV and 27–37 MeV proton peak intensities and fluences: Helios and IMP 8 observations, Astrophys. J., 653, 1531–1544, 2006. [NASA ADS] [CrossRef] (In the text)
 Lario, D., A. Aran, N. Agueda, and B. Sanahuja, Radial dependence of proton peak intensities and fluences in SEP events: Influence of the energetic particle transport parameters, Adv. Space Res., 40 (3), 289–294, 2007. [CrossRef] (In the text)
 Lee, M.A., Coupled hydromagnetic wave excitation and ion acceleration at interplanetary traveling shocks, J. Geophys. Res., 88, 6109–6119, 1983. [NASA ADS] [CrossRef] (In the text)
 Lee, M.A., Coupled hydromagnetic wave excitation and ion acceleration at an evolving coronal/interplanetary shock, Astrophys. J. Suppl. Ser., 158, 38–67, 2005. [NASA ADS] [CrossRef] (In the text)
 Li, G., G.P. Zank, and W.K.M. Rice, Energetic particle acceleration and transport at coronal mass ejectiondriven shocks, J. Geophys. Res., 108, 1082, 2013. [CrossRef] (In the text)
 Ng, C.K., and D.V. Reames, Focused interplanetary transport of approximately 1 MeV solar energetic protons through selfgenerated Alfvén waves, Astrophys. J., 424, 1032–1048, 1994. [NASA ADS] [CrossRef] (In the text)
 Ng, C.K., and D.V. Reames, Shock acceleration of solar energetic protons: The first 10 minutes, Astrophys. J., 686, L123–L126, 2008. [NASA ADS] [CrossRef] (In the text)
 Pomoell, J., and R. Vainio, Influence of solar wind heating formulations on the properties of shocks in the corona, Astrophys. J., 745, 151, 2012. [NASA ADS] [CrossRef] (In the text)
 Pomoell, J., R. Vainio, and R. Kissmann, MHD simulation of the evolution of shock structures in the solar corona: implications for coronal shock acceleration, Astrophys. Space Sci. Trans., 7, 387–394, 2011. [NASA ADS] [CrossRef] (In the text)
 Reames, D.V., Particle acceleration at the Sun and in the heliosphere, Space Sci. Res., 90, 413–491, 1999. (In the text)
 Reames, D.V., and C.K. Ng, Streaminglimited intensities of solar energetic particles, Astrophys. J., 504, 1002–1005, 1998. [CrossRef] (In the text)
 Rice, W.K.M., G.P. Zank, and G. Li, Particle acceleration and coronal mass ejection driven shocks: Shocks of arbitrary strength, J. Geophys. Res., 108, 1369, 2003. [CrossRef] (In the text)
 Vainio, R., On the generation of Alfvén waves by solar energetic particles, A&A, 406, 735–740, 2003. [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
 Vainio, R., and T. Laitinen, Monte Carlo simulations of coronal diffusive shock acceleration in selfgenerated turbulence, Astrophys. J., 658, 622–630, 2007. [NASA ADS] [CrossRef] (In the text)
 Vainio, R., and T. Laitinen, Simulations of coronal shock acceleration in selfgenerated turbulence, J. Atmos. Sol. Terr. Phys., 70, 467–474, 2008. [NASA ADS] [CrossRef] (In the text)
 Vainio, R., N. Agueda, A. Aran, and D. Lario, Modeling of solar energetic particles in interplanetary space. In: Space Weather, J., Lilensten, Editor, Astrophys. Space Sci. Lib, 344, Springer: Dordrecht, 27–38, 2007. [CrossRef] (In the text)
 Verkhoglyadova, O.P., G. Li, G.P. Zank, Q. Hu, and R.A. Mewaldt, Using the path code for modeling gradual SEP events in the inner heliosphere, Astrophys. J., 693, 894–900, 2009. [CrossRef] (In the text)
 Verkhoglyadova, O.P., G. Li, G.P. Zank, Q. Hu, C.M.S. Cohen, R.A. Mewaldt, G.M. Mason, D.K. Haggerty, T.T. von Rosenvinge, and M.D. Looper, Understanding large SEP events with the PATH code: Modeling of the 13 December 2006 SEP event, J. Geophys. Res., 115, A12103, 2010. [NASA ADS] [CrossRef] (In the text)
 Verkhoglyadova, O.P., G. Li, X. Ao, and G.P. Zank, Radial dependence of peak proton and iron fluxes in solar energetic particle events: Application of the PATH code, Astrophys. J., 757, 75, 2012. [CrossRef] (In the text)
 Zank, G.P., W.K.M. Rice, and C.C. Wu, Particle acceleration and coronal mass ejection driven shocks: A theoretical model, J. Geophys. Res., 105, 25079–25096, 2000. [NASA ADS] [CrossRef] (In the text)
Cite this article as: Vainio R, Pönni A, Battarbee M, Koskinen H.E.J, Afanasiev A, Laitinen T: A semianalytical foreshock model for energetic storm particle events inside 1 AU. J. Space Weather Space Clim., 2014, 4, A08.
All Tables
Simulation parameters. The first column gives the label of the run, the second column the shocknormal speed, the third the shocknormal angle, the fourth the initial position of the shock, the fifth the simulation time, the sixth the final distance of the shock, the seventh and eight columns the parametrisation of the seed population, the ninth and tenth columns give the shock’s gas compression ratio and magnetic compression ratio and the eleventh column its scatteringcentre compression ratio.
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 = E_{0}. 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 semianalytical 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 powerlaw 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 p_{c}. 

In the text 
Fig. 6.
The location of the foreshock boundary as obtained from the semianalytical 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 CSAsimulated mean free path of 2.84MeV protons and the solid curve the result from the semianalytical model. The vertical dotdashed 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 1MeV proton energy as a function of radial distance resulting from the simulations. 

In the text 