PIC Simulations of Microinstabilities and Waves at Near-Sun Solar Wind Perpendicular Shocks: Predictions for Parker Solar Probe and Solar Orbiter

Collisionless shocks are of fundamental interests in astrophysics and space physics. They have been proposed as primary mechanisms for energy dissipation (e.g., Richardson et al. 2008; Parks et al. 2012) and particle acceleration (e.g., Zank et al. 2006; Guo & Giacalone 2010). Observational studies suggest that large solar eruptions are often accompanied with shocks driven by corona mass ejections (CMEs). They are rich in various plasma waves (Wilson et al. 2007; Liu et al. 2018) and usually associated with type II radio bursts (Bale et al. 1999; Liu et al. 2009). Previous investigations reveal that there are two main mechanisms for high-frequency electromagnetic (EM) wave emissions (e.g., O or X modes) at collisionless shocks. One is synchrotron maser instability (SMI), which usually refers to an electromagnetic wave emission mechanism at relativistic magnetized shocks (Hoshino & Arons 1991; Plotnikov & Sironi 2019). This mechanism is in favor of a positive slope in an electron velocity distribution function (VDF) perpendicular to the ambient magnetic field, such as due to a loss cone or ring velocity distribution, which can become a free energy source for wave emissions. Essentially the same mechanism may operate as long as weakly relativistic anisotropic electrons exist (even down to a few keV) in the case of auroral electrons (Wu & Lee 1979). So it may indeed be a potential candidate for radio emissions at nonrelativistic CME-driven shocks associated with enhanced electron intensities from <40 keV to about 200 keV (Liu et al. 2008). The other is the nonlinear three-wave interaction (Pulupa et al. 2010; Gao et al. 2017). First, suprathermal and nonthermal electrons are produced in the shock. Then, these energetic electrons stimulate the growth of high-frequency electrostatic (ES) waves, such as Langmuir waves. These ES waves interact with each other in the nonuniform background solar wind or rippling shock front to produce the observed EM emissions due to the nonlinear wave interactions (Umeda 2010; Ganse et al. 2012). However, it is still unclear which mechanism and what wave mode plays a major role at CME-driven shocks near the Sun.

Electrostatic waves, such as (1) Langmuir waves and their higher harmonics (Bale et al. 1999; Thejappa & MacDowall 2019) and (2) electron cyclotron harmonic waves (Bernstein waves; Wilson et al. 2010; Goodrich et al. 2018), have been widely observed at shocks by spacecraft (ISEE 1, WIND, STEREO, CLUSTER, MMS, etc.) around 1 au. They are commonly believed to accelerate electrons and provide possible free energy sources of EM emissions at shocks, magnetopause, and solar flares (Graham et al. 2018; Horký et al. 2018; Henri et al. 2019). Despite this, there is a lack of in situ observations accurately determining shock microstructures and associated plasma waves excited in the near-Sun solar wind. Bale et al. (2016) extracted some solar wind parameters in the near-Sun conditions based on data from HELIOS and models, which give us some inspiration. The study of wave properties at shocks under the near-Sun solar wind condition could help understand ES wave excitations and EM wave emissions at young CME-driven shocks observed by Parker Solar Probe (PSP) or Solar Orbiter near the perihelion at later encounters. Near the Sun, we expect a strong magnetization and relatively low values of plasma β (Bale et al. 2016) that may affect characteristics of waves excited at the shock front. This is our motivation to do the work.

Umeda et al. (2012a) directly extract the field components of the wave from a 2D shock simulation (ambient B₀ in-plane case) with a relatively small box along the shock front (∼1c/ω_pi, where c/ω_pi is the ion inertial length). They identify whistler waves excited by modified two-stream instability (Matsukiyo & Scholer 2003) by using a dispersion relation analysis of fluctuating electromagnetic field components δB and δE. Based on this method, we directly extract high-resolution fluctuating electromagnetic field components from 2D large-scale shock simulations (ambient B₀ out-of-plane case, including the complete particle gyromotion perpendicular to the background B₀), and check the wave modes from the leading edge of the foot to the ramp. Furthermore, observed linear waves are confirmed by a linear theory tool (BO, a new version of PDRK; Xie 2019).

