THE RELATIONSHIP BETWEEN THE PARTICLE INJECTION RATE AND THE DISPERSION OF THE SCATTERING ANGULAR DISTRIBUTION

Gradual solar energetic particles with a power-law energy spectrum are generally thought to be accelerated by the first Fermi acceleration in interplanetary shocks (Axford et al. 1978; Krymsky 1977; Bell 1978; Blandford & Ostriker 1978). It is well known that a diffusive shock can efficiently accelerate particles by means of the accelerated particles scattering off the instability of Alfvén waves which are generated by the accelerated particles themselves (Lagage & Cesarsky 1983; Gosling et al. 1981; Ostrowski 1988; Cane et al. 1990; Lee & Ryan 1986; Pelletier et al. 2006; Li et al. 2009). The diffusive shock acceleration (DSA) is so efficient that the back-reaction of the accelerated particles on the shock dynamics cannot be neglected. Thus the theoretical challenge is how to efficiently model the dynamics of a full shock (Jones & Ellison 1991; Caprioli et al. 2010b; Zank et al. 2000; Li et al. 2003; Bykov et al. 2008; Lee 2005). The first step toward non-linear diffusive shock acceleration (NLDSA), in which particle acceleration and shock dynamics are calculated self-consistently, is represented by the two-fluid model proposed by Drüry & Völk (1981). This approach is indeed very useful for studying the dynamics of a shock modified by the presence of accelerated particles but it cannot provide any spectral information about the cosmic rays (CRs) that are produced. There are three basic approaches based on the kinetic study of NLDSA: (1) semi-analytic solutions of the stationary diffusion–convection equation coupled with the gas-dynamic equations (Blasi et al. 2007; Malkov et al. 2000; Bell 2004; Liu et al. 2013), (2) fully numerical simulations, in which a time-dependent diffusion–convection equation for the CR transport is solved in connection with coupled conservation laws describing the gas dynamics (Kang & Jones 2007; Zirakashvili & Aharonian 2010; Feng et al. 2012; Sun et al. 2012), and (3) in the stationary Monte Carlo simulations, where the particle population is calculated by using the particles scattering off the scatter centers frozen into the background magnetic fields (Ellison et al. 1996; Vladimirov et al. 2006; Bykov et al. 2008). Since the velocity distribution of superthermal particles in the Maxwellian tail is not isotropic in the shock frame, the diffusion–convection equation cannot directly follow the injection from the thermal "pool" into the superthermal CR population. Except for the Monte Carlo simulation, the other models require a statistical description of the transport of superthermal particles based on a priori assumption of the transparency function for thermal leakage (Blasi et al. 2005; Kang & Jones 2007; Vainio & Laitinen 2007). However, the dynamical Monte Carlo simulations can be used to model the shock dynamics time-dependently and also to eliminate uncertainty arising from an assumption of injection (Knerr et al. 1996; Wang & Yan 2011, 2012). Similar to the hybrid simulation, since the proton's mass is much larger than the electron's mass, the Monte Carlo model also treats the total plasma fluid as one species of proton fluid with a massless electronic fluid which just balances the electric charge state for maintaining a neutral fluid (Leroy et al. 1982). Owing to there being no distinction between thermal and non-thermal particles, the injection of particles is intrinsically defined by the prescribed scattering properties, and so it is not controlled with a free parameter (Caprioli et al. 2010b).

Wang & Yan (2011) have considered an anisotropic scattering angular distribution based on the previous dynamical Monte Carlo model (Knerr et al. 1996). This anisotropic scattering law includes four different dispersions of Gaussian distribution functions. As a result, a series of similar energy spectra showing small differences with respect to the power-law tail is obtained and the results show that the energy spectral index is affected by the prescribed scattering law. Specifically, the total shock's energy spectral index is less than 1 and shows an increasing function of the dispersion of the scattering angular distribution, but the subshock's energy spectral index is greater than 1 and shows a decreasing function of the dispersion of the scattering angular distribution.

Shen et al. (2007) present an improved method for calculating the Alfvén speed and shock strength near the Sun (>2 R_☉) in which as many observations as possible are used. For example, the SOHO, ACE, Wind, STEREO-A, and STEREO-B spacecrafts are used to study the eruption process of coronal mass ejections (CMEs; Wang 2006; Jin et al. 2009; Feng et al. 2013). They use decameter-hectometric (DH) type II radio bursts instead of a one-dimensional global density model to deduce the ambient plasma density, use LASCO and EIT images to deduce the heliocentric distance and shock speed at the time of DH type II bursts, and use the Wilcox Solar Observatory synoptic charts and CSSS models to extrapolate the coronal magnetic field. The plasma density, height and speed of shocks, and magnetic field strength are all deduced independently. Two events, a relatively slow CME on 2001 September 15 with a CME speed of 752 km s⁻¹and a very fast CME on 2000 June 15 with a CME speed of 1409 km s⁻¹, are chosen to describe the calculation process. The calculation results suggest that the first event had a stronger shock, with a Mach number larger than 3.43 and up to 4.18, but the second event had a relatively weak shock, with a Mach number less than 3.21 and as low as 1.90. In view of the NLDSA, we suggest that there are different degrees of anisotropic scattering that can exist in shocks driven by different speeds of CME flow. Here, the speed of CME flow should represent the upstream bulk flow speed, but it should not be equal to the shock evolutional velocity V_sh in the present simulation as mentioned below. In fact, the shock evolutional speed V_sh should represent the strength of the shock, which would be correlated with the relative difference in velocity (ΔU = U_u − U_d) between the upstream and downstream bulk flow speeds. Thus, it is necessary to further study whether the anisotropy in particle scattering would play an important role in the particle injection rate which proceeds from the thermal downstream to the precursor region upstream in the simulated system. The particle injection rate and acceleration efficiencies depend on many factors. For example, Ellison et al. (1990) have concluded that both the nonlinear shock structure and the mean free path would play a role in the particle injection rate. However, this paper will focus on the extent of anisotropy in particle scattering which would affect the energy injections and energy losses in simulations at a given "free escape boundary" (FEB) and scattering time. The previous simulated results indicated a slightly different energy spectrum in all cases (Wang & Yan 2011). To amplify this difference for further investigation of the energy injection rate in different anisotropic scattering cases, we consider a certain cross-field boundary (CFB) based on the original simulation box. In this scenario, the CFB could be helpful for preventing the particle from dissipating to infinity via cross-field diffusion and make it easy for injected particles from the downstream region to move into the upstream region via parallel-field diffusion. Thus the effect of this CFB could also amplify the difference in the injection rate and acceleration efficiency among the four cases using a different anisotropic scattering angular distribution. In addition, although the solar wind bulk flow is a one-dimensional fluid in this simulated system, the particle's thermal velocity still includes the three-dimensional components v_Lx, v_Ly, and v_Lz. The CFB is actually used to protect the velocity components v_Ly and v_Lz which have not yet been dissipated in the cross-field direction. Once particles reach the CFB, their velocity components V_Ly and V_Lz will be reflected by the CFB, resulting in the negative velocity components −V_Ly and −V_Lz along the corresponding opposite directions. For the sake of studying the particle injection rate, we just add the CFB to the original simulation system which uses a Gaussian scattering angular distributions with four standard deviation values.

