- Email: [email protected]
Simulations of coronal shock acceleration in self-generated turbulence
ARTICLE IN PRESS
Journal of Atmospheric and Solar-Terrestrial Physics 70 (2008) 467–474 www.elsevier.com/locate/jastp
Simulations of coronal shock acceleration in self-generated turbulence R. Vainioa,, T. Laitinenb a
Department of Physical Sciences, University of Helsinki, Finland b Department of Physics, University of Turku, Finland Accepted 27 August 2007 Available online 5 October 2007
Abstract We have recently presented a numerical model to simulate the acceleration and transport of energetic particles at coronal shock waves, which includes the back-reaction of the energetic particles on the Alfve´n waves responsible for their scattering in the ambient medium [Vainio, R., Laitinen, T., 2007. Monte Carlo simulations of coronal diffusive shock acceleration in self-generated turbulence. Astrophysical Journal 658, 622–630]. The model results imply that the acceleration of energetic protons to energies exceeding 100 MeV can occur in some minutes in parallel coronal shocks driven by fast CMEs. In this paper we give a detailed description of the numerical model and its foreseen extensions. We extend the parameter study presented in the first paper to cover (shorter and) longer time scales, and present the intensities observed in the interplanetary medium during the initial phases of an SEP event. We also study the scattering mean free path resulting from the self-consistent computations. Finally, the potential of the model to explain ion acceleration up to relativistic energies is discussed. r 2007 Elsevier Ltd. All rights reserved. Keywords: Shock waves; Sun; Particle emission; Turbulence
1. Introduction Solar energetic particle (SEP) acceleration is associated with two types of eruptive phenomena: flares and coronal mass ejections (CMEs). While impulsive flares seem to be able to accelerate ions to energies of some tens of MeVs, SEP events related to fast CMEs can have accelerated ions in the GeV range. The most plausible model to account for these large SEP events is diffusive shock acceleration Corresponding author. Tel.: +358 9 191 50676;
fax: +358 9 191 50610. E-mail addresses: rami.vainio@helsinki.fi (R. Vainio), timo.laitinen@utu.fi (T. Laitinen). 1364-6826/$ - see front matter r 2007 Elsevier Ltd. All rights reserved. doi:10.1016/j.jastp.2007.08.064
(Axford et al., 1977; Bell, 1978; Blandford and Ostriker, 1978; Krymsky, 1977) in coronal shocks driven by fast CMEs. In this model, particle acceleration occurs via the first-order Fermi mechanism: particles scattered by magnetic fluctuations around the shock feel the compression of the plasma and get accelerated to high energies by repeatedly crossing the shock front. In order to explain acceleration to high energies, the shock needs to propagate into a very turbulent upstream medium. The upstream turbulence, however, can be generated by the accelerated particles themselves through the streaming instability of Alfve´n waves (AWs) resonant with the accelerated ions (Bell, 1978; Lee, 2005). In steady state and
ARTICLE IN PRESS 468
R. Vainio, T. Laitinen / Journal of Atmospheric and Solar-Terrestrial Physics 70 (2008) 467–474
planar geometry, the shock-accelerated ion distribution and the power-spectrum of the waves can be solved analytically (Bell, 1978). Recent models of coronal diffusive shock acceleration (e.g., Rice et al., 2003; Lee, 2005) have started from an assumption that the acceleration process is adiabatic, i.e., that the dynamical time scale of shock propagation, r=V s , equals the time available for particle acceleration. This assumption, combined with the estimate of the scattering conditions from the steady-state theory, predicts that the particle energy spectrum at the shock is a power law with a cutoff at high energies determined by the steadystate acceleration rate and the available acceleration time. Particle escape from the shock to the far upstream region (i.e., observer) is modeled either by using an ad hoc free-escape boundary a few diffusion lengths upstream of the shock (Rice et al., 2003) or by a self-consistent steady-state mechanism based on adiabatic focusing in the upstream region (Lee, 2005). Following another approach, based on an investigation of the time evolution of the unstable waves, Vainio (2003) also derived an energy spectrum of the escaping energetic protons. He found that a relatively high number of particles needs to escape from the shock before the waves have grown appreciably in the vicinity of the shock wave. Motivated by this result, we developed a numerical model to simulate the coupled SEP acceleration and AW generation in a fully timedependent manner (Vainio and Laitinen, 2007). Our results indicated that while the acceleration process close to the shock is probably well approximated— at least in quasi-parallel regions of the shock—by the quasi-steady-state modeling, particle escape is a more dynamical process. In this paper we will give a detailed description of the numerical method to simulate coupled particle acceleration and turbulence evolution. We will extend the parameter study presented by Vainio and Laitinen (2007) to cover longer time evolution, and discuss the modeled wave intensities in a more detailed manner. We will also present results on particle acceleration beyond 100 MeV and discuss the possibility to achieve relativistic energies in our model. Finally, we will outline the next steps of our model development in order to reduce the effects of the current simplifications to the model results, and to study whether the remaining discrepancies between the model and the observations of SEP events can be reduced.
2. The model We consider the acceleration of protons in a single magnetic flux tube around a radial central magnetic field line and a super-radially expanding crosssectional area, A / B1 . Specifically, we take the field magnitude to be BðrÞ ¼ B0 ðr =rÞ2 ½1þ bf ðr =rÞ6 , where B0 ¼ 2:90 nT, r ¼ 1 AU, bf ¼ 1:9, and r is the solar radius. The electron density, nðrÞ, the solar–wind speed, uðrÞ, and the Alfve´n speed, pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi vA ðrÞ ¼ B= m0 m=e n, are calculated by assuming the conservation of mass in the solar wind, un=B ¼ const:, the constancy of the group speed of the AWs in the inertial frame, V ¼ u þ vA ¼ const:, and the values of V ¼ 400 km s1 and nðr Þ ¼ 10 cm3 . (m0 is the permeability of vacuum.) We consider a hydrogen plasma, i.e., we take the solar–wind mass per electron to equal the proton mass, m=e ¼ mp . The assumption of the constancy of V is of ad hoc nature and made out of convenience. In reality, of course, there is no physical process in the corona that would force V to be constant. However, despite of their great simplicity our assumptions lead to a reasonable model of a slow solar–wind stream (Vainio et al., 2003; Vainio and Laitinen, 2007). The AWs are transported along the magnetic field lines at the constant group speed V following the principles of geometric optics. Their generation is governed by the exponential growth caused by the streaming of resonant ions. We use a simplified wave–particle resonance condition, k ¼ r1 L , where k and rL ¼ gv=ocp are the AW wavenumber and the proton Larmor radius, with g, v and ocp ¼ eB=mp denoting the proton Lorentz factor, speed and angular cyclotron frequency.1 Then the growth rate can be approximated by (Vainio, 2003) pSp ðr; tÞ . (1) nvA R þ1 Here S p ðr; tÞ ¼ 2pp2 v 1 mF ðr; p; m; tÞ dm is the (wave-frame) streaming per unit momentum of protons to be evaluated at the resonant (waveframe) momentum p ¼ mp Vf cp =f , F is the distribution function of protons in the wave frame, f cp ¼ ocp =ð2pÞ, and f is the wave frequency in the inertial frame. We also introduce a frequency diffusion term, mimicking the effect of wave–wave
sðf ; r; tÞ ¼ p2 f cp
1 The full quasi-linear resonance condition reads k ¼ 1=ðrL mÞ, implying that particles with a certain momentum p resonate with waves at a range of wavenumbers, jkjXeB=p.
ARTICLE IN PRESS R. Vainio, T. Laitinen / Journal of Atmospheric and Solar-Terrestrial Physics 70 (2008) 467–474
interactions, to the AW transport equation, DP~ qP~ qP~ q qP~ þV ¼ sP~ þ Dff , (2) Dt qt qr qf qf with P~ ðV 2 =BvA ÞP. The diffusion coefficient is set 2=3 to Dff ¼ ðV =r Þf 8=3 f b , which makes the power spectrum erode at high frequencies and tend, at 1 AU, toward the observed form (e.g., Horbury et al., 1996) of Kolmogorov spectrum of turbulence, P / f 5=3 , at frequencies greater than f b ¼ 1 mHz. In addition to AWs, our simulation keeps track of a large number of protons propagating in the flux tube. The AW growth rate s is self-consistently computed from the streaming of the simulated protons. Proton propagation is governed by the conservation of first adiabatic invariant, p2? =B, and momentum p in the inertial frame. These effects are represented by the equations of motion (in the inertial frame) r_ ¼ mv;
v ¼ const.
1 m2 q ln B . v m_ ¼ qr 2 In addition to the adiabatic propagation, protons are subject to the effect of resonant AWs, which scatter them (elastically in the local wave frame) at rate fPðf ; r; tÞ , (3) B2 where again the resonance condition f ¼ f cp mp V =p, with p measured in the wave frame, is applied. Thus, AWs and SEPs form a coupled, nonlinear system. The agent responsible for the acceleration of SEPs is a shock front postulated to expand radially at a constant speed V s . We only consider particle propagation upstream of the shock front and treat the front as a boundary condition in the equation of motion of the protons. In each interaction with the shock, particles returning back upstream have gained a small amount of energy (see below). Thus, protons scattered to the shock many times will gain energies up to hundreds of MeVs and beyond. For further details of the model, see Vainio and Laitinen (2007). nðr; p; tÞ ¼ p2 f cp
3. The numerical method We solve the transport equation of waves and the equation of motion of the particles simultaneously. We use a finite difference scheme to solve the wave transport equation on a Lagrangian grid (taking
469
care of wave propagation), which has a logarithmic spacing in the initial radial distance, i.e., xi ¼ ri Vt with xi ¼ r0 expðiDx Þ and a constant Dx 51. Also the frequency grid is logarithmic, f j ¼ f 0 expðjDf Þ with a constant Df . Wave growth is proportional to the net flux of particles integrated over the time step of wave evolution, Dt: ~ i ; f j ; t þ DtÞ P~ ij ðt þ DtÞ Pðx Z tþDt 0 0 ~ sðxi ; f j ; t Þ dt ¼ Pij ðtÞ exp ¼ P~ ij ðtÞ exp
Z t
t tþDt
pSp ðxi ; t0 Þ 0 dt , p f cp ðri Þ nðri ÞvA ðri Þ 2
where the integral can be computed by adding up the net number of particles crossing the point x ¼ xi from below during the time step from t to t þ Dt. Each particle contributes to the balance at two grid levels of frequency, f j and f jþ1 , just below and above the resonant frequency f ¼ f cp ðri ÞmV =p, respectively, with weights wi;j ¼
ln f ln f j p2 f cp ðri Þ W l, Df Aðri Þnðri ÞvA ðri ÞDf
wi;jþ1 ¼
ln f jþ1 ln f p2 f cp ðri Þ Wl, Aðri Þnðri ÞvA ðri ÞDf Df
where the sign is positive [negative] for protons crossing the spatial grid level xi from below [above], and W l is the weight of the particle in the simulation. We use AðrÞ ¼ r2 B =BðrÞ and, thus, consider a flux tube which has a solid angle of 1 sr at the solar surface. After the wave growth, we update the effect of diffusion to the wave spectrum at each spatial grid level by using a standard Crank– Nicholson scheme (Morton and Mayers, 1994). Particles are traced allowing each of them to have its own time step fulfilling dtl ¼ 2nl Dt, where nl 2 N0 and l numbers the particles. The time step for each particle is selected so that it tries to remain constrained by nl dtl 2 ½0:02; 0:05. However, a still smaller time step is required if the particle would otherwise be able to cross two spatial grid levels at once, i.e., dtl oDxi =vl is required. The value of nl can be decreased after each time step dtl whenever the conditions call for this, but it can be increased only in the beginning of each wavetransport cycle, Dt. Particle propagation, focusing caused by the conservation of first adiabatic invariant, and scattering by the AWs are solved using a leapfrog method, i.e., first moving the particle by half a time
ARTICLE IN PRESS R. Vainio, T. Laitinen / Journal of Atmospheric and Solar-Terrestrial Physics 70 (2008) 467–474
470
step, rl ¼ rl þ vl ml dtl =2, then performing the changes to the velocity and, finally, moving the particle by another half time step rl rl þ vl ml dtl =2. The effects of focusing and scattering on the pitch-angle cosine ml and speed vl are updated at the intermediate position, at the solar fixed and wave frame, respectively, as follows: focusing:
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 þ 4dðml þ dÞ , ml 2d vl dtl q ln B d¼ . qr r 2 1 þ
l
Scattering: ðvl ; ml Þ LTðvl ; ml ; V Þ, qffiffiffiffiffiffiffiffiffiffiffiffiffi ml ml cos W þ 1 m2l sin W cos j, ðvl ; ml Þ
LTðvl ; ml ; V Þ,
where LTðv; m; uÞ denotes a Lorentz transformation to a system moving at speed u along the magnetic field, and the scattering angles W and j are obtained using a random generator. (See, e.g., Kocharov et al., 1998; Vainio et al., 2000, for a detailed description of the scattering process.) Specifically, W2 and j are picked for each scattering from an exponential distribution with a mean of 2nl dtl and a uniform distribution between 0 and 2p, respectively, and the scattering frequency nl ¼ nðrl ; pl ; tÞ is obtained using the simulated wave spectrum in Eq. (3) with logarithmic interpolation. The method for simulating focusing is the implicit Euler scheme, which is stable and always gives a value inside the allowed range, m 2 ½1; 1. In case the particle crosses a spatial grid level xi of the wave spectrum, the time-integrated growth rate is updated there, as described above. We inject low-energy particles at a constant rate to the simulation at the inner boundary, taken to coincide with the propagating shock front. (Thus, the inner boundary of our simulation box is moving at speed V s V with respect to the wave grid.) The velocity distribution of the injected particles (in the shock frame) is dN=dv ¼ N inj Hðv u1 Þ= v1 expððv u1 Þ=v1 Þ, where u1 ¼ V s uðrs Þ is the local bulk speed of the upstream plasma in the frame of the shock and v1 is a model parameter which describes the spread of the suprathermal injection, set to v1 ¼ 375 km s1 in this paper. The total number of injected particles inside the con-
sidered flux tube is N inj . The pitch-angle cosines of the injected particles are drawn from a fluxweighted (see Jones and Ellison, 1991) isotropic distribution in the shock frame, i.e., PðmÞ ¼ 2jmj. Particle acceleration takes place at the shock front: every time a simulated particle hits the shock front from upstream, we compute its speed, v0 , in the rest-frame of the downstream scattering centers, moving with speed V 2 ¼ V 1 =Rsc with respect to the shock toward the Sun. Here, V 1 ¼ V s V is the speed of the AWs with respect to the shock in the upstream region and Rsc is the scattering-center compression ratio of the shock, which we regard as a constant model parameter.2 After that, we turn the particle back toward the upstream region drawing its downstream wave-frame pitch-angle cosine m0 from an isotropic distribution weighted by the shock frame parallel speed, i.e., Pðm0 Þ ¼ 2
jm0 þ V 2 =v0 j ; ð1 V 2 =v0 Þ2
m0 2 ½1; V 2 =v0 ,
m0 ¼ 1 corresponding to the direction of the shockframe flow (i.e., toward the Sun). However, as the downstream bulk flow is away from the shock, not all particles can make it back to the shock and, thus, we have to decrease the weight of the particles in each shock encounter by 0 v V2 2 Wl W l Pret ; Pret ¼ 0 , v þ V2 where Pret is the probability of return from the downstream region calculated assuming that the particle distribution around the shock is quasiisotropic in the rest frame of the downstream scattering centers3 (Jones and Ellison, 1991). We also assume that the returning particles spend a negligible amount of time in the downstream region, which is equivalent to assuming that the scattering mean free path is very small there. By considering the Lorentz transformation between the upstream and downstream wave frames, it is easy to see that the energy of the returning particle in the upstream wave frame is increased by DE ¼ DVp0 ðm01 m02 Þ, where DV ¼ V 1 V 2 and m01½2 is the pitch-angle cosine of the particle (in the 2 Note that in the case of AWs, Rsc is not identical to the gas compression ratio, R ¼ n2 =n1 , of the shock; Rsc can be larger than R in a parallel shock wave, even substantially (Vainio and Schlickeiser, 1998). 3 This assumption holds by construction for the particles returning from downstream, and rather well for the incident particles at energies of interest in our study.
ARTICLE IN PRESS R. Vainio, T. Laitinen / Journal of Atmospheric and Solar-Terrestrial Physics 70 (2008) 467–474
downstream wave frame), measured before [after] the interaction with the shock. Thus, particles gain a small amount of energy (hDEi 43 DVp0 ) in each interaction with the shock. By returning all particles from the downstream region and decreasing their weight accordingly, we keep the statistics of the simulated particles relatively uniform over the whole energy spectrum. 4. Results and discussion Vainio and Laitinen (2007) presented a parameter study, which concentrated on finding out the dependence of the achieved maximum proton energy on selected model parameters, namely N inj , Rsc and V s . They considered cases with a simulation time of tmax ¼ 640 s, initial shock radius of r0 ¼ 1:5r , and an ambient wave spectrum solved from Eq. (2) using s ¼ 0, the assumption of a steady state, and the inner boundary condition fPðf ; r0 Þ ¼ P B2 ðr0 Þ, with P ¼ 6 107 , corresponding to a mean free path of l0 ¼ 1r of 100-keV protons at r ¼ r0 . In this paper, we will take a closer look to three model runs with r0 ¼ 1:5r , V s ¼ 1500 km s1 and Rsc ¼ 4.
471
We will consider (A) N inj ¼ 3 1035 sr1 , tmax ¼ 640 s, P ¼ 6 107 ; (B) N inj ¼ 3 1035 sr1 , tmax ¼ 160 s, P ¼ 1:2 107 ; and (C) N inj ¼ 1:69
1036 sr1 , tmax ¼ 4320 s, P ¼ 2:0 108 . By studying these cases we aim at establishing: (i) how the particle mean free path in our simulation model compares with analytical estimates from diffusive shock acceleration (Bell, 1978); (ii) whether the ambient mean free path has any effect on the acceleration efficiency of the shock; (iii) if the injection rate has a direct effect on the acceleration efficiency or is the total number of injected particles the decisive parameter. In Fig. 1, we present the mean free path and the particle intensities of run A in the end of the simulation as a function energy (for rs r0 þ V s tmax 2:88r , rs þ r and rs þ 10r ) and position (for E ¼ 0:29, 2.4 and 20 MeV) along with the predictions from the one-dimensional analytical steady-state theory of Bell (1978) for the 2.4-MeV mean free path and intensity as a function of distance and the mean free path and intensity at the shock as a function of energy. Following Vainio and Laitinen (2007) we present the differential intensities, IðE; rÞ, scaled with the flux-tube crosssectional area and normalized to 1 AU, i.e.,
107 106 105 104 103 102 101 101 100 10−1 10−2 10−3
10−3
10−2
10−1
100
101
10−1
100
101
102
E [MeV] Fig. 1. Particle intensity (scaled to 1 AU) and scattering mean free path as a function of distance from the shock (left column) and energy (right column) and in the end of the model run A. In the left column, the black curves denote the simulated values at E ¼ 0:29 MeV (dashed), 2.4 MeV (solid), and 20 MeV (dash-dotted) and the gray curves denote the analytical estimates for 2.4 MeV. In the right column, the black curves denote the simulated values at z ¼ 0 (solid), 1r (dashed) and 10r (dash-dotted) and the gray curves denote the analytical estimates for z ¼ 0.
ARTICLE IN PRESS R. Vainio, T. Laitinen / Journal of Atmospheric and Solar-Terrestrial Physics 70 (2008) 467–474
I ðE; rÞ ðA=A ÞIðE; rÞ. Such a scaling is necessary when comparing with one-dimensional theory (where the flux tube has a constant cross-section). Assuming vanishing wave and particle intensities far upstream, the analytical one-dimensional steadystate result is (Bell, 1978; Vainio and Laitinen, 2007) z0 , z0 þ z
ð4Þ
Vs V ðz þ z0 Þ, v
ð5Þ
I ðE; zÞ ¼ I ðE; rs Þ lðz; pÞ ¼ 3
where z ¼ r rs is the distance from the shock, z0 ¼ hðp=p0 Þb3 , b ¼ 3Rsc =ðRsc 1Þ ¼ 4, p0 is the average injection momentum and 3b N inj =tmax b p I ðE; rs Þ , ð6Þ A V s 4pp p0 V s Aðrs Þnðrs Þ vA ðrs Þ . ð7Þ h p2 b1 N inj =tmax f cp ðrs Þ Although not a perfect fit, the analytical theory provides a valid approximation to the simulation results at energies below the spectral cutoff and at distances close to the shock. Note, however, that at the lowest energies the simulated mean free path is clearly higher than the predicted value, which is most likely due to the erosion of the spectrum by frequency diffusion (not included in the analytical theory) to the region without any resonant particles. Further out from the shock ðz\r Þ and at large energies ðE\10 MeVÞ the discrepancies are large, as expected, because of the effects of spherical geometry and finite acceleration time, and the advance of the frequency diffusion at larger heliocentric distances. One can also see the coherently propagating escaping particles in the upper left panel of Fig. 1. These particles are the ones that escape to the ambient medium before the waves have grown enough to turbulently trap them in front of the shock wave. As shown by Vainio and Laitinen (2007), the time integrated flux of this particle population is very well approximated by the theoretical prediction of Vainio (2003), who noted that in order for the waves to have grown appreciably in a given position and frequency, i.e. that growth rate becomes large, R 0 the time-integrated dt sðf ; r; t0 Þb1, the time-integrated net flux of particles has to exceed a threshold value, i.e., Z t nvA dt0 pS p ðr; t0 Þb 2 . p f cp 0
Before this occurs, the particles can easily escape from the system if the ambient conditions are not very turbulent, i.e., if the initial mean free path satisfies l0 ðr; pÞ\ðq ln B=qrÞ1 r=2, as in our simulations. Simulation run B has the same amount of injected particles as run A, but the ambient mean free path is 5 times larger, and the simulation time is 4 times shorter. (Thus, the injection rate is 4 times greater than in run A.) The scaled spectrum at the shock, I ðE; rs Þ, is very similar in the two cases: both runs show good coincidence with the theoretical powerlaw energy spectrum up to a cutoff energy, which is very similar in both cases, E\10 MeV. Thus, while the normalization of the energy spectrum at the shock is proportional to the injection rate, N inj =tmax , it is the total number of injected protons that determines the maximum energy in our model, at least in the simulated cases. This is actually the prediction of the diffusive shock acceleration theory as well, which yields for the cutoff momentum (Vainio and Laitinen, 2007) ðRsc 1Þ=3 1 M 1 pc ðtmax Þ p0 C N inj , (8) Rsc 1 where M ¼ V s =V and C is a constant, which has a value of o2 1033 sr for the conditions in our model corona. In Fig. 2 we present the proton intensities of run C at a distance of z ¼ 0:3 AU upstream from the shock. (Note that these are local intensities not scaled to 1 AU.) The run is now much longer and we can follow the prompt phase of the event as it would be seen by spacecraft far upstream. This run had 104 I [1 / (cm2 sr s MeV)]
472
103
102
101
0
20
40 t [min]
60
80
Fig. 2. Particle intensity at z ¼ 0:3 AU as a function of time in the simulation run C at four energies: 6 MeV (dashed), 10 MeV (solid), 30 MeV (dot-dashed) and 50 MeV (dotted).
ARTICLE IN PRESS R. Vainio, T. Laitinen / Journal of Atmospheric and Solar-Terrestrial Physics 70 (2008) 467–474
a very long ambient mean free path (30 times that of run A, i.e., several AUs at r ¼ 0:3 AU), which leads to a very efficient decoupling of the promptly escaping particles from the shock and these particles appear as an impulsive injection from the Sun although the acceleration process was all due to a coronal shock wave with a constant speed, compression ratio and injection rate of low-energy seed particles. The very large value of the ambient mean free path leads to a very fast decay of the intensity, and for a more realistic interplanetary mean free path the prompt component would have a slower increase to the time of maximum and slower decay. Our simulations do, however, imply that a remotely observed time-intensity profile that shows two maxima does not necessarily mean that there are two phases of acceleration at the Sun. One can also see from Fig. 2 that the lowest energy channels show the typical velocity dispersion, i.e., low-energy particle intensities have later onset than the high-energy particle intensities. However, above 10 MeV, our model produces the opposite effect. This is because the acceleration process, for the parameters of the simulation run, takes a long time and the first particles above 10 MeV appear clearly later than at lower energies. Thus, typical velocity dispersion results with highest energies having the earliest onset imply that the acceleration process is much faster than in the present simulation run. This can be achieved most easily by increasing the injection rate. At first sight, there would be nothing preventing the simulation model to generate protons at relativistic energies: all we seem to need is to increase the number of injected particles (see Eq. (8)). Our model, however, relies on the quasilinear theory of wave–particle interactions. Thus, we should make sure that the value of scattering frequency n stays well below f cp , implying (see Eq. (3)) that fPðf ; rÞ 5p2 0:1 B2 ðrÞ would have to be fulfilled. For the model parameters V s ¼ 1500 km s1 , Rsc ¼ 4 and tmax ¼ 640 s, the simulations imply that the limit of the validity of quasi-linear approximation will be reached roughly when pc mp c. Thus, our simulation model may need more physics before it can safely be used to describe particle acceleration to the relativistic energies.
473
5. Outlook and conclusions We have employed a number of simplifications to our model in order to achieve computational efficiency. However, now that the feasibility of the approach has been demonstrated, it is crucial to study the influences of the effects we have neglected so far. The solar wind model, while conveniently allowing us to use a Lagrangian grid with constant spacings for transporting the wave power, limits our modeling to speed profiles, which do not have a maximum of the Alfve´n velocity in the corona, often appearing at height of a few solar radii in both semiempirical and analytical modeling of the solar wind. Such a feature may prove to be important: as the shock crosses the region of maximum vA its Mach number gets temporarily decreased, but at the same time the wave growth becomes most efficient (Vainio, 2003; Vainio and Laitinen, 2007). In order to model these effects, it is probably most reasonable to abandon the Lagrangian grid approach, and use a simple simulation scheme, such as upwind differencing, to model the wave transport in more general solar wind models. Our modeling of the spectral diffusion was based on a simple linear diffusion coefficient. However, in a more realistic approach reproducing the effects of wave–wave interactions, non-linear diffusion approach would be needed, such as that described in Vainio et al. (2003). The particle modeling also suffers from simplifications which must be accounted for. First, treating the shock simply as lower boundary condition with a constant scattering-center compression ratio neglects important effects due to wave transmission/reflection at the shock (Vainio and Schlickeiser, 1998). To account for these, particles and waves should be self-consistently traced in the downstream region, including modes propagating both toward and away from the Sun. Second, the simplified resonance condition used in these simulations essentially maps the particle to one resonant wave number, while the exact, pitch-angle dependent resonance condition spreads the resonance to a range of wavenumbers. This will influence the evolution of the wave spectra and, thus, also the spectrum of escaping particles: since fast particles can now generate waves that scatter the slower ones, the escaping spectrum will probably also depend on the temporal evolution rate of the maximum energy in the spectrum. Using the full quasi-linear
ARTICLE IN PRESS 474
R. Vainio, T. Laitinen / Journal of Atmospheric and Solar-Terrestrial Physics 70 (2008) 467–474
resonance condition, k ¼ 1=ðrL mÞ, will, however, require two major changes in the model: (i) the scattering frequency has to be taken as m-dependent, which typically requires much shorter time steps than an isotropic scattering process (Kocharov et al., 1998) and thus slows down the computational efficiency of the simulations by a considerable factor; and (ii) the AW growth rate has to be evaluated using the full quasi-linear theory (see, e.g., Vainio, 2003), rendering the simple book-keeping method of particles crossing grid cells invalid. Fortunately a similar method may be devised for the full resonance condition, since the growth rate in this case becomes proportional to the net flux of particles in momentum space, i.e., to Dmm qF =qm, where Dmm is the pitch-angle diffusion coefficient, evaluated at the resonant parallel momentum pk ¼ mocp =k. Thus, we may compute the timeintegrated growth rate now by a similar bookkeeping of particles crossing the resonant frequency levels f i through scattering. However, implementing these effects to the model will be considerably more tedious than the modifications to the solar–wind and wave-transport models. In conclusion, we have presented a numerical method to simulate particle acceleration and coupled AW evolution in CME-driven shocks in the solar corona. Our model is capable of accelerating protons to 100 MeV and beyond in a few minutes. The model is based on a simplified version of quasi-linear theory, and thus its range of validity is constrained by the requirement of linear wave amplitudes. This indicates that although relativistic energies (p mc) may be reached in the model within 10 min, acceleration faster than that may require more complicated approach on wave–particle and wave–wave interactions than employed in our model. The model produces a coherently escaping particle population with a hard energy spectrum ð/ p1 Þ, but this feature of the model may be especially sensitive to the simplifications made on modeling the wave–particle interactions. Further modifications to the model were outlined, by which these and other remaining simplifications of the model may be removed.
Acknowledgments We thank the Finnish Meteorological Institute for the use of their computing facilities. References Axford, W.I., Leer, E., Skadron, G., 1977. The acceleration of cosmic rays by shock waves. In: Proceedings of the 15th International Cosmic Ray Conference, vol. 11, pp. 132–137. Bell, A.R., 1978. The acceleration of cosmic rays in shock fronts. I. Monthly Notices of the Royal Astronomical Society 182, 147–156. Blandford, R.D., Ostriker, J.P., 1978. Particle acceleration by astrophysical shocks. Astrophysical Journal 221, L29–L32. Horbury, T.S., Balogh, A., Forsyth, R.J., Smith, E., 1996. The rate of turbulent evolution over the Sun’s poles. Astronomy and Astrophysics 316, 333–341. Jones, F.C., Ellison, D.C., 1991. The plasma physics of shock acceleration. Space Science Reviews 58, 259–346. Kocharov, L., Vainio, R., Kovaltsov, G.A., Torsti, J., 1998. Adiabatic deceleration of solar energetic particles as deduced from Monte Carlo simulations of interplanetary transport. Solar Physics 182, 195–215. Krymsky, G.F., 1977. A regular mechanism for the acceleration of charged particles on the front of a shock wave. Doklady Academie Nauk SSSR 243, 1306–1308 (Translation: Soviet Physics—Doklady 22, 327–328.). Lee, M.A., 2005. Coupled hydromagnetic wave excitation and ion acceleration at an evolving coronal/interplanetary shock. Astrophysical Journal Supplement 158, 38–67. Morton, K.W., Mayers, D.F., 1994. Numerical Solution of Partial Differential Equations. Cambridge University Press, Cambridge, UK. Rice, W.K.M., Zank, G.P., Li, G., 2003. Particle acceleration and coronal mass ejection driven shocks: shocks of arbitrary strength. Journal of Geophysical Research 108, CiteID 1369, 10.1029/2002JA009756. Vainio, R., 2003. On the generation of Alfve´n waves by solar energetic particles. Astronomy and Astrophysics 406, 735–740. Vainio, R., Laitinen, T., 2007. Monte Carlo simulations of coronal diffusive shock acceleration in self-generated turbulence. Astrophysical Journal 658, 622–630. Vainio, R., Schlickeiser, R., 1998. Alfve´n wave transmission and particle acceleration in parallel shock waves. Astronomy and Astrophysics 331, 793–799. Vainio, R., Kocharov, L., Laitinen, T., 2000. Interplanetary and interacting protons accelerated in a parallel shock wave. Astrophysical Journal 528, 1015–1025. Vainio, R., Laitinen, T., Fichtner, H., 2003. A simple analytical expression for the power spectrum of cascading Alfve´n waves in the solar wind. Astronomy and Astrophysics 407, 713–723.
Journal of Atmospheric and Solar-Terrestrial Physics 70 (2008) 467–474 www.elsevier.com/locate/jastp
Simulations of coronal shock acceleration in self-generated turbulence R. Vainioa,, T. Laitinenb a
Department of Physical Sciences, University of Helsinki, Finland b Department of Physics, University of Turku, Finland Accepted 27 August 2007 Available online 5 October 2007
Abstract We have recently presented a numerical model to simulate the acceleration and transport of energetic particles at coronal shock waves, which includes the back-reaction of the energetic particles on the Alfve´n waves responsible for their scattering in the ambient medium [Vainio, R., Laitinen, T., 2007. Monte Carlo simulations of coronal diffusive shock acceleration in self-generated turbulence. Astrophysical Journal 658, 622–630]. The model results imply that the acceleration of energetic protons to energies exceeding 100 MeV can occur in some minutes in parallel coronal shocks driven by fast CMEs. In this paper we give a detailed description of the numerical model and its foreseen extensions. We extend the parameter study presented in the first paper to cover (shorter and) longer time scales, and present the intensities observed in the interplanetary medium during the initial phases of an SEP event. We also study the scattering mean free path resulting from the self-consistent computations. Finally, the potential of the model to explain ion acceleration up to relativistic energies is discussed. r 2007 Elsevier Ltd. All rights reserved. Keywords: Shock waves; Sun; Particle emission; Turbulence
1. Introduction Solar energetic particle (SEP) acceleration is associated with two types of eruptive phenomena: flares and coronal mass ejections (CMEs). While impulsive flares seem to be able to accelerate ions to energies of some tens of MeVs, SEP events related to fast CMEs can have accelerated ions in the GeV range. The most plausible model to account for these large SEP events is diffusive shock acceleration Corresponding author. Tel.: +358 9 191 50676;
fax: +358 9 191 50610. E-mail addresses: rami.vainio@helsinki.fi (R. Vainio), timo.laitinen@utu.fi (T. Laitinen). 1364-6826/$ - see front matter r 2007 Elsevier Ltd. All rights reserved. doi:10.1016/j.jastp.2007.08.064
(Axford et al., 1977; Bell, 1978; Blandford and Ostriker, 1978; Krymsky, 1977) in coronal shocks driven by fast CMEs. In this model, particle acceleration occurs via the first-order Fermi mechanism: particles scattered by magnetic fluctuations around the shock feel the compression of the plasma and get accelerated to high energies by repeatedly crossing the shock front. In order to explain acceleration to high energies, the shock needs to propagate into a very turbulent upstream medium. The upstream turbulence, however, can be generated by the accelerated particles themselves through the streaming instability of Alfve´n waves (AWs) resonant with the accelerated ions (Bell, 1978; Lee, 2005). In steady state and
ARTICLE IN PRESS 468
R. Vainio, T. Laitinen / Journal of Atmospheric and Solar-Terrestrial Physics 70 (2008) 467–474
planar geometry, the shock-accelerated ion distribution and the power-spectrum of the waves can be solved analytically (Bell, 1978). Recent models of coronal diffusive shock acceleration (e.g., Rice et al., 2003; Lee, 2005) have started from an assumption that the acceleration process is adiabatic, i.e., that the dynamical time scale of shock propagation, r=V s , equals the time available for particle acceleration. This assumption, combined with the estimate of the scattering conditions from the steady-state theory, predicts that the particle energy spectrum at the shock is a power law with a cutoff at high energies determined by the steadystate acceleration rate and the available acceleration time. Particle escape from the shock to the far upstream region (i.e., observer) is modeled either by using an ad hoc free-escape boundary a few diffusion lengths upstream of the shock (Rice et al., 2003) or by a self-consistent steady-state mechanism based on adiabatic focusing in the upstream region (Lee, 2005). Following another approach, based on an investigation of the time evolution of the unstable waves, Vainio (2003) also derived an energy spectrum of the escaping energetic protons. He found that a relatively high number of particles needs to escape from the shock before the waves have grown appreciably in the vicinity of the shock wave. Motivated by this result, we developed a numerical model to simulate the coupled SEP acceleration and AW generation in a fully timedependent manner (Vainio and Laitinen, 2007). Our results indicated that while the acceleration process close to the shock is probably well approximated— at least in quasi-parallel regions of the shock—by the quasi-steady-state modeling, particle escape is a more dynamical process. In this paper we will give a detailed description of the numerical method to simulate coupled particle acceleration and turbulence evolution. We will extend the parameter study presented by Vainio and Laitinen (2007) to cover longer time evolution, and discuss the modeled wave intensities in a more detailed manner. We will also present results on particle acceleration beyond 100 MeV and discuss the possibility to achieve relativistic energies in our model. Finally, we will outline the next steps of our model development in order to reduce the effects of the current simplifications to the model results, and to study whether the remaining discrepancies between the model and the observations of SEP events can be reduced.
2. The model We consider the acceleration of protons in a single magnetic flux tube around a radial central magnetic field line and a super-radially expanding crosssectional area, A / B1 . Specifically, we take the field magnitude to be BðrÞ ¼ B0 ðr =rÞ2 ½1þ bf ðr =rÞ6 , where B0 ¼ 2:90 nT, r ¼ 1 AU, bf ¼ 1:9, and r is the solar radius. The electron density, nðrÞ, the solar–wind speed, uðrÞ, and the Alfve´n speed, pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi vA ðrÞ ¼ B= m0 m=e n, are calculated by assuming the conservation of mass in the solar wind, un=B ¼ const:, the constancy of the group speed of the AWs in the inertial frame, V ¼ u þ vA ¼ const:, and the values of V ¼ 400 km s1 and nðr Þ ¼ 10 cm3 . (m0 is the permeability of vacuum.) We consider a hydrogen plasma, i.e., we take the solar–wind mass per electron to equal the proton mass, m=e ¼ mp . The assumption of the constancy of V is of ad hoc nature and made out of convenience. In reality, of course, there is no physical process in the corona that would force V to be constant. However, despite of their great simplicity our assumptions lead to a reasonable model of a slow solar–wind stream (Vainio et al., 2003; Vainio and Laitinen, 2007). The AWs are transported along the magnetic field lines at the constant group speed V following the principles of geometric optics. Their generation is governed by the exponential growth caused by the streaming of resonant ions. We use a simplified wave–particle resonance condition, k ¼ r1 L , where k and rL ¼ gv=ocp are the AW wavenumber and the proton Larmor radius, with g, v and ocp ¼ eB=mp denoting the proton Lorentz factor, speed and angular cyclotron frequency.1 Then the growth rate can be approximated by (Vainio, 2003) pSp ðr; tÞ . (1) nvA R þ1 Here S p ðr; tÞ ¼ 2pp2 v 1 mF ðr; p; m; tÞ dm is the (wave-frame) streaming per unit momentum of protons to be evaluated at the resonant (waveframe) momentum p ¼ mp Vf cp =f , F is the distribution function of protons in the wave frame, f cp ¼ ocp =ð2pÞ, and f is the wave frequency in the inertial frame. We also introduce a frequency diffusion term, mimicking the effect of wave–wave
sðf ; r; tÞ ¼ p2 f cp
1 The full quasi-linear resonance condition reads k ¼ 1=ðrL mÞ, implying that particles with a certain momentum p resonate with waves at a range of wavenumbers, jkjXeB=p.
ARTICLE IN PRESS R. Vainio, T. Laitinen / Journal of Atmospheric and Solar-Terrestrial Physics 70 (2008) 467–474
interactions, to the AW transport equation, DP~ qP~ qP~ q qP~ þV ¼ sP~ þ Dff , (2) Dt qt qr qf qf with P~ ðV 2 =BvA ÞP. The diffusion coefficient is set 2=3 to Dff ¼ ðV =r Þf 8=3 f b , which makes the power spectrum erode at high frequencies and tend, at 1 AU, toward the observed form (e.g., Horbury et al., 1996) of Kolmogorov spectrum of turbulence, P / f 5=3 , at frequencies greater than f b ¼ 1 mHz. In addition to AWs, our simulation keeps track of a large number of protons propagating in the flux tube. The AW growth rate s is self-consistently computed from the streaming of the simulated protons. Proton propagation is governed by the conservation of first adiabatic invariant, p2? =B, and momentum p in the inertial frame. These effects are represented by the equations of motion (in the inertial frame) r_ ¼ mv;
v ¼ const.
1 m2 q ln B . v m_ ¼ qr 2 In addition to the adiabatic propagation, protons are subject to the effect of resonant AWs, which scatter them (elastically in the local wave frame) at rate fPðf ; r; tÞ , (3) B2 where again the resonance condition f ¼ f cp mp V =p, with p measured in the wave frame, is applied. Thus, AWs and SEPs form a coupled, nonlinear system. The agent responsible for the acceleration of SEPs is a shock front postulated to expand radially at a constant speed V s . We only consider particle propagation upstream of the shock front and treat the front as a boundary condition in the equation of motion of the protons. In each interaction with the shock, particles returning back upstream have gained a small amount of energy (see below). Thus, protons scattered to the shock many times will gain energies up to hundreds of MeVs and beyond. For further details of the model, see Vainio and Laitinen (2007). nðr; p; tÞ ¼ p2 f cp
3. The numerical method We solve the transport equation of waves and the equation of motion of the particles simultaneously. We use a finite difference scheme to solve the wave transport equation on a Lagrangian grid (taking
469
care of wave propagation), which has a logarithmic spacing in the initial radial distance, i.e., xi ¼ ri Vt with xi ¼ r0 expðiDx Þ and a constant Dx 51. Also the frequency grid is logarithmic, f j ¼ f 0 expðjDf Þ with a constant Df . Wave growth is proportional to the net flux of particles integrated over the time step of wave evolution, Dt: ~ i ; f j ; t þ DtÞ P~ ij ðt þ DtÞ Pðx Z tþDt 0 0 ~ sðxi ; f j ; t Þ dt ¼ Pij ðtÞ exp ¼ P~ ij ðtÞ exp
Z t
t tþDt
pSp ðxi ; t0 Þ 0 dt , p f cp ðri Þ nðri ÞvA ðri Þ 2
where the integral can be computed by adding up the net number of particles crossing the point x ¼ xi from below during the time step from t to t þ Dt. Each particle contributes to the balance at two grid levels of frequency, f j and f jþ1 , just below and above the resonant frequency f ¼ f cp ðri ÞmV =p, respectively, with weights wi;j ¼
ln f ln f j p2 f cp ðri Þ W l, Df Aðri Þnðri ÞvA ðri ÞDf
wi;jþ1 ¼
ln f jþ1 ln f p2 f cp ðri Þ Wl, Aðri Þnðri ÞvA ðri ÞDf Df
where the sign is positive [negative] for protons crossing the spatial grid level xi from below [above], and W l is the weight of the particle in the simulation. We use AðrÞ ¼ r2 B =BðrÞ and, thus, consider a flux tube which has a solid angle of 1 sr at the solar surface. After the wave growth, we update the effect of diffusion to the wave spectrum at each spatial grid level by using a standard Crank– Nicholson scheme (Morton and Mayers, 1994). Particles are traced allowing each of them to have its own time step fulfilling dtl ¼ 2nl Dt, where nl 2 N0 and l numbers the particles. The time step for each particle is selected so that it tries to remain constrained by nl dtl 2 ½0:02; 0:05. However, a still smaller time step is required if the particle would otherwise be able to cross two spatial grid levels at once, i.e., dtl oDxi =vl is required. The value of nl can be decreased after each time step dtl whenever the conditions call for this, but it can be increased only in the beginning of each wavetransport cycle, Dt. Particle propagation, focusing caused by the conservation of first adiabatic invariant, and scattering by the AWs are solved using a leapfrog method, i.e., first moving the particle by half a time
ARTICLE IN PRESS R. Vainio, T. Laitinen / Journal of Atmospheric and Solar-Terrestrial Physics 70 (2008) 467–474
470
step, rl ¼ rl þ vl ml dtl =2, then performing the changes to the velocity and, finally, moving the particle by another half time step rl rl þ vl ml dtl =2. The effects of focusing and scattering on the pitch-angle cosine ml and speed vl are updated at the intermediate position, at the solar fixed and wave frame, respectively, as follows: focusing:
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 þ 4dðml þ dÞ , ml 2d vl dtl q ln B d¼ . qr r 2 1 þ
l
Scattering: ðvl ; ml Þ LTðvl ; ml ; V Þ, qffiffiffiffiffiffiffiffiffiffiffiffiffi ml ml cos W þ 1 m2l sin W cos j, ðvl ; ml Þ
LTðvl ; ml ; V Þ,
where LTðv; m; uÞ denotes a Lorentz transformation to a system moving at speed u along the magnetic field, and the scattering angles W and j are obtained using a random generator. (See, e.g., Kocharov et al., 1998; Vainio et al., 2000, for a detailed description of the scattering process.) Specifically, W2 and j are picked for each scattering from an exponential distribution with a mean of 2nl dtl and a uniform distribution between 0 and 2p, respectively, and the scattering frequency nl ¼ nðrl ; pl ; tÞ is obtained using the simulated wave spectrum in Eq. (3) with logarithmic interpolation. The method for simulating focusing is the implicit Euler scheme, which is stable and always gives a value inside the allowed range, m 2 ½1; 1. In case the particle crosses a spatial grid level xi of the wave spectrum, the time-integrated growth rate is updated there, as described above. We inject low-energy particles at a constant rate to the simulation at the inner boundary, taken to coincide with the propagating shock front. (Thus, the inner boundary of our simulation box is moving at speed V s V with respect to the wave grid.) The velocity distribution of the injected particles (in the shock frame) is dN=dv ¼ N inj Hðv u1 Þ= v1 expððv u1 Þ=v1 Þ, where u1 ¼ V s uðrs Þ is the local bulk speed of the upstream plasma in the frame of the shock and v1 is a model parameter which describes the spread of the suprathermal injection, set to v1 ¼ 375 km s1 in this paper. The total number of injected particles inside the con-
sidered flux tube is N inj . The pitch-angle cosines of the injected particles are drawn from a fluxweighted (see Jones and Ellison, 1991) isotropic distribution in the shock frame, i.e., PðmÞ ¼ 2jmj. Particle acceleration takes place at the shock front: every time a simulated particle hits the shock front from upstream, we compute its speed, v0 , in the rest-frame of the downstream scattering centers, moving with speed V 2 ¼ V 1 =Rsc with respect to the shock toward the Sun. Here, V 1 ¼ V s V is the speed of the AWs with respect to the shock in the upstream region and Rsc is the scattering-center compression ratio of the shock, which we regard as a constant model parameter.2 After that, we turn the particle back toward the upstream region drawing its downstream wave-frame pitch-angle cosine m0 from an isotropic distribution weighted by the shock frame parallel speed, i.e., Pðm0 Þ ¼ 2
jm0 þ V 2 =v0 j ; ð1 V 2 =v0 Þ2
m0 2 ½1; V 2 =v0 ,
m0 ¼ 1 corresponding to the direction of the shockframe flow (i.e., toward the Sun). However, as the downstream bulk flow is away from the shock, not all particles can make it back to the shock and, thus, we have to decrease the weight of the particles in each shock encounter by 0 v V2 2 Wl W l Pret ; Pret ¼ 0 , v þ V2 where Pret is the probability of return from the downstream region calculated assuming that the particle distribution around the shock is quasiisotropic in the rest frame of the downstream scattering centers3 (Jones and Ellison, 1991). We also assume that the returning particles spend a negligible amount of time in the downstream region, which is equivalent to assuming that the scattering mean free path is very small there. By considering the Lorentz transformation between the upstream and downstream wave frames, it is easy to see that the energy of the returning particle in the upstream wave frame is increased by DE ¼ DVp0 ðm01 m02 Þ, where DV ¼ V 1 V 2 and m01½2 is the pitch-angle cosine of the particle (in the 2 Note that in the case of AWs, Rsc is not identical to the gas compression ratio, R ¼ n2 =n1 , of the shock; Rsc can be larger than R in a parallel shock wave, even substantially (Vainio and Schlickeiser, 1998). 3 This assumption holds by construction for the particles returning from downstream, and rather well for the incident particles at energies of interest in our study.
ARTICLE IN PRESS R. Vainio, T. Laitinen / Journal of Atmospheric and Solar-Terrestrial Physics 70 (2008) 467–474
downstream wave frame), measured before [after] the interaction with the shock. Thus, particles gain a small amount of energy (hDEi 43 DVp0 ) in each interaction with the shock. By returning all particles from the downstream region and decreasing their weight accordingly, we keep the statistics of the simulated particles relatively uniform over the whole energy spectrum. 4. Results and discussion Vainio and Laitinen (2007) presented a parameter study, which concentrated on finding out the dependence of the achieved maximum proton energy on selected model parameters, namely N inj , Rsc and V s . They considered cases with a simulation time of tmax ¼ 640 s, initial shock radius of r0 ¼ 1:5r , and an ambient wave spectrum solved from Eq. (2) using s ¼ 0, the assumption of a steady state, and the inner boundary condition fPðf ; r0 Þ ¼ P B2 ðr0 Þ, with P ¼ 6 107 , corresponding to a mean free path of l0 ¼ 1r of 100-keV protons at r ¼ r0 . In this paper, we will take a closer look to three model runs with r0 ¼ 1:5r , V s ¼ 1500 km s1 and Rsc ¼ 4.
471
We will consider (A) N inj ¼ 3 1035 sr1 , tmax ¼ 640 s, P ¼ 6 107 ; (B) N inj ¼ 3 1035 sr1 , tmax ¼ 160 s, P ¼ 1:2 107 ; and (C) N inj ¼ 1:69
1036 sr1 , tmax ¼ 4320 s, P ¼ 2:0 108 . By studying these cases we aim at establishing: (i) how the particle mean free path in our simulation model compares with analytical estimates from diffusive shock acceleration (Bell, 1978); (ii) whether the ambient mean free path has any effect on the acceleration efficiency of the shock; (iii) if the injection rate has a direct effect on the acceleration efficiency or is the total number of injected particles the decisive parameter. In Fig. 1, we present the mean free path and the particle intensities of run A in the end of the simulation as a function energy (for rs r0 þ V s tmax 2:88r , rs þ r and rs þ 10r ) and position (for E ¼ 0:29, 2.4 and 20 MeV) along with the predictions from the one-dimensional analytical steady-state theory of Bell (1978) for the 2.4-MeV mean free path and intensity as a function of distance and the mean free path and intensity at the shock as a function of energy. Following Vainio and Laitinen (2007) we present the differential intensities, IðE; rÞ, scaled with the flux-tube crosssectional area and normalized to 1 AU, i.e.,
107 106 105 104 103 102 101 101 100 10−1 10−2 10−3
10−3
10−2
10−1
100
101
10−1
100
101
102
E [MeV] Fig. 1. Particle intensity (scaled to 1 AU) and scattering mean free path as a function of distance from the shock (left column) and energy (right column) and in the end of the model run A. In the left column, the black curves denote the simulated values at E ¼ 0:29 MeV (dashed), 2.4 MeV (solid), and 20 MeV (dash-dotted) and the gray curves denote the analytical estimates for 2.4 MeV. In the right column, the black curves denote the simulated values at z ¼ 0 (solid), 1r (dashed) and 10r (dash-dotted) and the gray curves denote the analytical estimates for z ¼ 0.
ARTICLE IN PRESS R. Vainio, T. Laitinen / Journal of Atmospheric and Solar-Terrestrial Physics 70 (2008) 467–474
I ðE; rÞ ðA=A ÞIðE; rÞ. Such a scaling is necessary when comparing with one-dimensional theory (where the flux tube has a constant cross-section). Assuming vanishing wave and particle intensities far upstream, the analytical one-dimensional steadystate result is (Bell, 1978; Vainio and Laitinen, 2007) z0 , z0 þ z
ð4Þ
Vs V ðz þ z0 Þ, v
ð5Þ
I ðE; zÞ ¼ I ðE; rs Þ lðz; pÞ ¼ 3
where z ¼ r rs is the distance from the shock, z0 ¼ hðp=p0 Þb3 , b ¼ 3Rsc =ðRsc 1Þ ¼ 4, p0 is the average injection momentum and 3b N inj =tmax b p I ðE; rs Þ , ð6Þ A V s 4pp p0 V s Aðrs Þnðrs Þ vA ðrs Þ . ð7Þ h p2 b1 N inj =tmax f cp ðrs Þ Although not a perfect fit, the analytical theory provides a valid approximation to the simulation results at energies below the spectral cutoff and at distances close to the shock. Note, however, that at the lowest energies the simulated mean free path is clearly higher than the predicted value, which is most likely due to the erosion of the spectrum by frequency diffusion (not included in the analytical theory) to the region without any resonant particles. Further out from the shock ðz\r Þ and at large energies ðE\10 MeVÞ the discrepancies are large, as expected, because of the effects of spherical geometry and finite acceleration time, and the advance of the frequency diffusion at larger heliocentric distances. One can also see the coherently propagating escaping particles in the upper left panel of Fig. 1. These particles are the ones that escape to the ambient medium before the waves have grown enough to turbulently trap them in front of the shock wave. As shown by Vainio and Laitinen (2007), the time integrated flux of this particle population is very well approximated by the theoretical prediction of Vainio (2003), who noted that in order for the waves to have grown appreciably in a given position and frequency, i.e. that growth rate becomes large, R 0 the time-integrated dt sðf ; r; t0 Þb1, the time-integrated net flux of particles has to exceed a threshold value, i.e., Z t nvA dt0 pS p ðr; t0 Þb 2 . p f cp 0
Before this occurs, the particles can easily escape from the system if the ambient conditions are not very turbulent, i.e., if the initial mean free path satisfies l0 ðr; pÞ\ðq ln B=qrÞ1 r=2, as in our simulations. Simulation run B has the same amount of injected particles as run A, but the ambient mean free path is 5 times larger, and the simulation time is 4 times shorter. (Thus, the injection rate is 4 times greater than in run A.) The scaled spectrum at the shock, I ðE; rs Þ, is very similar in the two cases: both runs show good coincidence with the theoretical powerlaw energy spectrum up to a cutoff energy, which is very similar in both cases, E\10 MeV. Thus, while the normalization of the energy spectrum at the shock is proportional to the injection rate, N inj =tmax , it is the total number of injected protons that determines the maximum energy in our model, at least in the simulated cases. This is actually the prediction of the diffusive shock acceleration theory as well, which yields for the cutoff momentum (Vainio and Laitinen, 2007) ðRsc 1Þ=3 1 M 1 pc ðtmax Þ p0 C N inj , (8) Rsc 1 where M ¼ V s =V and C is a constant, which has a value of o2 1033 sr for the conditions in our model corona. In Fig. 2 we present the proton intensities of run C at a distance of z ¼ 0:3 AU upstream from the shock. (Note that these are local intensities not scaled to 1 AU.) The run is now much longer and we can follow the prompt phase of the event as it would be seen by spacecraft far upstream. This run had 104 I [1 / (cm2 sr s MeV)]
472
103
102
101
0
20
40 t [min]
60
80
Fig. 2. Particle intensity at z ¼ 0:3 AU as a function of time in the simulation run C at four energies: 6 MeV (dashed), 10 MeV (solid), 30 MeV (dot-dashed) and 50 MeV (dotted).
ARTICLE IN PRESS R. Vainio, T. Laitinen / Journal of Atmospheric and Solar-Terrestrial Physics 70 (2008) 467–474
a very long ambient mean free path (30 times that of run A, i.e., several AUs at r ¼ 0:3 AU), which leads to a very efficient decoupling of the promptly escaping particles from the shock and these particles appear as an impulsive injection from the Sun although the acceleration process was all due to a coronal shock wave with a constant speed, compression ratio and injection rate of low-energy seed particles. The very large value of the ambient mean free path leads to a very fast decay of the intensity, and for a more realistic interplanetary mean free path the prompt component would have a slower increase to the time of maximum and slower decay. Our simulations do, however, imply that a remotely observed time-intensity profile that shows two maxima does not necessarily mean that there are two phases of acceleration at the Sun. One can also see from Fig. 2 that the lowest energy channels show the typical velocity dispersion, i.e., low-energy particle intensities have later onset than the high-energy particle intensities. However, above 10 MeV, our model produces the opposite effect. This is because the acceleration process, for the parameters of the simulation run, takes a long time and the first particles above 10 MeV appear clearly later than at lower energies. Thus, typical velocity dispersion results with highest energies having the earliest onset imply that the acceleration process is much faster than in the present simulation run. This can be achieved most easily by increasing the injection rate. At first sight, there would be nothing preventing the simulation model to generate protons at relativistic energies: all we seem to need is to increase the number of injected particles (see Eq. (8)). Our model, however, relies on the quasilinear theory of wave–particle interactions. Thus, we should make sure that the value of scattering frequency n stays well below f cp , implying (see Eq. (3)) that fPðf ; rÞ 5p2 0:1 B2 ðrÞ would have to be fulfilled. For the model parameters V s ¼ 1500 km s1 , Rsc ¼ 4 and tmax ¼ 640 s, the simulations imply that the limit of the validity of quasi-linear approximation will be reached roughly when pc mp c. Thus, our simulation model may need more physics before it can safely be used to describe particle acceleration to the relativistic energies.
473
5. Outlook and conclusions We have employed a number of simplifications to our model in order to achieve computational efficiency. However, now that the feasibility of the approach has been demonstrated, it is crucial to study the influences of the effects we have neglected so far. The solar wind model, while conveniently allowing us to use a Lagrangian grid with constant spacings for transporting the wave power, limits our modeling to speed profiles, which do not have a maximum of the Alfve´n velocity in the corona, often appearing at height of a few solar radii in both semiempirical and analytical modeling of the solar wind. Such a feature may prove to be important: as the shock crosses the region of maximum vA its Mach number gets temporarily decreased, but at the same time the wave growth becomes most efficient (Vainio, 2003; Vainio and Laitinen, 2007). In order to model these effects, it is probably most reasonable to abandon the Lagrangian grid approach, and use a simple simulation scheme, such as upwind differencing, to model the wave transport in more general solar wind models. Our modeling of the spectral diffusion was based on a simple linear diffusion coefficient. However, in a more realistic approach reproducing the effects of wave–wave interactions, non-linear diffusion approach would be needed, such as that described in Vainio et al. (2003). The particle modeling also suffers from simplifications which must be accounted for. First, treating the shock simply as lower boundary condition with a constant scattering-center compression ratio neglects important effects due to wave transmission/reflection at the shock (Vainio and Schlickeiser, 1998). To account for these, particles and waves should be self-consistently traced in the downstream region, including modes propagating both toward and away from the Sun. Second, the simplified resonance condition used in these simulations essentially maps the particle to one resonant wave number, while the exact, pitch-angle dependent resonance condition spreads the resonance to a range of wavenumbers. This will influence the evolution of the wave spectra and, thus, also the spectrum of escaping particles: since fast particles can now generate waves that scatter the slower ones, the escaping spectrum will probably also depend on the temporal evolution rate of the maximum energy in the spectrum. Using the full quasi-linear
ARTICLE IN PRESS 474
R. Vainio, T. Laitinen / Journal of Atmospheric and Solar-Terrestrial Physics 70 (2008) 467–474
resonance condition, k ¼ 1=ðrL mÞ, will, however, require two major changes in the model: (i) the scattering frequency has to be taken as m-dependent, which typically requires much shorter time steps than an isotropic scattering process (Kocharov et al., 1998) and thus slows down the computational efficiency of the simulations by a considerable factor; and (ii) the AW growth rate has to be evaluated using the full quasi-linear theory (see, e.g., Vainio, 2003), rendering the simple book-keeping method of particles crossing grid cells invalid. Fortunately a similar method may be devised for the full resonance condition, since the growth rate in this case becomes proportional to the net flux of particles in momentum space, i.e., to Dmm qF =qm, where Dmm is the pitch-angle diffusion coefficient, evaluated at the resonant parallel momentum pk ¼ mocp =k. Thus, we may compute the timeintegrated growth rate now by a similar bookkeeping of particles crossing the resonant frequency levels f i through scattering. However, implementing these effects to the model will be considerably more tedious than the modifications to the solar–wind and wave-transport models. In conclusion, we have presented a numerical method to simulate particle acceleration and coupled AW evolution in CME-driven shocks in the solar corona. Our model is capable of accelerating protons to 100 MeV and beyond in a few minutes. The model is based on a simplified version of quasi-linear theory, and thus its range of validity is constrained by the requirement of linear wave amplitudes. This indicates that although relativistic energies (p mc) may be reached in the model within 10 min, acceleration faster than that may require more complicated approach on wave–particle and wave–wave interactions than employed in our model. The model produces a coherently escaping particle population with a hard energy spectrum ð/ p1 Þ, but this feature of the model may be especially sensitive to the simplifications made on modeling the wave–particle interactions. Further modifications to the model were outlined, by which these and other remaining simplifications of the model may be removed.
Acknowledgments We thank the Finnish Meteorological Institute for the use of their computing facilities. References Axford, W.I., Leer, E., Skadron, G., 1977. The acceleration of cosmic rays by shock waves. In: Proceedings of the 15th International Cosmic Ray Conference, vol. 11, pp. 132–137. Bell, A.R., 1978. The acceleration of cosmic rays in shock fronts. I. Monthly Notices of the Royal Astronomical Society 182, 147–156. Blandford, R.D., Ostriker, J.P., 1978. Particle acceleration by astrophysical shocks. Astrophysical Journal 221, L29–L32. Horbury, T.S., Balogh, A., Forsyth, R.J., Smith, E., 1996. The rate of turbulent evolution over the Sun’s poles. Astronomy and Astrophysics 316, 333–341. Jones, F.C., Ellison, D.C., 1991. The plasma physics of shock acceleration. Space Science Reviews 58, 259–346. Kocharov, L., Vainio, R., Kovaltsov, G.A., Torsti, J., 1998. Adiabatic deceleration of solar energetic particles as deduced from Monte Carlo simulations of interplanetary transport. Solar Physics 182, 195–215. Krymsky, G.F., 1977. A regular mechanism for the acceleration of charged particles on the front of a shock wave. Doklady Academie Nauk SSSR 243, 1306–1308 (Translation: Soviet Physics—Doklady 22, 327–328.). Lee, M.A., 2005. Coupled hydromagnetic wave excitation and ion acceleration at an evolving coronal/interplanetary shock. Astrophysical Journal Supplement 158, 38–67. Morton, K.W., Mayers, D.F., 1994. Numerical Solution of Partial Differential Equations. Cambridge University Press, Cambridge, UK. Rice, W.K.M., Zank, G.P., Li, G., 2003. Particle acceleration and coronal mass ejection driven shocks: shocks of arbitrary strength. Journal of Geophysical Research 108, CiteID 1369, 10.1029/2002JA009756. Vainio, R., 2003. On the generation of Alfve´n waves by solar energetic particles. Astronomy and Astrophysics 406, 735–740. Vainio, R., Laitinen, T., 2007. Monte Carlo simulations of coronal diffusive shock acceleration in self-generated turbulence. Astrophysical Journal 658, 622–630. Vainio, R., Schlickeiser, R., 1998. Alfve´n wave transmission and particle acceleration in parallel shock waves. Astronomy and Astrophysics 331, 793–799. Vainio, R., Kocharov, L., Laitinen, T., 2000. Interplanetary and interacting protons accelerated in a parallel shock wave. Astrophysical Journal 528, 1015–1025. Vainio, R., Laitinen, T., Fichtner, H., 2003. A simple analytical expression for the power spectrum of cascading Alfve´n waves in the solar wind. Astronomy and Astrophysics 407, 713–723.