The effect of cavity density on the formation of electrostatic shock in the lunar wake: 1‐D hybrid simulation

One‐dimensional hybrid simulations are carried out to study the plasma refilling process in the lunar wake. Previous theoretical and simulation studies have shown that ion‐ion acoustic (ⅡA) instability can be initiated and electrostatic shock can be formed under the condition Te≫Ti . We find that the time evolution of ⅡA instability and the formation of electrostatic shock strongly depend on initial cavity density. The initial position of the electrostatic shock is dependent on the ratio between initial cavity density and background solar wind density, i.e., the farther away the initial position, the lower is the ratio. When the initial cavity density is low enough, the density and electric field profile across the wake become much complex. Meanwhile, the back‐to‐back electrostatic shock is unstable in the case of lower cavity densities; at the late evolution stage, a new shock‐like structure can be formed at the central region of the lunar wake.


Introduction
Known as a non-magnetized non-conducting celestial body without an atmosphere, the Moon absorbs the solar wind particles that impact on its surface, while the Interplanetary Magnetic Field (IMF) remains unperturbed until it reaches the region on the lunar night side in the rarefied wake of non-equilibrium plasma. As a result, a nearly plasma-free cavity is formed just behind the obstacle (Moon). Therefore, a plasma refilling process is expected due to the pressure imbalance between the cavity and the ambient solar wind plasma. Since the thermal velocity of electrons is much higher than the solar wind speed, at the very initial stage fast electrons enter the whole region of the geometrical shadow, charging it negatively. As soon as the negative potential reaches ~ (where is the electron temperature) further charge separation stops, and electrons continue to fill the cavity along with solar protons that are accelerated to the ion-acoustic velocity (where is the proton mass) by the negative potential in the cavity (Michel, 1968). One consequence of the refilling process predicted by theory (Michel, 1967(Michel, , 1968 is the formation of a shock wave at the place where the collapse of the cavity is halted. However, such a shock-like structure has never been observed or modeled until recently. Israelevich and Ofman (2012) successfully used 1-D hybrid simulation to produce the shock wave trailing the moon. By setting the electron temperature to be larger than the ion temperature, they demonstrated that an ion-acoustic-like instability can be excited and acts as the main mechanism of shock dissipation. In their Case A (expansion along the magnetic field), a pair of electrostatic shocks can be formed at the center of the cavity. It is necessary to point out that in 1-D hybrid simulation, the quasi-neutral condition requires that the V ey ≡ V iy (here the subscripts i and e stand for ion/proton and electron); i.e., the relative motion between ion and electron in Y direction (along the simulation box/background IMF) is exactly equal to 0. Therefore, in a plasma with two ion populations, the relative motion between the l-th (l = 1, 2) component ions and electrons is smaller than that between the two ion components. On the other hand, as described by Gary and Omidi (1987), the threshold of ionion acoustic instability lies well below the threshold of ion-electron instability, such that the ion-acoustic-like instability in Israelevich and Ofman (2012) is essentially the ion-ion acoustic (IIA) instability, which is excited by the relative drift between two ion populations under the condition T e ＞ T i . , 1991), where V d is the relative drift speed between two ion populations, θ is the angle between V d and the wave vector k of the waves with the largest linear growth rate. Note that the speed of ion beams or the speed of cavity collapse in 1-D hybrid simulation is approximately , or the thermal velocity of ions under the electron temperature, such that the relative drift speed between the ions coming from both sides of the wake is approximately . Therefore, theoretically, the IIA instability should be stable under this condition. On the other hand, to avoid numerical difficulty, a tenuous but finite ion population should be pre-setup in the cavity, which is assumed to be nearly a vacuum. Those cavity ions might interact with the ion beams accelerated by the potential well and thus play a crucial role in the formation of electrostatic shocks and in determining the evolution of the structure of the lunar wake.
To address the above problems, we carried out a series of 1-D hybrid simulations. The outline of the rest of the paper is as follows: First, in Section 2, we give a short description of our simulation code and model. In Sections 3.1 and 3.2, respectively, we present simulations of IIA instability excited by two equal dense ion beams, and of plasma refilling. Finally, in Sections 4 and 5 we discuss the simulation results.

Simulation Model
The hybrid code used in this paper is adapted from Wang XY and Lin Y (2003); it has been successfully applied to study the ionbeam instability of foreshock region turbulence (Li HM et al., 2013). In the hybrid model the ions are represented as particles, neglecting collisions, while the electrons are described as a finitetemperature massless fluid that maintains the quasi-neutrality of the plasma (the isothermal electron fluid assumption is applied in the current paper). A second-order splitting scheme is applied to reduce the noise and artificial numerical heating. We also invoke uniform artificial resistivity to damp out small-scale fluctuations. Note that for all runs in the current paper, we carefully set up the initial condition such that only electrostatic waves/instabilities can be excited; therefore, the electric current is negligible, as is, also, the electric field caused by the artificial resistivity.

Interaction Between Two Equal-Dense Ion Beams
We start with simplest case: interaction between two ion beams under the condition . In this section, the length, time, and charge density are normalized to , , and , respectively, where is the proton inertial length, is the angular gyro-frequency of the proton, and is the initial proton number density. Since we are concerned with the interaction of the ion beams coming from both sides of the lunar wake, we set the two ion beam densities to be equal. The plasma of each ion component is set to be in all runs in this section to make sure that no electromagnetic instability (e.g. ion-ion right hand resonant/nonresonant instability) occurs. The ratio between electron and ion temperatures is 20. The relative drift speeds are , respectively, for the four runs in this section.
The charge density time evolutions of Runs 1-4 are shown in Figure 1. It can be seen that when is small (Run 1, panel (a)), the IIA mode is stable; the variations of charge density in space are simply the thermal fluctuation. In Run 2 (panel (b)), we increase to be ; ion-ion acoustic instability is excited, as indicated by the large amplitude density and also the electric field fluctuations (not shown here), while the magnetic field fluctuations remain within the thermal noise level. The relative drift speed is increased to in Run 3 (panel (c)); again, IIA instability can be initiated, though the growth rate is smaller than that of Run 2 and, at the final stage, the charge density amplitude is substantially lar-  JIn Y and Pang Y: Electrostatic shock in lunar wake ger than that in Run 2. Eventually, when we increase to be , the two-beam plasma system (panel (d)) is stable again, just as in Run 1. From Runs 1-4, one conclusion can be drawn: the relative speed of the two ion components is the crucial factor determining whether IIA instability can be excited or not. If the relative speed is neither too high nor too low, IIA instability cannot be excited. V d (Note: We recognize that the density ratio between the two ion components and the temperatures of two ion components play roles in determining the criteria for IIA instability (Gary and Omidi, 1987). However, this topic falls beyond the scope of the current research, which focuses on relative drift speed ; the effects of other variables should be studied in the future.)

Electrostatic Shock in the Lunar Wake
We have demonstrated that IIA instability cannot be excited in 1-D hybrid simulation if the relative drift speed of ions is exceed certain threshold V threshold . Although we don't know the exact number of V threshold , it must lie between and based the simulation results of Section 3.1 (Run 3 and 4). Also, as mentioned before, the relative speed of ion beams accelerated by the potential well in the Lunar wake can reach as high as . Therefore, it seems that the IIA instability cannot be initiated and serves as the main dissipation mechanism of electrostatic shock. To understand how the electrostatic shocks are formed in the wake, several simulations are carried out in this section.
First, to understand the formation mechanism of electrostatic shock in the wake, we reproduce the results of Case A of Israelevich and Ofman (2012). We will refer it as "Case A" hereafter. Typical parameters of the solar wind plasma are set to be the same as those in Case A, which are: n 0 = 5 cm -3 , , and . The initial cavity density is approximately n c = 0.5 cm -3 . The temperature ratio between electron and ions is 20 and the diameter of the Moon is 3500 km ( ). We have used 1600 particles for each cell to reduce statistical noise to acceptable levels (i.e., the random noise in velocity space is much less than the thermal distribution). The method of weighting is applied to simulate a "nearly-empty" cavity behind the Moon. For simplicity, a periodic boundary condition is applied and we utilize a large simulation box to make sure that such a choice of boundary conditions has no/negligible effect on the simulation results. The simulation box consists of cells and the total length is about . The time evolution of charge density (electron density) in Run A is shown in Figure 2. Note that we convert the time coordinate to distance along the -axis: and we assume that the velocity of the solar wind is . Our result is qualitatively the same as that of Case A of Israelevich and Ofman (2012). A pair of electrostatic shocks is formed around the wake axis ( ). To study the shock formation process in detail, profiles of charge density and ion distribution along the simulation box at different positions (time) are plotted in Figure 3. Each column represents a certain position (time). The positions (time) of those line-cut plots in Figure 3 are marked by the vertical dashed black line in Figure 2. Note that in this paper we apply the isothermal electron assumption to complete the kinetic equation set. Therefore, the charge density also can be viewed as a good proxy of the electric potential (unscaled) along the simulation box. The black (red) dots in the second and third rows represent the ions that were located outside (inside) the cavity at or and for simplicity we will refer to the "black (red)" ions as "solar wind ions (cavity ions)". The charge densities are shown in the first row; the green, black, and red lines in the first row are charge density integrated over all, solar wind, and cavity ions, respectively. The green and cyan dashed lines in the second and third rows are Y-component bulk velocity integrated over solar wind and cavity ions, respectively. Note that we evaluate separately the bulk velocity of solar wind ions from each side of the wake. The ion acoustic speeds are marked by the horizontal dashed magenta lines in the second and third rows. At the very initial stage (X = 2600 km or ), a piecewise-like charge density (electric potential) profile can be seen in panel (a) and, correspondingly, ion velocity jumps can be found across the strong density gradient region (strong electric field region). The fast ion beams ( ) have practically no interaction with the ambient ions (both cavity and solar wind ions) until the end of Run A. It is interesting to note that, although the two counterstreaming slow ion beams haven't run into each other, IIA instability has been initiated at the off-axis region (around ), which is due to the relative drift ( ) between the solar wind ions and cavity ions. The initiation of such instability is indicated by "ripples" in the distribution plot of cavity ions (see the rectangular domain surrounded by the blue solid lines in panel (b), Figure 3). At X = 5700 km or , the instability evolves into its nonlinear stage, as indicated by the vortex-like distribution of cavity ions (see the rectangular domain surrounded by the blue solid lines in panel (e), Figure 3). By comparing panel (b) and panel (e), one might notice that the bulk speed of the slow beam is slightly decreased to just above the ion acoustic speed due to the energy exchange between slow beam ions and cavity ions. Meanwhile, the two slow ion beams have met with each other and a high potential wall begins to build up around the center of the wake. It can be found that the contribution of cavity ions is dominant in the setup of the "central potential wall". At or , the central potential wall is fully developed and the density at the wake axis reaches as high as . The contribution t = 23Ω −1 P of solar wind ions becomes important in the building-up of the high potential (density) wall, but the contribution of cavity ions cannot yet be ignored. Also, the bulk speed of the slow ion beams around the wake axis decreases to become slightly smaller than the ion acoustic speed. Such decreases are probably due to (1) the interaction of cavity and solar wind ions and (2) the existence of the central potential wall. Furthermore, IIA instability is initiated by the relative drift between the slow ion beams that come from the two sides of the lunar wake (see the rectangular domain surrounded by the blue solid lines in panel (i), Figure 3). At X = 12000 km or , the two slow ion beams strongly interact with each other and a pair of electrostatic shocks begins to form. In simulation Run A, we study in detail the process of formation of electrostatic shocks in the lunar wake. Our results show that the cavity ions play an important role in decelerating the ion beams down to the threshold of IIA instability. A following question is: what will happen when we reduce the initial number density of cavity ions. Two more runs are carried out to answer the question. In Runs B and C we keep all parameters unchanged except for the initial density of cavity ions, which is changed to and , respectively. The result of simulation Run B is shown in Figure 4. In general, electrostatic shock can still be formed; however, it takes more time to evolve (farther away from the Moon). The shocks begin to form at about X = 12000 km or in Run B, later than the comparable time (X = 12000 km or ) in Run A. In Run C, we further reduce the initial cavity density by half; the charge density evolution is shown in panel (a), Figure 5. The time evolution/spatial distribution is greatly differ-   tions of those shocks are more off-center (farther away from the wake axis) than in the above two runs. On the other hand, the overall density structure across the wake becomes much more complex. Multiple "density clumps" can be seen among the central shocks and the outer boundary of the wake expansion cone; this will be discussed later. First of all, we focus on the physics processes happening in the central region. Panel (b) of Figure 5 shows in more detail the area marked by the black dashed rectangle in panel (a). The initial position of the electrostatic shocks is approximately ; they then expand toward and away from the wake axis. Four rarefaction waves (R waves) attached to each shock can also be seen. The inner two R waves (marked by black lines) and also the inner two shocks (marked by black dashed lines) expand toward the wake axis. At roughly X = 30000 km or , the two R waves run into each other. The consequence of the interaction of the two R waves is formation of a low-density "new wake" at the central region. Shortly after (at about or ) the two R waves run into each other, the expansions of the two inner shocks are halted and begin to move outwards (away from the wake axis). Eventually, a new pair of electrostatic shocks are formed at about X = 40000 km or .
To understand the difference between Run C and the first two runs (Runs A and B), it is necessary to check the evolution of the particle distribution, as shown in Figure 6. At approximately X = 5200 km or (see panels (a) and (e), the slow ion beams from both sides run into each other, which leads to enhancement of charge density at the center; in Run C, the contribution of solar wind ions to such an enhancement is dominant over that of cav- ity ions. Similar to Run A, IIA instability has been initiated due to the relative drift of solar wind and cavity ions (black and red dots). However, since in Run C the abundance of cavity ions is small, the slow ion beams hardly get decelerated by the initiation of IIA instability. The relative drift speed between the two solar wind ion beams is far above around the wake axis. Therefore, the interaction between the two slow ion beams is negligible and they simply penetrate each other. As a result of the lack of deceleration mechanism and tenuous cavity ions, the profile of charge density in the central wake becomes more flattened (0.4-0.5n 0 ) and no high potential wall can be built up (see panels (b) and (f). At roughly X = 20300 km or , when the ion beams with slower bulk velocity (close to ) begin to interact with each other, two pairs of electrostatic shock begin to form at around . Those shocks are fully developed by approximately X = 25500 km or ; the rarefaction waves are denoted by the ion holes in phase space -see the blue and red shaded areas in panels (d) and (h). The inner two shocks are moving toward each other and the inner two ion holes (rarefaction waves) merge into one at roughly X = 32800 km or . The scale of the new ion hole is about 1000 km or , much larger than that in Run A and B; since the density inside the ion hole is around , X = 5200 km X = 12000 km X = 20300 km X = 25500 km (e) (f) (g) (h) X = 32800 km X = 37500 km X = 44300 km X = 51600 km (i) (j) (k) (l) X = 32800 km X = 37500 km X = 44300 km X = 51600 km (m) (n) (o) (p) Figure 6. Line-cut plot of charge density and particle distribution of Run C. See the description in context for detail. we will refer it as the "new wake" hereafter. On the other hand, those fast beams going across the wake axis have interacted with the background solar wind ions due the decrease in the relative velocity between them, as indicated by the strong disturbances of phase space distribution and density at the left and right of panels (i) and (m) (see the blue shaded areas). By comparing panels (i), (m) with (j), (n), one can find that the expansions of the inner two shocks have halted and reversed, which is probably due to the low density of the new wake. The lower the density of the wake, the higher will be the shock potential wall. Therefore, the ions entering into the downstream of the shock from the wake are fewer than those entering into the "new wake" from the downstream of the shock. Another consequence is that the average density of the new wake increases to and a high potential wall at the wake axis begins to build up. Such a high potential wall continues to grow; eventually a new pair of shocks is formed at the central region -see panels (k) and (o). It is necessary to point out that by the end of the simulation ( ), the two ion holes that are closest to the wake axis are still unstable. One can see from panels (i) and (p) that a density hump (high potential region) marked by the black arrows begins to form in the center of each ion hole. Therefore, it is reasonable to anticipate that new high density region/electrostatic shock can be formed if the simulation is continued. The plasma refilling process is a long standing issue in efforts to understand the lunar wake. Due to the strong electric potential gradients at the boundaries between the solar wind and the nearly-empty wake, ions can be accelerated to super-sonic/superacoustic speed. From the thermodynamic equilibrium condition: we can estimate roughly the negative potential inside the wake by means of the equation , where and are the solar wind number density and wake number density, respectively (Israelevich and Ofman, 2012). Note that during the refilling, the ratio will getting smaller and smaller. It can be seen from the above equation that, if the initial cavity/wake density is large enough, i.e. the is small enough, the speed of accelerated ion beams will not exceed ion acoustic speed. Therefore, there will be no shock to form and IIA instability to be excited. The plasma refilling process will be less drastic and in a more MHD (MagnetoHydroDynamic) way: Slow expansion waves will bring solar wind ions into the wake and eventually make the wake density the same as the adjacent solar wind density.

Discussion
It is necessary to discuss the role of IIA instability in the form of electrostatic shock. Similar to Earth's bow shock, super-Alfvenic solar wind is just one necessary condition for the formation of bow shock; another necessary condition is Earth's Magnetosphere, or an obstacle. In Run B, for example, consider the IIA instability initiated in the off-axis region, as marked by the blue rectangle in Figure 3, panel (e): the IIA will cause fluctuations in the ion number density. Note that in our hybrid mode, electron number density is supposed to be the same as local ion number density while the electron temperature is assumed to be a constant. Therefore, the electron pressure is proportional to the ion num-ber density; there will thus be a local electric field pointing from the regions of high density toward those of low density. Those high density regions can be viewed as potential walls that can decelerate ions and will cause additional accumulation of ions and increases of local ion number density. Eventually, in the late nonlinear stage, two large potential walls are built (see panel (j) of Figure 3), and the electric fields located at the boundaries of the potential walls are strong enough to serve as the dissipation mechanism of shock, and shock-like structures are formed.
It is also interesting to point out that, in the final stage of Run C, multiple finite-width potential walls or high density regions are formed. Two adjacent potential walls will accelerate ions toward each other, the speed of those accelerated ion also can exceed the local ion acoustic speed. Therefore, two electrostatic shocks can be formed in between the two adjacent potential walls, or in another words, each finite-width potential wall/high density region in the wake can be viewed as a pair of back-to-back electrostatic shocks.
Another interesting question is What will happen if the initial cavity density is chosen to be exactly 0? This question cannot be answered by hybrid simulation, due to limitations of the numerical method itself. However, we still can make some postulations based on the results of simulation Runs A, B, and C. First of all, due to the lack of initial cavity ions, the IIA instability can be excited only by the interaction between the two ion beams coming from both sides of lunar wake. The IIA instability is unlikely to be excited at the central region of wake, as we have previously shown. The ion beam may possibly interact with those relatively slow ions after they have passed each other. There is a great chance that the initial IIA instability will be excited in the off-axis region, so that the initial strong potential wall will be formed there, at the off-axis region. The distance between the two initial potential walls will determine the late evolution: i.e., to form a pair of back-to-back electrostatic shocks as observed in Run B, or form multiple pairs of back-to-back shocks, as observed in Run C. If the distance is large enough (larger than the wave length of the linear IIA mode), then the region between the two initial potential walls will be unstable. New potential walls might be generated in late stages of the evolution.

Summary and Conclusions
T e ≫ T i 0.05n 0 In the above sections, we study the evolution of density structures and formation of electrostatic shocks in the lunar wake, via 1-D hybrid simulation under the condition . We have found that the amplitude of the initial cavity density plays a crucial role in determining how the density structures/electrostatic shocks evolve in the plasma refilling process of the lunar wake. Our simulation results show that the ion-ion acoustic instability between the solar wind ion beams (accelerated by the strong electric field located at the wake boundary) and the cavity ions is excited at the initial stage. If the initial cavity density is large enough ( ), the energy exchange between solar wind ions and cavity ions will be strong enough to substantially decelerate solar wind beam ions; meanwhile, large initial cavity densities (Runs A and B) also contribute significantly to the build-up of high density regions (high potential walls) at the center of the wake, which can addi- tionally decelerate solar wind ion beams. Eventually, the two decelerated solar ion beams will interact with each other, leading to formation of electrostatic shock. On the other hand, if the initial cavity ions (Run C) are sufficiently tenuous, fast ion beams simply penetrate each other; electrostatic shocks will be formed when the slower ion beams interact with each other. The initial positions of the shocks in Run C are more "off-axis" and farther away from the Moon than those in Runs A and B, which can be understood by the fact that the threshold and growth rate of ion-ion acoustic instability are strongly affected by the relative drift speed, density ratio, and temperature ratio between the two ion beams. However, because a full parametric survey of the growth rate and threshold of the ion-ion instability is beyond the scope of current paper, we simply point out that the threshold of the ionion acoustic mode can be greatly reduced if the densities of the two ion beams are unequal (Gary and Omidi, 1987), which probably leads to the deviation in initial positions of shocks. In Run C, the fast ions (above ), which have gone across the wake axis before the formation of electrostatic shocks, will interact with the slow-moving ions on the other side of the wake, and ion-ion acoustic instability will be excited and developed. Another interesting finding is that in 1-D simulation, the upstream regions of two "back-to-back" electrostatic shocks are unstable. The reflected ions from both shocks will form an ion hole at the central wake. Depending on the average density and scale, the ion hole can collapse into smaller ones or remain stable. It is necessary to point out that in 1-D simulation no oblique modes are allowed, which could have an important effect on the evolution of the plasma refilling process of the wake (Karimabadi et al., 1991). Meanwhile, fully kinetic simulation also shows that the electron two-stream instability can substantially disrupt the ion beams (e.g. Farrell et al., 1998;Birch and Chapman, 2001;Zhou M et al., 2014Zhou M et al., , 2018. Those issues need to be addressed in the future.