In Section 2, the basic simulation method is introduced with respect to the Gaussian scattering angular distributions for obtaining the energy injection and loss functions of time in each case. In Section 3, we present the data analysis for all cases with four types of scattering angle distributions. Section 4 includes a summary and the conclusions.

The Monte Carlo model is a general model, and although it is rather expensive computationally, it is important in many applications because it includes the dynamical effects of the NLDSA in simulations. Since the prescribed scattering law can be used to replace the calculation of the complex electromagnetic field in plasma simulations (Su et al. 2012; Giacalone 2004; Winske & Omidi 1996), we assume that the particles scatter elastically off the background scattering centers with their scattering angles following a Gaussian distribution in the local frame. In this scattering scenario, the assumption of elastic scattering requires that scattering centers are frozen into the background fluid (Ellison & Double 2004); simultaneously, the assumption of a constant scattering time (i.e., collision time) for all particles means that the particle's mean free path is proportional to its local velocity in the local frame with

$$\begin{equation} \lambda = V_{L}\cdot \tau , \end{equation} \tag{ 1 }$$

where τ is the average scattering time. These dependencies of the prescribed scattering law are similar to those in the hybrid simulation (λ∝v^1/2) (Ellison et al. 1993). In order to discuss the effect of the anisotropy of particle scattering on the particle injection rate in a simulation system, we chose a Gaussian distribution to represent the anisotropy in scattering angular distributions. Under the prescribed scattering law, the particle injection rate is closely correlated with those particles from the "thermal pool" in the downstream region that become superthermal particles (Ellison et al. 2005). This dynamical Monte Carlo code can be used to simulate a plane-parallel collisionless shock; in this scenario the large-scale electric and magnetic fields, which correspond to the parallel shock where the effects of fluctuations and waves are contained in the assumed scattering can be neglected.

In these simulations, the entire shock is simulated in a one-dimensional box as shown in Figure 1. The initial continuous inflow enters the box from the left boundary with a supersonic bulk velocity (U₀), and there is a stationary reflective wall at the right boundary of the box acting to form a piston shock moving from right to left. After a certain time, a steady compression region (i.e., a downstream region) will be formed in front of the reflective wall. The bulk velocity in the downstream region will become zero, since the particles dissipate in the downstream region and their large translational energy is converted into isotropic, random energy. To model the finite size of the system and the lack of sufficient scattering far upstream to turn particles around, the simulation includes the escape of energetic particles at an upstream FEB. In the absence of a better description of this phenomenon (i.e., the difficulty from the technical point of view is that it is not clear which particles actually do escape the system), so far, the best way to handle the escape flux is to impose a reasonable location for an FEB from a mathematical point of view (Mitchell et al. 1983; Caprioli et al. 2010a). As shown by the Ellison et al. (1990) model, the bow shock is estimated to moved at a speed in the range of 10–100 km s⁻¹. Given the time of ∼150 s for the shock to happen, one can figure out that the bow shock moved 1500–15,000 km. This distance was ∼40λ₀, and thus λ₀ = 40–400 km. The FEB distance in front of the shock was ∼ −93λ₀, or between 0.5 and 5R_e from the shock. Since the FEB parameterizes the size of the acceleration region, our model assumes that an FEB distance in front of the shock of ∼3R_e is plausible. We run the simulation with one size of FEB, which is equal to 3R_e. However, Knerr et al. (1996) performed two simulations; one has the same parameters as ours with the same size of FEB, and the other has the same parameters as ours with an FEB that is five times larger, 15R_e. If not for computational constraints on the size of the simulation grid, a simulation with an even larger FEB (say, 100 times larger) would be more appropriate for testing the highest energetic particles and the dynamical shock structure. However, the results of Knerr et al. (1996) demonstrate that a subshock structure persists in both of the two simulations with different sizes of FEB. It seems that the average shock flow profile for a simulation with a larger FEB shows a smoother shock structure than a simulation with a smaller FEB.

Zoom In Zoom Out Reset image size

