Formation of foreshock transients and associated secondary shocks
Xin An, Terry Z. Liu, Jacob Bortnik, Adnane Osmane, Vassilis Angelopoulos
DDraft version August 27, 2020
Typeset using L A TEX default style in AASTeX63
Formation of foreshock transients and associated secondary shocks
Xin An, Terry Z. Liu,
2, 3
Jacob Bortnik, Adnane Osmane, and Vassilis Angelopoulos Department of Atmospheric and Oceanic Sciences, University of California, Los Angeles, California 90095, USA. Department of Earth, Planetary and Space Sciences, University of California, Los Angeles, California 90095, USA. University Corporation for Atmospheric Research, Boulder, Colorado 80307, USA. Department of Physics, University of Helsinki, Helsinki 00014, Finland. (Received November 20, 2019; Revised July 28, 2020; Accepted August 7, 2020)
Submitted to ApJABSTRACTUpstream of shocks, the foreshock is filled with hot ions. When these ions are concentrated andthermalized around a discontinuity, a diamagnetic cavity bounded by compressional boundaries, re-ferred to as a foreshock transient, forms. Sometimes, the upstream compressional boundary can furthersteepen into a secondary shock, which has been observed to accelerate particles and contribute to theprimary shock acceleration. However, secondary shock formation conditions and processes are notfully understood. Using particle-in-cell simulations, we reveal how secondary shocks are formed. From1D simulations, we show that electric fields play a critical role in shaping the shock’s magnetic fieldstructure, as well as in coupling the energy of hot ions to that of the shock. We demonstrate that largerthermal speed and concentration ratio of hot ions favors the formation of a secondary shock. From amore realistic 2D simulation, we examine how a discontinuity interacts with foreshock ions leading tothe formation of a foreshock transient and a secondary shock. Our results imply that secondary shocksare more likely to occur at primary shocks with higher Mach number. With the secondary shock’spreviously proven ability to accelerate particles in cooperation with a planetary bow shock, it is evenmore appealing to consider them in particle acceleration of high Mach number astrophysical shocks.
Keywords:
Plasma astrophysics, Shocks, Planetary bow shocks INTRODUCTIONAlthough shocks are among the most powerful particle accelerators in the universe, how they accelerate particles toextreme energies, e.g., generate cosmic rays, is still not fully understood (see review in Treumann 2009). Recent in-situobservations at the Earth’s bow shock have shown that nonlinear structures in the foreshock can play an important rolein shock acceleration (Liu et al. 2016; Wilson III et al. 2016; Liu et al. 2017a,b, 2018; Turner et al. 2018; Liu et al. 2019).A region upstream of the shock, the foreshock, is magnetically connected to the shock and filled with particles comingfrom it (Eastwood et al. 2005). The foreshock is very dynamic and there are many nonlinear transient structureswith large plasma and field fluctuations referred to as foreshock transients. Because of the density and magnetic fieldvariations and/or plasma deflection, foreshock transients can result in dynamic pressure perturbations. When foreshocktransients convect towards and connect with the planetary bow shock, the bow shock surface can be disturbed. Suchperturbation can propagate to the magnetosheath, the magnetopause and thus the magnetosphere (e.g., Sibeck et al.1999; Turner et al. 2011). Hot flow anomalies (HFAs) (Schwartz et al. 1985; Thomsen et al. 1986, 1988; Omidi &Sibeck 2007) and foreshock bubbles (FBs) (Omidi et al. 2010; Turner et al. 2013) are two most significant types offoreshock transients due to their large spatial scale (from several foreshock ion gyroradii to almost the entire widthof the ion foreshock) and the significant plasma perturbations they entail. Foreshock transients have a lifetime of
Corresponding author: Xin An, Terry Z. [email protected], [email protected] a r X i v : . [ phy s i c s . s p ace - ph ] A ug An et al. several minutes. Their main characteristic is a hot, diamagnetic cavity surrounded by a compressional boundary or asecondary shock [Figs. 1(a) and 1(b)].Recent observations found that foreshock transients, especially HFAs and FBs, can accelerate particles and contributeto the primary shock acceleration (e.g., Wilson III et al. 2013; Liu et al. 2017a). As foreshock transients convect withthe upstream flow, particles enclosed within their boundary and the primary shock can experience Fermi acceleration(Liu et al. 2017b, 2018; Turner et al. 2018). Secondary shocks have also been observed to accelerate upstream particleson their own through the shock drift acceleration and even to form a secondary foreshock (Liu et al. 2016). Secondaryshocks can also capture and further energize primary shock-accelerated electrons through betatron acceleration (Liuet al. 2019). Recent observations also identified magnetic reconnection inside foreshock transients contributing tothe electron acceleration/heating (Liu et al. 2020). These additional accelerations by foreshock transients and theirsecondary shocks provide a possible solution to Fermi’s ‘injection problem’ (an unresolved seed population of energeticparticles for further acceleration (Treumann 2009)) and increase the acceleration efficiency of primary shocks. Todate foreshock transients have only been observed at planetary bow shocks where in-situ observations are available(e.g., Masters et al. 2009; Turner et al. 2013; Collinson et al. 2015). Whether foreshock transients exist at high Machnumber astrophysical shocks, such as supernova-driven shocks, is unknown or can only be inferred from simulations(e.g., Giacalone & Burgess 2010). Therefore, to infer their existence and contribution at other shock systems, it isimportant to fully understand and quantify the formation process of foreshock transients.Based on various simulations, conceptual models have been proposed for the formation of HFAs and FBs. Testparticle simulations by Burgess (1989) show that when convection electric field points towards the tangential discon-tinuity (TD) that intersects a shock, specularly reflected ions either traveling away upstream (quasi-parallel shock) ornot (quasi-perpendicular shock) can be trapped by the TD and channeling along it. Hybrid simulations (e.g., Thomaset al. 1991; Lin 2002; Omidi & Sibeck 2007) show that such scenario can lead to the formation of HFAs. In order tohave sufficient time to trap foreshock ions and form an HFA, TDs need to move along the bow shock surface very slowlyas shown in the statistical study by Schwartz et al. (2000). Additionally, hybrid simulations by Omidi et al. (2010)show that when field-aligned foreshock ions cross a rotational discontinuity (RD), an FB can form. To understandhow the interaction between foreshock ions and a TD/RD leads to the formation of an HFA or FB, Archer et al.(2015) and Liu et al. (2015) suggest that, because of the magnetic field direction change across the discontinuity, partof the foreshock ion’s parallel speed projects onto the perpendicular direction increasing the thermal energy in theperpendicular direction (thermalization). The decrease of the parallel speed also results in the increase of the foreshockion density (concentration), due to conservation of mass flux. The thermalization and concentration of foreshock ionsincrease their thermal pressure resulting in an expansion and formation of a foreshock bubble or a hot flow anomaly.However, such a model is too qualitative. For example, because foreshock ion gyradii are comparable to the systemsize, the concept of thermal pressure of foreshock ions is invalid. The role of electrons in the formation process isneglected. Under what condition there will be a secondary shock is unknown.The rest of the paper is organized as follows. In Section 2, we will study a simple initial value problem using particle-in-cell (PIC) simulations in one spatial dimension. In this configuration, we initialize a layer of hot ions and simulatethe plasma expansion perpendicular to the magnetic field and the resulting formation of foreshock transients. We willinvestigate the detailed roles of hot foreshock ions, cold ambient ions, electrons, and the associated electromagneticfield during the formation process. We will further reveal the critical parameters that determine the secondary shockformation. In Section 3, we will study a more realistic boundary value problem using a PIC simulation in two spatialdimensions. In this configuration, we continuously inject energetic ions at a boundary and simulate the interaction ofinjected ions with a rotational discontinuity and the formation of a foreshock transient and a secondary shock. Wewill compare the 2D simulation with 1D simulations and show effects of the additional spatial dimension. In Section4, we will summarize and discuss our results. EXPANSION PERPENDICULAR TO THE MAGNETIC FIELD: 1D PIC SIMULATIONSTo explore the foreshock transient formation, we carried out a series of PIC simulations using the osiris code(Fonseca et al. 2002; Hemker 2015), including 1 spatial ( x ) and 3 velocity components ( v x , v y , v z ). All the runs werein the rest frame of the upstream flow, so there is no background flow in the initial setup (see Appendix A for detailsof the simulation setup). A uniform background magnetic field B is oriented in the + z direction [Fig. 1(c)]. Theinitial plasma density n is uniform in the computation domain. In the foreshock region, when ions backstreamingfrom the primary shock encounter certain discontinuities in the upstream flow, they are concentrated and thermalized oreshock transients and secondary shocks (cid:54) x (cid:54) ρ h [Fig. 1(c)], where ρ h is the gyroradius corresponding to the initial thermalvelocity v T h of hot ions. Inside this layer, the densities of hot ions, ambient ions and electrons are 0 . n , 0 . n , and n , respectively; outside it, there are only ambient ions and electrons, each with density n . The hot ions have aninitial thermal velocity of v T h = 9 v A ; the ambient ions and electrons have initial thermal velocities of v T i = 0 . v A and v T e = 3 v A , respectively. Here the Alfv´en velocity is v A = B / √ πn m i , and the speed of light is c = 300 v A .The highly nonequilibrium state of the initial setup is consistent with the observed characteristics of different particlespecies at an early stage (less than an ion cyclotron period from the birth of the hot ion core) of foreshock transientformation (Turner et al. 2013; Omidi et al. 2010). Given the reduced ion-to-electron mass ratio m i /m e = 100, theordering of the typical gyroradii of the three species is ρ h = 30 ρ i = 300 ρ e . The wide gap between the gyroradii of thehot ions and those of the ambient ions and electrons, which captures the essential ordering in spacecraft observations,will have a significant effect on the evolution of the system. Later we also discuss runs that vary the initial thermalvelocity and the concentration ratio of hot ions. Figure 1.
A foreshock transient and its associated secondary shock. (a) A sketch of the foreshock transient upstream of the bowshock. On the right side of the discontinuity (gray board), foreshock ions (yellow arrows) move along magnetic field lines (bluearrows). Because of their large gyroradius, some foreshock ions gyrate across the discontinuity instead (curved yellow arrow),are trapped and become concentrated and thermalized (orange region), resulting in a fast expansion that forms a secondaryshock (purple surface). (b) A typical observation of a foreshock transient upstream of the Earth’s bow shock (Turner et al.2013; Wilson III et al. 2016; Liu et al. 2019). From top to bottom are magnetic field, ion bulk velocity in the solar wind restframe, plasma density, and ion energy flux spectrum in the Earth’s rest frame. Here the x direction is defined as the normaldirection of the secondary shock; the z direction is defined as the direction of maximum variation of the magnetic field; and the y direction completes the right-handed coordinate system (as sketched in (a)). The colors of the shaded regions correspond tothose shown in (a). (c) A sketch of the simulation setup. The coordinate system represents the geometry in (a) and (b). Formation process
Here we explain the details of the foreshock transient formation process. In our simulations, the concentrated hotions drive the system to form a diamagnetic cavity and a secondary shock [Fig. 2]. A fundamental question is, howis the energy of the hot ions transferred to the magnetic field and ambient plasma flows tied to understanding ofelectric fields, which control energy transfer between particles and fields. And, indeed, two types of electric fields wereidentified in the simulations: electrostatic fields in the − x direction [Fig. 2(c)] and induction electric fields in the + y direction [Fig. 2(e)]. As the energetic ions move out of the hot layer across the magnetic field in the first gyration(0 < t < τ ci , τ ci being the ion cyclotron period in B ), an electrostatic field E x directed from these displaced hot ionsto the initial hot layer (i.e., E x <
0) is generated. Correspondingly, the electrostatic potential decreases in the hot layerand is enhanced outside it [Fig. 2(c)]. This electrostatic field causes a drift u y,e = − cE x /B z of electrons, an electroncurrent in the − y direction [Figs. 2(d) and 1(c)]. The electron current and the hot ion (partial diamagnetic) current(due to the density gradient and partial gyration of hot ions at the interface), both flowing in the − y direction, cannotbe offset by the ambient ion (partial diamagnetic) current (flowing in the opposite, + y , direction). The resulting netcurrent, Hall current, (see Appendix B for the ion diamagnetic current and total current), reduces the magnetic fieldon one side (the cavity) and enhances it on the other side (the compressional boundary) [Fig. 2(a)]. Simultaneously,an induction electric field E y is generated by these magnetic field variations [Figs. 2(e) and 1(c)], in order to keep An et al. the electrons drifting across the magnetic field region extending in the + x direction, encompassing the hot ions thatgyrate out. Ambient ions, whose gyration is completely immersed in the interaction region of 2 ρ h due to their smallergyroradius, also gain momentum in the + x direction via the induction electric field, eventually forming a streamingflow with velocity u x,i = cE y /B z (see Fig. 2(f) and Supplemental Video 1 ). Their acceleration is similar to the ionpick-up process at a comet (e.g., Biermann et al. 1967; Gloeckler et al. 1986). As a result, the ambient plasma, aswell as the magnetic flux, are transported outward, and thus the density and magnetic field strength are depleted onone side, forming a cavity and piled up on the other side, forming the compressional boundary (see Figs. 2(a), 2(b)and Supplemental Video 2 ). This cycle repeats until the compressional boundary finally detaches from the hot ionsat 3 τ ci (cid:46) t < τ ci (see Supplemental Video 3 ). As explained below, the energy exchange between the fields and thehot ions ceases at the time of this detachment. The density of hot ions decreases in surrounding space, but remainsconcentrated in the cavity region. The compressional boundary, which continues to move at supermagnetosonic speed,steepens into a shock with the Mach number M f = 2 . M f , defined as the ratio of the upstream flow speed in theshock normal incidence frame to the local fast magnetosonic speed v f ). x / ρ h B z / B (a)0 2 4 6t / τ ci x / ρ h n e / n (b) −1 0 1 φ [x10 −2 ](c)0 2 4 6t / τ ci −0.4 0.0 0.4 J y,e [x10 −2 ](d) E y [x10 −3 ](e)0 2 4 6t / τ ci u x,i [x10 −2 ](f) Figure 2.
The formation of a foreshock transient and a secondary shock. (a)-(f) The spatio-temporal evolution of relevantfield and particle quantities. The position x is normalized to ρ h . The time t is normalized to τ ci . The arrows next to (a) and(b) indicate the original boundary between hot and ambient plasmas. (a) The magnetic field B z in the z direction. (b) Theelectron plasma density n e . The evolution of n e is identical to that of B z because their transport equations have the same formin 1D. (c) The electrostatic potential φ . The corresponding electric field, E x = − ∂φ/∂x , resides in the large electron densitygradient between the magnetic cavity and the compressional boundary. (d) The electron current J y,e in the y direction. (e) Theinduction electric field E y . (f) The fluid velocity of the ambient ions u x,i in the x direction. In panels (c)-(f) and hereafter, theelectric potential has the dimension m e c /e , the current density has the dimension n ec , the electric field has the dimension m e c e · c/ω pe , and the velocity has the dimension c . Version 1.0 of Supplemental Video 1 is archived on Zenodo https://doi.org/10.5281/zenodo.3951168. Version 1.0 of Supplemental Video 2 is archived on Zenodo https://doi.org/10.5281/zenodo.3951168. Version 1.0 of Supplemental Video 3 is archived on Zenodo https://doi.org/10.5281/zenodo.3951168. oreshock transients and secondary shocks < t < τ ci , a surge in the rate of work done by the electric field E y and E x on the particles is seen [Figs. 3(a) and 3(b)], which is associated with an energy decrease in the hot ions and anenergy increase in the ambient ions and electrons [Fig. 3(c)]. In this process, the hot ions transfer energy to the fieldsdominantly through the partial gyration against the induction electric field E y ( J y,h · E y < E x on hot ions is relatively small, because the spatial location of E x lags behind that of thehot ion current J x,h . The electron current J y,e is a generator ( J y,e · E y < E x tends to accelerate electrons [Fig. 3(b)], however, leading to a netelectron energy increase [Fig. 3(c)], which is consistent with a recent statistical study showing that electrons are almostalways heated/accelerated inside foreshock transients (Liu et al. 2017a). Accelerated by the induction electric field E y [Fig. 3(a)], ambient ions gain energy in the form of streaming in the + x direction [Fig. 2(f)]. In this way, the inductionelectric field mediates momentum and energy coupling between hot and ambient ions. We also note that in addition,the electrostatic field E x creates the electron current and thus plays a critical role in shaping the magnetic cavity andcompression. These electric fields have also been observed in other scenarios, such as solar wind-barium interactionin the Active Magnetospheric Particle Tracer Explorers (AMPTE) experiment (e.g., Papadopoulos et al. 1987) andlaser-produced plasma expansion in laboratory experiments (e.g., Bonde et al. 2015; Bondarenko et al. 2017). Becausethe decrease in the total kinetic energy of particles is matched by the increase in the magnetic field energy [Fig. 3(d)],global energy is conserved. The energy transfer between fields and particles almost vanishes after t = 4 τ ci , as hot ionsare detached from the outward-propagating compressional boundary (i.e., the region where E y is located). −2024−2024 < J y ⋅ E y > [ x − ] electronsambient ionshot ionstotal(a)−4−2024 < J x ⋅ E x > [ x − ] (b)020406080100120140 E p [ m e c ] (c)0 2 4 6time / τ ci −15−10−5051015 ∆ E [ m e c ] magnetic fieldparticle(d) Figure 3.
Energy exchange between fields and particles. (a),(b) The history of the rate of work done by E x and E y onelectrons (red), ambient ions (blue), hot ions (green), and all particles (black). The spatial average of any quantity Q is definedby (cid:104) Q (cid:105) = L x (cid:82) L x Q ( x ) dx , where L x is the domain size. (c) The kinetic energy of electrons (red), ambient ions (blue), and hotions (green) with respect to time. (d) The net change in total kinetic energy (solid) and magnetic field energy (dashed) withrespect to time. Compared with the magnetic field energy, the electric field energy is negligible. When swept over by the shock, ambient ions from the upstream are accelerated by the finite E y and begin to gyrate inlarge radii. These transmitted ions have gyrating velocities comparable to the sheath flow velocity, as shown in Figure4. The transmitted ions are non-gyrotropic (e.g., Sckopke et al. 1983, 1990; Burgess et al. 1989) and therefore causes An et al. magnetic oscillations along the background magnetic field [Fig. 2(a)], which is consistent with theoretical predictions(Gedalin 2015) and spacecraft observations (Pope et al. 2019). In fact, electric oscillations of both E x (inferred fromthe oscillations of electron current shown in Fig. 2(d)) and E y [Fig. 2(e)] are also present due to the non-gyrotropy oftransmitted ions. Due to gyrophase mixing, the electromagnetic oscillations gradually decrease in amplitude furtherdownstream from the shock ramp. The electromagnetic oscillations significantly disturb the ion flow and heat theambient ions in the sheath region (see Fig. 4), leading to dissipation of the shock structure. In the meantime, magneticflux is slowly transported from the sheath to the cavity, as seen in the return flow ( u x,i <
0) of ambient ions [Figs. 2(f)and 4] and electrons (not shown). Eventually the magnetic field in the cavity is expected to approach its ambientvalue B , i.e., the cavity will vanish. Figure 4.
The phase space portrait of ambient ions at t = 6 . τ ci . Colors denote the phase space density of ambient ions. Atthis instant, the electromagnetic oscillations are located at 9 < x/ρ h <
13, and the return flow from the sheath to the cavity islocated at 2 < x/ρ h <
5. The thin dashed line at v x = 0 indicates which part of the phase space is flowing towards/away fromthe center. Parameter scans
To determine the conditions necessary for secondary shock formation and to connect it to the Mach number of theprimary shock, we perform a parameter scan of the thermal velocity and the concentration ratio of hot ions [Fig. 5].Looking into the evolution of the system in each run (see Appendix C for the magnetic field structure in each run),we see that the large thermal velocities and high concentration ratios of hot ions favor secondary shock formation. Inthe simulations, at the lower limit of hot ion thermal velocity and concentration, weak magnetic bumps, rather thansecondary shocks, appear and propagate approximately at the fast magnetosonic speed of the ambient plasma. Aseither the thermal velocity or the concentration of hot ions is increased, shock structures emerge, and both the Machnumber and the magnetic compression ratio of these shocks increase [Figs. 5(a) and 5(b)]. The relation between theMach number and the magnetic compression ratio for the simulated secondary shocks agrees with that derived fromRankine-Hugoniot relation [Fig. 5(c)]. On the one hand, because the hot ion/foreshock ion speed is proportional tothe upstream flow speed, the ratio of hot ion thermal speed to the fast magnetosonic speed is proportional to the Machnumber of the primary shock (Burgess et al. 2012). The ratio of the foreshock ion density to the incident upstreamion density, on the other hand, also increases with the Mach number of the primary shock up to ∼ . ∼ oreshock transients and secondary shocks v T h / v A M f (a) 0.1 1.0024681012 v T h / v A η B s ho ck / B (b)1 10M f B s ho ck / B v Th /v A = 12, 9, 6, 3, 1(c) Figure 5.
Parameter scan examining the conditions under which secondary shocks form. (a) The Mach number of secondaryshocks and (b) the magnetic compression ratio across the secondary shock, as a function of the hot ion thermal velocity v Th = (12 v A , v A , v A , v A , v A ) for a given concentration ratio of hot ions η = (0 . , . , . , . , . EXPANSION IN TWO SPATIAL DIMENSIONS: INJECTION AT A DISCONTINUITYIn a more realistic configuration, foreshock ions interact with different types of discontinuities (e.g., tangentialdiscontinuities or rotational discontinuities). The resulting HFAs and FBs expand in both perpendicular and par-allel directions. To account for such a scenario, we perform a PIC simulation with two spatial dimensions. Thecomputational setup is shown in Figure 6(a). We use perfectly matched layers as absorbing boundary conditionsfor electromagnetic fields (Vay 2000) and also absorbing boundary conditions for particles. The orientation of thebackground magnetic field B changes from − ◦ to +30 ◦ with respect to the x -axis at x = 6 c/ω pi , which gives arotational discontinuity. The simulation is in the solar wind rest frame, so the discontinuity is not moving. Initially,both ambient ions and electrons have a uniform density n in the computation domain. To make computational costof the 2D simulation affordable, the ion-to-electron mass ratio is m i /m e = 25, and the normalized Alfv´en velocityis v A /c = 1 / v T i = 0 . v A and v T e = 1 . v A ,respectively. Energetic ions are continuously injected from a cathode on the boundary x = 0 (see the blue strip locatedat 51 c/ω pi (cid:54) y (cid:54) c/ω pi in Figure 6(a)). Injected ions stream along B with an initial beam velocity v bi = 4 . v A .The density of injected ions is 0 . n , and their mass is m bi = 100 m e . To maintain the continuous injection of ions,electrons are injected simultaneously with ions, otherwise the initially injected ions will create a space charge potential An et al. that prevents further injection of ions. Injected electrons have a density of 0 . n and an initial beam velocity of 1 . v A .Using this setup, a secondary shock and a cavity is developed as shown in Figure 6(b), and the structure is qualitativelyconsistent with that from the recent 3D global hybrid simulations (Wang et al. 2020). Below we examine the magneticfield structure and the associated current system of the secondary shock and the cavity, and discuss similarities anddifferences between the 2D and 1D results. ω pi ]020406080100 y [ c / ω p i ] Rotationaldiscontinuity B z(a) 0 20 40 60 80x [c/ ω pi ] n e / n (b) Figure 6.
Formation of foreshock transients in the 2D PIC simulation. (a) The computational setup. (b) The total electrondensity. The secondary shock, the sheath and the cavity are clearly visible. This and other snapshots hereafter from 2Dsimulations are taken at the end of the simulation t = 35 . ω − ci . Because the discontinuity spatial scale is far smaller than the gyroradius of injected ions and the timescale of transientformation is comparable to ion gyroperiod, injected ions are demagnetized [Figs. 7(c) and 8(c)], while injected electronsare nearly always magnetized ( ρ e ≈ . c/ω pi ) and move along the magnetic field lines across the discontinuity[Fig. 8(d)]. As injected ions cross the rotational discontinuity, they cannot change their velocities immediately so thatthe injected ion velocity remain in the x - y plane near the discontinuity ( x = 6 – 10 c/ω pi , y ≈ c/ω pi in Fig. 8(c)).This velocity projects to both perpendicular and parallel directions, resulting in a pitch angle about 60 ◦ with respectto B . As injected ions move further away from the discontinuity, they gradually gyrate to the z direction ( x = 15– 25 c/ω pi , y ≈ c/ω pi in Fig. 7(c)). The motions of demagnetized injected ions and magnetized electrons lead to alarge-scale Hall current in the z direction [Fig. 7(b)] and x - y plane [Fig. 8(b)] (also see the contribution from ambientplasma in Appendix E), which should be maintained by the continuous particle injection. Therefore, a vortex-likevector field of ( B x , B y ) is generated by the z component of the Hall current (the center of the vortex is located where J z of injected ions is maximized in Figures 7(b) and 7(c)). This gives a region of magnetic cavity on one side [Fig. 7(a)]and a region of enhanced magnetic field on the other side, which steepens into a secondary shock (see another currentat its surface carried by ambient plasma at x = 0 – 60 c/ω pi , y = 0 – 40 c/ω pi in Figs. 7(b) and 8(b)). The x and y components of the Hall current [Fig. 8(b)] generate a bipolar field of B z [Fig. 8(a)], which is located inside the cavityand can reach ± . B . There is no such a bipolar B z in 1D simulations, because the spatial variation of electriccurrent in the x - y plane was not resolved. Consequently, the profile of magnetic field strength is not exactly the sameas the density profile in 2D, which is more consistent with observations. oreshock transients and secondary shocks ω pi ]020406080100 y [ c / ω p i ] sqrt(B x2 +B y2 )/B (a) 0 20 40 60 80x [c/ ω pi ] −0.4 −0.2 0.0 0.2 0.4 J z [x10 −2 ] (total)(b) 0 20 40 60 80x [c/ ω pi ] −0.2 −0.1 0.0 0.1 0.2 J z [x10 −2 ] (injected ions)(c) Figure 7.
In-plane magnetic field and associated currents. (a) The vector field ( B x , B y ). The color scale represents thenormalized strength of the in-plane magnetic field (cid:112) B x + B y /B . (b) The total electric current in the z direction. (c) Thecurrent J z contributed by the injected ions. The electromagnetic fields caused by the perpendicular (gyrating) and parallel (streaming) motions of injectedions have distinctively different characteristics. On the one hand, driven by the gyrating motion of injected ions, acompressional boundary with enhanced plasma density and magnetic field is formed [Figs. 6(b) and 7(a)] and steepensinto a shock (the Mach number M f = 1 . ω − k (cid:107) v (cid:107) = − ω cb (e.g., Wilson 2016;Weidl et al. 2019), where v (cid:107) ≈ v bi sin 30 ◦ = 2 . v A is the parallel velocity of injected beam ions and ω cb is the cyclotronfrequency of beam ions. Such wave signatures have also been seen in the recent 3D global hybrid simulations by Wanget al. (2020). Because part of the free energy is released in the parallel expansion, the expansion speed of foreshocktransients would decrease in comparison with that of 1D simulations. CONCLUSIONS AND DISCUSSIONWe have shown the detailed process of diamagnetic cavity and secondary shock formation at foreshock transients. Ourresults provide clear evidence of the critical role electrostatic and induction electric fields play in this formation process,and reveal the energy transfer between different particle species and electromagnetic fields in the formation. Our studyalso demonstrates how a rotational discontinuity interacts with foreshock ions and leads to the formation of a foreshockbubble. Such a process is not simply increasing the thermal pressure of foreshock ions as explained in the previousmodels, but demagnetizing foreshock ions and generating the Hall current. Our ensemble of simulations indicates thatthe expansion speed of foreshock transients is proportional to the hot ion density ratio and thermal speed, suggestingthat foreshock transients with secondary shocks are more prevalent at high Mach number astrophysical shocks thanthose already observed at planetary bow shocks. Since these ion-kinetic structures have been shown to accelerateparticles cooperatively with primary shocks with high efficiency (Wilson III et al. 2016; Liu et al. 2017b, 2018; Turneret al. 2018; Liu et al. 2019), they could significantly contribute to high Mach number astrophysical shock acceleration,e.g., the generation of the cosmic rays. Therefore, they must be included in shock acceleration models in general.In our 2D simulation, we only consider a rotational discontinuity with a certain magnetic field configuration. Thebasic physical process, however, is general. Different types of discontinuity and magnetic field configurations affect thedetails of, e.g., how foreshock ions are demagnetized and how the corresponding Hall current changes the backgroundmagnetic field. In the future, these details will be examined by more advanced simulations. By collaborating with the0
An et al. y [ c / ω p i ] −0.4 −0.2 0.0 0.2 0.4 B z / B (a) sqrt(J x2 +J y2 ) [x10 −2 ] (total)(b)0 20 40 60 80x [c/ ω pi ]020406080100 y [ c / ω p i ] sqrt(J x2 +J y2 ) [x10 −2 ] (injected ions)(c) 0 20 40 60 80x [c/ ω pi ] sqrt(J x2 +J y2 ) [x10 −2 ] (injected electrons)(d) Figure 8.
Out-of-plane magnetic field and associated currents. (a) The color map of B z . (b) The vector field ( J x , J y ). Thecolor scale represents the normalized strength of the in-plane current (cid:112) J x + J y . (c), (d) The in-plane currents contributed byinjected ions and electrons, respectively. oreshock transients and secondary shocks osiris and for providingaccess to the osiris A. COMPUTATIONAL SETUPBelow we present the details of the computational configuration, how the parameters are scaled from the spaceenvironment to the numerical experiments, and the effects of including more spatial dimensions.A.1.
Configuration
We used the massively parallel, fully electromagnetic PIC code osiris (Fonseca et al. 2002; Hemker 2015) for oursimulations. The simulations have one dimension ( x ) in configuration space and three dimensions ( v x , v y , v z ) in velocityspace. The computational domain is − L x (cid:54) x (cid:54) L x [Fig. 9]. The specific size of the system L x is chosen based on thecondition of the hot ions as described below. The cell length ∆ x is 2 λ D , where λ D = v T e /ω pe is the initial electronDebye Length, v T e is the electron thermal velocity, and ω pe is the electron plasma frequency. Each cell contains500 particles per species. The time step is set as 0 . x /c to satisfy the CourantFriedrichsLewy condition in onedimension, where c is the speed of light. The ambient magnetic field B is along the + z direction. The electroncyclotron frequency ω ce is equal to ω pe /
30. Given the reduced ion-to-electron mass ratio m i /m e = 100, the Alfv´envelocity is v A = c/ − ρ h (cid:54) x (cid:54) ρ h [Fig. 9], where ρ h is the gyroradius corresponding to the initial thermal velocity v T h of hotions. Inside the hot layer, the densities of hot and ambient ions are ηn and (1 − η ) n , respectively, where η denotesthe fraction of hot ions. Outside the hot layer, the only positively charged species is ambient ions with density n .Electrons are initialized with a uniform density n in the entire domain. The initial thermal velocities of ambientelectrons and ions are v T e = 3 v A and v T i = 0 . v A , respectively. In the nominal run presented in the main manuscript,the hot ion concentration ratio is η = 0 .
2, and the initial thermal velocity of hot ions is v T h = 9 v A . To explore theconditions under which secondary shocks are formed, we complete 25 runs (including the nominal run), varying thehot ion thermal velocity v T h and the hot ion concentration ratio η . In each run, the hot ion thermal velocity is variedas v T h = (12 v A , v A , v A , v A , v A ) for a given hot ion concentration ratio in the sequence η = (0 . , . , . , . , . ρ h . For each hot ion thermal velocity in the sequence v T h = (12 v A , v A , v A , v A , v A ),2 An et al.
Figure 9.
A sketch of the simulation box. See the text for details of the simulation setup. the corresponding domain size is chosen as L x = (36 d i , d i , d i , d i , d i ), such that the full evolution of thesecondary shocks can be resolved within the computation domain. Here d i denotes the ion inertial length.The domain’s upper half 0 (cid:54) x (cid:54) L x (simulation box) and lower half − L x (cid:54) x (cid:54) x = 0. Although the additional image box doubles the computational cost, it allows us to use the periodicboundary condition for both fields and particles. When presenting the results, we focus on the simulation box.A.2. Scale-down numerical experiments
Table 1 shows how the parameters are scaled from Earth’s foreshock to the nominal simulation. The absolute valuesof the thermal velocities of different species and the alfven velocity in the numerical experiments are larger thanthose in the space environment, for the sake of saving computing time. But the key dimensionless parameters are thescale separations between the gyro-radii of electrons, ambient ions and hot ions, because they determine the chargeseparation and hence the electrostatic fields. We keep the ratio between the gyro-radii of different species on the sameorder as the measured values to capture the key ingredient for foreshock transient formation. [km/s] v A v Te v Ti v Th Earth’s foreshock 50 2000 50 400 − v A /c m i /m e ρ i /ρ e ρ h /ρ e Earth’s foreshock 1 / − /
300 100 10 300
Table 1.
The comparison of relevant dimensional and non-dimensional parameters between Earth’s foreshock and the nominalsimulation.
A.3.
Role of Alfv´en speed
To make the computational cost affordable, we have chosen a larger Alfv´en speed v A than the realistic value inthe simulations. Thus it is necessary to justify how v A affects the formation process of foreshock transients. Theelectrostatic field is caused by a small charge separation. This small charge separation is induced by the differencebetween the gyro-radii of hot ions and other species. To estimate the magnitude of the electrostatic field, we invokethe Gauss’s Law Eρ h ∼ π(cid:15)en h , where ρ h can be further replaced with v Th /ω ci , and ω ci is the ion cyclotron velocity. (cid:15) is a small dimensionlessparameter standing for the “small” charge separation ( (cid:15) (cid:28) n h = α · n ,v Th = β · v A . oreshock transients and secondary shocks E ∼ (cid:15)αβ · πen v A ω ci = (cid:15)αβ · ω pi ω ci v A c B = (cid:15)αβ · cv A B. Here ω ci = Be/m i c and ω pi = 4 πn e /m i . We have also made use of v A /c = ω ci /ω pi . It is more appropriate tonormalize the electrostatic field as (cid:18) cEB (cid:19) /v A ∼ (cid:15)αβ (cid:18) cv A (cid:19) , which is also the ratio of electric drift velocity to Alfv´en velocity. As v A decreases, the electrostatic field in fact increaseswith the scaling ( c/v A ) . Intuitively, as the strength of magnetic field decreases (lower v A ), the difference between thegyro-radii of hot ions and electrons becomes larger, inducing stronger charge separation and larger electrostatic field.Indeed, in Liu et al. (2017a), it is seen that lower interplanetrary magnetic (IMF) field has higher probability to formforeshock transients. Therefore lower v A favors the formation of foreshock transients. In this sense, the secondaryshocks in reality may form more easily than in our simulations. B. DIAMAGNETIC CURRENTTo show the contribution of each species to the total current that shapes the magnetic cavity and compression inforeshock transients, we plot the current J y by each species and the total current in Fig. 10. In the first and secondgyrations of the hot ions (0 < t/τ ci (cid:46) − y direction consistent with the sense of ion gyrationis seen [Fig. 10(a)]. The ambient ions, on the other hand, form a current in the + y direction [Fig. 10(b)] because moreambient ions are distributed outside the hot layer than inside it. As explained in the main manuscript, electrostaticfields E x cause an electron current in the − y direction to form [Fig. 10(c)]. The electron current and the hot ion(partial diamagnetic) current in the − y direction dominate the ambient ion (partial diamagnetic) current in the + y direction. This results in a net Hall current [Fig. 10(d)] that reduces the magnetic field on one side (the cavity) andenhances it on the other side (the compressional boundary). As the magnetic cavity is developed and continues toexpand ( τ ci (cid:38) − y directionat the leading edge of the cavity as some hot ions gyrate outward and is along the + y direction inside the cavity as someother hot ions gyrate inward. The gyroradii of ambient ions and electrons are much smaller than the characteristicspatial scale of the electrostatic fields. This gives rise to similar electric current profiles for ambient ions and electrons[Fig. 10(b) and 10(c)] caused by the drift − cE x /B z . The current associated with the magnetic perturbations of themagnetosonic waves is also seen [Fig. 10(b), 10(c) and 10(d)]. C. PARAMETER SCANIn Fig. 11, we show the detailed evolution of the magnetic field for 25 runs in the parameter scan. Correspondingly,Fig. 12 shows the spatial profile of the magnetic field at the end of each run. As the hot ion thermal velocity ( v T h ) andthe hot ion concentration ratio ( η ) increase, depletion of the magnetic flux in the foreshock cavity, as well as pile upof the magnetic flux in the compressional boundary, becomes more enhanced. At the same time, the Mach number ofsecondary shocks increases with v T h and η . For a higher Mach number of the primary shock, higher thermal velocityand higher concentration ratio of foreshock ions/hot ions will be produced (Burgess et al. 2012; Paschmann & Sckopke1983). Thus, the parameter scan strongly indicates that foreshock transients and associated secondary shocks are morelikely formed at the high Mach number shocks.The saturation time t s of the foreshock transients, i.e., the time when energy transfer between particles and fieldsvanishes, can be estimated as ρ eff /v shock . Here v shock = M f v f is the propagation speed of the secondary shock, and ρ eff is the effective gyroradius of hot ions. The parameter ρ eff characterizes the typical length over which the compressionalboundary is detached from the hot ions and takes the value Cρ h , where C is a constant that depends primarily on thehot ion distribution. For the Maxwellian distribution used in this study, C takes the approximate value 5 based on allthe runs we perform. Therefore, we have t s = Cω ci · v T h M f ( v T h , η ) · v f , where we have written ρ h as v T h /ω ci . The fast magnetosonic Mach number M f is a function of both v T h and η . For agiven hot ion thermal velocity, the saturation time ( t s ∝ /M f ) decreases as the hot ion concentration ratio increases.4 An et al. x / ρ h −0.1 0.0 0.1 J y,h [x10 −2 ](a) −0.2 0.0 0.2 J y,i [x10 −2 ](b)0 2 4 6t / τ ci x / ρ h −0.4 0.0 0.4 J y,e [x10 −2 ](c) 0 2 4 6t / τ ci −0.1 0.0 0.1 J y,total [x10 −2 ](d) Figure 10.
The spatio-temporal evolution of total current contributed by different species. (a) Electric current in the y direction contributed by hot ions. (b) Electric current in the y direction contributed by ambient ions. (c) Electric current in the y direction contributed by electrons. (d) Total electric current in the y direction contributed by all three species. To emphasizethe current patterns, the color scales differ in each panel. For a given hot ion concentration ratio, the saturation time ( t s ∝ v T h /M f ) slowly increases as the hot ion thermalvelocity increases, indicating the dependence M f ∝ v γT h ( γ < v T h = v A , hot ions generate a small magnetic perturbationeach time when they gyrate out of the hot layer. These magnetic perturbations propagate at the fast magnetosonicspeed v f . In these cases, hot ions do not have enough time to transfer their free energy to the magnetic field in eachinteraction, the remaining free energy leads to the consecutive generation of propagating magnetic perturbations. D. CORRELATION BETWEEN THE OCCURRENCE RATE OF FORESHOCK TRANSIENTS AND PRIMARYSHOCK’S MACH NUMBERThe statistical study by Liu et al. (2017a) shows that high solar wind speed and low interplanetary magnetic fieldstrength favor the formation of foreshock transients at Earth’s bow shock, implying that a large solar wind Alfv´enMach number is a favorable condition. To confirm this, we use the same database as in Liu et al. (2017a) and plotthe probability distribution of the solar wind Alfv´en Mach number during the entire time interval of the database(black) and during each observation of a foreshock transient (red), as shown in Fig. 13. We see that the latter oneindeed shows higher probability at large Alfv´en Mach number than the former one. This result indicates that a largeprimary shock Mach number favors the formation of foreshock transients. Further statistical study is needed to testthe relationship between the primary shock Mach number and secondary shock formation. E. ELECTRIC CURRENTS BY AMBIENT PARTICLES AND ELECTRIC FIELDS IN 2D PIC SIMULATIONSThe electric field E z is shown in Figure 14(c), which is an induction electric field coupled with the time-varying( B x , B y ) of the compressional boundary. This induction electric field has also been observed in 1D simulations inSection 2. Ambient elecctrons and ions are advected by this electric field to co-propagate with the shock, as shown inFigures 14(a) and 14(b).The in-plane electric field ( E x , E y ) is shown Figure 15(c). It has the most significant magnitude inside the cavitywhere strong cross-field current of injected ions is located [Figs. 7(c) and 8(c)]. Appreciable electric field magnitude oreshock transients and secondary shocks x / ρ h x / ρ h x / ρ h x / ρ h τ ci x / ρ h τ ci τ ci τ ci τ ci Figure 11.
The spatio-temporal evolution of the normalized magnetic field B z /B for the parameter scan. From top to bottom,each row corresponds to the initial thermal velocity of hot ions in the sequence v Th = (12 v A , v A , v A , v A , v A ). From left toright, each column corresponds to the hot ion concentration ratio in the sequence η = (0 . , . , . , . , . ρ h is different for each row. The total simulation time is different for each panel, which is chosen toresolve the full evolution of the system. (cid:113) E x + E y is also distributed in other regions, such as the shock surface and the region of fast magnetosonic waves.Due to the E × B drift caused by the in-plane electric field, ambient electrons and ions have currents J z similar to1D simulations, as shown in Figures 15(a), 15(b) and 15(d). However, because the in-plane electric field has bothelectrostatic and electromagnetic components, and they cannot be easily separated using the present simulation setup,the contribution of electron current due to the electrostatic field is difficult to diagnose. On the shock surface, on theother hand, it is clear to see the static electric field, which gives the surface current on the secondary shock as shownin Figures 15(d) and 14(d).6 An et al. B z / B B z / B B z / B B z / B ρ h B z / B ρ h ρ h ρ h ρ h Figure 12.
The spatial profiles of the magnetic field B z /B for the parameter scan. The position of each panel is arranged inthe same way as in Fig. 11. Each profile is taken at the end of the simulation. REFERENCES
Angelopoulos, V., Cruce, P., Drozdov, A., et al. 2019, SpaceScience Reviews, 215, 9Archer, M., Turner, D., Eastwood, J., Schwartz, S., &Horbury, T. 2015, Planetary and Space Science, 106, 56Biermann, L., Brosowski, B., & Schmidt, H. 1967, SolarPhysics, 1, 254Bondarenko, A., Schaeffer, D., Everson, E., et al. 2017,Nature Physics, 13, 573 Bonde, J., Vincena, S., & Gekelman, W. 2015, PhysicalReview E, 92, 051102Burgess, D. 1989, Journal of Geophysical Research: SpacePhysics, 94, 472Burgess, D., M¨obius, E., & Scholer, M. 2012, Space ScienceReviews, 173, 5Burgess, D., Wilkinson, W., & Schwartz, S. 1989, Journalof Geophysical Research: Space Physics, 94, 8783 oreshock transients and secondary shocks Figure 13.