GLOBAL EXPLICIT PARTICLE-IN-CELL SIMULATIONS OF THE NONSTATIONARY BOW SHOCK AND MAGNETOSPHERE

Collisionless shocks are of particular interest for space, plasma, and astrophysics. In the shock transition the bulk energy of the plasma is converted into thermal energy in the absence of particle collisions (Tidmann & Krall 1971; Lembège et al. 2004; Burgess et al. 2005). The collisionless shock has been studied for many decades. However, the cyclic reformation of the shock structure is still a major unresolved issue for collisionless shock physics (Winske et al. 2003; Hada et al. 2003; Scholer et al. 2003; Lee et al. 2004; Burgess & Scholer 2007; Umeda et al. 2008; Yang et al. 2009; Rekaa et al. 2014). The term "self-reformation" describes a process where the particles reflected by the shock ramp accumulate ahead of the shock and form a shock foot, which then grows and becomes a new ramp. The new ramp starts to reflect incident particles, and the process repeats itself. One striking point is that the ramp width can be very narrow, covering a few electron inertial lengths during the reformation cycle (Hada et al. 2003; Scholer et al. 2003; Mazelle et al. 2010). At such narrow shock ramps, the cross-shock electric field is large and the incident particles can be efficiently accelerated to very high energies (Lee et al. 1996; Zank et al. 1996; Yang et al. 2009). The shock front self-reformation was initially predicted from one-dimensional hybrid (Quest 1985) and particle-in-cell (PIC; Lembège & Dawson 1987) simulations. This problem has also been investigated by theoretical studies (Krasnoselskikh et al. 2002). In situ measurements of the terrestrial bow shock made by the Cluster mission (Moullard et al. 2006; Lobzin et al. 2007; Mazelle et al. 2010) clearly show that the shock front is strongly nonstationary. The shock front nonstationarity is also important for the bow shocks at other planets (e.g., Uranus and Mercury) and the heliospheric termination shock (Burlaga et al. 2008; Richardson et al. 2008; Tiu et al. 2011; Sundberg et al. 2013). Thus, the shock front nonstationarity is a widespread natural physical phenomenon. However, the effects of curved geometry and the solar wind condition on the bow shock front nonstationarity are still not well understood.

Although global magnetohydrodynamics (MHD) simulations can successfully describe the macrostructures (Ogino 1986; Wang et al. 2014) and macroinstabilities (Farrugia et al. 1998; Li et al. 2013) of the magnetosphere, the hybrid and PIC simulations are needed to solve the problems of microstructures on ion and electron scales and microinstabilities. During the past two decades, global hybrid simulations have been performed to study bow shocks of different planets (Swift 1995; Lin 2003; Omidi et al. 2005, 2006; Trávníček et al. 2007; Lu et al. 2015). Turc et al. (2015) studied the interaction of a magnetic cloud with a bow shock by using a global hybrid simulation with a spatial resolution of 1 ion inertial length. Studies on foreshock and hot flow anomalies also require at least hybrid simulations (Omidi & Sibeck 2007; Omidi et al. 2010, 2014) instead of MHD models. A new hybrid-Vlasov code has been developed to investigate the ion distributions in the Earth's foreshock and to explain the THEMIS observations (Kempf et al. 2015; Palmroth et al. 2015). The hybrid-Vlasov approach ensures a uniform sampling of the ion distribution function in all spatial and velocity dimensions, as the full three-dimensional ion velocity distribution function is propagated in each real-space cell. However, it has a very high computational cost. In all of the hybrid models above, an artificial/anomalous resistivity must be included to generate dissipation on electron scales. Lipatov & Zank (1999) find that the different values of the resistivity can lead to very different shock structures.