In the dynamic Monte Carlo model, the position of the FEB moves with the shock front at the shock velocity (V_sh) and remains a constant distance in front of the shock position (i.e., FEB = 3R_e). This distance is large enough for a majority of the injected particles to move between the foreshock region and the downstream region. The size of the foreshock is the distance from the shock to the FEB and thus sets a limit on the maximum energy a particle can obtain. Since the injected particles cross the shock and diffuse upstream, they negatively contribute to the bulk velocity, and the bulk velocity becomes smaller from the FEB to the shock position. If the length of the foreshock region remains constant, after enough time has elapsed to create a larger number of accelerated particles, the system will form a steady state with respect to the amount of energy entering and exiting the upstream region. In addition, to amplify the effects of the anisotropy in particle scattering, we also consider a CFB with a size of R_{y, z} ∼ 1.5R_e for preventing the cross-field diffusion to infinity. Once particles reach the CFB, they will move along the reflected directions with negative velocity components (i.e., v_Ly = −v_Ly, v_Lz = −v_Lz). Simultaneously, the CFB can ensure that particles have efficient diffusive processes in the one-dimensional system along the x direction. Thus we are able to compare the differences in the particle injection rate in each case.

The most comprehensive way to study dissipation and wave generation in collisionless shocks is through particle-in-cell (PIC) simulations, where particle motion and field fluctuations are obtained as solutions of the Newton–Lorentz and Maxwell equations. Research on PIC simulations has largely, but not exclusively, focused on perpendicular shocks (Baring 2007; Forslund 1985; Spitkovsky 2008; Niemiec et al. 2008). In the Monte Carlo method, the calculations of the particle motion and field are based on a computational grid with a prescribed scattering law that allows the particle to elastically scatter off the scattering centers frozen into the background fluid. The total box length in this simulation system is X_max = 300, and it is divided into n_x = 600 grids. The initial number of particles in each grid is n₀ = 650. In addition, we use a flux-weighted inflow to ensure that the particles entering into the box with the same density flow upstream with time. This inflow in the "pre-inflow box" (PIB) is entered in the left boundary of the simulation box. The total simulation time is T_max = 2400, and it is divided into the number of time steps N_t = 72000 with a time step dt = 1/30. The distance of the FEB from the shock front is set as X_feb = 90. The radii of the CFB is set as R_{y, z} = 50. These simulation codes consist of three substeps. (1) Individual particles move along the x, y, and z axes with their local velocities in each component being, respectively,

$$\begin{eqnarray} &&x=x_{0}+v_{x}\cdot t \end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray} &&y=y_{0}+v_{y}\cdot t \end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray} &&z=z_{0}+v_{z}\cdot t. \end{eqnarray} \tag{ 4 }$$

Since the magnetic field B₀ is parallel to the simulated shock's normal direction, the bulk fluid quantities describing the bulk fluid vary only in the x direction. (2) Summations of particle masses and velocities are collected on a background computational grid. In this substep, the statistical average bulk speed of each grid represents the velocity of each scattering center. Once the value of the bulk speed drops to zero, the position of the shock front is determined by the displacement of the corresponding grid, and that means that the shock position is moved with an evolutional velocity v_sh away from the stationary reflection wall. Simultaneously, the size of the downstream region is extended dynamically with a constant velocity v_sh. Similarly, a foreshock region or precursor with a bulk velocity gradient is formed by the "back pressure" of the backward diffused particles. The movement of the FEB is also parallel to the shock, moving with the same constant velocity v_sh. (3) An anisotropic scattering law is applied. According to the scattering rate (i.e., R_s = dt/τ, where R_s is the probability of the scattering events in time step dt, and τ is the average scattering time), the fractions of particles are chosen to scatter off the background scattering centers with their corresponding scattering angles obeying the given Gaussian distributions. Once the time step size is set to be less than the scattering time (dt < τ), this algorithm can accurately "resolve" the scattering and produce an exponential mean free path distribution. According to the pitch-angle scattering law in the steady-state Monte Carlo code, the tip of the particle's local-frame momentum vector p will undergo a random walk on the surface of a unit sphere. If the particle originally had a pitch angle θ, and after a time δt undergoes a change in direction of magnitude δθ, its new pitch angle θ' is related to the old one by (Ellison & Double 2004)

$$\begin{equation} \cos \theta ^{\prime } = \cos \theta \cos \delta \theta +\sin \theta \sin \delta \theta \cos \phi , \end{equation} \tag{ 5 }$$

where ϕ is the azimuthal angle of the momentum change dp measured relative to the plane defined by the original momentum p and B. Here ϕ is randomly chosen from a uniform distribution between -π and π, δθ is randomly chosen from a uniform distribution between 0 and δθ_max; the maximum value of δθ is determined by

$$\begin{equation} \delta \theta _{{\rm max}}=\sqrt{6\delta {t}/t_{c}}, \end{equation} \tag{ 6 }$$

where δt is the time between pitch angle scattering, and t_c is the collision time. The tip of the momentum vector randomly walks over the surface of a sphere with radius p = |P|.

In the present dynamical time-dependent Monte Carlo simulation, we define δθ and δϕ as the variations of the pitch angle θ and the location azimuth angle ϕ during the time step dt, respectively. The pitch angle δθ is changed to randomly vary from 0 to π, and the location azimuthal angle δϕ is changed to randomly vary from 0 to 2π in each time step dt, thus the new θ_N and new ϕ_N are changed as follows:

$$\begin{equation} \theta _{N} = \theta + \delta \theta \end{equation} \tag{ 7 }$$

$$\begin{equation} \phi _{N} = \phi + \delta \phi. \end{equation} \tag{ 8 }$$

Here, dt is the time step for the total evolution time. We randomly varied δθ and δϕ with a function f(δθ, δϕ) according to the given Gaussian distribution with four different dispersions of anisotropic scattering. In the present dynamical time-dependent code, δθ and δϕ are averaged over the collision time τ and their corresponding fraction of the collision events, which occurred in each computational grid during each time step dt. When the candidates with their local velocities scatter off the grid-based scattering centers, their corresponding variations of δθ and δϕ are randomly chosen from the given Gaussian function f(δθ, δϕ) with respect to the four different standard deviation values in the four cases. The scattered particles move along their paths until they experience new scattering. During the time step, if all the chosen particles have completed their scattering, the background bulk speed is subsequently changed. In turn, the varying background bulk speed will also change the particle's individual velocity in the local frame in the next time step. The entire simulation time consists of the number of time steps (N_t = 72, 000) involved in the above three substeps.

