Resonant micro-instabilities at quasi-parallel collisionless shocks: cause or consequence of shock (re)formation
RResonant micro-instabilities at quasi-parallel collisionless shocks:cause or consequence of shock (re)formation
Vladimir Zekovi´c a) Department of Astronomy, Faculty of Mathematics, University of BelgradeStudentski trg 16, 11000 Belgrade, Serbia (Dated: 5 March 2019)
A case of two interpenetrating, cold and quasi-neutral ion-electron plasmas is investigated with the multi-fluid approach. We consider that one plasma flows quasi-parallel to the lines of a background magnetic fieldembedded in another static plasma. If the flow turns super-Alfv´enic, we show that parallel R/L-modes andperpendicular X/O-modes become unstable and grow in amplitude. Within the linear theory, we find thatthe growth rate curve of an unstable mode has a maximum at some wavenumber specific to each mode. If weconsider a shock-like plasma configuration, we find that the fastest growing mode is the resonant one (with k ∼ r − gi ) which strongly interacts with ions. In Particle-In-Cell (PIC) simulations, we observe that a resonantwave with the same properties is excited during the early phases of shock formation. Once the wave becomesnon-linear, it efficiently scatters ions and triggers the initial shock formation. This implies that the typicalcompression ratio of ∼ I. INTRODUCTION
Although studied theoretically for many years and ob-served in our space surroundings, collisionless shocks stillshow a variety of physical processes that reveal complexphysics behind them. Some basic properties of shocks ingeneral, including non-linear steepening, are explained byMHD theory and simulations. Nevertheless, such mod-els could not account for particle acceleration to cosmic-ray (CR) energies and the formation of non-thermal par-ticle spectra. When the first-order Fermi accelerationmechanism, known as diffusive shock acceleration (DSA),was initially introduced , it was very successfulin explaining these important features of shocks. How-ever, the picture of how particles gain high enough ki-netic energy so that they can be injected into DSA wasnot completely clear. Later, particle-in-cell (PIC) sim-ulations showed how this mechanism may work in thecase of quasi-parallel shocks. From PIC and kinetic hy-brid simulations , we know that particles are at firstspecularly reflected and pre-accelerated by an upstreammotional electric field E = − v × B (where v is thevelocity of the flow and B is the background magneticfield). While particles drift along the shock surface, theygyrate around the mean magnetic field, and they gainenergy by E in every cycle. This process, known asshock drift acceleration (SDA) , is a dominant accelera- a) Electronic mail: [email protected] tion mechanism at quasi-perpendicular shocks. At quasi-parallel shocks, particles are energized up to a few times E sh = mv / in the upstream, whichfurther accelerates particles to energies much higher than ∼ E sh , through a DSA-like process.To understand how a return current of ions is formed,we need to consider the fine structure of a collisionlessshock. The emerging possibility that the steepening ofa wave, which is essential to the formation of a colli-sional shock, is not a required condition in the case of acollisionless shock, opens a question about the processesthat trigger and shape the shock structure in an unstableand turbulent plasma. These processes are highly non-linear, but the theory of corresponding micro-physics ismostly linear and gives a more qualitative than quanti-tative explanation. Kinetic simulations can provide uswith detailed insight into how this rather complex struc-ture forms and how it changes over time. The transitionregion at quasi-parallel shocks is shown to be much widerthan at quasi-perpendicular ones (which are moresimilar to MHD shocks). This shock interface shows avariety of different types of instabilities and turbulencesthat can grow, then saturate or disrupt, or even mutuallyoverlap. The type of these changes and their dynamicsdepend on the strength and level of magnetization of theshock itself. Recent papers review the efforts madeto explain such systems and how numerical simulations,observations, and theory agree. a r X i v : . [ a s t r o - ph . H E ] M a r In kinetic simulations of quasi-parallel shocks, filamen-tation instability can be induced by CRs , which firstexcite the transverse modes in the upstream. The shockitself could result from the non-linear steepening of thesemodes , which can occur quasi-periodically after ev-ery few gyrations; this is a phenomenon also known asthe shock reformation process . In case of a quasi-perpendicular shock, an incoming flow is reflected off thecompressed magnetic field. At weakly magnetized or un-magnetized shocks , the Weibel instability forms currentand density filaments in the overlapping region. Thesefilaments merge (through the second Weibel instability)and disintegrate, thus forming the shock transition.These different types of shock interfaces, together witha reforming shock barrier, reveal the complexity of thesystem. This was the motive for investigating whetherthere is some other type of wave-plasma interaction thatcan quickly put, from ground state, the two interpene-trating collisionless plasmas into a shock configuration.On the microscopic level of this interaction, particlesshould be effectively scattered, and their motion random-ized. To make the closure and fulfill the pre-accelerationconditions, we require this interaction to be able to formthe initial return current. We also require that the sametype of interaction later starts the reformation process.We assume that at the very beginning of the shock for-mation, a finely tuned resonant wave-particle interactionoccurs at the interface of a moving and steady plasma.We observe this in the PIC simulations presented here.As the interface propagates with the plasmas, this inter-action can contribute to the thermalization by scatteringthe particles resonantly. The scattering should occur onthe scales of the order of a particle gyroradius. At suchscales, most of the particles will become trapped by thewave. If there is no turbulence, the resonance acts toprevent the particles from leaking even from the regiona few gyroradii away from the interface. Still, some ofthe particles can be reflected at the edge of the interfacewith a velocity fast enough to escape the wave. This sce-nario is therefore restricted to super-critical shocks withmagnetic inclination (cid:46) ◦ .The properties of such an idealized resonant interac-tion were examined in the previous work . Therein,we assumed that the growth of circularly polarized EMwaves is seeded and governed at the interface, wherethe bulk plasma flow interpenetrates the steady plasma.Using EM N-body simulations, we studied the motionof test particles in the field of circularly polarized ion-cyclotron frequency ( ∼ ω ci ) waves. We showed that thescattering of test ions by this wave is the strongest, if thewavelength of such a mode is resonant ( λ ∼ πr gi ) withrespect to the longitudinal gyroradius ( r gi = v (cid:107) /ω ci = mv (cid:107) /qB ) of the ions inflowing with the velocity com-ponent v (cid:107) (parallel to B ). Not only are the particlevelocity vectors then randomized along the direction ofthe flow, but there is also a population of particles thatare resonantly backscattered (reflected) from the wave.We found that the initial return current of reflected ions can be quite strong ( (cid:38)
10% of the ions impinging on thewave are reflected). We also showed that if the wave-length is resonant, the flow disrupts even when the waveamplitude is as low as ≈ B in the non-linear regime.For shorter or longer wavelengths, the amplitude of afew times B or greater is required for a monochromaticwave to efficiently scatter and disrupt the flow.In this paper, we use these empirical results as trac-ers to search for the equivalent interaction in the linearplasma theory, which naturally favors the growth of res-onant modes. We also use these results to calibrate thederived equations so as to analytically show how suchresonant instability grows and triggers the formation ofa shock without the involvement of any other turbu-lent mechanisms. We here associate the resonant wavethat appears in PIC simulations at the plasma interfacewith the strong-beam instability. The upstream wave wemodel as a weak-beam instability. Both instabilities areinduced by the same type of beam-plasma interaction,and they differ only in the strength of the beams thatexcite them. The first refers to shock triggering, and thesecond refers to shock reformation. By using PIC sim-ulations, we examine the behavior of these waves in thenon-linear regime.We organized the paper as follows: the equations gov-erning the propagation of plasma waves and stability cri-terion are derived in Sec. II; the theoretical results ob-tained here are applied to shock waves and comparedto the results of numerical simulations in Sec. III; andfinally, the possible role of the resonant instabilities intriggering the collisionless shock formation along withtheir role in governing the shock reformation process isdiscussed in Sec. IV. II. CYCLOTRON MICRO-INSTABILITY
We start from the fundamental case of two cold andquasi-neutral ion-electron plasmas that collide in thepresence of a background magnetic field. We considerthat one plasma (of a finite extent) flows quasi-parallelto the magnetic field lines with the velocity v , throughanother static plasma. In the case of wave frequenciesmuch lower than ω ci , recent studies show that lin-early polarized Alfv´en waves can grow and become cou-pled with the flow-driven ( ω = kv ) instability. Suchinstabilities can draw the energy for their growth fromthe macroscopic kinetic energy pool of interpenetratingplasmas and the background magnetic field. The equa-tions show that in the frame of a flowing plasma, theinstability which starts to grow at the edge of the inter-face, advects further away with the stream. Thus, wehave the plasma flowing into the interface with the ve-locity − v , and the instability advecting away from theinterface at a velocity of − v / (1 + η ). Therefore, if in thenon-linear regime the wave is capable of triggering thescattering of the particles, it can easily provide a plasmaconfiguration which by its dynamics resembles a shockwave. However, the instability growth rate has linear de-pendence with k , and the solutions are restricted onlyto low-frequency waves. None of the resonant k is fa-vored within this theory, which thus fails to provide thesolution for the predecessor wave in the linear regime.To include frequencies up to the ion-cyclotron ω ci andelectron-cyclotron ω ce , each of the plasma species mustbe treated separately. Moreover, if we confine to thespace plasmas that can be considered as collisionless,then the Vlasov equation (or the fluid equations derivedfrom it) should be used. From the linear theory andsimulations, we know that various electrostatic and EMmicro-instabilities can grow within a configurationwhere the two plasmas collide with a relative drift veloc-ity v parallel or anti-parallel to an average (or uniform)background magnetic field. The most important are thecircularly polarized EM two-component modes, such asion/ion, electron/electron and ion/electron modes, whichhave already been studied in literature . The ion/ionresonant mode is the most interesting among them,since it is immune to the small changes of the plasma pa-rameters because of its large wavelength. Its growth ratehas a maximum near the wavenumbers at which the realfrequency part satisfies a cyclotron resonance condition ω r (cid:39) v · k ± ω ci . It has large application to solar windand terrestrial bow shock , where it appears as a right-hand resonant instability initiated by a return current ofions in the upstream region. However, these modes areassociated with relatively weak particle beams.We consider here the modes excited by a strong plasmabeam, whose density is comparable to target plasmadensity. The beam and target plasmas consist of bothspecies, ions and electrons. Note that in the case con-sidered here, the resonance we refer to is not given bythe cyclotron resonant condition ω r (cid:39) v · k ± ω ci asin . Since the beam is quite strong, all terms in thedispersion relation become significant, so the resonancedoes not occur in the same way. In this case, the reso-nance is given by the condition that the wavelength ofthe resonant mode λ res is of the order of r gi . We alsoassume that the relative drift velocity of the two plas-mas is much larger than the average thermal velocity ofthe plasma components ( v (cid:29) v T ). The equations canthen be written in the cold plasma limit (closed by theequation of state p = 0, where p corresponds to the hy-drodynamic pressure). A. Equations of cold interpenetrating plasmas
Fluid equations of cold plasma consisting of ions andelectrons are given by : m i (cid:18) ∂∂t + V i · ∇ (cid:19) V i = q i ( E + V i × B ) , (1) m e (cid:18) ∂∂t + V e · ∇ (cid:19) V e = q e ( E + V e × B ) , (2) where E and B are electric and magnetic field compo-nents, ρ i,e = n i,e m i,e is the density of ions (electrons),and q i,e their charge. For the flowing plasma, the veloc-ity of each species has the same constant component v added to its variable part v i,e , so that the total fluid ve-locity is V i,e = v + v i,e . For each variable, we considerperturbations of the form f ( r , t ) = f ( k , ω ) · e i ( kr − ωt ) .The dispersion matrix that constitutes a system of lin-earized equations (as derived in Sec. A in the Appendix)is D = P − n sin θ n sin θ cos θ S − n − iDn sin θ cos θ iD S − n cos θ ++ i(cid:15) ω σ f , (3)where θ is the angle between the wavevector k and thebackground magnetic field B ; and n is the refractive in-dex. The components P , S , D and matrix σ f are givenby Eqs. A16, A17, A18 and A19, respectively. Disper-sion relation is found by equating the determinant of thedispersion matrix to zero, |D| = 0.Two specific cases are considered in the following anal-ysis – solutions to the dispersion relation for wavevector k parallel (Sec. II B) and normal (Sec. II C) to the back-ground magnetic field B . B. Wavevector parallel to the magnetic field
Without loss of generality, the vector of the back-ground magnetic field is taken to be in the direction of the x -axis ( B = B ˆx ). In case of θ = 0, the unit wavevectoris ˆk = ˆx , and the dispersion matrix (3) reduces to D = P S − n − iD iD S − n −− (cid:88) α = i,e ω pαf ω k x v x ξω · k x v x ξω , (4)where ω pαf is the plasma frequency of the species α inthe flowing plasma, and ξ = 1 − v · k /ω is the modifica-tion parameter of the frequency due to the flow velocity.Dispersion relation is then found as P − (cid:88) α = i,e ω pαf ω k x v x ξω · (cid:18) k x v x ξω (cid:19) · (5) · (cid:2) ( S − n ) − D (cid:3) = 0 , whose roots are assorted according to the associatedeigenvectors of the electric field.The first root branch of this equation corresponds tothe dispersion relation of a longitudinal electrostatic os-cillation advected by the flow. After substituting P fromEq. (A16) and by making the use of the equivalence k x v x ≡ k · v , which is applicable in this parallel case,the first root becomes1 − (cid:88) α = i,e ω pαs ω − (cid:88) α = i,e η ω pαs ( ω − k · v ) = 0 , (6)where the density ratio of the flowing to the stationaryplasma is given by η = ω pif /ω pis = ω pef /ω pes . Thisis the standard cold plasma dispersion relation for elec-trostatic waves, with the inclusion of the velocity driftbetween the different plasma species . Depending onthe wave frequency, a longitudinal electrostatic wave canbe coupled with and advected by the flow. This modecan become unstable when the electrostatic wave andthe streaming plasma are close to being in resonance ω → k · v . Energy is thus transferred from the elec-trons to the wave, leading to the exponential growth ofthe wave amplitude. As shown in , the waves saturateby trapping the electrons, and then collapse, making theelectrons become rapidly thermalized.The second root branch of Eq. (5) is given by n = S ± D = R ( L ) ,R ( L )= 1 −− (cid:88) α = i,e ω pαs (cid:20) η ξ ξω ( ξω ± ω cα ) + 1 ω ( ω ± ω cα ) (cid:21) , (7)which is essentially the same as in the case of a two-species, three-component cold plasma consisting of a ten-uous beam that interpenetrates a dense core and thethird component of the other species.These right (R) and left (L) circularly polarized wavesgiven by Eq. (7) can both turn into instabilities. Forcomparable densities of the flowing and stationary plas-mas, two types of instabilities can occur. One of themis resonant and it has λ ∼ πr gi . The other type is sub-resonant in the sense that its wavelength is of the order ofthe longitudinal gyroradius of streaming ions ( λ ∼ r gi ).These two modes are very similar to EM ion/ion modes,as given in . We here solve the equations of each of themodes analytically (see Sec. B in the Appendix), whichis only possible when considering some specific range offrequencies. Real ( ω r ) and imaginary ( γ ) parts of thefrequency ω = ω r + iγ are thus derived as solutions tothe quartic equation.In the range of ion-cyclotron frequencies, condition ξω ∼ ω ci (cid:28) ω ce implies the interaction between ions ofthe flowing and stationary plasmas. Eq. (7) is thereforeapproximated by R ( L ) ≈ ω pis (cid:20) η ξ ω ci ( ω ci ± ξω ) + 1 ω ci ( ω ci ± ω ) (cid:21) . (8) For frequencies much lower than ω ci , two of the fourroots of this relation merge to form the solutions whichare the same as in Ref. . Eq. (8) therefore representsthe generalization of the case where at low frequenciesan Alfv´en wave couples with the flow driven ( ω = kv )mode .These two roots, given by the quartic pairs ω , or ω , , are thus both assigned to the R or L-mode, respec-tively, depending on whether “+” or “ − ” sign is chosenin Eq. (B6) in the Appendix. Real and imaginary partsof the complex frequency ω ... = ω r + iγ correspond tothe real frequency and growth rate of the wave, respec-tively. The two roots of the L-mode equation are plottedin Fig. 1. One of them grows with γ ≈ . ω ci , whilethe other one decreases at the same rate. - - - Ω r Ω ci k × r gi - - ΓΩ ci FIG. 1. The L-wave real frequency and growth rate de-pendences of the wavenumber for plasmas with the densityratio η = 3. The real part ( ω r ) of each of the roots isplotted as a continuous line and the imaginary part ( γ ) asa dashed line; both expressed in relative units of ω ci , whilethe wavenumber is given in units of r − gi . The curves of ω and ω are given in dark blue and red color, respectively.The position of the maximum is at the resonant wavenumber k L max = 2 π/λ res ≈ /r gi . Earlier, we determined by test particle simulations that a beam of ions interacts with the L-wave only in avery narrow band of wavelengths around λ res ≈ πr gi ,where r gi = m i v /q i B and v is the velocity of thebeam. From Fig. 1, we observe that both roots havea maximum at some wavenumber k L max . If we imposethe condition that the maximum must be located ex-actly at k L max = 2 π/λ res = 2 /r gi in the equations, theythen give the density ratio of the two plasmas. There-fore, when the wave is resonant, this ratio is η ≈
3. Orto put it another way, if the flowing plasma density isthree times the density of the stationary plasma, the L-mode (that emerges from the wave-plasma interaction)is then resonant. Near the beginning of the overlappingregion, where this resonant wave grows, there is an over-all density ratio of η + 1 ≈
4. The linear equations thusshow that the resonant wavenumber k L max is related tothe shock-like density ratio at the interface of the twointerpenetrating plasmas. Although this connection maybe a coincidence, it still implies that the resonant or sub-resonant modes should grow the fastest. - - - Ω r Ω ci k × r gi - - ΓΩ ci FIG. 2. Same as in Fig. 1, only for the R-wave. The curvesof ω and ω are given in green and yellow color, respectively.The growth of the wave is seeded at the sub-resonant wave-lengths k R max = 2 π/λ sub ≈ /r gi ≈ k L max . The position ofthe resonance k L max is marked by a vertical red line. There is a pair of solutions to the R-mode that alsohas the non-zero complex part of the frequency, which isshown in Fig. 2. It corresponds to a sub-resonant (withthe ions gyroradii) wave-mode. For comparable densitiesof the two colliding plasmas, the wavenumber of the max-imum growth of the sub-resonant wave k R max is roughlyproportional to the maximum of the resonant mode as k R max ∼ η · k L max . The absolute positions of these maximaand their position relative to each other do not depend onthe Afv´enic Mach number ( M A = v /v A ) of the flow, ifit is higher than ∼
10 (for M A (cid:46)
5, the growth rates aremodified). It is shown by simulations that ions pass-ing through such a sub-resonant R-wave do not interactsignificantly with it; they instead gain a small wigglingmotion around the equilibrium position along their flow.Even if the wavelength is set to resonant, interaction stillremains very weak. For the wave amplitudes of ∼ B ,polarization of the wave is therefore a decisive property;it determines whether or not interaction with ions willemerge.In Fig. 3 we see that for most of the values of η , thegrowing modes of L and R-waves at λ res have nearly thesame growth rate. They have slightly different real fre-quencies and, therefore, different phase speeds. In therange of η ∼ [0 . − γ ≈ . ω ci , and the difference betweenthe real frequencies of the two modes is in the range∆ ω r /ω ci ∼ [0 . − E L = E ( γt ) e i ( k · x − ω r t ) √ i , E R = E ( γt ) e i ( k · x − ω r t ) √ i e − i ∆ ω r t , E L + R = E ( γt ) e i ( k · x − ω r t ) √
01 + e − i ∆ ω r t i (1 − e − i ∆ ω r t ) . After multiplying the last expression by e − i ∆ ω r t/ /e − i ∆ ω r t/ , we apply the Euler’s formulato get: E L + R = 2 E ( γt ) e i [ k · x − ( ω r + ∆ ω r t ] · (9) · √ (cid:0) ∆ ω r t (cid:1) − sin (cid:0) ∆ ω r t (cid:1) . In the linear regime, both waves grow at the same rate γ . Eq. (9) implies that the superposed wave should growas linearly polarized. Its plane of polarization shouldrotate (the Faraday rotation). We estimate the rotationperiod ∆ T of the plane of polarization, by equating thephase of the rotation to 2 π . From the previous equationwe get: ∆ ω r τ rot = 2 π, τ rot = 2 · π ∆ ω r . For the resonant wavenumbers (∆ ω r ≈ ω ci ), the rota-tion period is τ rot ≈ · T ci , where T ci is the period oftime required for an ion to make one orbit around themagnetic field line. Due to a fast growth rate, the waveremains in the linear regime only for a few cycles, untilthe amplitude grows enough ( ∼ B ) to trigger the reso-nant scattering of ions . We think that in such a non-linear regime, the magnetic field line bends to a helicoidalshape, and the vector of the electric field starts to rotate.Since the R-wave only weakly interacts with ions , thiswill eventually become the point of separation of the Land R modes, where polarization of the wave turns froma rotating linear to a left circular.We emphasize here that in case of the cold and weak( η (cid:46) .
1) ion beam, the maximum of the growth rate forthe R and L modes shifts toward the lower wavenumbers.For the R-mode, the position of the maximum is thennear the point where k · r gi ≈
1. The obtained real fre-quency and the growth rate have the same values as thosestudied in for the ion/ion right-hand resonant insta-bility. The cyclotron resonant condition ω r (cid:39) v · k ± ω ci then applies too. For a fixed Alfv´enic Mach number,lowering the density of the flow makes the L-wave insta-bility turn off much earlier than the R-wave instability. - - - Ω r Ω ci k × r gi - - ΓΩ ci Η = - - - Ω r Ω ci k × r gi - - ΓΩ ci Η = - - - Ω r Ω ci k × r gi - - ΓΩ ci Η = - - - Ω r Ω ci k × r gi - - ΓΩ ci Η = - - - Ω r Ω ci k × r gi - - ΓΩ ci Η = - - - Ω r Ω ci k × r gi - - ΓΩ ci Η = - - - Ω r Ω ci k × r gi - - ΓΩ ci Η = - - - Ω r Ω ci k × r gi - - ΓΩ ci Η = - - - Ω r Ω ci k × r gi - - ΓΩ ci Η = - - - Ω r Ω ci k × r gi - - ΓΩ ci Η = - - - Ω r Ω ci k × r gi - - ΓΩ ci Η = - - - Ω r Ω ci k × r gi - - ΓΩ ci Η = - - - Ω r Ω ci k × r gi - - ΓΩ ci Η = - - - Ω r Ω ci k × r gi - - ΓΩ ci Η = - - - Ω r Ω ci k × r gi - - ΓΩ ci Η = - - - Ω r Ω ci k × r gi - - ΓΩ ci Η = - - - Ω r Ω ci k × r gi - - ΓΩ ci Η = - - - Ω r Ω ci k × r gi - - ΓΩ ci Η = - - - Ω r Ω ci k × r gi - - ΓΩ ci Η = - - - Ω r Ω ci k × r gi - - ΓΩ ci Η = - - - Ω r Ω ci k × r gi - - ΓΩ ci Η = - - - Ω r Ω ci k × r gi - - ΓΩ ci Η = - - - Ω r Ω ci k × r gi - - ΓΩ ci Η = - - - Ω r Ω ci k × r gi - - ΓΩ ci Η = - - - Ω r Ω ci k × r gi - - ΓΩ ci Η = - - - Ω r Ω ci k × r gi - - ΓΩ ci Η = FIG. 3. The L-wave (blue and red) and the R-wave (green and yellow) real frequency and growth rate dependences onwavenumber. Graphs of the L and R modes are sorted in pairs (L, R), for which the product of their density ratios is η L · η R = 1. The real part ( ω r ) of each of the roots is plotted as a continuous line and the imaginary part ( γ ) as a dashed line;both expressed in relative units of ω ci , while the wavenumber is given in units of r − gi . The L-mode instability of such a tenuous beam has thesame value and shape of ω r and γ dependences, as in thecase of the ion/ion non-resonant instability . How-ever, here, this mode propagates in the direction of the flow. This is probably the reason this instability is leftcircularly polarized, instead of being right circularly po-larized and propagating opposite to the flow, as studiedin .In Fig. 3, the R and L modes are shown together inpairs for the different density ratios of the two plas-mas. Besides the double peak that emerges at ratios η ∼ . − . C. Wavevector normal to the magnetic field
Without loss of generality, we take the wavevector tobe in the direction of z -axis. We find (see Sec. C inthe Appendix) that the plasma flowing in the directionparallel to the magnetic field lines has no influence on ex-traordinary EM (X-mode) waves. However, it can makeordinary EM (O-mode) waves subjected to the Weibelinstability. On the other hand, the plasma flowing in thedirection perpendicular to the magnetic field lines andparallel to k = k ˆz has no influence on the O-mode, butit can make the X-mode become unstable.If the flow is perpendicular to B and parallel to thewavevector ( v = v ˆz ), the polarization matrix given byEq. (3) takes the form: D = P − n S − n − iD iD S + (10)+ (cid:88) α = i,e ω pαf ω · − i ω cα kv ξ ω − ω cα i ω cα kv ξ ω − ω cα − ωkv ξ ω − ω cα ( ξ + 1) . The X-mode dispersion relation is given by the root:( S − n ) S M − D M = 0 , n = RL − SM S − DM D − M D S − M S , where S M and D M are derived in Sec. C in the Appendix.For frequencies around the ion-cyclotron or lower, thedispersion relation remains complicated and we solve itnumerically.Similarly to parallel modes of propagation, the modi-fied X-mode also becomes unstable if the flow is super-Alfv´enic. The difference here is that the wavenumber ofthe maximum growth of the X-mode instability is some-what lower ( ∼ r − gi ) than that of the R and L modes. Themaximum growth rate is also at least two times smallerfor this mode. The instability of the perpendicular flowis presented for different η in Fig. 4. Contrary to the caseof the parallel modes, the position of the maximum forthe X-mode does not depend on the velocity of the flow.Nevertheless, its growth rate does.Finally, if the plasma flow is inclined to the backgroundmagnetic field, the resulting instability is a compositeone. Its component parallel to B is composed of the Ror L mode, and the component which is normal to B consists of the X-mode. III. PARTICLE-IN-CELL SIMULATIONS OF SHOCKWAVES
In this section we use the kinetic simulations to resolvethe whole process of a quasi-parallel collisionless shockformation in a fully non-linear case. We emphasize thepotential role of resonant micro-instabilities in a shocktriggering. Once the shock is developed, we show howthe same type of instability is excited in the upstreamby a return plasma beam. We show how this upstreammicro-instability, once it becomes non-linear, governs theprocess of shock reformation.
A. Instability in the non-linear regime
We first consider the wave-plasma interaction thateventually leads to the triggering of a shock wave. Tounderstand what happens in the initial phases of thenon-linear regime, we combine the linear theory and theresults of non-linear test particle simulations . In thisestimation, we consider the interaction in a laboratoryframe, which is the frame of the ISM. We assume that theexcited wave in the non-linear regime does not steepensignificantly. For the monochromatic transverse wave, itis shown in that the energy will not be significantlytransferred from fundamental to higher harmonics. Thismeans that the waves that deviate from the linear dis-persion law are not much distorted due to non-linearity.Static plasma (ISM) is reflected off the moving ejecta,so initially the two interpenetrating plasmas have thesame densities ( η = 1). If the velocity of the flowingplasma is v (cid:38) v A , the fastest growing mode is the reso-nant one (see Fig. 3), which has λ = πr gi ( k = 4 /r gi ).For a very short period of time, the instability shouldgrow as a linearly polarized wave that propagates with aphase velocity v ph ≈ v in the direction of the flow. Inthe instability frame (the center-of-mass frame), the flow-ing plasma is moving with relative velocity ∆ f v ≈ v .For this relative flow, the wavelength of the instability isresonant: λ ≈ π m i v q i B = π m i ∆ f vq i B = π ∆ f r gi . Therefore, we expect that ions are strongly scattered once this wave reaches the amplitude ∼ B . Becauseshocks are strongly non-linear, with this estimate we canonly discuss the interaction during the period until theparticles become scattered by the instability. B. Simulation setup
In a supernova explosion, the stellar material is blownaway from the center of the star, and initially, it freelypropagates through ISM. Due to a much larger density - - Ω r Ω ci k × r gi - - - ΓΩ ci Η = - - Ω r Ω ci k × r gi - - - ΓΩ ci Η = - - Ω r Ω ci k × r gi - - - ΓΩ ci Η = - - Ω r Ω ci k × r gi - - - ΓΩ ci Η = - - Ω r Ω ci k × r gi - - - ΓΩ ci Η = - - Ω r Ω ci k × r gi - - - ΓΩ ci Η =
FIG. 4. Similar to Fig. 3, only for the X-wave instability. The curves of the growing and decreasing waves are given in darkblue and red color, respectively. The position of the maximum is near the resonant wavenumber, at k X max = 2 π/λ res ≈ /r gi . and amplified magnetic field, the expanding ejecta be-haves like a physical barrier, reflecting the ISM parti-cles it encounters. This process is frequently presentedin 2D or 3D PIC simulations , where the initiatedplasma flow is simply reflected off the wall that replacesthe ejecta. In the region where two counter-streamingbeams overlap, instabilities develop and make the plasmabecome turbulent, forming a shock wave. The simulationframe thus corresponds to the downstream (or ejecta)frame.We numerically solve kinetic plasma equations by us-ing the PIC code TRISTAN-MP . We initiate the shockby reflecting a cold stream of particles (ions + electronscoming from the right side of a simulation box) froma conducting wall that is on the left. In one case, weconsider the background magnetic field lines that are ex-actly parallel to the flow. In the other case, the linesare inclined relative to the direction of the flow withthe inclination angle θ = 15 ◦ . Simulation setup hereis similar to the case presented in . We use a fixedsize simulation box of a rectangular shape in the x − y plane, with periodic boundary conditions in the y direc-tion. The computational domain is 256 cells wide (along y ) and 50000 cells long (along x ), with the skin depthof 10 cells and each cell initially containing 4 particles(two electrons and two ions). In physical units, the sizeof the domain is ∼ × c/ω pe . The end time ofthe simulation is T = 4500 ω − pe . The bulk Lorentz fac-tor of the electron-ion flow is given by γ , and the beamthermal spread is ∆ γ = 10 − . Magnetization is de-fined as the ratio of magnetic to kinetic energy density σ = B / πγ n i m i c = ω ci /ω pi . We choose the constantvalues for electron magnetization σ e = ω ce /ω pe = 0 . m i /m e = 16, which implies σ = 0 . v A ≈ . c . All the simulation runs are restrictedto low M A < C. Resonant instability and shock triggering
Here, we present the resonant mechanism of a shocktriggering, which we observe in PIC simulations. Ini-tially, a beam of particles is reflected off the conduct-ing wall, and the two beams start to interpenetrate eachother. Right at the first stage of the collision, we ob-serve that the wave which is resonant with the velocityof the inflowing plasma is excited inside the overlappingregion (see Fig. 5). Once the wave grows to the ampli-tude ∼ B , we find that ions become strongly scatteredby this wave. In the phase spectrum, we observe that thetwo counter-streaming beams start to scatter and mergeperiodically at localized regions. We see that the period-icity of this merging is related to the wavelength of thegrowing instability. Such behavior in the phase spectrumis expected, if the interaction is of the resonant type (seeFig. 2 in Ref. ). The wave grows as linearly polarized,and once it couples non-linearly with the plasma, its po-larization switches to a left circular. We also find thatits wavelength depends on the ion mass and the velocityof the flow as λ ∼ m i v /q i B , which reveals the Alfv´enicnature of the wave.By previous properties, we find that this wave stronglymatches the strong-beam resonant instability that weanalyzed in Sec. II. In Fig. 5, we see that the shocktriggering occurs in the same way as we discussed inSec. III A. The resonant wave therefore disrupts the flowby capturing the ions. We found that in this stage,the wave remains nearly monochromatic, and its shapeis only slightly distorted. This is in part expected forthe transverse monochromatic wave that enters the non-linear stage .Immediately afterwards, ions that are not captured(the part of the beam slightly in front of the wave) arethus cut off by the resonant mode. This sparse bunch ofions continues to propagate and excites the wave of anopposite polarization ahead of the forming shock. Thiswave grows and advects towards the shock. In this stage, x [ c / pe ]1.00.50.00.51.0 i x , i x [ c / pe ]42024 B / B x [ c / pe ]0246 d e n s it y [ n ] x [ c / pe ]1.00.50.00.51.0 i x , i x [ c / pe ]42024 B / B x [ c / pe ]0246 d e n s it y [ n ] x [ c / pe ]1.00.50.00.51.0 i x , i x [ c / pe ]42024 B / B x [ c / pe ]0246 d e n s it y [ n ] x [ c / pe ]1.00.50.00.51.0 i x , i x [ c / pe ]42024 B / B x [ c / pe ]0246 d e n s it y [ n ] x [ c / pe ]1.00.50.00.51.0 i x , i x [ c / pe ]42024 B / B x [ c / pe ]0246 d e n s it y [ n ] x [ c / pe ]1.00.50.00.51.0 i x , i x [ c / pe ]42024 B / B x [ c / pe ]0246 d e n s it y [ n ] x [ c / pe ]1.00.50.00.51.0 i x , i x [ c / pe ]42024 B / B x [ c / pe ]0246 d e n s it y [ n ] x [ c / pe ]1.00.50.00.51.0 i x , i x [ c / pe ]42024 B / B x [ c / pe ]0246 d e n s it y [ n ] x [ c / pe ]1.00.50.00.51.0 i x , i x [ c / pe ]42024 B / B x [ c / pe ]0246 d e n s it y [ n ] x [ c / pe ]1.00.50.00.51.0 i x , i x [ c / pe ]42024 B / B x [ c / pe ]0246 d e n s it y [ n ] x [ c / pe ]1.00.50.00.51.0 i x , i x [ c / pe ]42024 B / B x [ c / pe ]0246 d e n s it y [ n ] FIG. 5. The p x - x phase space, transverse fields B y and B z , and density profile captured at the very first phases of shockformation t = 45 − ω − pe . On the horizontal axis, the x -coordinates are given in units of c/ω pe . we observe a sharp change in the phase angle of the mag-netic field, which clearly separates the two waves. Thereformation is triggered once the upstream wave growsenough so that its amplitude can be compressed to atleast ∼ − B at the shock. The shock reformationis further governed by the upstream wave, which com-pletely overtakes the processes of particle thermalization and shock transmission.For Mach numbers (cid:46)
5, the shock ramp forms quitegradually. In Fig. 6 (the case where M A ≈ E x profilepresented in Fig. 8, which indicates the overall slowing0down of the plasma flow. For M A ≈
3, this peak ap-pears at the junction of the upstream and downstreamwaves, which thus clearly separates the upstream anddownstream flows. The resonant wave that triggered thisshock did not change its properties significantly while itwas entering the non-linear stage. The weak upstreamwave was excited later by the escape ions. The down-stream structure actually corresponds to the imprint ofthe resonant wave, which entered the non-linear regime,but is still not affected by the shock reformation. Theupstream oscillations in E x profile for M A ≈ ∼ B . We also showed that the finely tuned resonantwavelength ( λ res ≈ πr gi ) and the left circular polariza-tion of the instability are the properties we get whenthe two plasmas are subjected to shock-like conditions.The downstream structure thus shows a great similaritywith the resonant instability that we studied in Sec. II.Although, the linear theory does not apply in this non-linear stage, it may point to why the wavelength andpolarization change.In the case of an inclined magnetic field (see Fig. 7), wesee that the wavevector of the upstream mode is alignedwith the lines of the inclined magnetic field. The resonantwave, however, propagates in the direction of the flow. Ifwe consider the case of super-Alfv´enic flow with a non-relativistic velocity ( c (cid:29) v (cid:29) v A ), we find from theEq. (8) that the approximate analytical expression forthe growth rate in a close region around the maximumis: γ = ηη + 1 κ − κ + 1 , (11)where γ M = γ/ω ci is the relative growth rate, and κ = k · v /ω ci is the relative wavenumber in the direction ofthe flow ( η is a density ratio, as before). Since the relativewavenumber is defined by the scalar product k · v , it willbe the highest for the wave modes that are aligned withthe direction of the flow. Therefore, from Eq. (11), wefind that the resonant mode, which propagates in thatdirection, will grow at the fastest rate. Consequently,this means that the wavevector of this mode should bealigned with the flow, rather than being aligned with B . We observe such behavior of the resonant wave inthe simulations. It is observed during the linear stage(before a shock is triggered) as well as in the non-linearstage, when shock reformation is ongoing. In the case of M A ≈ p x − x diagram,it can be seen that the flow is dispersed and compressedby a factor ∼ B with themean velocity of ∼ − . c . Once the plasma reaches thewall, it reflects with a velocity of ∼ . c . It thus doublesits density up to a ratio ∼
6. Because the two counter-streaming plasmas are sub-Alfv´enic in the downstream,the wave mode there changes. The wave decreases itsamplitude and allows the plasma to diffuse freely towardsthe upstream. This region spreads between the wall onthe left and the point on the right, where the compressionis r ∼
4. It gradually expands across the downstream,creating the observed density profile. We observe thisduring the whole simulation run.
D. Upstream non-resonant instability and shockreformation
In the simulations, we observe that the ions not cap-tured by the resonant wave during the shock triggering,escape and form the initial pulse of the return current.We find that this current excites a circularly polarizedEM wave in the upstream, which is the first to trig-ger shock reformation (once it grows to the amplitudeof ∼ B ). The initial bunch of escape ions continues topropagate and pre-heats the plasma ahead of the shock.It thus forms a weak precursor in the upstream. In thenext stage, shock reformation completely overtakes theprocess of particle scattering (thermalization) and fur-ther mediates the shock . The return current is formedby specular reflection of ions from the reforming shockbarrier, as described in the minimal model for ion injec-tion . The upstream wave is then coupled with the cur-rent of reflected ions to create a self-sustained structure.Depending on the effective Mach number, the upstreammode can be both – resonant or non-resonant , wherethe latter is also known as Bell’s instability . We alsofind that the weak-beam instability that we studied inSec. II has the same property. Lower Mach numbers ofthe current of reflected ions cause the upstream insta-bility to be resonant, and higher Mach numbers cause itto be non-resonant. As the instability accelerates par-ticles through the process of DSA , the return currentpenetrates deeper into the upstream. By this process,the instability spreads farther upstream, and the parti-cles are accelerated to higher energies. The non-thermalpeak appears in the particle spectra, shortly after theshock initially forms and the return current initiates theupstream mode.In Fig. 6 (the left column), we observe the upstreaminstability as a right circularly polarized wave that is anegative helicity wave. In the region farther away fromthe reforming shock barrier (800 − c/ω pe ), the waveis of a smaller amplitude ( < B ), and it does not perturb1 x [ c / pe ]1.00.50.00.51.0 i x , i l og f i ( p ) x [ c / pe ]505 B / B x [ c / pe ]02468 d e n s it y [ n ] x [ c / pe ]01020 y [ c / p e ] p f ( p ) T e = 0.480 m e c T p = 0.640 m e c x [ c / pe ]0.40.20.00.20.4 i x , i l og f i ( p ) x [ c / pe ]505 B / B x [ c / pe ]02468 d e n s it y [ n ] x [ c / pe ]01020 y [ c / p e ] p f ( p ) T e = 0.128 m e c T p = 0.128 m e c FIG. 6. From top to bottom: p x - x phase space, transversely averaged fields B y and B z , density profile, B z field in 2D, andspectra of particles. Two cases are shown for the beam velocity: v = 0 . c ( M A ≈
6) at time t = 3500 ω − pe = 875 ω − pi on theleft, and v = 0 . c ( M A ≈
3) at t = 5400 ω − pe = 1350 ω − pi on the right. the upstream flow, nor the return current. In this region,the upstream wave can approximately be considered likeit is in the linear stage. Therefore, we model this waveas a weak-beam, cold plasma instability that we alsoanalyzed here in Sec. II. As Eq. (8) implies, the instabil-ity is kinematically bound to the upstream plasma be-cause the upstream flow is much denser than the returnbeam of ions. By its polarization and wavelength, theupstream instability is non-resonant with the upstreamflow. However, reflected ions “feel” it as a left circu-larly polarized wave – it therefore interacts resonantlywith the return current. The cyclotron resonant condi-tion ω r (cid:39) v ri · k − ω ci applies to the return beam, where v ri is the velocity of returning ions relative to the upstreamflow.Here, we apply the results from Sec. II to estimate theamount η ri of reflected ions that drive the weak-beaminstability, which has the same properties that we observefor the upstream wave. We assume that ions are reflectedfrom the periodically reforming shock barrier and thatthe reflection occurs in the downstream frame (as shownin ). The mean velocity of injected ions in the ISM (lab)frame is then v ri ∼ . v sh , which corresponds to theinjection energy E inj ∼ E sh in the shock frame (whichwe observe). From Fig. 6, we see that the wavelength of the upstream mode is λ up ≈ c/ω pe . We inspect Eq. (8)to find the beam density η ri , for which the growth rate γ ( k ) has a maximum at k · r gi = k · v ri ω ci = 2 πλ up / (cid:104) cω pe (cid:105) · v ri v A (cid:114) m i m e . In the case M A ≈ .
34 (Fig. 6), for ions which areinjected with a velocity v ri and excite a wave with k · r gi ≈ .
4, we obtain that η ri ≈ . , where injectionfraction is ∼ . M A = 20 and λ up ∼ c/ω pi , we find that η ri ∼ .
01. If we use the same duty cycle for ion injection ∼
25% (meaning that the probability of a particle beingcaptured by a reforming barrier is
P ∼ . ,we find that the number of SDA cycles needed to obtain η ri = (1 − P ) N is N ∼ .
3. This implies the injectionenergy E inj ∼ E sh . In our case of M A ≈ .
34, the num-ber of cycles is
N ∼ .
4, the same as in . It means thatin the case of a quasi-parallel shock, specularly reflectedions are energized via SDA. Only those ions that expe-rience N cycles of SDA are able to overcome the shockbarrier and drive the upstream instability.2 x [c/ω pe ] ω [ c / ω p e ] −3−2−101230 200 400 600 800 1000 1200 x [c/ω pe ] ω [ c / ω p e ] x [c/ω pe ] ω [ c / ω p e ] x [c/ω pe ] ω [ c / ω p e ] x [c/ω pe ] ω [ c / ω p e ] FIG. 7. Two-dimensional plot of the transverse magneticfield component B z in the units normalized to B , obtainedfor γ = 15 at times t = { , , , , } ω − pe ,from bottom to top, respectively. Magnetic field lines areinclined to an angle θ = 15 ◦ . Therefore, we relate the weak-beam, cold plasma insta-bility with the upstream wave, which we observe in thesimulations. In the non-linear stage, this instability trig-gers shock reformation. While being transmitted throughthe shock, the upstream instability changes its polariza-tion and wavelength to form the (non-linear) resonantwave in the downstream. The amplitude of the instabil-ity grows to the value of (1 . − . B just in front of theshock interface. There, the upstream flow slows downand becomes compressed. The transverse component ofthe magnetic field also compresses, and the amplitude ofthe instability increases to r · (1 . − . B ≈ (6 − B ,which is commonly observed in kinetic simulations .At the shock transition, due to longitudinal compression,the instability wavelength also shrinks by the same ratio.We observe that in the region of the reforming shockbarrier, the beam of returning ions is periodically dis-rupted by the upstream mode in the longitudinal direc-tion. One part of the returning ions slows down and be-comes trapped locally along the path, at places where theamplitude of the upstream instability reaches its maxi-mum. We find that this deceleration occurs because of the wave resonance with the velocity of returning ions.Due to a slightly different velocity of the upstream modewith regard to the upstream flow, we also find that theflow is non-resonantly perturbed by the instability. In theregion where the wave becomes non-linear, ions fall downinto a periodic potential well (that is imposed by the up-stream wave), and form clumps. These clumps of ionsare bound to the instability, which advects them to theshock. While approaching the shock region, the plasmaclumps gradually slow down and become compressed inexactly the same way as the upstream instability. Den-sity spikes with η ∼ . η ∼ −
10, as weobserve in the density profile shown in Fig. 9. The spikesperiodically reform ahead of the shock, at the regionswhere the growing instability becomes compressed. Thisis the typical signature of the shock reformation processthat is commonly observed in kinetic simulations .We also notice that the upstream instability may beamplitude modulated, as shown in Fig. 6. This modu-lation directly reflects the upstream plasma density bymaking the particles bunch more or less in synchronic-ity with the change of the wave envelope. We observethis in the density profiles shown in Fig. 9, where wefound that the amplitude of the spikes is modulated inthe same way as the upstream instability. This meansthat the shock might reform with the periodicity of thecarrier wave (the upstream instability), which we see inthe upper-half of the density graph. However, in thelower-half of the graph, we observe the regions wherespikes are of a low amplitude or absent. These regionscorrespond to the minima in the amplitude of the mod-ulation waveform. Contrary to them, the regions withamplified spikes correspond to the maximum amplitudeof the modulating wave. As a result of the collectiveeffect, within the observed time interval, the shock re-forms on two time scales: shorter – with the length ofone oscillation of upstream instability, and longer – withthe length similar to the wavelength of the downstreammode. There is a shock-pause in profiles that reside inthe middle of Fig. 9 (around ω ci t ∼
350 400 450 500 550 600 650 x [c/ω pe ] y [ c / ω p e ] x [c/ω pe ] −200000200004000060000 E / E
350 400 450 500 550 600 650 x [c/ω pe ] −10−50510 ω / ω
250 300 350 400 450 x [c/ω pe ] y [ c / ω p e ] x [c/ω pe ] −5000050001000015000 E / E
250 300 350 400 450 x [c/ω pe ] −4−3−2−101234 ω / ω FIG. 8. From top to bottom: density field in 2D, transversely averaged fields E x , B y and B z . The two cases of M A ≈ M A ≈ E x profiles are associated with the collective slowing down of theflow: at the first maximum of the downstream instability (right) and at every maximum of the downstream and compressedupstream instabilities (left) of the existing and reformed shocks, respectively. IV. RESULTS AND DISCUSSION
We analyzed the interaction of the two interpenetrat-ing, cold, quasi-neutral and magnetized ion-electron plas-mas of comparable densities. We considered the plasmaflowing quasi-parallel to the lines of the background mag-netic field B , with a constant velocity v through an-other static plasma. The linear fluid equations showedthat, if v becomes super-Alfv´enic, the parallel R and Lwaves as well as the perpendicular X and O waves allbecome unstable. We showed that the instability growsonly within a certain range of wavelengths, where thegrowth rate has a maximum. For comparable densitiesof the two plasmas, this maximum resides in the range ofthe resonant wavelengths. To understand what happensin the first stages of the non-linear regime, we startedwith the case of two colliding plasmas that have equaldensities. We combined the linear theory with the re-sults of the test particle simulations and we found thatdue to resonance, the instability (which grows to the am-plitude ∼ B ) becomes able to capture the ions on itspropagation path and, thus, trigger shock formation.We ran PIC simulations of quasi-parallel(sub)relativistic collisionless shocks, considering dif-ferent ion masses and shock velocities. At the very firstphases of a shock formation, we detected a resonantwave of the amplitude ∼ B and a wavelength ∼ λ res . Aswe explained in Sec. III, this wave largely matches theresonant instability that we studied in Sec. II. Duringthis very short transient, we observed that ions arestrongly scattered by the wave. Therefore, it seems thatinitially the resonant instability disrupts the flow andtriggers thermalization, and then shortly afterwards, a shock is formed in the non-linear stage.Ions initially not captured by the resonant mode (thosethat reside in part of the flow which is not encompassedby the grown instability), thus manage to escape andform the initial return current. This sparse bunch of ionplasma excites a wave of opposite polarization ahead ofthe forming shock. Because this wave quickly grows toa non-linear instability that advects toward the shock,it is the first one to trigger the shock reformation pro-cess. Reformation quickly takes over, and it further bal-ances the coupling between the upstream wave and thereturn current. Polarization of the upstream wave (whichis commonly observed in kinetic simulations) is such thatit is non-resonant with respect to the upstream flow. It isknown as the cosmic-ray-induced streaming instability ,whose non-resonant mode is also referred to as Bell’s in-stability , or hybrid non-resonant instability . At thesame time, we find that polarization of this wave is res-onant with respect to the gyro-motion of the outflowingions. Based on this, we related the upstream wave tothe R-mode micro-instability (positive helicity whistlerwave) that is studied here in the weak-beam case and isalready known as EM ion/ion right-hand resonant insta-bility .In Sec. III, we explained how this self-sustaining insta-bility governs the shock reformation process in the caseof a quasi-parallel shock. The inflowing upstream plasmais non-resonantly perturbed by the instability, and refor-mation is triggered in the region where the wave grows toa few B in amplitude. The return beam of ions becomesperiodically disrupted and slowed down by the resonantmechanism. Nevertheless, one part of the returning ionscontinues to propagate, thus expanding the acceleration4 c / ω pi ω c i · t
75 100 125 150
FIG. 9. The density profile captured at different times showsan ongoing shock reformation process. On the horizontal axis,the x -coordinates are given in units of the ion skin depth( c/ω pi ); and on the vertical axis, the time is given normalizedto the ion gyro-period T ci = 2 π/ω ci . region for the DSA mechanism. We used the propertiesof the upstream wave that we observed at a later stage,when the wave-beam coupling is stabilized. In Sec. III,we showed that in the region further away from the shock,the weak-beam case can be applied to the return cur-rent. In the calculation explained in Sec. III D, where weused the linear theory presented in Sec. II, we obtainedthe value of ∼ .
8% for ion injection. For comparison,in the minimal model , the fraction of injected ions isfound to be ∼ .
6% of the inflowing ions.The wave observed at the shock is due to advectionand compression of the upstream instability. The shockis reformed at the first peak of the compressed upstreamwave. Afterwards, the wave slows down and couples to adownstream structure, which is of opposite polarizationand of a significantly different wavelength. We found thatthe wavelength of the downstream instability is highlyresonant ( ∼ πr gi ). In Sec. II, using linear equationswe related this resonance to the case of flow compres-sion of ∼
4. Therefore, the compression ratio similarto the one from Rankine-Hugoniot conditions can alsoarise due to resonant wave-plasma coupling during shock(re)formation. Nevertheless, because the shock is heav-ily non-linear, the results of the linear theory cannot be applied with complete certainty. On the other hand, wefound that for Mach numbers (cid:46)
5, the shock forms quitegradually. The upstream wave remains quasi-linear for alonger period of time, which postpones the beginning ofshock reformation. We found that the instability, whichtriggers the initial particle scattering, does not change itswaveform throughout shock formation. The peak in E x profile (which is clearly associated with the overall slow-ing down of the flow), and the boost in the wave phase(which also appears at the same point), thus clearly dis-tinguish the upstream and the downstream wave in thestage before the beginning of shock reformation.As in , we observe that shock reformation happenson the scales of the instability wavelength. However, weshowed that the upstream wave is amplitude modulated.The envelope of the modulated wave changes with thewaveform of the downstream mode. Therefore, we thinkit is very likely that the downstream instability modu-lates the intensity of the return current, which in returnmodulates the amplitude of the upstream mode. Con-sequently, this leads to the amplification or suppressionof the whole process on a larger scale, as the envelopeincreases or decreases with the waveform of the down-stream instability.The composite picture we get here, shows certain simi-larities to the model presented in . Firstly, the Alfv´enictype wave is excited in the upstream. Secondly, the insta-bility constrains the escape of ions by resonantly trappingand advecting them back to the shock. However, pertur-bation of the upstream flow by the wave and, therefore,the shock reformation process are not accounted therein.Kinetic simulations show that a return current isformed by the population of reflected ions that escape thereforming barrier. Even more important is that the ma-jor ion population cannot thermally leak , not evenfrom the near downstream. The more energetic, heav-ier nuclei can, however, diffuse from the far downstream and become re-accelerated while crossing the shock, asexplained in . V. SUMMARY
In addition to the existing models, in this paperwe want to disclose the significance of resonant micro-instabilities and their possible role in the triggering,transmission and reformation of a shock wave. We give asummary of our results found in the theory and observedin PIC simulations:1. From the linear theory, we found that resonant cir-cularly polarized micro-instabilities can be driven by twointerpenetrating, cold, quasi-neutral and magnetized ion-electron plasmas of comparable densities.2. During the very first phases of a collisionless shockformation, we found that a resonant wave is excited inPIC simulations. By its properties, we related this waveto resonant micro-instability, which is driven by the flowitself (a strong beam). Because of resonance, we found5that the wave strongly scatters the particles and triggersthermalization. Shortly afterwards, shock conditions areformed by non-linear processes.3. Ions initially not captured by a resonant wave, es-cape and form the initial return current. This currentexcites the wave ahead of the forming shock, which isthe first to trigger shock reformation process. We mod-eled the wave by the same type of micro-instability thatwe found in the weak-beam case. The amount of in-jected ions that we got from the equations finely agreeswith the minimal model for ion injection from quasi-periodic shock barrier.4. The downstream wave is likely due to the advectedupstream wave, which changes its wavelength and polar-ization at the shock interface and becomes resonant withthe upstream flow. The linear theory does not apply be-cause of the high non-linearity of the wave. Still, ourequations indicate that the flow compression of ∼ ACKNOWLEDGMENTS
Gratitude for this work I owe to my advisor BojanArbutina. For valuable discussions and suggestions, Ialso thank Dejan Uroˇsevi´c and Marko Pavlovi´c. For hergreat help in preparing this paper, I thank Aleksandra´Ciprijanovi´c. Special thanks to both reviewers for theirvery constructive comments. PIC simulations were runon PARADOX-IV supercomputing facility at Scientific Computing Laboratory of the Institute of Physics Bel-grade, and also on cluster “Jason” of Automated Rea-soning Group (ARGO) at the Department of ComputerScience, Faculty of Mathematics, University of Belgrade.The results of PIC simulations were visualized by “Iseult”- a GUI written by Patrick Crumley. This research hasbeen supported by the Ministry of Education, Scienceand Technological Development of the Republic of Ser-bia under project No. 176005.
Appendix A: Equations of cold interpenetrating plasmas
In the regime of small oscillations, the first order per-turbation equations for the flowing plasma in the fre-quency domain become:( − iω + i k · v ) v i = q i m i ( E + v i × B + v × B ) , (A1)( − iω + i k · v ) v e = q e m e ( E + v e × B + v × B ) . (A2)In the zeroth order, Eqs. (1)–(2) relate the equilibriumquantities as E + v × B , (A3)and show us that the flow of particles through the con-stant magnetic field induces a constant electric field,which then opposes the magnetic force imposed on themoving charges. The flow remains undisturbed and theparticles continue to drift through the background mag-netic field.For the stationary plasma ( v = 0), Eqs. (A1) and(A2) become: − iω v i = q i m i ( E + v i × B ) , (A4) − iω v e = q e m e ( E + v e × B ) . (A5)For simplicity, from now on, index “1” is dropped fromthe symbols for variable EM field components E and B .Perturbation velocities of stationary plasma particles v si,e are obtained from Eqs. (A4)–(A5) and can be ex-pressed as functions of fields: v si,e = q i,e m i,e − iω iωω − ω ci,e − ω ci,e ω − ω ci,e ω ci,e ω − ω ci,e iωω − ω ci,e · E , (A6)where cyclotron frequency of the ions/electrons is givenby ω ci,e = q i,e B /m i,e .6Similarly, the perturbation velocities v fi,e of the flowingplasma particles are obtained from Eqs. (A1)–(A2). Byusing Faraday’s law of induction ω B = k × E , the com-ponent of these equations, which contains a perturbedmagnetic field, can be written as v × B = − ( v · k ) ω E + ( v · E ) ω k = − ( v · k ) ω E + kv T0 ω E . (A7)By substituting this into the equations of the flowingplasma, the perturbation velocities are found to be v fi,e = q i,e m i,e − iξω iξωξ ω − ω ci,e − ω ci,e ξ ω − ω ci,e ω ci,e ξ ω − ω ci,e iξωξ ω − ω ci,e ·· (cid:18) ξ + kv T0 ω (cid:19) E , (A8)where the term that contains the tensor product is de-fined by kv T0 = k x v x k x v y k x v z k y v x k y v y k y v z k z v x k z v y k z v z . (A9)The frequency of the wave is now modified by the col-lective motion of plasma particles, where modification ismade through the parameter ξ = 1 − v · k ω . (A10)The Current is then found by summing up the contri-butions from all of the species of the flowing and station-ary plasmas: j = β = f,s (cid:88) α = i,e n βα q α v βα + (cid:88) α = i,e n f α q α v , (A11)where n f i,e stands for the density perturbation of thedifferent species of the flowing plasma. It satisfies thelinearized continuity equation − iξωn f i,e + in i,e k · v fi,e = 0 , n f i,e = n fi,e ξω k · v fi,e . Using standard procedure for obtaining the disper-sion relation, Amp`ere and Faraday laws are combined toget a linearized equation of the plasma-wave coupling k × ( k × E ) + ω c E + iωµ j = 0 . (A12)After substituting Eqs. (A6) and (A8) for velocities inEq. (A11), and by using that current, Eq. (A12) can berewritten as k × ( k × E )+ ω c K · E = 0 . (A13)Polarization matrix K can further be separated into twomatrices: K = κ + i(cid:15) ω σ f . (A14)The first matrix has the form as in the flowless case: κ = P S − iD iD S , (A15)but its parameters are modified by the flow through ξ : P = 1 − β = f,s (cid:88) α = i,e ω pαβ ω , (A16) S = 1 − (cid:88) α = i,e ξ ω pαf ξ ω − ω cα + ω pαs ω − ω cα , (A17) D = (cid:88) α = i,e ω pαf ω ξω cα ξ ω − ω cα + ω pαs ω ω cα ω − ω cα . (A18)The second matrix fully depends on v and vanishes com-pletely if the flow speed is zero: σ f = (cid:15) (cid:88) α = i,e ω pαf · (cid:18) v k T ξω (cid:19) − iξω iξωξ ω − ω cα − ω cα ξ ω − ω cα ω cα ξ ω − ω cα iξωξ ω − ω cα kv T0 ω ++ v k T ω − iξω iξωξ ω − ω cα − ω cα ξ ω − ω cα ω cα ξ ω − ω cα iξωξ ω − ω cα . (A19)Taking the wavevector to be in the x − z plane, thedispersion matrix that constitutes a system of equationsgiven by (A13) is7 D = P − n sin θ n sin θ cos θ S − n − iDn sin θ cos θ iD S − n cos θ ++ i(cid:15) ω σ f . (A20) Appendix B: Instability of R and L waves
To derive ω r and γ , the two frequency ranges are con-sidered: cases when ξω ∼ ω ci (cid:28) ω ce and ω ci (cid:28) ω ce ∼ ξω .Eq. (7) is at first rearranged as R ( L ) = 1 − η ξω · (cid:34) ω pis ξω ± ω ci + ω pes ξω ± ω ce (cid:35) − ω · (cid:34) ω pis ω ± ω ci + ω pes ω ± ω ce (cid:35) . (B1)It is then approximated by the use of previous conditionsfor the ion-cyclotron and electron-cyclotron frequencies.
1. Ion-cyclotron frequencies
Condition ξω ∼ ω ci (cid:28) ω ce implies the interactionbetween ions of the flowing and stationary plasmas.Eq. (B1) is therefore approximated by R ( L ) ≈ − η ξω · (cid:34) ω pis ξω ± ω ci ± ω pes ω ce (cid:35) − ω · (cid:34) ω pis ω ± ω ci ± ω pes ω ce (cid:35) . (B2)Equality ω pes /ω ce = − ω pis /ω ci is then used to get: R ( L ) ≈ ω pis (cid:20) η ξ ω ci ( ω ci ± ξω ) + 1 ω ci ( ω ci ± ω ) (cid:21) , (B3)where ξ is given by Eq. (A10).This is the cold plasma relation as given in . Inthe range of frequencies much lower than ω ci , these twocircularly polarized modes merge to form an Alfv´en wavethat is modified by the flow. Therefore, in approximation ω (cid:28) ω ci , Eq. (B3) reduces to n ≈ c v (cid:20) η (cid:16) − v c n (cid:17) (cid:21) , (B4)and after the substitution of n = k c /ω , it becomes ω k (cid:18) n f n s + v c (cid:19) − n f n s v ωk + n f n s v − v = 0 . (B5)This is the same dispersion relation as in Ref. , exceptfor the term V /c that here appears within the firstbracket and is caused by the term ω /k E of Eq. (A12),which is neglected in in the expression of the Amp`erelaw.Equation (B1) therefore represents the generalizationof the case where at low frequencies Alfv´en wave coupleswith the flow driven ( ω = kv ) mode . Eq. (B3) canbe solved analytically for ion/ion interactions. By thesubstitution of R ( L ) = n R /L = k c /ω in Eq. (B3),quartic polynomial equation is obtained: a i ω + b i ω + c i ω + d i ω + e i = 0 , (B6) a i = 1 ω ci v c ,b i = ± ω ci (cid:20) η + v c (cid:18) ∓ k · v ω ci (cid:19)(cid:21) ,c i = (cid:18) η + v c (cid:19) (cid:18) ∓ k · v ω ci (cid:19) ∓ η k · v ω ci − k v ω ci ,d i = (cid:18) η k · v ± k v ω ci (cid:19) (cid:18) − ± k · v ω ci (cid:19) ,e i = η ( k · v ) − k v (cid:18) ∓ k · v ω ci (cid:19) . General solutions to this equation are found as ω , = − b i a i − S i ± (cid:114) − S i − p i + q i S i , (B7) ω , = − b i a i + S i ± (cid:114) − S i − p i − q i S i , (B8)where coefficients and determinants are given by p i = 8 a i c i − b i a i ,q i = b i − a i b i c i + 8 a i d i a i ,S i = 12 (cid:115) − p i + 13 a i (cid:18) Q i + ∆ Q i (cid:19) ,Q i = (cid:115) ∆ + (cid:112) ∆ − , ∆ = c i − b i d i + 12 a i e i , ∆ = 2 c i − b i c i d i + 27 b i e i + 27 a i d i − a i c i e i .
2. Electron-cyclotron frequencies
When frequencies around the electron-cyclotron orhigher ( ξω ∼ ω ce (cid:29) ω ci ) are considered, and the con-dition in which the electron plasma frequency is muchhigher than the ion one ( ω pe (cid:29) ω pi ) is assumed, Eq. (B1)can be approximated by the expression: R ( L ) ≈ − ω pes ω ( ω ± ω ce ) − η ξ ω pes ξω ( ξω ± ω ce ) . (B9)This quartic polynomial equation can be rewritten inthe form similar to the case of the ion-cyclotron frequen-cies: a e ω + b e ω + c e ω + d e ω + e e = 0 , (B10) a e = 1 ω ce v e c ,b e = 1 ω ce (cid:20) − v e c (cid:18) k · v ω ce ∓ (cid:19)(cid:21) ,c e = v e c (cid:18) ∓ k · v ω ce (cid:19) − (cid:18) η + k v e ω ce (cid:19) ,d e = ( k · v ∓ ω ce ) (cid:18) η + k v e ω ce (cid:19) ∓ k v e ω ce ,e e = ± η k · v ω ce − k v e (cid:18) ∓ k · v ω ce (cid:19) . Solutions to this equation are found by the same pat-tern as given in Eqs. (B7, B8) using the coefficients withthe subscript e instead of i . Contrary to the case ofions, instability due to the bulk flow of electrons doesnot arise here, and they are therefore not subjected tofurther analysis. Appendix C: Instability of O and X waves
We distinguish the two cases for the waves that prop-agate normal to B . If the flow is parallel to thewavevector ( v = v ˆz ), the polarization matrix given byEq. (A20) takes the form: D = P − n S − n − iD iD S + (C1)+ (cid:88) α = i,e ω pαf ω · − i ω cα kv ξ ω − ω cα i ω cα kv ξ ω − ω cα − ωkv ξ ω − ω cα ( ξ + 1) . If the flow is perpendicular to the wavevector and par-allel to B ( v = v ˆx ), the scalar product k · v vanishesand ξ = 1, which reduces the Eq. (A20) to D = P − n S − n − iD iD S + (C2)+ (cid:88) α = i,e ω pαf ω · − k v ω − ω cα i ω cα kv ω − ω cα − ωkv ω − ω cα − i ω cα kv ω − ω cα − ωkv ω − ω cα .
1. Ordinary electromagnetic wave
From previous relations, it can be seen that the O-mode is modified by the flow only if v (cid:107) B , whichis described by the first root of Eq. (C2). For simplic-ity, we analyze this wave in the center-of-momentum ref-erence frame. This means that we have the two op-posite streaming plasmas, which satisfy the condition ρ v = − ρ v . Plasma frequencies and velocities are thusrelated as ω p = ηω p and v = − v /η , respectively. Thenon-diagonal terms in Eq. (C2) are therefore cancelled,and the first root simplifies: n = P − (cid:88) α = i,e ω pα ω k v ω − ω cα (cid:18) η (cid:19) . (C3)If we consider unmagnetized plasmas ( ω cα = 0),Eq. (C3) reduces to a relation that has a purely com-plex solution. It represents the growth rate of Weibelinstability , which is commonly observed in shock sim-ulations . In the case of magnetized flows, if ω ∼ ω ci ,Eq. (C3) simplifies to a biquadratic equation ω − ω ( ω p + k c + ω ci ) − (C4) − ω pi k v + ω ci ( ω p + k c ) = 0 , where ω p = (cid:80) i,e (1+ η ) ω pi and ω pi = (1+1 /η ) ω pi . Thesolution for the unstable wave mode, we find as ω = − ω pi v c
11 + ω p k c − (cid:18) η (cid:19) M . (C5)For B = 0, this solution necessarily becomes complex,and takes the known form of the growth rate of Weibelinstability. However, in the case of magnetized plasmas,for ω ∼ ω ci , Eq. (C5) implies that the stability criterionnow depends on the Alfv´enic Mach number. For a given M A , there is a cutoff wavenumber k cut ≈ r − gi (cid:112) η m i /m e ,at which the wave becomes unstable.9
2. Extraordinary electromagnetic wave
If the flow is along the magnetic field lines, the secondroot of Eq. (C2), as given in the center-of-momentumframe, is ( S − n ) S − D = 0 , n = RLS .
This mode is a simple elliptically polarized EM wave,with the electric field vector in a direction perpendicularto the magnetic field lines ( E ⊥ B ).If the flow is directed parallel ( v = v ˆz ) to thewavevector, the X-mode becomes modified. Its disper-sion relation is given by the second root of Eq. (C1):( S − n ) S M − D M = 0 , n = RL − SM S − DM D − M D S − M S , where S M and D M are expressed as S M = S − (cid:88) α = i,e ω pαf ω ωkv ξ ω − ω cα ( ξ + 1) ,D M = D + (cid:88) α = i,e ω pαf ω ω cα kv ξ ω − ω cα . Amato, E. and Blasi, P., Mon. Not. R. Astron. Soc. , 1591(2009). Anderson, D., Fedele, R., and Lisak, M., Am. J. Phys. , 1262(2001). Axford, W. I., Leer, E., and Skadron, G., in , Vol. 11 (Bulgarian Academy of Sciences,Conference Papers, Plovdiv, Bulgaria, 1977) pp. 132–137. Bell, A. R., Mon. Not. R. Astron. Soc. , 147 (1978). Bell, A. R., Mon. Not. R. Astron. Soc. , 550 (2004). Belmont, G., Grappin, R., Mottez, F., Pantellini, F., and Pel-letier, G.,
Collisionless Plasmas in Astrophysics (WILEY-VCH,2014). Blandford, R. D. and Ostriker, J. P., Astrophys. J. , L29(1978). Caprioli, D., Astrophys. J. Lett. , L38 (2015). Caprioli, D., in
Cosmic Ray Origin - Beyond the Standard Mod-els, San Vito di Cadore, 2016 , Vol. 297-299, edited by O. Tibolla,R. Blandford, S. Kaufmann, M. Persic, K. Mannheim, H. Voelk,and A. D. Angelis (Elsevier, Nucl. Part. Phys. Proc., Radarweg29, 1043 NX Amsterdam, The Netherlands, 2018) pp. 226–233. Caprioli, D., Pop, A.-R., and Spitkovsky, A., Astrophys. J. Lett. , L28 (2015). Caprioli, D. and Spitkovsky, A., Astrophys. J. Lett. , L20(2013). Caprioli, D. and Spitkovsky, A., Astrophys. J. , 91 (2014). Caprioli, D. and Spitkovsky, A., Astrophys. J. , 46 (2014). Caprioli, D., Zhang, H., and Spitkovsky, A., J. Plasma Phys. , 1 (2018). Chen, G. and Armstrong, T. P., in , Vol. 5 (Max-Planck-Institut f¨ur extrater-restrische Physik, Munich, West Germany, 1975) pp. 1814–1819. Dieckmann, M. E., Drury, L. O., and Shukla, P. K., New J.Phys. , 40 (2006). Dokuchaev, V. P., Sov. Phys. J. Exp. Theor. Phys. , 292(1961). Drury, L. O., Rep. Prog. Phys. , 973 (1983). Fitzpatrick, R.,
Plasma Physics: An Introduction (CRC Press,2014). Gary, S. P., Space Sci. Rev. , 373 (1991). Gary, S. P.,
Theory of Space Plasma Microinstabilities (Cam-bridge University Press, 1993). Gary, S. P., Smith, C. W., Lee, M. A., Goldstein, M. L., andForslund, D. W., Phys. Fluids , 1852 (1984). Guo, F. and Giacalone, J., Astrophys. J. (2013). Gurnett, D. A. and Bhattacharjee, A.,
Introduction to PlasmaPhysics: With Space and Laboratory Applications (CambridgeUniversity Press, 2005). Kan, J. R., Lyu, L. H., and Mandt, M. E., Space Sci. Rev. ,201 (1991). Krymskii, G. F., Sov. Phys.Dokl. , 327 (1977). Lapuerta, V. and Ahedo, E., Phys. Plasmas , 1513 (2002). Lee, R. E., Chapman, S. C., and Dendy, R. O., Astrophys. J. , 187 (2004). Malkov, M. A., Phys. Rev. E , 4911 (1998). Marcowith, A., Bret, A., Bykov, A., Dieckman, M., Drury, L.,Lembege, B., Lemoine, M., Morlino, G., Murphy, G., Pelletier,G., Plotnikov, I., Reville, B., Riquelme, M., Sironi, L., and Novo,A. S., Rep. Prog. Phys. , 1 (2016). Matsukiyo, S. and Scholer, M., Adv. Space Res. , 57 (2006). Sironi, L. and Spitkovsky, A., Astrophys. J. , 75 (2011). Spitkovsky, A., in
Astrophysical Sources of High Energy Particlesand Radiation , Vol. 801 (AIP Conference Proceedings, Melville,NY, 2005) pp. 345–350. Spitkovsky, A., Astrophys. J. Lett. , L39 (2008). Vedenov, A. A., Velikhov, E. P., and Sagdeev, R. Z., Nucl. Fusion , 82 (1961). Vranjes, J., Phys. Plasmas , 5 (2015). Wang, X. Y. and Lin, Y., Phys. Plasmas , 3528 (2003). Weibel, E. S., Phys. Rev. Lett. , 83 (1959). Zekovi´c, V. and Arbutina, B., in