In the global PIC simulation (Buneman et al. 1995; Cai et al. 2015; Peng et al. 2015b), which has higher numerical noise but lower computational costs than the hybrid-Vlasov simulations (Kempf et al. 2015; Palmroth et al. 2015), the accessibility to both ion and electron scales is automatically and self-consistently included. Buneman et al. (1995) created a three-dimensional electromagnetic full PIC code (TRISTAN) to simulate the solar wind–magnetosphere interaction. In recent years, Cai et al. (2015) parallelized the TRISTAN code and obtained a global structure of the magnetosphere with a resolution of ∼0.1 ion inertial length. This setup, however, has difficulty retrieving the nonstationary shock front as in previous local simulations (Lembège & Dawson 1987; Hada et al. 2003; Yang et al. 2009). Peng et al. (2015a, 2015b) recently improved the three-dimensional implicit PIC code (iPic3D) using a dipole magnetic field immersed in the flow of the plasma that showed the formation of a magnetosphere. The highest spatial and temporal resolutions used in their work are 0.05d_i and $0.15{\omega }_{{pi}}^{-1}$, where d_i and ${\omega }_{{pi}}^{-1}$ are the ion inertial length and the ion plasma frequency, respectively. In their simulation with the highest resolutions, the total size of the simulation box is 20d_i × 20d_i × 20d_i, and the dipole becomes very small. Therefore, the resulting bow shock is blurred and is on a scale of several ion inertial lengths, which is roughly the same as the shock front rippling scale observed in local PIC simulations (Hellinger et al. 2007; Lembège et al. 2009). Hence, it is difficult to retrieve the features of shock front nonstationarity. However, 3D global simulations with a high-performance configuration (high resolutions, large particle number per cell, large simulation domain, etc.) are still a very computation-consuming task. It is thought that the details of the bow shock and the magnetosphere can eventually be well obtained from the iPic3D code (Peng et al. 2015a, 2015b), and we expect to see the fascinating results in the near future.

In this paper, we use a two-dimensional (2D) explicit full PIC code to achieve a relatively high spatial and temporal resolution for simulating the nonstationary bow shock. A self-reforming curved bow shock is generated by the interaction of the solar wind with a 2D dipole field. The paper is organized as follows. In Section 2, we describe briefly the simulation model. The results from our global simulations are presented in Section 3, followed by a discussion of the results and conclusions in Section 4.

Previous high-resolution hybrid simulations (Δx = 0.1d_i) have shown that the artificial/anomalous resistivity η employed in hybrid simulations can strongly affect the shock front self-reformation process in the same plasma and Mach number conditions (M_A ∼ 2.1). Lembège et al. (2009) find that (1) for the high-resistivity run ($\eta ={10}^{-2}{\mu }_{0}{v}_{A}^{2}{\omega }_{{ci}}^{-1}$), the steepening of the ramp is restricted and cannot initiate a self-reformation even if the spatial resolution is relatively high; and (2) for the low-resistivity run ($\eta ={10}^{-4}{\mu }_{0}{v}_{A}^{2}{\omega }_{{ci}}^{-1}$), the self-reformation is partial in the sense that the foot amplitude is increasing but stays relatively weak and never reaches an amplitude comparable to that of the old ramp. Instead, the new ramp crashes down and restarts, reflecting new incoming ions. They conclude that using a high resistivity value and a low spatial resolution will stop almost simultaneously the self-reformation and the emission of nonlinear whistler waves. But for very high Mach number shocks (M_A ∼ 23), the self-reformation can appear in hybrid simulations with a high resistivity value ($\eta ={10}^{-2}{\mu }_{0}{v}_{A}^{2}{\omega }_{{ci}}^{-1}$) and a relatively high spatial resolution (Tiu et al. 2011). Hence, the choice of different resistivity can lead to different shock front structures, as mentioned by Lipatov & Zank (1999). In fact, it is difficult to tell how much artificial resistivity we should choose. In addition, the use of high-resolution hybrid simulations can be questionable in terms of computing costs when compared to 1D/2D PIC simulations performed with a reasonable mass ratio. Let us note that the question of accessibility to a small scale (lower than ion scale) is expressed differently in full PIC and hybrid simulations. In PIC simulations, the accessibility to both ion and electron scales is automatically and self-consistently included. Furthermore, the calculations of the electric field in hybrid and full PIC codes are different. The former usually use the generalized Ohm's law with an artificial resistivity, and the latter is based on the complete Maxwell equations without an artificial resistivity.