These simulations are all based on a one-dimensional parallel component of the Earth's bow shock. All the simulated parameters are listed in Table 1. The observed shock data are taken from the Earth's bow shock. An upstream supersonic flow U₀ with an initial Maxwellian thermal velocity V₀ in the local frame and an inflow in a PIB are both moving along the one-dimensional simulation box from the left to the right. The parallel magnetic field B₀ is along the x axis direction. An FEB with a constant length X_feb = 90 is in front of the shock position. The CFB is set as R_{y, z} = 50. The simulation box dynamically consists of three regions: upstream, precursor, and downstream. The bulk fluid speed in the upstream region is U = U₀, the bulk fluid speed in the downstream region is U = 0, and the bulk fluid speed with a gradient of velocity in the precursor region is U₀ > U > 0. The modified flow velocity's spatial profile in the shock deviates from the form of a familiar step-function in the test-particle acceleration scenarios, with the energetic particles pushing against the upstream flow and decelerating it far ahead of the shock discontinuity. Accordingly, an obvious upstream shock precursor is produced with a declining flow velocity profile from the FEB to the shock position (Baring 2007; Berezhko et al. 1996; Bykov & Toptygin 2005). To obtain detailed information on all the particles in the simulation processes at any instant of time, we need a matrix large enough to record the velocities, positions, and the elapsed time of all particles as well as the indices and the bulk speeds of all the grids. Then we can obtain the energy spectra from the downstream, precursor, and upstream regions. The mass of the particles that escape and their momentum and energy losses via the FEB can also be obtained. By calculating the particle injection from the thermal "pool" to the superthermal population in the downstream region and the energy losses via the FEB in the precursor region, we can compare the particle injection rate and energy losses in each case and find how the anisotropy in particle scattering affects the particle injection rate in each case.

To examine the relationships between the particle injection rate and the anisotropy in particle scattering, we perform these Monte Carlo simulations with four different dispersion of the pitch-angle scattering distributions using a new simulation system with a CFB. The simulated cases are presented by a Gaussian function with a standard deviation σ and an average value (i.e., the expected value in mathematical terms) μ = 0 involving four cases: (1) Case A: σ = π/4. (2) Case B: σ = π/2. (3) Case C: σ = π. (4) Case D: isotropic distribution.

3.1. Shock Structures

We present the entire shock evolution with the velocity profiles of the time sequences in each case as shown in Figure 2. The continuous inflow with a supersonic velocity U₀ moves from an initial boundary (X = 0) of the upstream region to the downstream region in the simulation box with time. The total bulk speed profiles consist of three regions with respect to time: the upstream region U = U₀, the precursor region 0 < U < U₀, and the downstream region U = 0. The profiles of the bulk speed are distinct in two positions of the FEB and the shock front with time. From Cases A, B, and C to Case D, the precursor explicitly shows an increasing slope for the bulk speed, and the shock's position X_sh also shows an increasing displacement increment in the x axis at the end of the simulation. This means the shock evolves with an increasing velocity V_sh from Cases A, B, and C to Case D. The simulated results of the dynamical shock in these four cases are listed in Table 2. In addition, by adding a CFB in the new simulation system, we also obtain the differences in the shock front position ΔX_sh compared with the previous simulations (Wang & Yan 2011) with an increasing value of (ΔX_sh)_A = −3.5, (ΔX_sh)_B = +6, (ΔX_sh)_C = +6, and (ΔX_sh)_D = +17.5 from Cases A, B, and C to Case D, respectively. It is obvious that the CFB enhances this difference in the shock structure in the four cases with different dispersion of particle scattering angular distributions.

Zoom In Zoom Out Reset image size

According to the relationship between the upstream and the downstream regions, we are able to calculate the total shock compression ratio r_tot in the shock frame in each case as follows:

$$\begin{equation} r_{{\rm tot}}=\frac{U_{0}+|V_{{\rm sh}}|}{|V_{{\rm sh}}|}. \end{equation} \tag{ 9 }$$

The different shock evolutional velocity V_sh in the different cases will probably lead to a different dynamical shock structure. To show this difference, we present the subtle velocity profiles at the end of the simulation in each case in Figure 3. Evidently, the fluctuation of the velocity between the V_sub and V_d has an obviously increasing value from Cases A, B, and C to Case D. The specific structure in each plot consists of three main parts: the precursor, the subshock, and the downstream. The smooth large-scale precursor is between the FEB and the position of the subshock, X_sub, where the bulk velocity gradually drops from U₀ to V_sub. The sharp short-scale subshock spans three grid-lengths involving a deep drop of the bulk speed abruptly from v_sub to v_d, where the scale of the three grid-lengths is approximately the same as that of the thermal mean free path of the thermalized particles in the downstream region. Thus the subshock's velocity can be defined by the value of V_sub in each case. The velocity V_d represents the downstream bulk speed at the shock position at the end of the simulation. The bottom solid line denotes the backward shock evolutional velocity V_sh with an increasing value from Cases A, B, and C to Case D. For a comparison, in the same shock frame, we can calculate the subshock's compression ratio r_sub as follows:

$$\begin{equation} r_{{\rm sub}}=\frac{V_{{\rm sub}}+|V_{{\rm sh}}|}{|V_{{\rm sh}}|}. \end{equation} \tag{ 10 }$$

Zoom In Zoom Out Reset image size

3.2. Particle Injection Rate

