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 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 a 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. The 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 non-relativistic CME-driven shocks associated with enhanced electron intensities from $<40keV$ to about $200keV$ (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 of 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 $\beta$ (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 2-D shock simulation (ambient $B_{0}$ in-plane case) with a relatively small box along the shock front ($\sim 1c/\omega_{pi}$, where $c/\omega_{pi}$ is the ion inertial length). They identify whistler waves excited by modified two-stream instability (MTSI) (Matsukiyo & Scholer, 2003) by using a dispersion relation analysis of fluctuating electromagnetic field components $\delta B$ and $\delta E$. Based on this method, we directly extract high resolution fluctuating electromagnetic field components from 2-D large scale shock simulations (ambient $B_{0}$ out-of-plane case, including the complete particle gyro-motion perpendicular to the background $B_{0}$), 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 paper 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 Parker Solar Probe 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 which 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 paper, collisionless shocks were generated by the so-called injection method (Matsukiyo & Scholer, 2012; Yang et al., 2015), in which particles are injected from the one side of the simulation boundary (at $X=0$) at super-Alfvénic speed $V_{inj}=6V_{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 cell is $n_{x}\times n_{y}=25,600\times 2,000$. The spatial resolution $\Delta X=\Delta y=0.0025c/\omega_{pi}$. The box size along the shock surface is $5c/\omega_{pi}$. The ion-to-electron mass ratio $m_{i}/m_{e}$ is 100 and the particle number per cell is 50. In the 2-D simulation, the ambient magnetic field $B_{0}$ is along $Z$ and strictly perpendicular to the $X-Y$ simulation plane (i.e., the shock normal angle $\theta_{Bn}$) is nearly $90^{\circ}$. This $B_{0}$ 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}$. Such as the magnetic field $B\sim 591nT$, the proton density $N_{p}\sim 1932cm^{-3}$, the proton temperature $T_{p}\sim 45.3eV$, and the Alfvénic velocity $V_{A}\sim 290km/s$. $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\ Rs$ can often exceed $1500\sim 2200km/s$ (Zhao et al., 2019), and a fraction of them have extremely high speeds ($\sim 3300km/s$) (Liu et al., 2013, 2019). Their corresponding average $\overline{M}_{A}$ is $5\sim 9$ and the extreme ones can reach $10\sim 15$. Such Mach numbers are much larger than that 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\sim 9$. Plasma beta values: $\beta_{e}=0.3$ and $\beta_{i}=0.1$ are employed for the real solar wind condition. We adopted the upstream magnetization $\omega_{pe}/\Omega_{ce}\sim 7.8$, on the order of 10 for the near-Sun solar wind conditions. In addition, we examined wave properties at 2-D shocks with similar setups for different Mach numbers: $M_{A}=5\sim 6$ and $>10$ for slower and faster shocks, respectively. A 3-D shock simulation is also carried out for comparing the ES wave property in higher dimensions. In this paper, we focus on the moderate-Mach-number 2-D shock with universal significance. The other cases will be discussed in section 4.

Figure 1a 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\Omega_{ci}^{-1}$ in this case. In order to study wave properties at the shock front with a relatively high resolution in $\kappa-\omega$ space, the electromagnetic field components are sampled in a long period: $t=3.5\sim 4.5\Omega_{ci}^{-1}$, and a large box including the whole shock front. Figure 1b-d represent the snapshots of the fluctuations of $\delta E_{x}$, $\delta E_{y}$ and $\delta B_{z}$ at a typical time $t=3.5\Omega_{ci}^{-1}$ within the sampling period. The zoomed $\overline{B}_{z}$ profile consisting typical shock structures: upstream, foot, ramp, overshoot, and downstream is shown in Figure 1d for reference. The shock ramp is located at about $X=56.4$ 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 is sampled for the wave analysis. Figure 1e-f illustrate a slice of the sampled shock profile $\delta E_{x}$ at $Y=L_{y}/2$ in regions A and B, respectively. In order to investigate the dispersion relation of waves, the 3-D fluctuation data in the $X-Y-t$ space are converted to the shock ramp rest frame. A 3-D Fourier transform is applied after Hanning windowing to compensate for the nonperiodicity of the data, in both $X$ and $t$ directions (Figure 1g-h). Similar process are also carried out for other fluctuation components $\delta E_{y}$ and $\delta B_{z}$.

Figure 2a-d shows corresponding phase space plots ($X-V_{x,y}$) of particles at $t=3.5\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 Figure 1b-d. In region A (Figure 2a-b), some electrons are trapped and accelerated by the excited ES waves (Figure 1b-c) at the leading edge of the foot (Umeda et al., 2009; Amano & Hoshino, 2009; Yang et al., 2018). They immediately form a ring-like velocity distribution relative to the hot drifting Maxwellian core. Figure 2c-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 $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 in consistent with previous simulations (Scholer et al., 2003; Yang et al., 2009). Figure 2e 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. Figure 2f-h represent the 2-D density profiles of electrons, reflected ions and incident ions. Combined with corresponding field fluctuations (Figure 1b-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.

Figure 3 shows the 3-D view of the wave dispersion relation diagram of $\delta E_{x}$, $\delta E_{y}$, and $\delta B_{z}$ in the $\kappa_{x}-\kappa_{y}-\omega$ space. Top and bottom panels denote results from regions A and B, respectively. In the simulation, the electron plasma frequency $\omega_{pe}$ and the Debye length $\lambda_{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 is averaged over the whole sampling period in the shock rest frame. More details are shown in Table 1.

Firstly, 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 Figure 3a-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 $+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 braches on the $+\kappa_{y}$ and $-\kappa_{x}$ side (Figure 4a-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 gyro-motion 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 Figure 3a-b. After the “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_{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 Figure 3a-b). (3) the third wave is an electromagnetic mode (marked by “EM-1” in Figure 3c). This high frequency EM wave is emitted facing the upstream and is visible in both $\delta E_{\perp}$ and $\delta B_{z}$ diagrams. In this case, the $\delta 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.

Secondly, the waves in region B (i.e., in the middle of the foot) are investigated. Figure 2e shows that the values of $n_{e}$ and $B$ are $\geq 2$ times of their upstream values $n_{e0}$ and $B_{0}$ in this region. The local ratios $\omega_{pe}/\omega_{pe0}$ and $\Omega_{ce}/\Omega_{ce0}$ are about 1.414 and 2, respectively. The frequency ratio $\omega_{pe}/\Omega_{ce}$ is about 5.5 and lower than its upstream value ($\omega_{pe0}/\Omega_{ce0}\approx 7.8$). This means that the magnetization of the plasma becomes greater. We keep using the upstream values $\omega_{pe0}$ amd $\lambda_{De0}^{-1}$ as units of $\omega$ and $\kappa$ for region B in bottom panels of Figure 3. Corresponding dispersion relations from the linear theory are shown in Figure 4c-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 3d, “ES-3” wave modes are excited by ECDI on $\pm\kappa_{x}$ directions. The wave frequency in the $\pm\kappa_{x}$ directions are comparable to each other as predicted by the linear theory (Figure 4c). (2) in Figure 3e, “ES-4” wave modes have more harmonic branches in the $+\kappa_{y}$ direction. This is because the old reflected ions are gyrating back towards the ramp and they have a large bulk velocity ($>10V_{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 1-D simulations (Muschietti & Lembège, 2006, 2013) which only discuss ECDI along the $x-$axis. Our work has taken a step forward on this basis, and find 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 3f) associated with the “ES-4” waves. Such low frequency EM wave has a modulating effect on the magnetic field and plasma density profiles (refer Figure 1d and Figure 2f-h). The analysis of such EM wave amplification might requires a theory considering nonuniform plasmas and gradient background magnetic fields. Hence, it will not be further doscussed in this paper. It is worthy noting that the X mode can exist well in region B due to the local relativistic electrons, and it has a wave vector towards upstream as in region A.

This paper presents PIC simulations of a perpendicular shock in the near-Sun solar wind condition, where the magnetization is relatively high, $\beta$ 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 electron cyclotron drift instability in the $-X$ and $+Y$ directions. The instability is triggered by the coupling between the incident electrons and gyro-reflected ions. The wave vector 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 this 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 a X mode and propagates upstream. Other potential candidates of X mode emission mechanisms are also need to be considered for quasi-perpendicular shocks. For an instance, a combination of wave growth due to electron cyclotron maser instability and nonlinear wave-coupling processes is suggested for plasmas with $\omega_{pe}/\Omega_{ce}\sim 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 begin 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 $\pm X$ and $\pm Y$ directions. These multi-directional 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 1-D 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 is strengthen as it approaches 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 $>3000\ km/s$ (i.e., $Ma>10$) (Liu et al., 2019). Normally this is not the case and the shock speed is as slow as that used in this paper. At slower shocks (Ma$\leq 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 3-D shock simulation is also carried out for comparing ES wave excitations (not shown here). Briefly, ES waves sampled in the $B_{0}$ out-of-plane from the 3-D simulation is similar to that observed in this paper. However, ES waves sampled in the $B_{0}$ in-plane propagate almost along the shock normal. This is because the main ion bulk velocity components is not included in this cross section. More completely, a 4-D FFT is required for the wave analysis in 3-D 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\sim 6$ years (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). Parker Solar Probe and Solar Orbiter are likely to observe high frequency EM emissions (e.g., X modes) toward upstream accompanied with 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.