For that reason, a 2D full particle code is used to simulate the bow shock formation due to the solar wind–magnetosphere interaction. Yang et al. (2015) have already described the simulation code used for modeling of planar shocks in considerable detail; hence, only a brief description is given here. The simulation plane corresponds to the x–z plane with x along the solar wind flow direction (Sun–planet line) and z pointing along the dipole axis (as in the GSM coordinate system). Consequently, the x–z simulation plane corresponds to the noon-midnight meridian plane with z pointing northward. The solar wind Alfvén Mach number ranges from 4 to 8, and the value of β_i (the ratio of the plasma pressure to the magnetic pressure) is set from 0.01 to 2. Based on previous local simulation results (Lembège & Savoini 1992; Hada et al. 2003; Hellinger et al. 2007; Yang et al. 2013; Hao et al. 2014), the shock is expected to be in the supercritical regime (Lembège & Dawson 1987), and the shock front can be nonstationary. It is generally believed that the distance between the center of the dipole and the magnetopause at the subsolar point is the shortest when the interplanetary magnetic field (IMF) is southward. To save computational time and decrease the simulation box, the IMF initially lies in the x–z plane, points southward (−z), and makes a 90° angle with the x-axis. A reduced mass ratio m_i/m_e = 20 is chosen, which is a little larger than in global PIC simulations (Cai et al. 2015). The size of the simulation box is L_x × L_z ≈ 100d_i × 100d_i, and the grid consists of 4096 × 4096 cells. For reference, the temporal-spatial resolution used in previous hybrid and PIC simulations is summarized in Table 1. Local hybrid and PIC simulations show that the shock front self-reformation can be seen in high-resolution cases, i.e., Δx ≤ 0.05d_i (Lembège et al. 2009). The grid spacing used in our simulations is Δx = Δz = 0.025d_i = 0.11d_e = 0.68λ_De, where d_e and λ_De are the electron inertial length and the Debye length, respectively. The chosen temporal resolution is ${\rm{\Delta }}t=0.001{{\rm{\Omega }}}_{{ci}}^{-1}=0.015{\omega }_{{pi}}^{-1}$. The initial setups of the dipole and IMF fields are similar to those used in previous simulations (Omidi et al. 2005; Peng et al. 2015b). The particles are evenly distributed from the solar wind to the nightside of the dipole to speed up the magnetosphere formation. In our explicit code, the electric field and the magnetic flux function instead of the magnetic strength are updated at each time step. This will help us directly obtain the accurate magnetic field lines by plotting the contour of the magnetic flux function on the 2D simulation plane.

Table 1.  Summary of Simulation Setups of Hybrid and PIC Global Simulations of the Magnetosphere

Simulations	Particles for Each Species per Cell	Spatial Resolution Δx	Temporal Resolution Δt
1. Global hybrid simulations
Omidi et al. (2005)	⋯	0.5di	$0.0025{{\rm{\Omega }}}_{{ci}}^{-1}$
Omidi et al. (2010)	15	1di	⋯
Trávníček et al. (2005)	⋯	0.2–0.25di	$0.005{{\rm{\Omega }}}_{{ci}}^{-1}$
Richer et al. (2012)	⋯	3di	$0.05{{\rm{\Omega }}}_{{ci}}^{-1}$
Lu et al. (2015)	50 (private communication)	2di	$0.05{{\rm{\Omega }}}_{{ci}}^{-1}$
2. Global PIC simulations
Baraka & Ben-Jaffel (2007) (TRISTAN)	0.8	10λDe	⋯
Cai et al. (2015) (TRISTAN)	8	0.1di	$0.005{{\rm{\Omega }}}_{{ci}}^{-1}$
Vapirev et al. (2014) (iPic3D)	250	0.078di	$0.125{{\rm{\Omega }}}_{{ci}}^{-1}$
Peng et al. (2015a, 2015b) (iPic3D)	25	0.05–0.12di	0.1–0.15 ${{\rm{\Omega }}}_{{ci}}^{-1}$
3. Local PIC simulations (reforming shocks)
Savoini & Lembège (1994) (2D)	4	0.79λDe	$0.0012{{\rm{\Omega }}}_{{ci}}^{-1}$
Lee et al. (2004) (1D)	200	1λDe	⋯
Lembège et al. (2009) (2D)	4	0.0167di	$7.5\times {10}^{-5}{{\rm{\Omega }}}_{{ci}}^{-1}$
Matsukiyo & Scholer (2011) (1D)	100	0.89λDe	$1.2\times {10}^{-5}{{\rm{\Omega }}}_{{ci}}^{-1}$
4. This work	4	0.025di (0.68λDe)	$0.001{{\rm{\Omega }}}_{{ci}}^{-1}$

Download table as:  ASCIITypeset image