We have monitored the energy of all the particles over time in different regions with respect to all cases. Figure 4 shows all the types of energy functions with time. The term E_tot is the summation of energy for all the particles in the complete simulation system over time. E_box is the energy summation of the actual particles in the simulation box over time. E_pib is the energy summation of the continuous particles entering the simulation box from the PIB over time. E_dow1 is the energy summation of the particles in the downstream region over time. E_dow2 is the energy summation of the particles whose local velocity is larger than the initial bulk velocity U₀ in the downstream region over time. E_dow3 is the summation of the initial kinetic energy of the particles injected from the "thermal pool" into the superthermal populations in the downstream region over time. E_feb is the summation of energy for the particles in the precursor region over the time. E_out is the summation of energy for the particles which escaped from the FEB over time. E_inj is equal to the E_dow2 minus E_dow3 over time. Clearly, the total energy E_tot in the simulation system at any instant in time is not equal to the actual box energy E_box at any instant in time in each plot. As is evident from the real-time functions in Figure 4, the non-linear divergence between the curves for E_box and E_tot is produced with a decreasing value from Cases A, B, and C to Case D. Also, the energy loss function E_out is produced with a decreasing value from Cases A, B, and C to Case D. Simultaneously, the energy injection E_inj shows an increasing value from Cases A, B, and C to Case D.

Zoom In Zoom Out Reset image size

As shown in Table 3, all the listed results of the particle injection and losses in each case are calculated at the end of the simulation. M_loss, P_loss and E_loss are the mass loss, the momentum loss, and the energy loss of the particles that escaped via the FEB, respectively. E_feb, E_inj, E_tot, and E_dow1, with units of initial box energy E₀, are all the energy values in their respective statistical volumes at the end of the simulation. R_inj represents the rate of energy injection E_inj with total downstream energy E_dow1 at the end of the simulation. R_loss represents the rate of energy loss E_loss with total energy in the system E_tot at the end of the simulation. These correlations are presented as follows:

$$\begin{equation} E_{{\rm inj}}=E_{{\rm dow2}}-E_{{\rm dow3}} \end{equation} \tag{ 11 }$$

$$\begin{equation} R_{{\rm inj}}=E_{{\rm inj}}/E_{{\rm tot}} \end{equation} \tag{ 12 }$$

$$\begin{equation} E_{{\rm loss}}=E_{{\rm out}} \end{equation} \tag{ 13 }$$

$$\begin{equation} R_{{\rm loss}}=E_{{\rm out}}/E_{{\rm tot}}. \end{equation} \tag{ 14 }$$

Table 2. The Results of the Shock Simulation

Item	Case A	Case B	Case C	Case D
Xsh	199.5	165.5	146	106.5
XFEB	109.5	75.5	56	16.5
Vsub	0.1075	0.1460	0.1757	0.2525
Vd	+0.0207	−0.0024	+0.0144	+0.0045
Vsh	−0.0419	−0.0560	−0.0642	−0.0806
rtot	8.1642	6.3532	5.6753	4.7209
rsub	2.3337	2.8114	3.1799	4.1328
Γtot	0.7094	0.7802	0.8208	0.9031
Γsub	1.6247	1.3281	1.1881	0.9788
VLmax	11.4115	14.2978	17.2347	21.6285

Download table as:  ASCIITypeset image

Table 3. The Particle Injection and Losses

Item	Case A	Case B	Case C	Case D
Mloss	1037	338	182	38
Ploss	0.0352	0.0189	0.0123	0.0025
Eloss	0.7468	0.5861	0.4397	0.0904
Etot	3.3534	3.4056	3.3574	3.4025
Efeb	0.8393	0.5881	0.5310	0.3397
Edow1	2.1451	2.5612	2.6359	2.6903
Einj	0.1025	0.1912	0.2873	0.3955
Rinj	3.02%	5.62%	8.45%	11.63%
Rloss	22.27%	17.21%	13.10%	2.66%

Notes. The units of mass, momentum, and energy are normalized to the proton mass M_p, initial total momentum P₀, and initial box energy E₀, respectively.

Download table as:  ASCIITypeset image

For a comparison, the mass loss, momentum loss, energy loss, and energy injection functions with time are calculated in Figure 5. Since the simulated systems are based on the computed calculations, the resulting energy losses are inevitable.

Zoom In Zoom Out Reset image size

Figure 5 shows that the mass loss, momentum loss, and the energy loss functions have a decreasing value at any instant of time from Cases A, B, and C to Case D. Among these loss functions, the energy loss function shows decreasing values of (E_loss)_A = 0.7468, (E_loss)_B = 0.5861, (E_loss)_C = 0.4397, and (E_loss)_D = 0.0904 at the end of the simulation from Cases A, B, and C to Case D. By contrast, the energy injection function shows increasing values of (E_inj)_A = 0.1025, (E_inj)_B = 0.1912, (E_inj)_C = 0.2873, and (E_inj)_D = 0.3955 at the end of the simulation from Cases A, B, and C to Case D. The energy injection rate also shows increasing values of (R_inj)_A = 3.02%, (R_inj)_B = 5.62%, (R_inj)_C = 8.45%, and (R_inj)_D = 11.63% at the end of the simulation from Cases A, B, and C to Case D. In the isotropic Case D, the energy injection rate is consistent with the observation that more than 10% of the incoming solar wind energy flux goes into superthermal particles without including escaping particles (Ellison et al. 1990). Because of the existence of energy loss in the simulated system, the shock compression ratios are naturally affected according to the Rankine–Hugoniot (RH) conditions. Therefore, the difference in the energy injection or loss, which is caused by the anisotropy of particle scatters, can directly affect all aspects of the simulated shock including the subtle shock structures, compression ratios, maximum particle energies, energy spectral shapes, etc. It is this self-consistent injection mechanism and computer grid techniques which enable the energy injection and loss functions to be calculated, further simplifying the analysis of energy in the DSA theory.

3.3. Maximum Energy

We select some individual particles from the downstream region at the end of the simulation to obtain the plots in the coordinates of phase, space, and time. The trajectories of the selected particles are shown in Figure 6.

Zoom In Zoom Out Reset image size