This Letter is organized as follows. We present a description of the simulation model and setup in Section 2. The shock front wave analysis is presented in Section 3, and we conclude with a summary in Section 4 and discuss the implications of our results for PSP and Solar Orbiter.

We carry out shock simulations using an open-source electromagnetic particle-in-cell (PIC) code named EPOCH (Arber et al. 2015) with a normalization that is the same as that used in our previous works (Yang et al. 2016, 2018; Lembege et al. 2020) to simulate the waves excited at moderate-Mach-number, perpendicular shocks. In this Letter, collisionless shocks were generated by the so-called injection method (Matsukiyo & Scholer 2012; Yang et al. 2015), in which particles are injected from one side of the simulation boundary (at X = 0) at super-Alfvénic speed ${V}_{\mathrm{inj}}=6{V}_{{\rm{A}}}$ in the +X direction and specularly reflected at the other side of the simulation boundary (X = L_x). The shock propagates in the −X direction in the present downstream rest frame. The periodic boundary condition is applied in the Y direction. The number of grid cells is n_x × n_y = 25,600 × 2000. The spatial resolution ΔX = Δy = 0.0025c/ω_pi. The box size along the shock surface is 5c/ω_pi. The ion-to-electron mass ratio m_i/m_e is 100, and the particle number per cell is 50. In the 2D simulation, the ambient magnetic field B₀ is along Z and strictly perpendicular to the X − Y simulation plane (i.e., the shock normal angle θ_Bn) is nearly 90°. This B₀ configuration is similar to that in previous simulations (Amano & Hoshino 2009; Matsumoto et al. 2013). Based on fitting methods similar to Bale et al. (2016), plasma parameters in the near-Sun solar wind (at about 10R_s) can be estimated by using PSP data observed around its perihelion in the first three encounters at about 36R_s. The magnetic field $B\sim 591\,{nT}$, the proton density ${N}_{p}\sim 1932\,{\mathrm{cm}}^{-3}$, the proton temperature T_p ∼ 45.3 eV, and the Alfvénic velocity ${V}_{{\rm{A}}}\sim 290\,\mathrm{km}\,{{\rm{s}}}^{-1}$. T_e/T_i = 3 is adopted based on previous observations from HELIOS (Liu et al. 2005) and PSP (Maksimovic et al. 2020). The speed of fast CME-driven shocks at 10 R_s can often exceed $1500\sim 2200\,\mathrm{km}\,{{\rm{s}}}^{-1}$ (Zhao et al. 2019), and a fraction of them have extremely high speeds ($\sim 3300\,\mathrm{km}\,{{\rm{s}}}^{-1}$; Liu et al. 2013, 2019). Their corresponding average ${\overline{M}}_{{\rm{A}}}$ is 5 ∼ 9 and the extreme ones can reach 10 ∼ 15. Such Mach numbers are much larger than those observed at interplanetary (IP) shocks at 1 au (Wilson et al. 2007; Liu et al. 2018). In this simulation, the Alfvénic Mach number of the shock is about 7 ∼ 9. Plasma beta values β_e = 0.3 and β_i = 0.1 are employed for the real solar wind condition. We adopted the upstream magnetization ω_pe/Ω_ce ∼ 7.8, on the order of 10 for the near-Sun solar wind conditions. In addition, we examined wave properties at 2D shocks with similar setups for different Mach numbers: M_A = 5 ∼ 6 and >10 for slower and faster shocks, respectively. A 3D shock simulation is also carried out for comparing the ES wave property in higher dimensions. In this Letter, we focus on the moderate-Mach-number 2D shock with universal significance. The other cases will be discussed in Section 4.