We present the overall structure of the bow shock and magnetosphere in the global PIC simulation. Our study is separated into three parts: (1) the time evolution of the bow shock and the entire magnetosphere, (2) the shock front self-reformation at different shock normal angles along the curved shock front, and (3) the impact of the plasma β value and upstream solar wind Mach number M_A on the nonstationarity of the bow shock.

3.1. Formation of the Bow Shock and Magnetosphere

Figure 1 shows the macroscopic and microscopic evolutions of the normalized ion number density log₁₀(N_i/N₀ + 1) in the meridional plane (i.e., in the simulation x–z plane). Initially, the magnetosphere undergoes an expanding stage. Figures 1(a)–(c) (right column) show the first stage of a bow shock formation. At $t=0.1{{\rm{\Omega }}}_{{ci}}^{-1}$ (cycle 100), the high ion density is located at the cusp region and behind the dipole center. At a later time ($t=1-3{{\rm{\Omega }}}_{{ci}}^{-1}$, cycles 1000–3000), more ions are accumulated in front of the dipole and the cusp region. A fraction of incident ions is reflected at the newborn magnetopause back toward the Sun. Instead of a perfect conductive wall, the solar wind ions are reflected by an elastic wall "magnetopause." The physical mechanism is similar to that in 1D local simulations of a collisionless shock by the reflection wall method. Figures 1(d)–(f) show the second stage of the bow shock formation. At $t=4{{\rm{\Omega }}}_{{ci}}^{-1}$ (cycle 4000), the magnetotail stretches in the x direction. The plasma sheet and magnetotail lobe are easily distinguished. The plasma mantle (light blue) around the lobe (dark blue) is also stretching, accompanied by the plasma sheet. At $t=6{{\rm{\Omega }}}_{{ci}}^{-1}$ (cycle 6000), the reflected ions ahead of the magnetopause are convected back together with the new incoming solar wind ions, and they are compressed in front of the magnetopause. Then a dense magnetosheath is formed. The mature shock is first generated at the subsolar point. At $t=8{{\rm{\Omega }}}_{{ci}}^{-1}$ (cycle 8000), the bow shock is almost completely formed and the magnetotail continues to stretch. Figures 1(g)–(i) (cycles start from 11000) show that both the bow shock and the entire magnetosphere have now reached a quasi-static state with minimal changes until the end of the simulation. In this low-β_i (=0.01) and high Mach number (M_A = 8) case, a nonstationary bow shock front with ion scale ripples of the order of ion inertial length d_i is observed.

Zoom In Zoom Out Reset image size

An overview of the bulk velocity and electromagnetic field components of the bow shock and the magnetosphere at $t=14{{\rm{\Omega }}}_{{ci}}^{-1}$ is shown in Figure 2. The magnetic field line (in black) is superimposed on the contours of the ion bulk velocity profile V_z (Figure 2(c)). In the magnetosheath (downstream region of the bow shock), Figures 2(a)–(f) show that the bow shock reduces the super-Alfvénic solar wind speed and the magnetopause deflects the downstream plasma flow. The sheath plasma diverts around the magnetosphere. The bulk velocity components V_x of both ions (Figure 2(a)) and electrons (Figure 2(d)) become sub-Alfvénic. One striking point is that the velocity component V_z (let alone the total bulk speed) can still remain super-Alfvénic in the sheath region due to the deflection motion. The shock only decreases the upstream inflow speed along the shock normal. Furthermore, the velocity moments of ions (Figure 2(b)) and electrons (Figure 2(e)) at the magnetopause and the magnetotail plasma sheet are in opposite directions, and they contribute to the magnetopause and cross-tail currents. In the magnetotail, the plasma flow and electromagnetic fields are quite similar to those obtained in local magnetic reconnection simulations (Daughton et al. 2006; Fu et al. 2006; Lu et al. 2010; Huang et al. 2014; Guo et al. 2015). Magnetic reconnection takes place in the magnetotail because of the southward magnetic field used in the simulation.

Zoom In Zoom Out Reset image size

3.2. A Baseline Case of a Nonstationary Bow Shock

To examine the impact of shock normal angle θ_Bn on the shock microstructure and the shock reformation process in our simulation, we show in Figure 3 the time evolution (stack plots) of the magnetic field components B_x, B_y, and −B_z for the shock profiles A, B, C, and D measured at different locations (from perpendicular shock to oblique shock regions) of the bow shock. The spatial ranges of measurements are marked by the red lines in Figure 2(g). In all panels, the same scales are used for the stack plots.