Among these selected particles in each case, one of these trajectories clearly shows the full acceleration processes of the maximum energy particle which undergoes multiple crossings with the shock front. The maximum value of the local velocity shows increasing values of (VL_max)_A = 11.4115, (VL_max)_B = 14.2978, (VL_max)_C = 17.2347, and (VL_max)_D = 21.6285 from Cases A, B, and C to Case D in each plot and in Table 3. Consequently, the cutoff energy at the "power-law" tail in the energy spectrum is given with an increasing value of (E_max)_A = 1.23 MeV, (E_max)_B = 1.93 MeV, (E_max)_C = 2.80 MeV and (E_max)_D = 4.41 MeV from Cases A, B, and C to Case D, respectively. For the particles that escaped, because their energies are higher than the cutoff energy, they are not available in the system due to their eventual escape via the FEB. Since the FEB is a constant distance (i.e., X_feb = 90) in front of the shock and maintains parallel movement of the shock front in each case, once an accelerated particle diffuses beyond the position of the FEB, this particle will be removed from the system. Table 3 shows the numbers of particles that escape at the end of the simulation with decreasing mass losses of (M_loss)_A = 1037, (M_loss)_B = 338, (M_loss)_C = 182, and (M_loss)_D = 38 from Cases A, B, and C to Case D. Also, the statistical data describing energy exhibit an energy loss rate with a decreasing value of (R_loss)_A = 22.27%, (R_loss)_B = 17.21%, (R_loss)_C = 13.10%, and (R_loss)_D = 2.66% from Cases A, B, and C to Case D, at the end of the simulation.

Except for the particles with maximum energy, there are also common energetic particles which are shown in Figure 6, some of which have finite energy accelerations from the multiple crossings with the shock and some of which do not have additional energy gains owing to their lack of probability for crossing back into the shock front. Provided that the cutoff energy of the simulated system is not affected by the anisotropy in particle scattering, the maximum particle energy in each case should be identical. Actually, this difference of the maximum particle energy in each case should be contributed by the different energy injection and losses caused by the anisotropy in particle scattering.

3.4. Heating, Acceleration, and Spectrum

As shown in Figure 7, the four energy spectra with power-law tails represent the four cases, averaged over the precursor region, at the end of the simulation. The thin solid curve with a narrow peak is the initial Maxwellian distribution in the shock frame. The four extended energy spectra all consist of two very different parts: a low energy part and a high energy part. The low energy part in the left side of the initial spectrum, which ranges from the low energy to the central peak, shows an irregular fluctuation in each case. The high energy part in the right side of the initial spectrum, which ranges from the central peak energy to the cutoff energy, shows a smooth power-law tail in each case. The irregular fluctuation indicates that the supersonic upstream fluid slows down in the precursor region via the back-reaction of the accelerated particles injected via the subshock. The power-law tail implies that the injected particles from the thermal pool in the downstream region scatter into the precursor region crossing the shock front for multiple energy gains and become superthermal particles.

Zoom In Zoom Out Reset image size

Looking closely at the extended curves in Figure 8, the low energy part in each case has a point that clearly joins with the high energy part in the energy spectra calculated in the precursor region combined with the downstream region for the presented simulations. We can see that the obvious joining points exist in the energy spectra calculated in the regions including the precursor, but do not exist in the pure downstream region in Figure 9. By comparing the spectral shapes in the precursor and downstream regions, we can verify that the injection rate is a vital factor in the reproduction of energetic particles in the standard DSA model. As shown in the extended curves in Figure 9, which are only calculated in the downstream region in the present simulations, the shape of the energy spectrum is not similar to that in the previous simulations performed by Wang & Yan (2011) and Knerr et al. (1996). Although this difference is not so obvious, the differences of the energy spectral indices in the corresponding group of simulations are clearly visible as shown in Figure 12. Compared with Figures 7 and 8, there are no irregular fluctuations in the energy spectra; we guess that the main difference between the energy spectra in the precursor region and the spectra without a precursor region is that the precursor region consists of very cold upstream particles and a minority of the superthermal particles that are injected from the downstream region via a subshock; the downstream region consists of thermal particles, which undergo a thermalization processes, and some of the accelerated particles. Since the thermalization processes hides the disconnection between the cold particles and the superthermal particles in the downstream region, there is no joining point in the energy spectrum in the volume of the downstream region. By constrast, there is no thermalization processes in the precursor region, in which the incoming cold plasma from the upstream region encounters the superthermal particles injected from the downstream region via a subshock, so a distinct joining point is revealed between the cold plasma and superthermal particles in the energy spectrum. Simultaneously, the position of the joining point shows an increasing energy value from Cases A, B, and C to Case D. Consequently, the corresponding cutoff energy at the power-law tail in the precursor region also shows an increasing value from Cases A, B, and C to Case D. Also this position of the joining point should be correlated with the average thermal velocity in the downstream region. As shown in Figure 10, the four thermal velocity functions are averaged over the downstream region with time. Each curve denotes the evolution of the average thermal velocity with time and becomes constant after a certain duration (i.e., t = 500). Eventually, the average thermal velocity V_th shows an increasing value from Cases A, B, and C to Case D, at any instant of time. As expected, the energy injection from the thermal pool in the downstream region shows an increasing value from Cases A, B, and C to Case D. Therefore, as shown in Figure 7, the energy spectrum in the precursor region shows an increasingly hard spectral slope as the dispersion value σ of the Gaussian scattering angular distribution increases. This correlation of the energy spectrum averaged over the precursor region with the prescribed scattering law is consistent with the energy spectrum averaged over the downstream region.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

According to DSA theory, the spectral index can be deduced by the following equations:

$$\begin{equation} dJ/dE\propto E^{-\Gamma } \end{equation} \tag{ 15 }$$

$$\begin{equation} \Gamma _{{\rm tot}} = (r_{{\rm tot}}+2) /[2\times (r_{{\rm tot}}-1)], \end{equation} \tag{ 16 }$$

$$\begin{equation} \Gamma _{{\rm sub}} = (r_{{\rm sub}}+2) /[2\times (r_{{\rm sub}}-1)], \end{equation} \tag{ 17 }$$

