Quasi-adiabatic and stochastic heating and particle acceleration at quasi-perpendicular shocks
aa r X i v : . [ phy s i c s . s p ace - ph ] S e p Draft version September 15, 2020
Typeset using L A TEX twocolumn style in AASTeX63
Quasi-adiabatic and stochastic heating and particle acceleration at quasi-perpendicular shocks
Krzysztof Stasiewicz
1, 2 and Bengt Eliasson Department of Physics and Astronomy, University of Zielona G´ora, Poland Space Research Centre, Polish Academy of Sciences, Warsaw, Poland SUPA, Department of Physics, University of Strathclyde, Glasgow, G4 0NG, United Kingdom (Accepted for publication in the Astrophysical Journal, September 10, 2020)
ABSTRACTBased on Magnetospheric Multiscale (MMS) observations from the Earth’s bow shock, we haveidentified two plasma heating processes that operate at quasi-perpendicular shocks. Ions are subjectto stochastic heating in a process controlled by the heating function χ j = m j q − j B − div( E ⊥ ) forparticles with mass m j and charge q j in the electric and magnetic fields E and B . Test particlesimulations are employed to identify the parameter ranges for bulk heating and stochastic accelerationof particles in the tail of the distribution function. The simulation results are used to show that ionheating and acceleration in the studied bow shock crossings is accomplished by waves at frequencies(1–10) f cp (proton gyrofrequency) for the bulk heating, and f > f cp for the tail acceleration. Whenelectrons are not in the stochastic heating regime, | χ e | <
1, they undergo a quasi-adiabatic heatingprocess characterized by the isotropic temperature relation
T /B = ( T /B )( B /B ) / . This is obtainedwhen the energy gain from the conservation of the magnetic moment is redistributed to the parallelenergy component through the scattering by waves. The results reported in this paper may also beapplicable to particle heating and acceleration at astrophysical shocks. Keywords: acceleration of particles – shock waves – solar wind – turbulence – chaos INTRODUCTIONElectron and ion acceleration and heating at col-lisionless shocks is an important problem in as-trophysics and space physics, which has been ad-dressed over many years by a number of au-thors (Bell 1978; Lee & Fisk 1982; Wu et al. 1984;Goodrich & Scudder 1984; Blandford & Eichler1987; Balikhin & Gedalin 1994; Gedalin et al. 1995;Treumann 2009; Burgess et al. 2012; See et al. 2013;Mozer & Sundqvist 2013; Krasnoselskikh et al. 2013;Guo et al. 2014; Wilson III et al. 2014; Park et al. 2015;Cohen et al. 2019; Xu et al. 2020). In the above citedpublications one can find a long list of waves, instabil-ities, and processes that could play a role in ion andelectron heating/acceleration – however, the question ofthe exact heating mechanisms working at shocks is stillin an inconclusive state. [email protected]@strath.ac.uk
New generation of high-resolution space instrumentswith simultaneous 4-spacecraft measurements on theMagnetospheric Multiscale (MMS) mission (Burch et al.2016) opened unprecedented possibility for testing heat-ing and acceleration mechanisms that operate at colli-sionless shocks in reality, and not only in theory. Forexample, MMS offers the capability of computing gra-dients of plasma parameters and of electric and mag-netic fields on spacecraft separation distances of ∼ Shock compression of the density N and the mag-netic field B → diamagnetic current → lower hybrid K. Stasiewicz & B. Eliasson drift (LHD) instability → electron E × B drift → elec-tron cyclotron drift (ECD) instability → heating: quasi-adiabatic ( χ j < ), stochastic ( χ j > ) .Stochastic heating is a single particle mechanismwhere large electric field E gradients due to spacecharges destabilize individual particle motions in a mag-netic field B , rendering the trajectories chaotic in thesense of a positive Lyapunov exponent for initiallynearby states. The stochastic heating function of parti-cle species j ( j = e for electrons and p for protons) is(Stasiewicz 2020a) χ j ( t, r ) = m j q j B div( E ⊥ ) (1)where m j and q j are the particle mass and charge. Theparallel (to the magnetic field) electric field E k is hereexcluded since it does not directly contribute to thestochasticity, leaving only the perpendicular field E ⊥ in(1). Stochastic heating typically occurs when | χ j | & | χ j | < χ j can beregarded as a measure of demagnetization. Particles aremagnetized (adiabatic) for | χ j | <
1, and demagnetized(subject to non-adiabatic heating) for | χ j | & χ p is typi-cally in the range 10 −
100 in the bow shock and themagnetosheath, which implies that the ions are stronglydemagnetized and can be subject to stochastic heatingprocesses in these regions. In Section 4, test-particlesimulations are carried out for a range of parametersprimarily relevant for stochastic heating of protons.In order to also demagnetize and stochastically heatelectrons we need χ e > χ p > m p /m e =1836, which requires either very strong E -gradients orlow B -fields, or both, as implied by Eq. (1). Elec-tron heating at perpendicular shocks based on χ e (withthe divergence reduced to ∂E x /∂x ) has been referredto as the kinematic mechanism (Balikhin et al. 1993;Gedalin et al. 1995; See et al. 2013). The required gra-dient of the electric field is associated with a macroscopicelectric field in the direction normal to the shock. Un-fortunately, in perpendicular shocks the observed thick-ness of the shock ramp and measured values of the nor-mal electric field do not allow χ e > χ e are mostly below the stochastic threshold for elec-trons, because of the increasing values of B ≈ − χ e ∝ B − . In Section 3, we demonstrate that such situa-tions instead lead to quasi-adiabatic electron heating ,characterized by electron heating on the compression ofthe magnetic field, combined with scattering by waves,leading to the isotropic temperature relation T /B =( T /B )( B /B ) / . This is to our knowledge a novelconcept identified and explained for the first time in Sec-tion 3. MULTIPLE CROSSINGS AND WAVES INSHOCKSWe analyse recent MMS measurements from 2020-01-03 obtained by the 3-axis electric field (Lindqvist et al.2016; Ergun et al. 2016; Torbert et al. 2016) and mag-netic field vectors measured by the Fluxgate Magne-tometer (Russell et al. 2016), and the number density,velocity, and temperature of both ions and electronsfrom the Fast Plasma Investigation (Pollock et al. 2016).Figure 1 shows multiple crossings of shocks caused bythe oscillatory movements of the bow shock with an am-plitude of 6–10 km/s estimated from the time shifts ofthe density signals. The speed is with respect to theMMS spacecraft moving at 1.9 km/s earthward. Thespacecraft position at time 14:30 was (13.5, 10.3, -1.8)R E GSE (geocentric solar ecliptic) coordinates, and theaverage inter-spacecraft distance was 22 km (minimum10 km). Shown in Fig. 1 are: the electron number den-sity N , the magnetic field B , the ion and electron tem-peratures, and the ratio T ⊥ /B (not to scale). Notably inFig. 1c, the parallel and perpendicular electron tempera-tures are almost equal, indicating that an isotropizationprocess takes place.Complementary plasma parameters for the shockcrossings are displayed in Fig. 2, which shows (a) the an-gle between the magnetic field vector and the geocentricradial direction to the spacecraft (a proxy for the shocknormal), (b) the Alfv´en Mach number, M A = V i /V A ,i.e, the ratio between the ion bulk speed and the Alfv´enspeed, and (c) the ion and electron plasma beta, i.e.,the ratio between the thermal energy density of parti-cles and the magnetic field energy density. It is seenthat all shocks have quasi-perpendicular configurations,the plasma beta is β i ∼ , β e ∼
1, and the Alfv´en Machnumber 6-8, outside the shocks.The ion and electron temperatures indicate rapidheating at the shock ramp, and the repetitive events uasi-adiabatic and stochastic heating N [ c m - ] (a) B [ n T ] (b) T i [ e V ] || (c)T /B 14:30 15:00 15:30 16:00 T e [ e V ] || (d)T /B Figure 1.
A series of shock crossings caused by the os-cillatory movement of the bow shock. Panel (a) shows theelectron number density N , panel (b) the magnetic field B ,panel (c) shows T i ⊥ , T i k , and the ratio T i ⊥ /B (not to scale),and panel (d) shows the same parameters as (c) but for elec-trons. Note the different behaviors of the ratio T ⊥ /B for theions and electrons in the shock ramps, with humps for theions and dips for the electrons. ( B , R ) deg (a) M A (b) 14:30 15:00 15:30 16:00 be t a i e (c) Figure 2.
Complementary plasma parameters for Fig. 1: (a)the angle between the magnetic field vector and the geocen-tric radial direction (a proxy to the shock normal direction),(b) the Alfv´en Mach number, and (c) ion and electron β parameter. offer a great opportunity to study heating processes op-erating at quasi-perpendicular shocks. The ratio T ⊥ /B derived from measurements is an excellent indicator ofthe heating processes. A flat ratio across the shockwould indicate adiabatic perpendicular heating comingfrom the conservation of the magnetic moment. Thisshould be accompanied by unchanged parallel tempera- Figure 3.
Time-frequency spectrogram of the perpendicularelectric field for the first shock ramp in Fig. 1. Over-plottedare the proton cyclotron frequency f cp , the lower hybrid fre-quency f lh , the electron and ion temperatures (eV), and theelectron cyclotron frequency f ce . Waves between f cp − f lh are attributed to the LHD and above f lh to the MTS (modi-fied two-stream) instabilities, and for frequencies around f ce and above to the ECD instability. Note the vertical stria-tions that start from below 1 Hz (LHD instability) and gothrough the MTS and ECD instabilities up to 4 kHz, indi-cating co-location and common origin of these instabilities. ture. We see that the ion ratio T i ⊥ /B has humps andthe parallel temperature is smaller, which is indicativefor non-adiabatic perpendicular heating and less efficientparallel heating. The electron ratio T e ⊥ /B instead hasdips and the temperature is nearly isotropic.Both particle species in Figure 1 manifest non-adiabatic behavior but of a different character. To studythe processes and wave modes involved in the heating weshow in Fig. 3 the time-frequency power spectrogram ofthe perpendicular electric field sampled at the rate 8192s − for the first shock crossing. Over-plotted are the pro-ton cyclotron frequency f cp , the lower hybrid frequency f lh , the electron temperature T e ⊥ , the ion temperature T i ⊥ , and the electron cyclotron frequency f ce . Near thepeak of the shock at around 14:24 UT, the mean valuesof the proton plasma frequency is f pp ≈ f pe ≈
42 kHz. The lower-hybrid frequency is f lh = [ f − pp + ( f cp f ce ) − ] − / ≈ f lh are related to the LHD instabil-ity, which has maximum growth rate at k ⊥ r e ∼ K. Stasiewicz & B. Eliasson lation results of Daughton (2003) indicate that they havelonger wavelengths k ⊥ ( r e r p ) / ∼ k ⊥ is thewavenumber perpendicular to the magnetic field, and r e = v T e /ω ce is the electron thermal Larmor radius,and v T e = ( T e ⊥ /m e ) / is the electron thermal speed, r p is the proton thermal Larmor radius. For r e ≈ V i ⊥ ≈
250 km/s with re-spect to the spacecraft, the LHD waves at frequencies f < f lh and wavelengths λ M ∼ πr e corresponding tothe maximum growth rate will be upshifted in frequencyby ∆ f = V i ⊥ /λ M and observed up to 50 Hz in Figure 3.Notably, there are other wave modes in this frequencyrange, as for example magnetosonic whistlers at ∼ f ce and aboveare associated with the ECD instability and coupledwith ion-acoustic (IA) waves (Wilson III et al. 2010;Breneman et al. 2013; Goodrich et al. 2018; Stasiewicz2020a). Because of short wavelengths ∼ r e they couldbe Doppler downshifted by up to hundreds Hz.The LHD and ECD instabilities are cross-field currentdriven instabilities caused by relative electron-ion drifts.Diamagnetic drift of protons V d = T p ( m p ω cp L N ) − = v T p ( r p /L N ) caused by gradients of the density leads tothe LHD instability when the ratio between the scaleof the density gradient L N = N |∇ N | − and the pro-ton thermal gyroradius r p obeys the condition L N /r p < ( m p /m e ) / (Huba et al. 1978; Drake et al. 1983). Thedensity compression starts at the foot of the shock andis strongly amplified in the ramp, which can be seen inFigure 3 in the profile f lh , which is proportional to B ,but also representative for N .In the nonlinear stage, the LHD waves produce largeamplitude electric fields resulting in efficient E × B drifts of the electrons. Due to the large ion gyrora-dius compared to the wavelength of the LHD waves,the ions do not experience significant E × B , andhence there is a net current set up by the electronsonly. When the differential drift velocity betweenthe electrons and ions exceeds the ion thermal ve-locity, the modified two-stream (MTS) instability cantake place resulting in waves at frequencies above f lh (Krall & Liewer 1971; Lashmore-Davies & Martin 1973;Gary et al. 1987; Umeda et al. 2014). Below, we willnot differentiate between the MTS and LHD instabil-ities since they belong to the same dispersion surface(Silveira et al. 2002; Yoon & Lui 2004), but we will usethe term LHD instability in the sense of a generalizedcross-field current driven instability in the lower-hybridfrequency range. When the relative electron-ion drift speed becomesa significant fraction of the electron thermal speed, V E = | E × B | /B ∼ v T e , the ECD instability is ini-tiated, which creates even larger electric fields on spa-tial scales of r e . (Forslund et al. 1972; Lashmore-Davies1971; Muschietti & Lemb´ege 2013). The ECD instabil-ity takes place near cyclotron resonances ( ω − k ⊥ V de ) = nω ce , where ω ce = eB/m e is the angular electron cy-clotron frequency, k ⊥ = 2 π/λ is the perpendicularwave number, λ is the wavelength, and n is an inte-ger (Janhunen et al. 2018). Here, V de ≈ V E is theperpendicular electron drift velocity in the rest frameof the ions. This resonance condition can be written k ⊥ V E ≈ nω ce and expressed as k ⊥ r e ≈ nv T e /V E . For r e = 1 km and n = 1 their wavelengths are λ ≈ πr e n V E v T e ≈ . V E v T e . (2)This means that the ECD waves with n = 1 and electricdrift velocities V E > v T e (see Figure 4b) have wave-lengths large enough to enable correct gradient compu-tations in the calculations of div( E ) by the MMS space-craft constellation. Contribution of shorter waves with n > χ may be underestimated bythe gradient computation procedure. The ECD wavesresonate/couple with structures created by the LHD in-stability, k ⊥ r e ∼
1, when nv T e /V E = 1. The n =1 ECDmode can be naturally excited in drift channels createdby the LHD instability when V E = v T e . There is smoothtransition and co-location of LHD and ECD waves, seenin Fig. 3 which is possibly related to the matching con-dition between these two instabilities.The scale of the density gradient L N shown in Fig. 4ais computed directly from 4-point measurements usingthe method of Harvey (1998). As a verification, we showalso the gradient scale L B = B |∇ B | − for the magneticfield. They coincide in the shock proper, as expected forfast magnetosonic structures.In the pioneering work on the LHD instability,Krall & Liewer (1971) used an expression for the elec-tron drift in the form V de ∝ ( N − ∂ x N − B − ∂ x B + .... ),which implied that the current due to the ∇ B termfrom fluid integration would cancel the diamagnetic cur-rent due to the ∇ N drift, in case when the gradientscale lengths are the same ( L N = L B ) and both gradi-ents point in the same direction. Because these scalesare the same at the bow shock, several influential au-thors (Lemons & Gary 1978; Zhou et al. 1983; Wu et al.1984), and many others, claimed that there would beno LHD instability at the bow shock. This erroneousconclusion affected many researchers afterwards and hasled to a 40-year delay in the identification of the LHD uasi-adiabatic and stochastic heating Figure 4.
Diagnostic parameters for the case in Fig. 3: (a)The gradient scales of the plasma density L N (blue) and ofthe magnetic field L B (red) are normalized with thermal pro-ton gyroradius r p (mean = 91 km). Regions with L N /r p . E × B drift speed V E (blue), the electron thermal speed v Te (red), and perpendicular ion speed V i ⊥ (magenta). Regionswith V E ∼ v Te indicate presence of the ECD instability. (c)The stochastic heating function χ p for protons derived fromthe data with Eq. (1) for the electric field 0.25–512 Hz. instability as the prevailing ion heating mechanism incompressional shock waves, and as a possible trigger forthe ECD instability (Stasiewicz 2020a,b). As a mat-ter of fact, in a homogeneous plasma with a spatiallyvarying magnetic field, the ∇ B drift cancels with otherterms due to the gyration of particles in the magneticfield, and therefore does not contribute to macroscopiccurrents as explained in section 7.4 of the textbook byGoldston & Rutherford (1995). Thus, the diamagneticdrift current is not canceled by the magnetic gradientdrift term, and the LHD instability can be excited atthe bow shock.Figure 4b shows the computed E × B drift speed whichis increased and comparable to the electron thermalspeed in the ion and electron heating regions in Fig. 3.The drift velocity was computed in the frequency range0 - 512 Hz, because for frequencies larger than the elec-tron cyclotron frequency the drift approximation is notvalid. For comparison we also show the plot of the mea-sured perpendicular ion speed V i ⊥ in magenta color. Ithas the same value as the computed E × B drift speed inthe solar wind up to 14:23, but deviates strongly insidethe shock. Large difference between the electron drift V E and the measured perpendicular drift of ions V i ⊥ wouldinduce sequentially the LHD and ECD instabilities asmentioned in the Introduction. The diagnostic param-eters support the interpretation of the waves shown in Fig. 3 as caused by the LHD and ECD instabilities,and that these waves are spatially co-located, which isseen as striations extending over the whole frequencyspectrum.Stochastic heating is controlled by the stochasticheating function (1) computed directly from 4-pointmeasurements using the method of Harvey (1998).Goodrich et al. (2018) raised concerns that the axialdouble probe (ADP) on MMS, which uses rigid axialbooms shorter than the wire booms of the spin-planedouble probe experiment (SDP), produces a larger am-plitude response for short, tens of meter (Debye length)waves such as the ion-acoustic (IA) wave. This instru-mental difference may affect the computations of thedivergence of the electric field and the resulting valueof χ . Therefore, the computations of div( E ) are madein the despun spacecraft coordinates (DSL), which sep-arates E z provided by the ADP, from ( E x , E y ) providedby the SDP. This enables removal of the highest fre-quency components above f ce , which may contain suchshort waves around the ion plasma frequency f pp , fromthe analysis. We have therefore removed the highestfrequency components before computing the gradients,and χ p shown in Fig. 4c is computed for the frequencyrange 0.25–512 Hz. Frequencies below 0.25 Hz are re-moved to avoid spurious effects at the spin frequencyand its harmonics. Not using the ADP E z componentat all produces χ ca 20% smaller.The computed χ p shown in Fig. 4c indicates that theions are demagnetized and likely to undergo stochasticheating as seen in detail in Fig. 3. On the other handthe value of χ p ∼
100 corresponds to χ e ∼ .
06, which istoo small to demagnetize the electrons and subject themto stochastic heating. The computed contribution to χ e from short n > χ are the same as in measurements of theelectric field, i.e., 1 mV/m, or ca 10% for large amplitudefields (Lindqvist et al. 2016).Figure 3 is representative for all 9 shock crossingsshown in Fig. 1. It shows that the ion heating is ob-served at the foot of the shock, earlier than the elec-tron heating, and correlates well with the power of LHDwaves. Electron temperature correlates well with the in-creased wave activity in the LH/ECD frequency rangeand maximizes in the region of the most intense ECDwaves around 14:24. On the other hand, it appears tocorrelate also with the magnetic field strength, repre-sented here by the electron gyrofrequency f ce = ω ce / π ,which suggests that an adiabatic behavior T e ⊥ ∝ B should be also considered here. However, this appar- K. Stasiewicz & B. Eliasson T / B [ e V / n T ] model measured Figure 5.
Comparison of the measured ratio T e ⊥ /B (red)and the modeled dependence expressed by Eq. (4) (blue) forthe second event in Fig. 1. It shows excellent agreement inthe ramp of the shock, where (4) is applicable. ent correlation is not exact, as seen in the T e ⊥ /B ratioshown in Fig. 1, with humps for the ions and dips for theelectrons at the shock ramps. This will be explained inthe next section as quasi-adiabatic electron heating in-volving the compression of the magnetic field combinedwith isotropization by the scattering on waves. QUASI-ADIABATIC ELECTRON HEATINGThe computed value of the heating function (1) shownin Fig. 4c indicates that the stochastic heating may notbe available for electrons in the analyzed shock cross-ings. The behavior of the ratio T e ⊥ /B and the isotropyof the electron temperature, discussed in Section 2, sug-gests a different kind of heating process. Let us assumethat the electrons obtain perpendicular energy from theconservation of the magnetic moment (they are mag-netized, consistent with χ e < T ⊥ /B = const , the differential temperature increase is dT ⊥ = T ⊥ B − dB . If the energy gain from 2 degrees offreedom (2 dT ⊥ ) is redistributed by pitch angle scatter-ing to 3 degrees of freedom (3 dT ) the conservation ofenergy implies 3 dT = 2 T B − dB (3)for T = T ⊥ = T k . This can be easily integrated to give TB = T B (cid:18) B B (cid:19) / (4)which predicts a dip of T /B where B has a maximum.Figure 5 shows a detailed comparison of Eq. (4)with the measured ratio for the second shock struc-tures of Fig. 1. The model is in excellent agreementwith measurements, which supports the validity of theabove described type of non-adiabatic heating, hence-forth referred to as quasi-adiabatic heating . All shock crossing in Fig. 1 show similar signatures of electronquasi-adiabatic heating with χ e <
1. The agreement be-tween the two curves in Fig. 5 indicates an outstandingquality of the particle measurements by the FPI instru-ment (Pollock et al. 2016), which is able to reproducethe subtle effects of Eq. (4) at sharp gradients of theshock ramps seen in Fig. 1.It comes as a surprise from this analysis that thestrong ECD waves (Fig. 3) with electric field ampli-tudes of ∼
150 mV/m and short wavelengths of ∼ r e do not directly heat electrons. Such expectations havealso been expressed by Mozer & Sundqvist (2013), whonoted that the wave potential of the ECD waves signif-icantly exceeds the thermal energy of the electrons, sothat some amount of heating would be anticipated. Theelectron reluctance to stochastic heating appears to berelated to the dependence χ e ∝ B − , and high values of B in the shock ramp, which keeps χ e < ∼ − SIMULATIONS OF NON-ADIABATICSTOCHASTIC HEATINGFor a possible stochastic heating of particle speciesof mass m , charge q we have available wave frequenciesfrom dc to 4096 Hz shown in Fig. 3, and spatial scalesranging from above ∼ r e ∼ B is in the z direction, and a macroscopicconvection electric field E y drives particles into an elec-trostatic wave with amplitude E x propagating in the x -direction. We keep the magnetic field constant to sep-arate purely stochastic heating from the quasi-adiabaticheating discussed above. Thus, in the Doppler frame ofthe satellite, the drifting plasma is characterized by theconvecting electric field and the time-dependent waveelectric field. The governing equations are m dv x dt = qE x cos( k x x − ωt ) + qv y B , (5) m dv y dt = qE y − qv x B , (6) dxdt = v x , dydt = v y . (7) uasi-adiabatic and stochastic heating v x = e V x + ωk x , x = e X + ωk x t (8)we obtain the system m d e V x dt = qE x cos( k x e X ) + qv y B , (9) m dv y dt = q e E y − q e V x B , (10) d e Xdt = e V x , dydt = v y , (11)where the shifted convection electric field is e E y = E y − ωk x B . (12)In this frame, the electric field is time-independent andgoverned by the electrostatic potential e Φ( e X, y ) = − E x k x sin( k x e X ) − e E y y. (13)On the other hand, by a change of frame into that ofthe E × B -drift velocity, v x = V x + E y B , x = X + E y B t (14)we obtain instead m dV x dt = qE x cos( k x X − e ωt ) + qv y B , (15) m dv y dt = − qV x B , (16) dXdt = V x , dydt = v y , (17)where the convecting electric field has been eliminatedand absorbed into the frequency in the plasma frame e ω = ω − k x E y B . (18)These two different approaches indicate that the modelof waves in the plasma frame with the wave frequency (18) that absorbs the convection field is equivalent to thestatic wave structures superposed with the convection(Stasiewicz 2007).Without loss of generality we choose to simulate thesystem (15)–(17). A suitable normalization of variables(Karney 1979; Fukuyama et al. 1977; McChesney et al.1987) with time normalized by ω − c , space by k − x andvelocity by ω c /k x with ω c = qB /m being the angularcyclotron frequency, gives the system of dimensionless,primed variables, dv ′ x dt ′ = χ cos( x ′ − Ω t ′ ) + v ′ y , (19) dv ′ y dt ′ = − v ′ x , (20) dx ′ dt ′ = v ′ x , dy ′ dt ′ = v ′ y , (21)in which there are only two parameters, the normalizedwave frequency in the plasma frame,Ω = e ωω c , (22)and the stochastic heating parameter, equivalent to (1), χ = mk x E qB = k x ω c E B , (23)which represents the normalized wave amplitude. Animportant third parameter is the initial velocity ofthe particles, since stochastic motion takes place onlyin restricted regions in phase space (Karney 1979;Fukuyama et al. 1977; McChesney et al. 1987). For astatistical description of the particles, the initial con-dition can be described by a Maxwellian distributionfunction F = N πv T exp (cid:18) − ( v x + v y )2 v T (cid:19) , (24)where v T = ( T /m ) / is the initial thermal speed and T is the initial temperature. In the normalized variableswith F = N ( k x /ω c ) F ′ it is written F ′ = 12 πv ′ x exp (cid:18) − ( v ′ x + v ′ y )2 v ′ x (cid:19) (25)where v ′ x = k x r c is the normalized thermal speed and r c = v T /ω c is the thermal Larmor radius. The valueof v ′ x determines the initial temperature in the velocitydistribution function, which due to the normalization, isin fact proportional to the ratio of the gyroradius to thewavelength λ = 2 π/k x .We carry out a set of test particle simulations for M = 10 000 particles, which are Maxwell distributed K. Stasiewicz & B. Eliasson
Figure 6.
A color plot of stochastic heating showing the difference T ′ − T ′ between the normalized kinetic temperature T ′ = k x T /mω c at the end of the simulation and the initial value T ′ = ( v ′ x ) = k x T /mω c after 3 cyclotron periods for chargedparticles in an electrostatic wave with normalized amplitude χ = 60. Here f c = ω c / π , T ′ = v ′ x is the normalized initialtemperature and v ′ x = k x v T /ω c with the thermal speed v T = ( T /m ) / . The insets show distribution functions in ( x , v x )space at t = 0 and 3 f − c for different values of Ω and v ′ x . Bulk heating takes place for Ω .
10, while for Ω &
10 there issignificant heating only for thermal velocity comparable to the phase velocity, or v ′ x ∼ Ω in the normalized variables (dashedline) leading to a distribution function having a high energy tail of particles. in velocity and uniformly distributed in space. The sys-tem (19)–(21) is advanced in time using a St¨ormer-Verletscheme (Press et al. 2007). The input variables for thesimulations are: the normalized wave frequency Ω in therange 10 − to 10 , and the initial normalized thermal ve-locity v ′ x = k x r c spanning 10 − to 10 . The normalizedamplitude of the electrostatic wave is set to χ = 60, con-sistent with the observations in Fig. 4c. The simulationis run for a relatively short time of 3 cyclotron periodsof the particles, motivated by the observations of rapidion heating within a few cyclotron periods, see Fig. 3.The kinetic temperatures resulting from the stochasticheating is calculated as T = 12 M M X k =1 m ( v x,k + v y,k ) . (26)Simulations are carried out for different values of Ωand v ′ x to produce the color plot in Fig. 6, which showsthe difference T ′ − T ′ between the normalized kinetictemperature T ′ = k x T /mω c at the end of the simulationand the initial value T ′ = ( v ′ x ) = k x T /mω c .The most interesting regions are the ones with redcolor, representing a temperature increase of ∼ (20 − mω c /k x . The frequency region 1 . Ω .
10 repre-sents bulk heating where the cold population is signifi-cantly heated to a temperature of ∼ (20 − mω c /k x .The bulk heating region would expand to higher val-ues of Ω for larger χ as well as to lower Ω for longertimes. The inset plots for Ω = 3 . v ′ x = 10 − and 1 show that theparticles are bulk heated and spread almost uniformlyin velocity space up to a maximum speed ∼ ω c /k x ,and the distributions achieve a kinetic temperature of ∼ mω c /k x after 3 cyclotron periods. This is rel-evant for the heating of protons by the low-frequencywaves observed in Fig. 3. Somewhat similar cases butfor Ω < χ ∼ mω c /k x after 3 cyclotron periods.On the other hand, for Ω significantly larger than 10only particles with a high enough initial thermal veloc-ity comparable to the phase velocity, or v ′ x ∼ Ω in thenormalized variables (dashed line in Fig. 6) are furtheraccelerated, leading a warm component with extendedenergy tails in the distribution function. Such cases of uasi-adiabatic and stochastic heating v ′ x = 10 andΩ = 20 the normalized temperature increases a factor2 within 3 cyclotron periods, which may be relevant forwaves below the lower hybrid frequency seen in Fig. 3.Below a threshold initial temperature, the distributionis not affected by the wave, and there is a gap in theheating for low initial temperatures, seen in the lowerright blue-colored region of Fig. 6 including the phasespace plots for v ′ x = 10 − and Ω = 20. In this region,the particles oscillate in the wave field without beingheated. Finally, for Ω ≪ DISCUSSIONWith the simulation results shown in Section 4 we arenow in the position to assess which of the broad spec-trum of waves in Fig. 3 are likely to provide stochasticheating of protons at the bow shock.Figure 7 shows the function χ p from Fig. 4c de-composed into discrete frequency dyads with orthogo-nal wavelets (Mallat 1999). The signal is divided intodiscrete frequency layers (dyads) that form 2 − n f N hi-erarchy starting from the Nyquist frequency ( f N is halfof the sampling frequency). Orthogonality means thatthe time integral of the product of any pair of the fre-quency dyads is zero, and the decomposition is exact,i.e., the sum of all components gives the original sig-nal. The y-labels are dyad numbers with the unit am-plitude corresponding to χ p = 70. We see that χ p inthe frequency channels from 0.25 Hz and above havesufficient amplitude, and correlate well with ion heatingseen in Fig. 3 in the time interval 14:23-14:24. On theother hand, in Fig. 6 we see that bulk proton heatingoccurs for f ≈ (1–10) f cp , while the stochastic acceler-ation of suprathermal particles can be done by waves f > f cp . Full kinetic simulations of the LHD instabil-ity (Daughton 2003) show that the instability developsat longer wavelengths k ⊥ ( r e r p ) / ≈ k ⊥ r p ≈
10 in our case, andhas lower frequencies f cp < f . f cp , with signifi-cant magnetic component (Gary 1993; Daughton 2003;Huang et al. 2009). This puts these waves in the bulkheating region of Fig. 6. They could be Doppler up-shifted by 5 Hz, so the possible frequency range for wavesthat could heat bulk protons is most likely 0.25-8 Hz in d y ad no , un i t =
64 Hz 32 Hz 16 Hz 8.0 Hz 4.0 Hz 2.0 Hz 1.000 Hz 0.500 Hz 0.250 Hz
Figure 7.
Decomposition of χ p from Fig. 4c in the range0.25-64 Hz shows that heating of protons shown in Fig. 3can be associated with waves at f ≥ .
25 Hz. One plot unitcorresponds to χ p = 70. the satellite frame of Figs. 3 and 7. Please note that k ⊥ ≡ k x throughout this paper.It is the result of the simulations, that the LHD wavesat lower frequencies f cp < f . f cp , and longer wave-lengths k ⊥ r p .
30 are found here to be responsiblefor the bulk proton heating. Incidentally, they also ap-pear to play a key role in the heating of plasma in themagnetotail and at the magnetopause (Zhou et al. 2014;Graham et al. 2019; Ergun et al. 2019).Waves at dyads 16 Hz and above in Figure 7 maycorrespond to shorter LHD wavelengths k ⊥ r e ∼ f lh (Davidson et al. 1977;Drake et al. 1983), which are Doppler upshifted tohigher frequencies, and also to the modified two-streaminstability which could be triggered by LHDI when theelectron drift velocity exceeds the ion thermal veloc-ity (Lashmore-Davies & Martin 1973; Gary et al. 1987;Umeda et al. 2014). This means that they may haveΩ ∼ f lh /f cp ∼
40 and k x r c = k ⊥ r p .
90 in Fig. 6.They can be associated with the heating of suprather-mal ions in the region Ω >
10, around the dashed line k x r c ∼ Ω in Fig. 6. CONCLUSIONSThis research has shown that there are two majorheating mechanisms in quasi-perpendicular shocks, asimplied from the analysis of 9 crossings of the bow shockby the MMS spacecraft.In this particular event, the electrons do not reach thestochastic heating level with | χ e | < K. Stasiewicz & B. Eliasson isotropization by LH/ECD waves excited by the den-sity compression. The quasi-adiabatic heating processis supported by the observed isotropic temperature re-lation
T /B = ( T /B )( B /B ) / which predicts a dipin the electron temperature-to-magnetic field ratio whenthe magnetic field increases.The ions instead undergo rapid non-adiabatic stochas-tic heating by electric field gradients perpendicular tothe magnetic field ( χ p ≈ T /B ratio in-stead shows an increase in the shock region. Test par-ticle simulations show that efficient stochastic heatingwithin 3 cyclotron periods takes place for a range of pa-rameters in space (Ω , v ′ x , χ ) where Ω = e ω/ω c is thewave frequency in plasma frame normalized by the cy-clotron frequency, v ′ x = k x r c is the normalized thermalspeed proportional to the ratio between the Larmor ra-dius to the wavelength, and χ = ( k x /ω c )( E x /B ) isthe stochasticity parameter representing the normalizedwave amplitude. We have identified in this space therange of the ion bulk heating and the range for acceler-ation of suprathermal ion tails (Figure 6).It is found that in the analyzed cases the ion bulk heat-ing is most likely accomplished by waves in the frequencyrange 0.25-8 Hz in the spacecraft frame, or (2-10) f cp inthe plasma frame, with k ⊥ r p .
30, i.e., with λ >
20 km.Waves at frequencies larger than 8 Hz in the spacecraftframe ( f > f cp in plasma frame) and with shorterwavelengths can provide acceleration of the tail of theion distribution function, producing diffuse energetic ionpopulation observed at shocks. The agreement between the theoretical and numericalresults with the MMS observations gives a more com-plete picture of the heating processes involved in theEarth’s quasi-perpendicular bow shock.The chain of the physical processes described in thispaper is initiated by a single event – namely – the com-pression of the plasma density N and the magnetic field B . The induced diamagnetic current triggers consecu-tively three cross-field current driven instabilities: LHD → MTS → ECD, which produce stochastic bulk heat-ing and acceleration of ions and electrons, in addition toa common quasi-adiabatic heating of electrons on com-pressions of B . Thus, the presented model has universalapplicability, and the processes described could occur inother types of collisionless shock waves in space, associ-ated with the density compression. The results may alsobe applicable to theories and models of particle heatingand acceleration in astrophysical shocks.ACKNOWLEDGMENTSThe authors thank members of the MMS missionfor making available the data. MMS science datais made available through the Science Data Cen-ter at the Laboratory for Atmospheric and SpacePhysics (LASP) at the University of Colorado, Boul-der: https://lasp.colorado.edu/mms/sdc/public/. B.E.acknowledges support from the EPSRC (UK), grantsEP/R004773/1 and EP/M009386/1. Software:
Data analysis was performed us-ing the IRFU-Matlab analysis package available athttps://github.com/irfu/irfu-matlab.REFERENCES