Zoom In Zoom Out Reset image size

Figures 3(a)–(c) show the magnetic field components of the shock profile A, which is a nearly perpendicular shock (${\theta }_{{Bn}}\approx 90^\circ$) measured at z = 51d_i in the bow shock (see Figure 2(g)). The upstream southward IMF B field is in the −z direction, and thus the main compressed magnetic field component in the shock transition is B_z. The Mach number of the incident solar wind along the shock normal is V_SW × sinθ_Bn = 8V_A. Figure 3(c) shows that the shock is in the supercritical regime and the front is undergoing self-reformation. For example, at about $t=10.8{{\rm{\Omega }}}_{{ci}}^{-1}$, a foot (at x = 7d_i) is formed ahead of the shock ramp (x = 8.5d_i). The newborn foot propagates together with the injected solar wind toward the right-hand side. Finally, the foot grows and becomes a new ramp at about $t=12.5{{\rm{\Omega }}}_{{ci}}^{-1}$, and the process repeats itself. Different reformation cycles are marked by red arrows. The observed period of the reformation cycle is consistent with previous local hybrid and PIC simulations (Lembège & Savoini 1992; Matsukiyo et al. 2003; Yang et al. 2009; Yuan et al. 2009).

Figures 3(d)–(f) display similar plots for the shock profile B, which is a quasi-perpendicular shock (${\theta }_{{Bn}}\approx 70^\circ$) measured at z = 66d_i (see Figure 2(g)). First, the magnetic field component B_x measured in the southern part of the bow shock is almost positive (e.g., the profile A), and that measured in the northern part is negative (profiles B, C, and D). It depends on the sampling locations (refer to Figure 2(g)) because the IMF is curved inside the magnetosheath. Second, the Rankine–Hugoniot equations specify only the change in magnetic field and plasma parameters from one side of the shock to the other and do not specify how these parameters change inside the shock. Thus, while the upstream magnetic field, downstream magnetic field, and shock normal $\hat{n}$ must all be in the same plane, the magnetic field within the shock layer could deviate from the plane and become non-coplanar. The magnetic field in fact frequently deviates from this plane in observations of Earth's bow shock (Friedman et al. 1990), in numerical simulations of collisionless shocks (Thomsen et al. 1987), and in theoretical models that include at least two fluids (Gedalin 1996). In our simulation, the component B_y is the non-coplanar magnetic field component (Figure 3(e)). Third, the Mach number of the incident solar wind along the shock normal is V_SW × sinθ_Bn = 7.5V_A, which is smaller than that in profile A, but the shock is still in the supercritical regime. Hence, the shock front is still nonstationary, and the self-reformation process takes place. However, the maximum amplitude of the profile $| {B}_{x}|$ at the overshoot is lower than the profile for A.

Figures 3(g)–(i) and (j)–(l) display the stack plots for shock profiles B and C, respectively. These two shocks are sampled at the oblique shock region (z = 81d_i and 96d_i) of the bow shock. Their corresponding shock normal angles θ_Bn are about 50° and 30°, respectively. The Mach numbers of the incident solar wind along the shock normal for these two shocks are V_SW × sinθ_Bn = 6.1V_A and 4V_A, respectively. Due to the decreased upstream solar wind speed in the shock normal direction, shock compression at the wings of the bow shock is weaker than that at the subsolar point, as expected.