where dJ/dE in Equation (15) is the energy flux and the Γ is the spectral index. Equations (16) and (17) represent two types of spectral indices, Γ_tot and Γ_sub, which correspond to the total compression ratio r_tot and the subshock's compression ratio r_sub, respectively. The two types of spectral indices, Γ_tot and Γ_sub, in each case are deduced from the compression ratios and are listed in Table 2. The total shock energy spectral index shows an increasing value of (Γ_tot)_A = 0.7094, (Γ_tot)_B = 0.7802, (Γ_tot)_C = 0.8208, and (Γ_tot)_D = 0.9031 from Cases A, B, and C to Case D. However, the subshock's spectral index shows a decreasing value of (Γ_sub)_A = 1.6247, (Γ_sub)_B = 1.3281, (Γ_sub)_C = 1.1881, and (Γ_sub)_D = 0.9788 from Cases A, B, and C to D. Our simulation results also suggest that the subshock's spectral index value in Case D, being less than 1, may accompany high acceleration efficiencies caused by the effect of the CFB which makes it easier to re-cross the shock front for further energy gains. These simulated results are further compared with the SEP observations in two CME-driven shock events with different speeds of the CME on 2001 September 15 and 2000 June 15. We find that the slow CME drove a strong Mach number shock with an efficient energy injection rate caused by isotropic scattering while the fast CME drove a relatively weak Mach number shock with a low efficient energy injection rate caused by anisotropic scattering.

As shown in Figure 11, the energy injection rate of all cases is denoted by triangles on the solid line, indicating that the energy injection rate is increasing as the dispersion of the scattering angular distribution increases from Cases A, B, and C to Case D. However, the energy loss rate of all cases denoted by circles on the dashed line indicates that the energy loss rate is decreasing as the dispersion of the scattering angular distribution increases from Cases A, B, and C to Case D. This anti-correlation between the energy injection and energy loss rate verified that the fast CME drove a weak Mach number shock and the slow CME drove a strong Mach number shock. From the analysis of the energy injection and energy loss rate, we can suggest that anisotropic scattering exists in the fast CME shock which leads to a relatively large amount of energy escaping via the FEB. Isotropic scattering exists in the slow CME shock which leads to a relatively efficient energy injection from the thermal downstream region to the precursor region. Generally, isotropic scattering ensures that the incoming solar wind energy is conserved at the given FEB position in the system, but anisotropic scattering would not ensure that energy is conserved in a situation with the same FEB position and evolution time of the system. If energy conservation is required for anisotropic scattering, then the scattering process would take a long time and a large FEB in the shock system. Compared with isotropic scattering in the non-relativistic diffusive shock, the case of anisotropic scattering can accelerate the particles with a relatively small energy injection rate and inefficient diffusion, which would exist in the shock driven by a fast CME.

Zoom In Zoom Out Reset image size

Simultaneously, we compared the spectral indices with the previous simulation without a CFB. Although this difference in the spectral slope is not so obvious, the differences in the spectral indices in the corresponding group of simulations are clearly visible in Figure 12. As shown in Figure 12, the empty circles and triangles represent the previous spectral indices without CFB, and the filled circles and triangles represent the present spectral indices with CFB. In the previous simulation, we discussed that the subshock spectral index has decreased as the energy loss decreased from Cases A, B, and C to Case D. Also, the total shock's spectral index has increased as the energy loss decreased from Cases A, B, and C to Case D. Here, we simply discuss whether the differences in the spectral indices depend on the existence of a CFB. Solid circles on the solid line, which denote the present simulations with CFB, compared with empty circles on the dashed line, which denote the previous simulations without CFB, show that the subshock's spectral index has an enhanced deviation higher than the standard value of 1 from anisotropic scattering to isotropic scattering. Simultaneously, filled triangles on the solid line compared with empty triangles on the dashed line also show that the total shock's spectral index has an intensive variation approximating a standard value of 1 from anisotropic to isotropic scattering. By calculating the difference of the spectral index between Case A and Case D (i.e., δΓ = (Γ)_D − (Γ)_A), we find that the subshock spectral indices have a larger difference between Case A and Case D with a value of δΓ_sub = 0.6459 in the present simulation for the case of adding CFB than the value of δΓ_sub = 0.4633 in the previous simulation without adding CFB; the total shock spectral indices also show a larger difference between Case A and Case D with a value of δΓ_tot = −0.1938 in the present simulation instead of the value of δΓ_tot = −0.1500 in the previous simulation. The enhanced deviation of the spectral indices indicates that the cross-field diffusion boundary could effectively prevent the particles from being dissipated in the perpendicular direction of the magnetic field. The particles that escape could have escaped only from upstream via the FEB. Since these particles can only escape from the FEB, a shock driven by a fast CME could lead to greater energy loss from anisotropic scattering with a lower injection rate; a shock driven by a slow CME could lead to lower energy loss from isotropic scattering with a large injection rate. This significant difference in the spectral index verified the observation of a radio II burst and EIT images from the SOHO spacecraft (i.e., the anti-correlation of the strength of the shock with the speed of the CME). According to DSA theory, if the energy loss is limited to a minimum, the simulation models based on the computer will more closely fit the realistic physical situation. All of the values of the total shock's spectral index are less than 1 (i.e., Γ_tot < 1), and the solid lines denote the total shock's spectral index with an increasing value from Cases A, B, and C to Case D as the dispersion of the scattering angular distribution increases. As seen from Cases A, B, and C to Case D, both the subshock's spectral index and the total shock's spectral index approximate the standard value of 1 (i.e., Γ ∼ 1) as the dispersion of the scattering angular distribution increases. As predicted, the RH jump conditions allow us to derive the relation of the compression ratio with the Mach number as r = (γ_a + 1)/(γ_a − 1 + 2/M²). For a non-relativistic shock with an adiabatic index γ_a = 5/3, if the Mach number M ≫ 1, then the maximum compression ratio should be 4. According to the RH conditions, the total shock compression ratio should be less than the standard value of 4, and the total shock's corresponding spectral index should be less than the standard value of 1 for a non-relativistic shock (Pelletier 2001). Simultaneously, we can see that if the dispersion of the scattering angular distribution increases from anisotropy to isotropy, the subshock's spectral index will closely approximate or even exceed the standard value of 1. This study presents the relationship between the energy injection or energy loss and anisotropy in particle scattering to verify the observation of an anti-correlation between the strength of a CME-driven shock and the speed of a CME. The present study also suggests that different degrees of anisotropic scatter will occur in shocks driven by different speeds of CMEs. These relationships will also be very helpful to improving simulation models by using the best choice of the prescribed scattering law. By using the prescribed scattering law instead of assuming a transparent function in the thermal leakage mechanism, a dynamical Monte Carlo model can calculate the energy injection and loss rate in a fully self-consistent and time-dependent manner.