Figure 1(a) shows an overview of the time-evolving shock magnetic field ${\overline{B}}_{z}$ averaged along Y in the simulation. The shock becomes mature and has reached a fully evolved state after t > 2.5Ω_ci⁻¹ in this case. In order to study wave properties at the shock front with a relatively high resolution in κ − ω space, the electromagnetic field components are sampled in a long period: $t=3.5\sim 4.5{{\rm{\Omega }}}_{{ci}}^{-1}$, and a large box including the whole shock front. Figures 1(b)–(d) represent the snapshots of the fluctuations of δE_x, δE_y, and δB_z at a typical time $t=3.5{{\rm{\Omega }}}_{{ci}}^{-1}$ within the sampling period. The zoomed ${\overline{B}}_{z}$ profile consisting of typical shock structures (upstream, foot, ramp, overshoot, and downstream) is shown in Figure 1(d) for reference. The shock ramp is located at about X = 56.4 c/ω_pi and marked by a vertical black line. From left to right, fluctuation profiles reveal different wave properties. Two regions are selected for wave analysis (marked by red and blue horizontal lines: A and B). In region A, ES waves are dominant and propagate in nearly the same direction. In region B, ES waves have different propagation directions and pass through each other. Long-wavelength EM waves appear and gradually increase from the middle of the foot to the ramp. Let us take the fluctuation $\delta {E}_{x}(x,y,t)$ as an example to show how the data are sampled for the wave analysis. Figures 1(e)–(f) illustrate a slice of the sampled shock profile δE_x at Y = L_y/2 in regions A and B, respectively. In order to investigate the dispersion relation of waves, the 3D fluctuation data in the X − Y − t space are converted to the shock ramp rest frame. A 3D fast Fourier transform (FFT) is applied after Hanning windowing to compensate for the nonperiodicity of the data, in both the X and t directions (Figures 1(g)–(h)). Similar processes are also carried out for the other fluctuation components, δE_y and δB_z.

Zoom In Zoom Out Reset image size

Figures 2(a)–(d) show corresponding phase-space plots ($X-{V}_{x,y}$) of particles at $t=3.5{{\rm{\Omega }}}_{{ci}}^{-1}$. The solar wind is coming from the left-hand side, and the shock ramp is marked by a black vertical line as in Figures 1(b)–(d). In region A (Figures 2(a)–(b)), some electrons are trapped and accelerated by the excited ES waves (Figures 1(b)–(c)) at the leading edge of the foot (Amano & Hoshino 2009; Umeda et al. 2009; Yang et al. 2018). They immediately form a ring-like velocity distribution relative to the hot drifting Maxwellian core. Figures 2(c)–(d) show that a fraction of incident ions are reflected at the ramp, and the others are directly transmitted to the downstream. At the same time, the electron bulk velocity shifts in both the X and Y directions to keep the quasi-neutrality. The ion bulk velocity component V_iy can be larger than V_ix in region A due to their gyromotion and the shock acceleration along the shock surface. This overall picture is consistent with previous simulations (Scholer et al. 2003; Yang et al. 2009). Figure 2(e) shows Y-averaged number density profiles of electrons (green), incident ions (blue), and reflected ions (red). Here, we only focus on the wave analysis region and quickly separate the reflected ions before the ramp as in Otsuka et al. (2019). From the foot to the ramp, the percentage of reflected ions increases. Figures 2(f)–(h) represent the 2D density profiles of electrons, reflected ions, and incident ions. Combined with corresponding field fluctuations (Figures 1(b)–(d)), we find that the coupling between reflected ions and electrons is strong in region A. The coupling between incident ions and electrons along the shock normal becomes noticeable in region B.

Zoom In Zoom Out Reset image size

Figure 3 shows the 3D view of the wave dispersion relation diagrams of δE_x, δE_y, and δB_z in the κ_x − κ_y − ω space. Top and bottom panels denote results from regions A and B, respectively. In the simulation, the electron plasma frequency ω_pe and the Debye length λ_De as units are measured in the far upstream undisturbed solar wind. By using a unified tool for plasma waves and instabilities analysis, BO (Xie 2019), corresponding dispersion relations calculated by the linear theory are shown in Figure 4 to assist in identifying wave modes observed in the simulation. The input plasma parameters for the linear analysis are directly issued from the simulation data, which are averaged over the whole sampling period in the shock rest frame. More details are shown in Table 1.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 1.  Setups of Species Parameters for the Linear Analysis