In contrast to the 1D shock model, the downstream region of a 2D bow shock is more complicated, because the incident and transmitted solar wind particles at the perpendicular region can convect to the downstream region of the quasi-perpendicular and oblique shocks. The bulk velocities of ions and electrons in the downstream have already been shown in Figures 2(c) and (f). To show the single particle motions in the magnetosheath, we have traced the particle trajectories from the simulation. Figures 4(a) and (b) show the trajectories of solar wind ions and electrons injected at the nearly perpendicular shock (i.e., subsolar point). The locations of the bow shock (black curve) and the magnetosphere (black dashed) are also shown for reference. These locations are obtained by tracing their visual outlines, shown in Figure 2. Initially, the ions and electrons drift together in the solar wind toward the bow shock (black dots). Later on, the incident particles reach the bow shock (blue dots). Then a fraction of the incident ions are reflected by the shock, leading to a self-reformation of the nonstationary shock front, and almost all of the electrons are directly transmitted across the shock (green dots). At a later time, both ions and electrons diffuse into the downstream region and the particles are convected to the downstream regions of quasi-perpendicular and oblique shock regions (yellow dots). Finally, the particles become more dispersed and a fraction of electrons leak out from the magnetosheath at the wings of the bow shock, that is, at the oblique shock region (red dots). Figures 4(c) and (d) display similar plots for ions and electrons with different initial locations (the particles injected at the northern part of the bow shock). A large number of ions and electrons can be reflected at the oblique shock front. At a later time, a fraction of ions can enter the magnetosphere on the nightside. Some electrons can be trapped by the dipole field, and these trapped electrons include bounce motion on the magnetic field on the dayside. These particles could affect the ring current in the inner magnetosphere, but the analysis requires at least a 3D model, which is beyond the scope of this article. Other electrons are convected to the southern part of the bow shock.

Zoom In Zoom Out Reset image size

3.3. Impact of $\beta$ and ${M}_{A}$ on the Shock Nonstationarity

In this part, we study the impact of the plasma β_i and the solar wind Mach number M_A on the shock front nonstationarity. It is generally thought that the quasi-perpendicular shock front self-reformation and local instabilities can be affected by the ion beta and the upstream inflow speed (Hellinger et al. 2002; Scholer & Matsukiyo 2004; Yang et al. 2013). The results of the baseline case (Run 1: β_i = 0.01 and M_A = 8) have already been shown in Sections 3.1 and 3.2 above. Two comparative cases have been carried out: (1) Run 2, β_i = 2 and M_A = 8, and (2) Run 3, β_i = 2 and M_A = 5. The other setups are kept unchanged. In Run 2, we increase the ion beta value only. In Run 3, we increase the ion beta but decrease the solar wind Mach number. Figure 5 shows the main magnetic field component B_z from Runs 1, 2, and 3, respectively. By comparing Figures 5(a) and (b), we find that the ripples on the bow shock surface become blurred and the spatial scale of the rippling wave length becomes larger. However, the bow shock is still nonstationary. The magnetosheath at the subsolar point becomes thicker in Run 2. In addition, the magnetotail reconnection is visible in Runs 1 and 2. In contrast, Figure 5(c) shows the low Mach number case (Run 3). In this case, the shock front is quasi-stationary, and the magnetotail reconnection disappears. To double-check the shock front nonstationarity in different runs, the main magnetic field component B_z is sampled at different locations of the bow shock. Figures 6(e)–(h) show the stack plots of the shock profiles measured at different locations of the bow shock obtained in Run 3 at different times. The corresponding plots for Run 1 are also shown for reference in Figures 6(a)–(d). It should be clear that the bow shock with a high plasma beta value and a low solar wind speed is quasi-stationary. This behavior is consistent with previous 1D and 2D local simulations of planar shocks (Scholer & Matsukiyo 2004; Yang et al. 2013).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In this paper, we have presented a 2D global PIC model to study the interaction between the solar wind and magnetosphere. We demonstrate the capability of the model focusing on the kinetic effects associated with the bow shock nonstationarity and magnetic reconnections. By sampling the shock profiles at different times and different locations of the bow shock and tracing the particles, we have shown the following results:

• 1.  
  We identify different stages in the macroscopic evolution of the bow shock and the entire formation of the magnetosphere in the meridian plane. The macrostructures reach a quasi-stable state after several ion cyclotron cycles, consistent with previous 3D implicit PIC simulations (Peng et al. 2015b). Furthermore, a relatively high temporal-spatial resolution employed in our simulation provides an opportunity to examine the microstructure of the bow shock and the magnetic reconnection.
• 2.  
  In the southward IMF condition, the shock around the subsolar point is quasi-perpendicular, and the shock at the wings is oblique. The angle ${\theta }_{{Bn}}$ and solar wind speed along the shock normal decrease with the distance away from the subsolar point. At quasi-perpendicular regions of the shock, a self-reformation of the shock front is found and the cyclic period of the reformation is similar to those observed in 1D and 2D local simulations. At the oblique region, the shock becomes weak due to the slow solar wind speed in the shock normal direction.
