1 - Finnish Meteorological Institute, PO Box 503, 00101 Helsinki, Finland
2 - Department of Physical Sciences, PO Box 64, 00014 University of Helsinki, Finland

Abstract
We consider particle acceleration at the scatter-free limit in quasi-planar shock waves propagating in inhomogeneous magnetic fields. It is shown that both non-constant magnetic field intensity and field-line curvature may lead to efficient acceleration of particles at shocks propagating through the structure. Shocks propagating towards increasing magnetic field intensity trap energetic particles, and as the field increases at the shock front the particles, by conserving their magnetic moment $\mu_B=E_\perp/B$ , increase their perpendicular energy by the ratio of maximum field magnitude to the field magnitude at the point of injection, $E_{\perp,\rm max} = E_{\perp, \rm inj} B_{\rm max}/B_{\rm inj}$ . This may result in energy gains by factor of 100 in the solar corona. In addition, shocks propagating in curved magnetic fields may trap particles and accelerate them to high energies on field lines on which the shock-normal angle gradually increases toward 90 $^\circ$ . Suitable field-line geometries should be common in many astrophysical objects, such as stellar coronae and quasi-perpendicular parts of supernova shocks.

Key words: acceleration of particles - ISM: cosmic rays - shock waves - Sun: particle emission

Acceleration of ions in coronal mass ejection (CME) driven shock waves is the most commonly accepted and best developed theory of the genesis of gradual solar energetic particle (SEP) events (e.g., Reames 1999; Lee 1997). The theory is based on the diffusive shock acceleration mechanism (Krymsky 1977; Axford et al. 1977; Blandford & Ostriker 1978; Bell 1978). This mechanism requires turbulent conditions ahead of the shock in order to be rapid enough to account for the observed particle energy spectra and intensity-time profiles. In principle, the required turbulence can be generated by the accelerated particles themselves (Lee 2005) in the largest SEP events, in which the SEP intensities are streaming limited (Reames & Ng 1998). As shown by Vainio (2003), however, the wave-generation by the accelerated particles in smaller SEP events is not expected to be important. Likewise, at relativistic energies the particle fluxes seem too small to be able to generate waves in reasonable time scales to account for their diffusive acceleration at coronal shocks (Vainio 2003).

As noted by Vainio & Laitinen (2001), fluctuating magnetic fields required by coronal particle acceleration might result also from plasma motions at the solar surface: high-frequency Alfvén waves propagating upwards from the base of the corona, responsible for coronal heating on open field lines, will also scatter energetic particles. Small gradual SEP events, with maximum proton energies in the range of a few tens of MeVs, have been shown to be consistent with diffusive shock acceleration in such coronal turbulence (Vainio & Laitinen 2001; Vainio & Khan 2004). According to SEP observations, however, the ion mean free path in the solar wind is far too large to be consistent with the values predicted by propagating such coronal Alfvén waves to the solar wind using WKB theory (Vainio & Laitinen 2001). But the highly turbulent magnetic fields may be limited to a narrow layer around the Sun, if an active turbulent cascade operates in the solar wind to dissipate the magnetic flucutations (Vainio et al. 2003). Detailed numerical modeling of diffusive shock acceleration in such a turbulent coronal layer was performed by Kocharov et al. (2005). They confirmed that in order to accelerate particles to the highest energies the mean free path of the accelerated particles has to be very short. It is unclear, however, if enough turbulence at high frequencies is available in the corona to provide the required efficient scattering - the ambient mean free paths could be even orders of magnitude larger than the values required for the diffusive acceleration mechanism. In this case, particle acceleration during the initial phases of the eruption would have to be essentially scatter free.

Theories for scatter-free shock acceleration exist as well, but they are, with the exception of the shock surfing mechanism (Zank et al. 1996; Lee et al. 1996), restricted to rather small energy gain factors, (see, e.g., Webb et al. 1983). The main idea of diffusive shock acceleration is to enable particles to interact with the shock many times to overcome this problem. In this paper, we investigate the effects of inhomogeneous upstream magnetic fields on scatter-free acceleration at shocks. Our study indicates that under quite plausible circumstances large-scale magnetic field inhomogeneities may lead to multiple shock encounters and, thus, to efficient particle acceleration as well. We will consider two basic types of large scale inhomogeneities ahead of planar shocks: (i) gradients of field intensity and (ii) curvature of field lines (Fig. 1). Although we will treat them separately in our initial analysis, a realistic magnetic environment is expected to contain both types of inhomogeneities simultaneously. This will be addressed by numerical modeling.

$\begin{figure} \par\includegraphics[width=8.1cm,clip]{4754fi01.eps}\end{figure}$	Figure 1: Two basic types of inhomogeneities on upstream magnetic field lines (solid curves): non-constant field magnitude ( upper panel, left) and field-line curvature ( lower panel, left). A fast-mode shock (dashed line) is propagating upwards at speed $V_{\rm s}$ . The figures on the right show the changes in field magnitudes in the two cases.
Open with DEXTER

In the following section we first examine the theoretical basis for the particle acceleration in inhomogeneous magnetic fields. Implications and applicability of the theory are then examined by numerical test particle simulations in Sects. 3 and 4. Details of the simulations are given in Appendix A.

2 Shock acceleration in inhomogeneous fields

2.1 Effect of field-intensity gradient

Consider a magnetic field with a non-zero component of the field-magnitude gradient along the lines of force and a planar fast-mode shock wave propagating at constant speed $V_{\rm s}$ relative to the upstream plasma. Assume that the compression ratio, $r=\rho_2/\rho_1$ , of the shock is constant. Consider particles on a straight upstream field line making a constant angle $\theta _{\rm n}$ with the shock normal and assume that the shock is propagating toward increasing magnetic field intensity. In the local de Hoffmann-Teller frame (HTF), the upstream plasma flows toward the shock along the field line at speed $V_{\rm s}/\! \cos\theta_{\rm n}$ while the magnetic field intensity increases as a function of time. The field intensity at the shock is compressed,