Run	B/B0	N/N0	Species	n%	${\beta }_{//}$	β⊥	${V}_{{dx}}/{V}_{{\rm{A}}}$	Vdy/VA	Vdz/VA	Vring/VA	M
1 (Region A)	1.0	1.0	Inc. H+	83.5	0.1	0.25	9.00	−0.16	0	0	0
			Ref. H+	16.5	0.1	0.75	−1.60	14.8	0	0	0
			Core e−	90.0	0.3	0.75	7.25	2.30	0	0	1
			Ring e−	10.0	0.3	0.1	7.25	2.30	0	25.0	1
2 (Region B)	2.0	2.0	Inc. H+	55.0	0.1	0.5	7.82	−0.83	0	0	0
			New ref. H+	25.4	0.1	0.5	−3.50	5.0	0	0	0
			Old ref. H+	19.6	0.1	0.5	3.63	12.5	0	0	0
			Core e−	90.0	0.3	2.0	4.12	3.26	0	0	1
			Ring e−	10.0	0.3	0.1	4.12	3.26	0	28.0	1

Note. (1) The species can be treated as either magnetized (M = 1) or unmagnetized (M = 0). (2) J = 8 is used for the J-pole Padé expansion. (3) The electron VDF is contributed by a hot drifting Maxwellian core and a relatively cool Maxwellian ring. The ring velocity distribution and its drift across the field are modeled as that used by Umeda et al. (2012b) and Umeda & Nakamura (2018), respectively. The ion VDF is a superposition of two or three drifting Maxwellian subpopulations. As shown in Figure 2, "Inc." and "Ref." are abbreviations of "incident" and "reflected," respectively. "New ref." represents ions that are newly reflected at the ramp and moving toward upstream. "Old ref." refers to the ions that are reflected at earlier time, and they are gyrating back toward the downstream at this time. (4) ${V}_{{dx},y,z}$ and V_ring indicate the drift velocity components and the radius of modeled ring distributions. (5) We consider the nearly perpendicular wave modes. Wave normal angles (WNAs) θ_kB = 89.9° and 89° are employed for regions A and B, respectively. (6) All parameters are issued from the simulation and averaged over the sampling period in the shock ramp rest frame.

Download table as:  ASCIITypeset image

First, we study the waves in region A (the leading edge of the foot). In region A, three main wave modes are observed: (1) the first is an ES wave marked by "ES-1" in Figures 3(a)–(b). This wave is rapidly excited by the relative drift between incident electrons and gyro-reflected ions. As shown in Figure 2, the main bulk velocity component of reflected ions is in the +Y and −X directions at the beginning of the foot. So the reflected ion beam is strongly coupled with the Doppler-shifted electron cyclotron harmonic branches on the +κ_y and −κ_x sides (Figures 4(a)–(b)). This coupling process is discussed as electron cyclotron drift instability (ECDI) by Muschietti & Lembège (2013). They assumed that the reflected ion beam is straight along the shock normal and the gyromotion is not considered (i.e., only focus on the X direction). From the simulation above, we realize that the ion bulk velocity in the Y direction can be large, and even plays a more important role. This is a new point. (2) The second is also an ES wave marked by "ES-2" in Figures 3(a)–(b). After "ES-1" is excited, a fraction of incident electrons are trapped and accelerated by the "ES-1" wave. One possible acceleration mechanism of these electrons is the shock surfing acceleration where the electric field or potential plays an important role (Zank et al. 1996; Hoshino & Shimada 2002; Amano & Hoshino 2009; Matsumoto et al. 2017). The accelerated electrons immediately form a ring-like velocity distribution. This VDF has a weak positive slope in the velocity space perpendicular to the magnetic field. Free energy is quickly released through the coupling between the ring and the hot Maxwellian core. The growth rate of "ES-2" peaks at about $\kappa =\pm 0.23{\lambda }_{\mathrm{De}}^{-1}$ in the linear theory around the Doppler-shifted upper hybrid frequency (Figure 4). It is consistent with the simulation results (purple regions in Figures 3(a)–(b)). (3) The third wave is an electromagnetic mode (marked by "EM-1" in Figure 3(c)). This high-frequency EM wave is emitted facing the upstream and is visible in both the δE_⊥ and δB_z diagrams. In this case, the δE_z only has a background noise (not shown here). Hence, the EM mode is an extraordinary electromagnetic mode (X mode). This X-mode emission cannot be described in Figure 4 because kinetic relativistic effects are not included in the present version of our linear solver yet. One possible scenario is the synchrotron maser radiation excited by relativistic electrons with a ring-type velocity distribution at strongly magnetized plasmas (Hoshino & Arons 1991; Plotnikov & Sironi 2019) under a Doppler-shifted condition. In summary, this EM wave is locally emitted at the shock front and propagates toward the upstream.

Second, the waves in region B (i.e., in the middle of the foot) are investigated. Figure 2(e) shows that the values of n_e and B are ≥2 times their upstream values n_e0 and B₀ in this region. The local ratios ${\omega }_{{pe}}/{\omega }_{{pe}0}$ and Ω_ce/Ω_ce0 are about 1.414 and 2, respectively. The frequency ratio ω_pe/Ω_ce is about 5.5 and lower than its upstream value (ω_pe0/Ω_ce0 ≈ 7.8). This means that the magnetization of the plasma becomes greater. We keep using the upstream values ${\omega }_{{pe}0}$ and ${\lambda }_{{De}0}^{-1}$ as units of ω and κ for region B in the bottom panels of Figure 3. Corresponding dispersion relations from the linear theory are shown in Figures 4(c)–(d). In region B, both reflected and incident ion beams couple with the electron cyclotron harmonic branches in all directions of the simulation plane. The main excited waves are as follows: (1) in Figure 3(d), "ES-3" wave modes are excited by ECDI on ±κ_x directions. The wave frequency in the ±κ_x directions are comparable to each other as predicted by the linear theory (Figure 4(c)). (2) In Figure 3(e), "ES-4" wave modes have more harmonic branches in the +κ_y direction. This is because the old reflected ions are gyrating back toward the ramp and they have a large bulk velocity ($\gt 10{V}_{{\rm{A}}}$) in the +Y direction. In contrast, the bulk velocity of incident ions in the  −Y direction caused by the deflection ahead of the ramp is relatively low. In summary, the middle of the foot could be a zoom of ECDI, which is triggered in different directions. This is interesting and brand new relative to previous 1D simulations (Muschietti & Lembège 2006, 2013) that only discuss ECDI along the X-axis. Our work has taken a step forward on this basis, and finds that the ES waves can be excited along the shock surface as well as along the shock normal. In addition, we also find some low-frequency EM waves (marked by "EM-2" in Figure 3(f)) associated with the "ES-4" waves. Such a low-frequency EM wave has a modulating effect on the magnetic field and plasma density profiles (refer to Figures 1(d) and 2(f)–(h)). The analysis of such EM wave amplification might require a theory considering nonuniform plasmas and gradient background magnetic fields. Hence, it will not be further discussed in this Letter. It is worth noting that the X mode can exist well in region B due to the local relativistic electrons, and it has a wavevector toward the upstream as in region A.

This Letter presents PIC simulations of a perpendicular shock in the near-Sun solar wind condition, where the magnetization is relatively high, β is relatively low, and the average shock Mach number is below 10, corresponding to a young, fast CME-driven shock propagating in the pristine solar wind. The simulated electromagnetic fluctuations show that the shock foot can be segmented into two parts by different wave features.

1. At the leading edge of the foot, ES waves are excited by ECDI in the −X and +Y directions. The instability is triggered by the coupling between the incident electrons and gyro-reflected ions. The wavevector of this ES wave is oblique to the shock normal, and is mainly along the shock surface. This is because the bulk velocity of reflected ions changes from the −X direction to the +Y direction during the gyro-reflection process.

2. In the same region, a fraction of incident electrons can be trapped and accelerated by the above exited ES waves. A secondary instability occurs around the Doppler-shifted upper hybrid frequency due to the coupling between these accelerated ring-like electrons and the hot Maxwellian core electrons. In addition, some weak relativistic electrons in the ring VDF could lead to a high-frequency EM emission probably induced by SMI. The emitted EM wave is an X mode and propagates upstream. Other potential candidates of X-mode emission mechanisms are also need to be considered for quasi-perpendicular shocks. For instance, a combination of wave growth due to electron cyclotron maser instability and nonlinear wave-coupling processes is suggested for plasmas with ω_pe/Ω_ce ∼ 10 in the outer corona (Ni et al. 2020).

3. From the middle of the foot all the way to the ramp, the incident ions begin to couple with the electrons. In this region, the incident ion beam begins to deflect in the  −Y direction when approaching the ramp, and the old gyro-back-reflected ions have a large bulk velocity in the +Y direction. Therefore, ES waves can be excited by ECDI in both the ±X and ±Y directions. These multidirectional ES waves propagate across each other and form fibrous-like pattern in electric field and particle number density profiles. This is brand new and different from early 1D simulations and theoretical models of shock front ES wave excitations.

4. The long-wavelength low-frequency EM waves associated with ES harmonics in the +Y direction are strengthened as they approach the ramp. In addition, the high-frequency X mode is also observed in this region.

Furthermore, we carried out two additional shock simulations with lower and higher Mach numbers as mentioned in Section 2. Preliminary results indicate that the ECDI is robust and can be observed in all cases. Buneman instability appears only at extremely fast CME-driven shocks with a speed $\gt 3000\ \,\mathrm{km}\,{{\rm{s}}}^{-1}$ (i.e., M_A > 10; Liu et al. 2019). Normally this is not the case and the shock speed is as slow as that used in this Letter. At slower shocks (M_A ≤ 6), the reflected ion beam is modulated by the shock front self-reformation (Amano & Hoshino 2009) and results in a intermittent ES wave excitation.

A 3D shock simulation is also carried out for comparing ES wave excitations (not shown here). Briefly, ES waves sampled in the B₀ out-of-plane from the 3D simulation is similar to that observed in this Letter. However, ES waves sampled in the B₀ in-plane propagate almost along the shock normal. This is because the main ion bulk velocity components are not included in this cross section. More completely, a 4D FFT is required for the wave analysis in 3D shock cases. This is still infusible under most of the current storage conditions.

As we all know, the Sun is entering the solar maximum in the next 5 ∼ 6 yr (solar cycle 25). More solar eruptions accompanied with fast CME-driven shocks are expected. The perihelion of PSP will be closer to the Sun (from $27.8\ {R}_{s}$ to less than $10\ {R}_{s}$; Bale et al. 2016; Fox et al. 2016; Kasper et al. 2016). PSP and Solar Orbiter are likely to observe high-frequency EM emissions (e.g., X modes) toward upstream accompanied by local energetic electrons and electrostatic waves induced by ECDI or BI along the shock normal or along the shock surface at fast CME-driven shocks.

We thank the referee for helpful comments. We acknowledge the NASA Parker Solar probe Mission, and FIELDS and SWEAP instruments for use of data. The authors are grateful to Dr. L. Muschietti from UC Berkeley for helpful discussions on PDRK/BO benchmarks. The computations are performed by Numerical Forecast Modeling R&D and VR System of State Key Laboratory of Space Weather, and HPC of Chinese Meridian Project. This work is supported by the NSFC (41574140, 41674168, 41774179), the Specialized Research Fund for State Key Laboratories of China, Youth Innovation Promotion Association of the CAS (2017188), the Open Research Program Key laboratory of Geospace Environment CAS (GE2017-01), the Open Research Program Key laboratory of Polar Science, MNR (KP202005), Beijing NSF (1192018), Beijing Municipal Science and Technology Commission (Z191100004319001, Z191100004319003), Beijing Outstanding Talent Training Foundation (2017000097607G049), and the Strategic Priority Research Program of CAS (XDA14040404). S.M. acknowledges partial support by Grant-in-Aid for Scientific Research (C) No.19K03953 and (B) No.17H02966 from JSPS. J.G. is supported by the Shandong Provincial National Natural Science Foundation (ZR2017MD012).