• 3.  
  Different from 1D and 2D local simulations, particles in the downstream region are more complicated at the bow shock. By tracing the ions and electrons in the simulation, we find that the solar wind ions injected at the subsolar point will become diffuse and can be convected to the cusp region and the downstream of the oblique region of the bow shock. Electrons are more diffusive than ions in the downstream region. The ions injected at the northern part of the bow shock will fill in the northern part of the magnetosheath. A fraction of electrons injected at the same place will be convected to the southern part of the magnetosheath. The particles injected and reflected by the quasi-perpendicular shock region can affect the local velocity distributions and electromagnetic fluctuations in the transition and downstream regions of the oblique shock at the wings of the bow shock. This feature makes the oblique shock regions become a little different from those observed in local planar shock simulations.
• 4.  
  The impact of the upstream plasma beta β_i and solar wind Mach number M_A on the bow shock and the magnetosphere show that the shock front becomes quasi-stationary in the high beta and low Mach number cases. In addition, the magnetotail reconnection disappears in this quasi-stationary case due to the low compression of the magnetotail under a low solar wind pressure condition. If we keep the Mach number of the incident solar wind unchanged and only increase the beta value, the bow shock is still nonstationary. Contrasting the low beta case, the ripple scale becomes larger in the high beta case.
• 5.  
  To confirm the impact of the solar wind condition on the bow shock front nonstationarity, we also studied the time evolution of shock profiles sampled at different locations of the bow shock. The results support the conclusions in point 4.

It is worth noting that there are two limitations in PIC simulations: (1) the unrealistic mass ratio m_i/m_e and (2) the unrealistic ratio of electron plasma frequency to cyclotron frequency $\left(\frac{{\omega }_{{pe}}}{{{\rm{\Omega }}}_{{ce}}}=\tfrac{c}{{v}_{A}}\sqrt{\tfrac{{m}_{e}}{{m}_{i}}}\right)$. Krasnoselskikh et al. (2013) note that the electric field fluctuations are overestimated in PIC simulations with low values of ω_pe/Ω_ce. The impact of these two parameters on the shock front self-reformation has been investigated by 1D PIC simulations (Matsukiyo et al. 2003). They conclude the following: (1) In small ion-to-electron mass ratio runs, the reformation is due to the accumulation of gyrating reflected ions. Furthermore, at the extremely small mass ratio, the Buneman instability is generated in the foot. In the realistic mass ratio run, however, the modified two-stream instability excited in the foot leads to the reformation. Hence, the self-reformation also occurs in the realistic mass ratio run, but the associated instability is changed. Of course, a higher mass ratio allows for an easier separation between ion and electron scales but requires relatively large computer capacities, in particular for 2D or 3D simulations. (2) The self-reformation is not a low ω_pe/Ω_ce process but occurs also in ${({\omega }_{{pe}}/{{\rm{\Omega }}}_{{ce}})}^{2}\gg 1$. Nevertheless, they do mention that in the solar wind at Earth's orbit the quantity of ω_pe/Ω_ce is 100–200. However, in most simulations the value of ω_pe/Ω_ce is assumed to be of the order of 1, i.e., simulations are done for shocks in the strongly magnetized condition because of computational constraints of PIC codes. In summary, recent 1D simulations have shown that the shock front self-reformation can occur even in high mass ratio and high plasma-to-cyclotron frequency ratio conditions. The impact of such parameters on self-reformation in 2D and 3D simulations is still open due to the enormous computation burden.

Future work will include the background turbulence in the solar wind because hybrid simulations show that the turbulence can affect the bow shock and magnetosphere structures (Guo & Giacalone 2012; Karimabadi et al. 2014). Furthermore, it is interesting to see the impact of the third population of ${{\rm{O}}}^{+}$ ion outflows (Seki et al. 2001; Lennartsson et al. 2004; Wiltberger et al. 2010) on the magnetotail reconnection by using the global PIC code. This has implications for understanding how planets begin to experience a runway greenhouse effect (Zhang et al. 2012).

The authors are grateful to D.S. Cai, X.Y. Wang, Y. Lin, and G. Lapenta for helpful discussions on numerical methods and global simulations. This research was supported by NSFC under Grant Nos. 41574140 and 41374173, the Recruitment Program of Global Experts of China, and the Specialized Research Fund for State Key Laboratories of China.