$\begin{displaymath}% r_B=\frac{B_2}{B_1} = \sqrt{\sin^2\theta_{\rm n} r^2\frac{(M^2-1)^2}{(M^2-r)^2}+\cos^2\theta_{\rm n}}, \end{displaymath}$

(1)

$\displaystyle % \frac{{\rm d}E}{{\rm d}t}$	=	$\displaystyle \mu_B\frac{\partial B}{\partial t}$	(4)
$\displaystyle m\frac{{\rm d}v_\parallel}{{\rm d}t}$	=	$\displaystyle -\mu_B\frac{\partial B}{\partial z}$	(5)
$\displaystyle \frac {{\rm d}z}{{\rm d}t}$	=	$\displaystyle v_\parallel$	(6)

$\begin{displaymath}% \frac{{\rm d}N}{{\rm d}E}=\int_0^t {\rm d}t_0 \int_0^\infty... ...E_0\frac{B(t)}{B(t_0)}\right)\frac{{\rm d}\dot N}{{\rm d}E_0}, \end{displaymath}$

(8)

$\displaystyle % \frac{{\rm d}N}{{\rm d}E}$	=	$\displaystyle \int_0^t {\rm d}t_0 \:\delta\left(E-E_0\frac{B(t)}{B(t_0)}\right)\dot N_0$
	=	$\displaystyle \frac{B_i}{E}\frac{\dot N_0(t_i(E))}{\vert\dot B(t_i(E))\vert} =\... ...{\rm d}N_0/{\rm d}z\vert _{z_i(E)}}{\vert\partial B/\partial z\vert _{z_i(E)}},$	(9)

$\begin{displaymath}% \frac{{\rm d}N}{{\rm d}E}=\frac {L(z_i(E/E_0))}{E}\frac{{\rm d}N_0}{{\rm d}z} \end{displaymath}$

(11)

2.2 Effect of field-line curvature

Let us next consider the effect of upstream field line curvature on particle acceleration at shocks. We will again consider a planar shock in the local HTF of a single upstream field line. This time, however, we take the upstream field magnitude to be constant and $\theta _{\rm n}$ to be time-dependent. Thus, the HTF is non-inertial, and we need to consider the effects of the corresponding inertial force in the equations of motion. We will neglect the drift effects caused by the perpendicular component of the inertial force. The inertial drift can be shown to be small

$\begin{displaymath}% \frac{V_{\rm HT,drift}}{V_{\rm HT}} = \frac{\partial \theta... ...rm n}/ \partial t}{\omega _{{\rm gyro}} \sin \theta _{\rm n}}, \end{displaymath}$

(13)

$\begin{displaymath}% \frac{V_{\rm HT,drift}}{V_{\rm HT}} \approx \frac{1}{23 ~ 900} \frac{1}{\sqrt{ (x/R)^2 - (x/R)^4} }, \end{displaymath}$

(14)