Zoom In Zoom Out Reset image size

In summary, we performed a dynamical Monte Carlo simulation using Gaussian scattering angular distributions based on the Matlab platform by monitoring the particles' mass, momentum, and energy at any instant in time. The specific energy injection and loss functions with time are presented. We examined the correlation of the prescribed Gaussian scattering angular distributions with the energy injection and loss functions. Simultaneously, this dependence is further enhanced by considering the CFB in the simulated system. The role of the CFB demonstrates that the isotropic scattering law produces a strong subshock and the anisotropic scattering law produces a weak subshock. This behavior is consistent with observations of a fast CME driving a low Mach number shock due to a large injection rate and a slow CME driving a high Mach number shock due to a relatively small injection rate. Thus one would think that this visible difference in the spectral indices between the two groups of simulations would be produced by the factor of injection rate in the different CFB conditions. As is well known, the DSA is basically an efficient acceleration mechanism for solar CRs, which would probably dominate in the SEP production. This acceleration mechanism is also expressed by the standard model and first-order Fermi acceleration mechanism. We apply this model in the parallel part of the Earth's bow shock to test the acceleration efficiency of SEP production. For a total shock, its global structure is complicated, and this complexity is related to the obliqueness, different parameters, strength of the magnetic field, etc. Here, besides focusing on the parallel shock, in which the shock normal is parallel to the magnetic field direction, we will also discuss other structures of the Earth's bow shock according to the solar wind magnetic field direction with the shock normal direction. We apply these simulation results of the Earth's bow shock to discuss how the ejection rate can affect the acceleration efficiency in similar CME-driven shocks. Combined with observations of the CME-driven shock, we can explain the correlation of the shock strength with the DSA efficiency driven by different speeds of the CME flow in the standard DSA model. As far as the large-scale structure of the Earth's bow shock or CME-driven shock is concerned, we could discuss the efficiency of the shock acceleration by dividing the global shock model into three components: the perpendicular shock model, the oblique shock model, and the parallel shock model. Thus we can see Earth's total bow shock consists of three parts according to the different angles between the solar wind magnetic field direction and the shock normal direction: the parallel shock, the perpendicular shock, and the oblique shock. Combined with observations from spacecraft, we find that most of the energetic ions seem to be easy to detect from the parallel shock regions of the bow shock because the particles obtain enough additional energy from the Fermi acceleration mechanism through multiple crossings in the shock front. However, the particles in the perpendicular shock experience a crossing only once in the shock front just like the particles running in the induced electric field with relatively less efficient acceleration. The other shock model is the oblique shock model, which combines the parallel shock and the perpendicular shock. In this model, the particle acceleration efficiency depends on which component is dominant in terms of the tilt of the solar wind magnetic field direction versus the shock normal direction. For the CME-driven shock, there are similar plasma behaviors for the solar wind interacting with the solid mass during shock formation.

In conclusion, the relationship between the energy injection or energy loss and anisotropy in particle scattering confirms that the particle injection rate dominates viscous subshock modification and particle acceleration efficiency. As expected, the maximum energy of accelerated particles is closely correlated with the particle injection rate from the thermal pool to the superthermal population. Therefore, we find that the particle injection rate increases as the standard deviation value of the scattering angular distribution increases. This agrees with the observations of SEP fluxes of the shock driven by different CME speeds. In these multiple scattering angular distribution scenarios, the prescribed scattering law affects the energy injection and energy loss. This dynamical Monte Carlo model uses a self-consistent energy injection mechanism instead of assuming a thermal leakage function, which is used in other numerical models. Consequently, we ensure that there is one possibility where a fast CME flow could lead to an anisotropic scattering distribution in the corresponding shock driven by a fast CME, which is equal to anisotropic pitch-angle scattering in our simulations, and then leads to further energy loss via the FEB and weak Mach numbers for the shock. By contrast, another possibility exists where a slow CME flow can lead to an isotropic scattering distribution in the corresponding shock driven by a slow CME, which is equal to isotropic pitch-angle scattering in our simulations, which causes less energy loss via the FEB and strong Mach numbers for the shock. In the CME-driven shock, observations from spacecraft also show that the parallel shock region has a privileged acceleration efficiency on the global shock model. In essence, the efficiency of the particle acceleration in all kinds of complicated shock structures will be dependent on the particle ejection rate in the shock acceleration model. The particle ejection rate decreases with shock model from the parallel shock to the oblique shock to the perpendicular shock model according to the tilt of the magnetic field versus the shock normal direction. Thus the particle injection rate in DSA is a vital factor for the production of SEP in the interplanetary shock.

This work is funded in part by a KLAS of CAS grant Y300101001 and the National Basic Research Program of MOST (grant Nos. 2011CB811401 and 2009CB824800). In addition, this work is also supported by the National Natural Science Foundation of China through grant Nos. 11173041 and 11173042. We thank Prof. Hongbo Hu, who is working on cosmic-ray physics at the Institute of High Energy Physics, CAS, for providing much help during this work. Simultaneously, we thank the referee for raising a lot of helpful points.