In the upstream rest frame, the projected speed of the shock along the magnetic field is $V'_{\rm s} =V_{\rm s}/\! \cos\theta_{\rm n}(t)$ . Thus, the equations of motion, including the parallel component of the inertial force in the local HTF, are

		$\displaystyle \frac{{\rm d}E}{{\rm d}t} = mv_\parallel \frac{{\rm d}V'_{\rm s}}{{\rm d}t}$	(15)
		$\displaystyle \frac{{\rm d}v_\parallel}{{\rm d}t} = \frac{{\rm d}V'_{\rm s}}{{\rm d}t}-\mu_B\frac{\partial B}{\partial z}$	(16)
		$\displaystyle \frac {{\rm d}z}{{\rm d}t} = v_\parallel.$	(17)

We note, finally, that if the shock is propagating into a closed loop, i.e., if both ends of an upstream field line are connected to the shock, particles performing bounce motion between the two ends of the field line will be Fermi-accelerated inside the collapsing magnetic trap by the conservation of the longitudinal adiabatic invariant. However, as this mechanism will give the particles more and more parallel momentum without increasing the perpendicular, it will also be limited to rather small energy gains, because the particles will soon find themselves in the loss-cones of the magnetic mirrors provided by the shock. This mechanism is, therefore, much less efficient than the two other mechanisms described above.

$\begin{figure} \par\includegraphics[width=8.5cm,clip]{4754fi02.eps}\end{figure}$	Figure 2: Setup for line current runs. At t=0 the shock plane is situated at x=0 (dashed line) and moves from right to left with velocity $V_{\rm s}$ . 10 keV protons are injected to the field line marked with a diamond. Magnitude of $\vec{B}$ is shown with the contour lines, along with two field lines.
Open with DEXTER

3 Numerical simulations

We have numerically simulated shock acceleration in three magnetic field configurations where the two scatter-free acceleration mechanisms can be studied independently as well as operating together. The test particle method and approximations we employ are similar to the ones used by Decker (1988) and more recently by Giacalone (2005a). We model the shock as a planar discontinuity propagating with a constant speed $V_{\rm s}$ and gas compression ratio r in a given upstream magnetic field. The electromagnetic magnetic fields behind the shock are calculated analytically under the approximation of infinite Mach number of the shock. Particles are treated as test particles and their trajectories are computed using a Boris-Buneman leapfrog integrator (see Birdsall & Langdon 1985, Ch. 4), modified to reduce the timestep when a particle is interacting with the shock front. The code has been tested against analytical test-particle calculations at shocks by Webb et al. (1983) and by replacing the solver with an embedded Runge-Kutta Prince-Dormand 8-9 method. Details of the numerical model are given in Appendix A.

In all setups presented we have used three shock speeds, $V_{\rm s} \in [400,700,1000]$ km s^-1, and two compression ratios, $r \in [3,4]$ , as shock parameters. Protons are injected in front of the shock so that their guiding centers reside in the given positions, $r_{\rm inj} = r_{\rm gc}(t=0)$ . The exact positions and velocities are then calculated using randomized gyro phase angle and isotropic distribution in pitch (= $\cos{\alpha }$ , where $\alpha$ is the pitch angle), all calculated in the shock normal incidence frame (SNIF). The isotropic pitch distribution should not be taken as an assumption of upstream plasma distribution but rather as an approach "put in all kinds of particles and see what happens''.

The results of the simulations with different shock parameters and magnetic configurations can be classified in many ways. We have found it useful to use two values, namely the number of high energy particles (HEP) produced and the maximum energy reached, which can be always calculated. A particle is considered to be a HEP when its final energy $E_{\rm final} \geq 100$ keV. We give the number of HEPs as percentage of the total number of injected particles. The maximum energy reached $E_{\rm max}$ is calculated as the average of the energies of 25 highest energy particles observed. Both values are computed separately for particles escaping to upstream and downstream regions.

The reader is cautioned that we define the pitch angle of the particle as the angle between the exact velocity vector of the particle and the magnetic field, $\cos{\alpha } = \vec{v} \cdot \vec{b} / v$ . This definition means that $0 \leq \alpha \leq 90^{\circ}$ for particles moving parallel and $90^{\circ} \leq \alpha \leq 180^{\circ}$ for particles moving antiparallel relative to the magnetic field. When analyzing the figures presented in the following sections, the direction of $\vec{B}$ should be checked in each case from the figure showing the corresponding setup.

$\begin{figure} \par\includegraphics[width=8.5cm,clip]{4754fi03.eps}\end{figure}$	Figure 3: Initial pitch angle vs. energy gain ratio $\Delta E / E_{0}$ for a line current run with $V_{\rm s}=400$ km s^-1 and r=4. Each dot is one simulated particle.
Open with DEXTER

3.1 Magnetic field of a line current

In the first set of runs we consider the field of a line current, which is a simple example of magnetic configuration where only the effects due to the curvature of field lines are important for acceleration. Initially at t=0 the shock is situated so that the current $\vec{I}$ resides at the shock plane. Protons with initial energy E₀=10 keV are injected to position $\vec{r}_{\rm inj} =(-10\:{\rm km},R_{\rm curv},0)$ , i.e. 10 km in front of the shock (see Fig. 2). The magnitude of the current I was chosen so that the radius of curvature of $\vec{B}$ is much larger than the radius of gyration, $R_{\rm curv} \gg r_{\rm gyro}$ , even for the resulting HEPs. The curvature radii we have used are $R_{\rm curv} \in [2500,5000,10~000]$ km. In each case B = 0.1 mT at the injection field line. Particles were removed from the simulation when their distance to the shock front $d \geq 4000$ km in the downstream region.

The shock compresses the downstream plasma by the factor r, which means, in the approximation that the field lines are frozen into the plasma, that an initially circular field line turns into an ellipse with major axis $R_{\rm curv}$ and minor axis $R_{\rm curv} / r$ . The tangential component of $\vec{B}$ will also be enhanced by a factor of r, as indicated by the shock jump conditions (A.1-A.4). This in turn means that there will be a second maximum of B in the downstream region (visible in Fig. 2), in addition to the maximum at the shock front. Particles can get trapped between these maxima.

Figure 3, which shows initial pitch angle vs. energy gain ratio $\Delta E / E_{0}$ for the case with $V_{\rm s}=400$ km s^-1 and r=4, is an example of what is typically observed in the line current runs. As indicated by the figure the accelerated protons can be grouped to two classes, namely to Fermi-accelerated and HT-resonant particles. The former propagate initially away from the shock and gain energy by bouncing between the two approaching mirror points located at the points where the shock plane intersects the injection field line. The latter is a population with initial pitch angle $85^{\circ} \la \alpha \la 90^{\circ}$ , propagating slowly away from the shock. However, since the field line curves, the shock catches up with these particles and gives a (small) boost to their parallel speed. The abovementioned procedure is then repeated for a few times until the shock turns completely perpendicular and passes over the injection field line.

Table 1: Results for line current runs as functions of the curvature radius of the field line used, shock speed $V_{\rm s}$ and compression ratio r. Shown are the relative amount of HEP escaping to the downstream region (ds), maximum energy reached $E_{\rm max}$ (in $\log{}_{10}({\rm eV})$ ) and the spectral index $\sigma$ with error limits, ${\rm d}N/{\rm d}E \propto E^{\sigma }$ . In this setup there are no particles escaping to the upstream region.

Figure 3 shows actually two Fermi-accelerated populations. The extra population with $115^{\circ} \la \alpha _{0} \la 140^{\circ}$ has reflected from the previously mentioned downstream maximum of B and moved back to upstream, where these particles have then been Fermi-accelerated as usual. Particles with $90^{\circ} \la \alpha _{0} \la 115^{\circ}$ did not have enough parallel speed to catch the shock again early enough - this population got trapped between the shock and the downstream maximum. However, the formation of this second Fermi-accelerated population is sensitive to initial conditions (especially to the injection position) and shock parameters and is not observed in other runs. A summary of the results from the line current runs is shown in Table 1 as a function of shock parameters.

3.2 Dipole magnetic field

3.2.1 Polar region runs

In this setup we consider the effects of magnetic field gradients alone. The shock propagates towards the positive pole of a dipolar magnetic field. Protons with E₀=10 keV are injected to the straight field line emerging from the pole, which is at an angle $\theta _{\rm n}$ relative to the shock normal (see Fig. 4). In these runs we have used five different oblique angles $\theta _{\rm n} \in [30^{\circ},45^{\circ},60^{\circ},70^{\circ},80^{\circ}]$ with same shock parameters as in the line current runs. The injection positions were chosen so that with each $\theta _{\rm n}$ used the protons' distance to the dipole origin is 8000 km and 300 km to the shock front. In these runs particles were removed if they were $d_{1} \geq 1500$ km to the downstream of the shock or when their distance to the dipole origin $d_{2} \leq 1200$ km. The main purpose of the latter condition is to prevent particles from entering too strong magnetic field, where the particle mover starts to fail due to the value of $\Delta t$ used in the computations. However, as all magnetic flux tubes found in nature have a finite size, this artificial boundary can also be understood to represent, e.g., the solar surface. In the following results we take the particles hitting this surface to be escaping to the upstream region.

$\begin{figure} \par\includegraphics[width=8.5cm,clip]{4754fi04.eps}\end{figure}$	Figure 4: Setup for runs where the shock (dashed line) moves through a flux tube connected to the positive pole of a dipole field. Protons are injected to the "straight'' field line marked with a diamond, which is at an angle $\theta _{\rm n}$ relative to the shock normal. Magnitude of $\vec{B}$ is shown with the contour lines along with some field lines.
Open with DEXTER

In general terms, a certain population of injected particles gets trapped between the shock and the upstream maximum of $\vec{B}$ near the surface of the dipole. In the SNIF, these particles gain parallel energy when they reflect from these maxima. At the same time the injection flux tube contracts and trapped particles are pushed towards the surface, where $\nabla B$ is stronger and $\alpha \rightarrow 90^{\circ}$ . All this then means that the gradient dominated acceleration can be viewed as a form of adiabatic compression, and that the final energy of the particles should scale as $E_{\rm final} \sim E_{0} ~ (B_{\rm final} / B_{\rm initial})$ . In this setup $E_{\rm final} \sim 10^{6.47}$ eV in ideal case, since $B \sim 1/r^{3}$ above the poles of a dipolar field. The corresponding energy gain ratio is $\Delta E / E_{0} \sim 296$ .

$\begin{figure} \par\includegraphics[width=8.5cm,clip]{4754fi05.eps}\end{figure}$	Figure 5: Initial pitch angle vs. energy gain ratio $\Delta E / E_{0}$ for a polar region run with $V_{\rm s}=400$ km s^-1, $\theta _{\rm n} = 45^{\circ }$ and r=4.
Open with DEXTER

However, this picture is somewhat changed by the fact that the gradients affect the particles' trajectories during their encounter with the shock. The gradient force $- \mu \nabla B$ tends to push particles towards the downstream region while efficient acceleration requires that the particles must be able to run away from the shock to upstream region. If the shock is too oblique, the low energy particles are forced to the downstream region, encountering the shock only once. The end result resembles shock drift acceleration (SDA) in homogenous field and is not effective at producing HEPs.

Table 2: Results for the polar region runs as functions of the shock velocity, compression ratio and oblique angle. Shown are the no. HEP escaping to upstream (us) and downstream (ds) regions and maximum energies reached (in $\log{}_{10}({\rm eV})$ ), respectively. Omitted results show no acceleration.

Figure 5 is a typical example of what is observed in the polar region runs. Since there is a mirror point at the shock at all times, a significant portion of injected particles get reflected and accelerated. The figure also shows some horizontal stripes, which are due to the bouncing motion of particles. As in the line current runs, the particles whose parallel speed is initially close to zero (in the SNIF) are accelerated most. This is because these particles are trapped near the shock front and are gradually pushed towards the dipole surface, encountering the shock very many times.

Summary of the results in the polar region runs is given in Table 2. As mentioned before, effective acceleration stops when the projected shock speed $V_{\rm s}'$ is large enough. Fastest shocks with $V_{\rm s} = 1000$ km s^-1 can accelerate particles only up to oblique angle $\theta _{\rm n} \sim 60^{\circ}$ , while the slower ones get to $\theta _{\rm n} \sim 70 ^{\circ}$ . Also a considerable portion of HEPs produced are pushed to the dipole surface when $\theta _{\rm n} \geq 45^{\circ}$ with most shock parameters. These particles also have somewhat higher energies than the ones escaping to downstream region.

3.2.2 Rotated dipole

This final setup is an example of the general situation where the effects of both curvature and gradients are important. The dipole has been rotated $135^{\circ }$ relative to the shock normal. E₀=10 keV protons are injected to position $\vec{r}_{\rm inj} = (-100,10~000,0)$ km, which is 100 km in front of the shock plane at t=0. The injection field line has been chosen so that the shock will eventually turn completely perpendicular (see Fig. 6). We have used the same shock parameters and similar conditions for removing the particles as in the previous setups.

$\begin{figure} \par\includegraphics[width=8.5cm,clip]{4754fi06.eps}\end{figure}$	Figure 6: Setup for runs where the dipole has been rotated $135^{\circ }$ relative to the shock normal. The shock (dashed line) moves from right to left. Protons are injected to a field line (marked with a diamond) having both non-zero gradient and curvature. Magnitude of $\vec{B}$ is shown with contour lines along with some field lines.
Open with DEXTER

In these runs we observe that in the initial phase, when the shock is quasi-parallel, trapped particles bounce between the dipole and the shock. During this time interval the acceleration is mainly gradient-driven adiabatic compression, which is effective to $\theta _{\rm n} \sim 70 ^{\circ}$ depending on the shock velocity and magnitude of $\nabla B$ along the field line. The curvature dominated acceleration starts to operate when the shock is quasi-perpendicular, roughly $60^{\circ} \la \theta _{\rm n} \la 90^{\circ}$ , during which time the acceleration is very rapid (see the line current runs). It should be noted here that, as in the polar region runs, the gradients of the magnetic field intensity start to prevent acceleration when the shock is nearly perpendicular, reducing the effectiveness of HT-resonance.

$\begin{figure} \par\includegraphics[width=8.5cm,clip]{4754fi07.eps}\end{figure}$	Figure 7: Initial pitch angle vs. energy gain ratio $\Delta E / E_{0}$ for a rotated dipole run with $V_{\rm s}=700$ km s^-1 and r=4.
Open with DEXTER

The Fig. 7 is a typical example of results obtained from rotated dipole runs. Most of the particles moving initially towards the shock simply transmit to the downstream region, without gaining any energy. This is because the shock is not a very good mirror in the early phases, since $\vec{B}$ is almost parallel to the shock normal. Again the particles with initial $v_{\vert\vert} \sim 0$ have been trapped for the longest time and gained the highest energies. Also visible are also some stripes due to the finite number of bounces, especially for the population with $120^{\circ} \la \alpha _{0} \la 170^{\circ}$ .

Table 3: Results for the rotated dipole runs as functions of the shock velocity, compression ratio and oblique angle. Shown are the no. HEP escaping to downstream (ds) region and maximum energy reached (in $\log{}_{10}(eV)$ ). In this setup no particles were removed in the upstream region.

Like in the line current runs, the shock traps the particles until it turns nearly perpendicular. The HEPs are then rapidly released as a high-energy beam. Table 3 shows the results for the rotated dipole runs. The maximum energies are also somewhat higher than in the line current setup, because the gradients preaccelerate the particles before they end up to the quasi-perpendicular region. The number of HEPs produced is also considerably higher. However, in contrast to the polar region runs, there were no HEPs escaping to upstream region with the same dipole surface distance d₂=1200 km.

4 Discussion

4.1 Time evolution

In the curvature dominated situations, there is a clear division to accelerated and non-accelerated particles. The former population consists of particles trapped to a region just upstream of the shock with parallel speed $v_{\vert\vert} \sim 0$ in the SNIF at all times. These particles have a good chance of encountering the shock multiple times if the field line they are following shows negative curvature, i.e. the oblique angle $\theta _{\rm n}$ increases with time. In the region of space where the shock is quasi-parallel, the resulting acceleration is modest at best. For shocks with $\theta _{\rm n} \sim 60^{\circ}$ the energy gain ratio for a single encounter is in the range $0 \leq \Delta E / E_{0} \leq 0.15$ for the reflected particles and reaches its maximum $\Delta E / E_{0} \sim 9$ for nonrelativistic particles just before the shock turns completely perpendicular (see the treatment by Webb et al. 1983).

$\begin{figure} \par\includegraphics[width=8.5cm,clip]{4754fi08.eps}\end{figure}$	Figure 8: Energy and position of a HT-resonant proton with initial pitch angle $\alpha = 90^{\circ }$ in the line current run $V_{\rm s}=400$ km s^-1, r=4, $R_{\rm curv} = 10~000$ km, all values are measured in the SNIF. Energy is plotted vs. x-position using the secondary axes ( right). Arrows mark initial values. See Figs. 2 and 3 for reference.
Open with DEXTER

$\begin{figure} \par\includegraphics[width=8.5cm,clip]{4754fi09.eps}\end{figure}$	Figure 9: Energy and position of a Fermi-accelerated proton with initial pitch angle $\alpha = 45^{\circ }$ in the line current run $V_{\rm s}=400$ km s^-1, r=4, $R_{\rm curv} = 10~000$ km, all values are measured in the SNIF. Energy is plotted vs. x-position using the secondary axes (right). Arrows mark initial values. See Figs. 2 and 3 for reference.
Open with DEXTER

The main reason for this kind of behaviour is that the increase in the energy of the reflected particle after a shock encounter depends on the projected speed of the field line (followed by the guiding center) along the shock front, $\Delta E \propto V_{\rm s} \tan{\theta _{\rm n} }$ , which grows very rapidly when $\theta _{\rm n} \rightarrow 90^{\circ}$ . Thus, in curvature dominated situations the acceleration process as a whole also tends to be very rapid, since the shock passes over quasi-perpendicular flux tubes in a relatively short time. Therefore, there should be little or no HEPs until the shock has transported the trapped particles to the quasi-perpendicular region, where they start to get accelerated. All the HEPs produced are then released to the downstream region in a short time interval.

For instance, in the line current runs (Fig. 2) the region where $\theta _{\rm n} > 80^{\circ}$ corresponds to a 1110 km wide area in the direction of shock normal. It takes roughly between 1 to 3 s for the shocks used in this study to travel this distance, and also to produce the observed HEPs. Figure 8 shows an example of a HT-resonant proton in the run $V_{\rm s}=400$ km s^-1, r=4, $R_{\rm curv} = 10~000$ km. The proton received a small boost to parallel velocity in the beginning of simulation run, making it travel near the shock for a long distance. It then hit the shock again in the quasi-perpendicular region, where it received its final energy in a short time and was finally transmitted to downstream region x > 0. Figure 9 shows an example of a Fermi-accelerated particle in the same run. Due to its higher parallel speed, this particle hit the shock for the first time when it was below the line current in y < 0 half plane. It was also transmitted to the downstream in the quasi-perpendicular region, but because of a smaller pitch it received a smaller final boost to its energy.

The line current setup can also be compared with recent test-particle simulations of Giacalone (2005a), who considered particle acceleration in a turbulent upstream medium. In Giacalone's model the wave lengths of the upstream turbulence extend to values much larger than the gyro radii of the accelerated particles, which in the quasi-perpendicular case creates loop structures in the upstream medium just as our line current model. Giacalone (2005a) found that perpendicular shocks were efficient in injecting particles to the acceleration process if such large-amplitude, large-wave-length fluctuations were present in the upstream medium (see also Giacalone 2005b). Another important finding of Giacalone (2005a) was that such upstream turbulence also makes quasi-parallel shocks accelerate particles faster than predicted by the diffusive shock acceleration theory. One possible explanation for this would be that the meandering upstream field lines overtaken by the shock form sites for HT resonance, which can operate much faster than the diffusive mechanism, which only relies on the fluctuations at resonance scales to provide chances for the particles to interact with shock.

The time evolution is completely different in the gradient dominated situations, where the energies of particles escaping to downstream region increase at a more or less constant rate. This is simply because the particles observed in the downstream region at a later time have, on average, encountered the shock more times than the particles observed at an earlier instant of time. Alternatively, we can say that the total adiabatic compression steadily increases over time. Figure 10 shows an example of an accelerated proton in the polar region run $V_{\rm s}=400$ km s^-1, r=4, $\theta _{\rm n} = 45^{\circ }$ . The proton has gained a small amount of energy during each bounce and was finally transmitted to downstream region x > 0with energy $\sim$ 350 keV.

$\begin{figure} \par\includegraphics[width=8.5cm,clip]{4754fi10.eps}\end{figure}$

Figure 10: Energy and position of a proton with initial pitch angle $\alpha = 100^{\circ }$ in the polar region run $V_{\rm s}=400$ km s^-1, r=4, $\theta _{\rm n} = 45^{\circ }$ , all values are measured in the SNIF. Energy is plotted vs. x-position using the secondary axes ( top and right), where the x-position (top axis, in km) has been shifted to right in order to clarify the figure. Arrows mark initial values. See Figs. 4 and 5 for reference.

Open with DEXTER

In most simulation runs we observe stripes in the phase space plots in one form or another, due to the bouncing motion of the particles, which are mainly a byproduct of the injection method we have used. We have injected all the particles in front of the shock at t=0 instead of injecting them at a constant rate along the field line, which would be a closer-to-realistic situation where the shock picks up particles from the upstream population present. In reality the contributions of all HEPs originating from different injection sites add up, after which the phase space stripes should disappear. However, in curvature-dominated situations, where the shock passes over loop structures, there should be a large number of counter-streaming high energy protons present, which in principle could be detected via satellite observations. The main difficulty is that suitable field configurations may be scarce in regions accessible to spacecraft. One possibility to detect these counter-streaming, accelerated protons would be in situations where an interplanetary shock driven by a fast CME is overtaking the closed magnetic field of a slower, previously ejected CME.

4.2 Number of accelerated particles

In addition to the differences in time evolution, shocks in gradient and curvature dominated configurations also produce very different amount of HEPs. When considering the initial pitch angle distribution, adiabatic compression operates on all particles not lying in the loss cone initially. In the polar region runs (see Fig. 5 for an example), protons with initial pitch angles $\alpha \la 25^{\circ}$ were in the loss cone of the shock and protons with $\alpha \ga 170^{\circ}$ hit the dipole surface. Between these limits, all other particles had a chance of gaining energy. The probability of getting accelerated is then determined solely by the initial phase angle of gyro motion, which is the only remaining quantity it can depend on. An example of this phase angle dependence is in Fig. 11, where we have plotted the final energy of particles as function of initial pitch and gyro phase angle for the polar region run $V_{\rm s}=400$ km s^-1, r=4, $\theta _{\rm n} = 45^{\circ }$ . Recall that the zero value for phase angle is in principle arbitrary but, of course, same for all particles in a given run.

$\begin{figure} \par\includegraphics[width=8.5cm,clip]{4754fi11.eps}\end{figure}$	Figure 11: Contour plot of initial phase and pitch angles vs. final energy for the polar region run $V_{\rm s}=400$ km s^-1, r = 4, $\theta _{\rm n} = 45^{\circ }$ . Particles with pitch $0 \leq \cos{\alpha } \leq 1$ are initially moving towards the shock.
Open with DEXTER

The large number of accelerated particles can be understood by considering a particle which is moving away from the shock along a magnetic field line in a setup similar to the polar region runs (Fig. 4). Let the magnetic field intensity scale as $B(s) = B_{0} (\frac{s_{0}}{s})^{n}$ along the field line. Assuming the conservation of magnetic moment, the pitch angle of particle approaches the value $90^{\circ}$ as it moves towards increasing B, $\sin{\alpha(s)} = \sin{\alpha _{0}} (s_{0} / s)^{n/2}$ . At the same time the distance between the shock and multipole origin shortens at a rate $V'_{\rm s} = V_{\rm s} / \cos{\theta _{\rm n} }$ . For sufficiently small $V'_{\rm s}$ (compared to initial parallel speed), most particles can move some distance away from the shock. Unless the particle resides in the loss cone of the upstream field, e.g. hits the dipole surface in our simulations, it will always encounter the shock again with higher pitch angle $\sin {\alpha } > \sin {\alpha _{0}}$ than initially. Since the condition of reflection off the shock is also proportional to $\sin {\alpha }$ , most particles initially moving towards the upstream region have a good chance of gaining energy. Exact treatment would require a solution of particle trajectory in upstream region which is difficult to obtain.

The catch here concerns the cases when the shock is nearly perpendicular. In order for, say, a 10 keV proton to run away from the shock, it needs to gain enough energy to make $v_{\vert\vert} > V'_{\rm s}$ during the first encounter with the shock. Since $V'_{\rm s}$ does not have an upper limit, there's always a limiting value for $\theta _{\rm n}$ after which all incident particles will be forced to the downstream region. After the $\theta _{\rm limit}$ has been reached, the acceleration process will resemble SDA in perpendicular geometry, which means that no HEPs are produced (according to our definition). In other words, the incident particles escape to downstream region, crossing the shock only once. This is visible in our polar region runs (Table 2), where $\theta _{\rm limit} < 70^{\circ}$ for $V_{\rm s} = 1000$ km s^-1 and $\theta _{\rm limit} < 80^{\circ}$ for other shock speeds.

This effect can be understood by considering the gradient force, $- \mu \nabla B$ . Since the undisturbed upstream field shows gradients, these persist even in the shocked regions, i.e., there is usually a gradient force pushing particles towards the downstream region on both sides of the shock front. There is of course another gradient force at the shock front due to the jump in tangential component of magnetic field, which is aligned towards the upstream region, but this affects particles only when they are crossing the shock. Under suitable conditions, the strong downstream gradient force can prevent particles from reaching the shock again during a single gyro orbit.

Regarding the cases where the curvature effects dominate, the phase angle dependence is less dramatic (see Fig. 12). For the Fermi-accelerated populations (compare to Fig. 3), the distribution of HEPs in phase angle space is almost isotropic. The HT resonant population with $0 \la \cos{\alpha } \la 0.2$ shows that there are small "hot spots'' in the velocity space, where the acceleration is most efficient. The phase angle dependence is evident also from Fig. 3.

$\begin{figure} \par\includegraphics[width=8.5cm,clip]{4754fi12.eps}\end{figure}$	Figure 12: Contour plot of initial phase and pitch angles vs. final energy for the line current run $V_{\rm s}=400$ km s^-1, r = 4, $R_{\rm curv} = 10~000$ km. Particles with pitch $0 \leq \cos{\alpha } \leq 1$ are initially moving away the shock.
Open with DEXTER

When analyzing the figures like 11 and 12, it is important to remember that the number of particles (typically 15 000) we have used in the simulations is not enough to produce smooth velocity space distributions. We feel that most "small area'' fluctuations in these figures should be regarded as particle noise rather than any real physical phenomena. A more detailed analysis requires that the simulations be rerun with a lot more particles and is beyond the scope of this paper.

4.3 Energies reached

Both types of field inhomogeneities are able to produce HEPs with energies of a few MeVs. The energies reached in curvature dominated geometry are slightly smaller than in gradient dominated one, but of course the number of HEPs are also quite different. What is perhaps a bit surprising is the dependence on the radius of curvature in the line current runs, since $R_{\rm curv} \gg r_{\rm gyro}$ even for all observed HEPs (see Table 1). However, the amount of time available for acceleration is, of course, very important and depends on the radius of curvature. If more time is available, particles can encounter the shock a larger number of times before ending up to the perpendicular region with higher energy.

In the polar region runs, the dependence on the oblique angle is what one would expect on the basis of SDA, i.e. the amount of energy given to particle during one shock encounter is higher when $\theta _{\rm n}$ is higher (see Table 2). If the shock is too oblique, $\theta _{\rm n} > \theta _{\rm limit}$ , the acceleration is cut off.

$\begin{figure} \par\includegraphics[width=8.5cm,clip]{4754fi13.eps}\end{figure}$	Figure 13: Energy and position of a proton with initial pitch angle $\alpha = 110^{\circ }$ in the rotated dipole run $V_{\rm s}=700$ km s^-1, r=4, all values are measured in the SNIF. Energy is plotted vs. x-position using the secondary axis (right) for energy values. Arrows mark initial values. See Figs. 6 and 7 for reference.
Open with DEXTER

The most interesting case (and also the most complicated one) is the general situation where the effects of curvature and gradients combine, like in the rotated dipole runs we have presented. A comparison of the results (Table 3) with results of other runs shows that the highest energies with all shock speeds and compression ratios have been produced in this geometry. Figure 13 shows an example of an accelerated proton in this geometry and the reason for the highest energies. Compare this figure with the setup for the rotated dipole runs (Fig. 6, particles with $\alpha > 90^{\circ}$ are initially moving away from the shock) to see the field geometry at each position.

At first glance it might seem that only the perpendicular region ( $y \sim 3000$ km) is important. However, in terms of energy gain ratio the acceleration here is not very efficient, $\Delta E/E_{0} \sim 3.5$ . This should be compared to the results of SDA. In perpendicular geometry $\Delta E/E_{0} \approx (r-1) \sin{}^{2} \alpha = 3$ at maximum with r=4 as in the example run. The observed energy gain in the perpendicular region is therefore only slightly higher than the one predicted by analytical theory. This is also in agreement with the results from the polar region runs, i.e. when $\theta _{\rm n} \ga 80^{\circ}$ in the presence of field gradients the acceleration process is similar to SDA in perpendicular geometry. What is more important is the acceleration during the whole time when the shock is transporting the particles to the perpendicular region. Here the field gradients have kept the particles close to the shock front at all times, since a proton with initial pitch angle $\alpha = 110^{\circ }$ should run further than 800 km away from the shock, for example. The most efficient acceleration, in terms of the energy gain ratio, has by far occurred during the first three bounces.

4.4 Energy spectra

We have calculated the energy spectra of accelerated particles only for the line current runs. In realistic situation the shock picks up particles continuously from the injection field line, whereas we have used a single injection point. Therefore a spectral index calculated from our results is not, in general, a meaningful quantity, since the simulations do not take into account the effect of particles injected at other positions.

However, the line current runs are less sensitive to this effect because of the time evolution of particle energies. Since the shock has to first transport the HT-resonant particles to the quasi-perpendicular region and because the Fermi-acceleration is of minor importance, we can expect the resulting energy spectra to be nearly invariant of the injection position. For the particles injected at the quasi-perpendicular region, the acceleration is mainly SDA at perpendicular geometry because the projected speed of the shock is too large for low-energy particles to be in HT-resonance with the shock. This results in an energy gain $\Delta E / E_{0} \approx 3$ at most, so these particles should affect only the low energy part of the resulting spectrum.

One should also note that in inhomogenous magnetic fields the energy spectra will not be time independent, since the energy of particles escaping to the downstream region increases as a function of time. The spectra ${\rm d}N/{\rm d}E \propto E^{\sigma }$ presented for the line current runs have been calculated from all HEPs, regardless of the time when the particles were removed from the simulation (Table 1). The error limits have been obtained from the least squares fitting. All spectra are linear in a logarithmic scale above the energy $E \ga 100$ keV with spectral indices in the range $\sigma \in [-3,-2]$ . Higher shock speeds and compression ratios produce harder spectra, reflecting the fact that the shocks with these properties are able to accelerate particles to higher energies. As a comparison, the spectra from the polar region runs are often flat, $\sigma \approx 0$ .

5 Conclusions and summary

We have shown that shocks propagating through inhomogenous upstream magnetic field configurations are able to produce high energy protons with energies in the range $\sim$ 1-10 MeV by scatter-free acceleration mechanisms. Magnetic geometries can be categorized to be gradient and curvature dominated, which both lead to situations where particles encounter the shock front multiple times. Accelerating protons interacting with the shock gain energy through (i) adiabatic compression, (ii) Fermi-acceleration and (iii) HT-resonance, which mainly operate in gradient (i) and curvature (ii & iii) dominated geometries.

The three acceleration mechanisms lead to different time evolution of HEPs in the upstream region. Fermi acceleration produces counter-streaming protons in loops being overtaken by shocks, where the energy of both trapped and escaping particles steadily increases with time. Also the adiabatic compression of the HEPs trapped in front of a shock propagating to an increasing field intensity operates continuously, but as the compression increases the perpendicular energy, the angular distributions are more of a pancake shape in this case. In contrast, HT-resonant particles trapped in front of the shock in a curved field are only modestly accelerated until they are convected to the quasi-perpendicular region, where they are quickly accelerated and released to the downstream region as a high energy beam.

Another noteworthy difference between the gradient and curvature driven mechanisms in our simulations is the number of high energy particles they produce. In the SNIF, only a small fraction of incident particles are HT-resonant and get accelerated, i.e., the ones with small parallel speed and pitch $0 \la \vert \cos{\alpha } \vert \la 0.2$ . Adiabatic compression, in contrast, picks up particles from the whole velocity phase space, excluding the loss cone of shock, resulting in a much larger number of HEPs.

We have studied the adiabatic compression and the HT resonance at non-relativistic energies. Both mechanisms could, in principle, be operative at relativistic energies as well. Adiabatic compression relies on the conservation of the first adiabatic invariant, which at the relativistic energies is $p_ \perp^2/B$ . Thus, the final energy in the relativistic case will be proportional to the square root of the magnetic field compression between the point of injection and the point of escape. The relativistic description of the HT resonance, however, is more complicated. Because the final velocity of the particle in the upstream frame is close to the projected shock speed at the time of escape, the limits of this mechanism would require detailed modeling of the trapping process in the non-inertial HT frame. This is beyond the scope of the present analysis, but will be an interesting subject for a future study.

HT resonance should be operative in the quasi-perpendicular parts of the supernova shocks, for example, because the ambient interstellar magnetic field is expected to have large-amplitude, long-wavelength MHD fluctuations. As the shock crosses these fluctuations, the upstream field consists of loops connected to the shock from both ends, which is an ideal setting for the mechanism. Since MHD turbulence is believed to be ubiquitous in the universe, the mechanism can be applied practicelly anywhere in connection with quasi-perpendicular shocks. Adiabatic compression, on the other hand, can be operative in any environment with shocks and magnetic field intensity gradients. The only requirement is that the shock has to propagate towards an increasing magnetic field. Plausible locations include stellar coronae and the magnetospheres of compact stars and giant planets.

Appandix A: Numerical model

The shock is treated as a planar discontinuity moving in upstream plasma rest frame with a constant velocity $\vec{V}_{\rm s}$ , which is perpendicular to the shock plane. We use an infinite-Mach-number approximation meaning that the shock is assumed to propagate with speed much larger than the Alfvén speed of the upstream plasma, $V_{\rm s} \gg V_{\rm A,n}$ , so that the tangential component of the downstream velocity field remains zero in the SNIF and the tangential magnetic field is compressed by the mass compression ratio. The compression ratio is taken as constant so that $1 \leq r \leq 4$ , where r=1 refers to the situation where there is no shock. Upstream plasma is also assumed to be at rest. With these assumptions the relevant ideal MHD shock jump conditions for the plasma flow and magnetic field in the SNIF are

$\begin{displaymath}% \frac{\partial \vec{B}_{1(2)}}{\partial t} = \nabla \times ( \vec{V}_{1(2)} \times \vec{B}_{1(2)} ), \end{displaymath}$

(A.5)

		$\displaystyle \vec{B}_{2}(t,x,y,z) = \vec{B}_{2}(t',x \rightarrow 0^{+},y,z),$	(A.6)
		t' = t - x / V₂,	(A.7)

The approximation used corresponds to a situation where the downstream plasma is suddenly jerked into motion with super-Alfvénic velocity $\vec{V}_{\rm s}$ , i.e. the shock is a kind of blast wave. The downstream solution is not valid for situations where the shock is driven by an obstacle (a CME, planet in the solar wind etc.), since in these cases the plasma must flow around the obstacle and $\vec{V}_{2}$ has to be obtained from, e.g., MHD simulation. The MHD solution also does not describe any microscopic structure of the shock (turbulence etc.), which may scatter escaping particles and return them back to the upstream, where they can be accelerated more. However, in the scatter-free limit the correct downstream conditions are really not needed, since the produced HEPs are accelerated in the upstream region. It should also be noted that a more general approach to jump conditions requires the specification of upstream velocity and mass density distribution.

The trajectories of the injected protons in the known electric and magnetic fields are computed using a second-order accurate leapfrog algorithm, commonly called the Boris-Buneman algorithm (for more information, see Ch. 4 of Birdsall & Langdon 1985), which has excellent energy conservation properties. All particles are injected on the same field line, i.e. the guiding center position of the injected particles is specified along with initial energy and pitch angle. Exact position and velocity vectors are then calculated from the single-particle theory using field values at guiding center position. The phase angle of the gyro motion is uniformly randomized to interval $\phi \in [0,2 \pi ]$ . The direction of $\phi = 0$ is, of course, same for all particles injected on a given position, but generally depends on the magnetic field configuration used.

Trajectories are computed until some specified end conditions are fulfilled. The ones we have used here are the distance to the dipole origin $d \leq d_{\rm min}$ , distance to the shock front in the downstream region $x \geq x_{\rm max}$ , and also a similar condition for time, $t \geq t_{\rm max}$ . When the trajectory computation for a given particle is finished, it's final state is stored. Various statistics are then calculated from the stored values, including the ones presented in this paper.

Particle quantities, specifically the energy, are computed in a local coordinate system where the particle drifts vanish. In this coordinate system the particle moves with a velocity $\vec{v}^{*} = \vec{v} - \vec{V}_{\rm drifts}$ . This system is not the same as the guiding center frame, since $\vec{v}^*$ has a nonvanishing component parallel to the magnetic field. Kinetic energy E^* calculated in this frame does not depend on the phase angle of gyro motion. The energy in the SNIF is $\langle E\rangle = E^{*} + \frac{1}{2}mV_{\rm drifts}^{2}$ , which corresponds to a gyro-averaged energy calculated directly from $\vec{v}$ . The advantage of this approach is that the computed energies do not fluctuate due to gyro motion, as can be seen, e.g., from Fig. 8.

		V_2n / V_1n = 1 / r,	(A.1)
		V_2t = V_1t = 0,	(A.2)
		B_2n / B_1n = 1,	(A.3)
		B_2t / B_1t = r,	(A.4)

Particle acceleration at shocks propagating in inhomogeneous magnetic fields

1 Introduction