Ion dynamics in laser-produced collisionless perpendicular shock: one-dimensional particle-in-cell simulation

1.   Introduction

Collisionless shock is a nonlinear physical process commonly observed in space and astrophysical plasma [1–4], and it is believed to generate high-energy charged particles through diffusive shock acceleration (DSA) [5–7]. In DSA, particles can be accelerated to very high energy and form a power-law spectrum through scattering back and forth between the upstream and downstream regions of the shock [8, 9]. However, the acceleration of thermal particles is inefficient in this mechanism. In other words, particles must be pre-accelerated to a certain level of energy before entering the DSA process. This pre-acceleration process has remained unclear until now, and is known as the 'injection problem' of shocks [10, 11].

A quasi-perpendicular shock has a shock normal angle (the angle between the shock normal and the upstream magnetic field) larger than 45°. When a quasi-perpendicular shock has a sufficiently high Mach number, it can reflect part of the upstream particles, and then they get accelerated by the motional electric field associated with the shock [12]. Depending on the process by which particles are reflected, two mechanisms called shock surfing acceleration (SSA) [13–15] and shock drift acceleration (SDA) [16–19] have been proposed to accelerate the reflected particles, and to solve the injection problem of shocks. The ions undergoing SSA are reflected by the shock potential, while ions undergoing SDA experience a drift along the shock surface because of the large gradient of the magnetic field.

Laboratory experiments provide a good platform to explore the injection problem. Magnetized collisionless shocks have been produced by laser facilities using a piston-driven method in recent years [20]. The key steps of shock formation and early-stage ion energization have been observed in such experiments [21, 22]. Related simulation works have also been presented to study the formation of these piston-driven shocks [23, 24]. However, the exact pre-acceleration mechanism of ions in this kind of shock has never been discussed in detail. Here, we present a one-dimensional (1D) particle-in-cell simulation (PIC) of a supercritical perpendicular shock (Alfvén Mach number: $M_{\mathrm{A}}~6.85-9.35$ ), in which typical parameters related to experiments are used. This is a widely used method to study laser-produced shocks, the validity of which has been recognized by its consistency with the experimental results [21, 22]. We first investigate the formation process of the shock and then analyze the contributions of SDA and SSA to ion acceleration. Finally, we present a twice-reflected ion which is accelerated to relatively high energy.

2.   Simulation method and setup

A 1D electromagnetic PIC simulation code is used to simulate the piston-driven perpendicular shock. The code is modified from an original two-dimensional PIC simulation code, which has been successfully employed to study collisionless magnetic reconnection, plasma waves, and shocks [25–29]. In the simulation, every particle has one spatial and three velocity components (1D3V), whose motion is calculated by Newton's equation using a numerical method called Boris pusher [30]:

A uniform plasma condition is assumed in both the y and z directions ($\frac{\partial }{\partial y}=\frac{\partial }{\partial z}=0$ ), which is reasonable for shock experiments as the variation of physical values in the y–z plane can be ignored compared with the variation in the x direction. Then, the electromagnetic field is integrated over time by the Maxwell equations:

where ${\epsilon }_0$ is the vacuum dielectric constant, ${\mu }_0$ is the permeability of a vacuum, ${\mathbf{e}}_y,$${\mathbf{e}}_z$ are the unit vectors in y and z directions, and the current density J in each grid is calculated by:

Δx is the grid length, $\sum\limits_p$ means a sum over all the particles, q_p is the charge of each particle, and S is a three-point interpolation function based on the relative position of the particle () to the grid () [30].

The size of the simulation domain is $L_x=10.56\mathrm{mm},$ with each cell $2.64\mu \mathrm{m}$ in length. The ions consist of a 1:1 mixture of ${\mathrm{C}}^{5+}$ and ${\mathrm{H}}^{+}.$ At the beginning of the simulation, the entire domain is divided into three parts: the ambient plasma is placed in the region of $2.67\mathrm{mm}<x<10.56\mathrm{mm},$ which is embedded in a magnetic field of $B_0=6\mathrm{T}$ in $+y$ direction, with the electron number density $n_{\mathrm{e}0}=3\times {10}^{18}{\mathrm{cm}}^{-3}$ and temperature $T_{\mathrm{C}}=T_{\mathrm{H}}=T_{\mathrm{e}}=30\mathrm{eV}.$ The piston plasma is placed in the region of $0\mathrm{mm}<x<2.64\mathrm{mm},$ moving toward the ambient plasma at a bulk speed of $453\mathrm{km}{\mathrm{s}}^{-1},$ with the electron number density $n_{\mathrm{e}1}=3\times {10}^{19}{\mathrm{cm}}^{-3}$ and temperature $T_{\mathrm{C}}=T_{\mathrm{H}}=T_{\mathrm{e}}=800\mathrm{eV}.$ The plasma density, temperature, bulk velocity, and magnetic field are set to decrease linearly in the region of $2.64\mathrm{mm}<x<2.67\mathrm{mm}$ to avoid numerical problems. A reduced light speed of $3.8\times {10}^6\mathrm{m}{\mathrm{s}}^{-1}$ is used, and the mass ratio between the three species is $m_{\mathrm{C}}:m_{\mathrm{H}}:m_{\mathrm{e}}=12:1:0.04,$ which is already enough to make electrons and ions distinct in the shock. Boundary conditions are open for fields and reflecting for particles at both x boundaries, while piston particles are continuously injected from the left boundary. The simulation is pushed at $3.42\times {10}^{-5}\mathrm{ns}$ time intervals, and lasts for $11.97\mathrm{ns}.$ Since the ions are highly collisionless in laser-produced shocks [20–22], the Coulomb collision module is not employed in the simulation.

Based on the number density of the ambient plasma, the proton inertial length is calculated as ${d_{{\rm{p}}0}} = \sqrt {{m_{\rm{H}}}/{\mu _0}{n_{{\rm{e}}0}}{e^2}} \approx 0.132\:{\rm{mm}}$. The simulation domain is thus about long, with each cell about in length. The upstream Alfvén speed is ${V_{{\rm{A}}0}} = {B_0}/\sqrt {{\mu _0}({m_{\rm{H}}}{n_{{\rm{H0}}}} + {m_{\rm{C}}}{n_{{\rm{C0}}}})} \approx 51.29\:{\rm{km}}\:{{\rm{s}}^{ - 1}}$. Hence, the piston moves with a speed of about in direction at the beginning, and the light speed is set as The gyroperiod of protons in the ambient plasma is which means the simulation lasts for about The Debye length in the upstream region is which is larger than the grid length, and the Courant condition is satisfied by These parameters are similar to those applied in shocks produced by the Shengguang Ⅱ laser facility [31].

3.   Results

3.1   Shock formation

The evolution of the magnetic field $B_y$ is shown in figure 1. The piston plasma in the left part of the simulation domain moves in the $+x$ direction with a super-Alfvénic velocity and causes a magnetic field compression ahead of it. This compression is marked out by the black dashed line, indicating that the piston front moves with an approximately uniform speed of $V_p\approx 7.55V_{\mathrm{A}0}.$ This is lower than the speed of the piston plasma given at the beginning due to the existence of the magnetic field. As the piston plasma moves forward, another peak of the magnetic field starts to be formed upstream at ${\mathrm{\Omega }}_{\mathrm{cp}0}t\approx 1.4.$ This new structure is considered as the shock precursor, which clearly separates from the piston and becomes a well-formed shock at ${\mathrm{\Omega }}_{\mathrm{cp}0}t\approx 2.6.$ The position of the maximum magnetic field of the shock (the shock overshoot) is marked out by the dotted line in the figure, indicating that the shock velocity varies from about $6.85V_{\mathrm{A}0}$ to $9.35V_{\mathrm{A}0}.$

Figure  1.  Temporal evolution of the magnetic field $B_y.$ The black dashed line and the orange dotted line mark out the positions of the piston front and the shock overshoot, respectively. The whole simulation domain is thus divided into three parts: the upstream region, the downstream region, and the piston.

DownLoad: Full-Size Img PowerPoint

Figure 2 shows the shock formation process from the perspective of the phase space distribution of the protons. The first and second rows are the ambient protons in the and phase space, respectively. The third row is the piston protons in the phase space, and the fourth row is the profile of the magnetic field and electron density. At the beginning (), there is an electric field in the direction at the piston front because of the pressure gradient. A small part of the ambient protons is accelerated by this electric field and then gyrates in the upstream region (population Ⅰ). Part of the piston protons also penetrate into the ambient plasma (population Ⅳ), which behave very similarly to population Ⅰ later on. These piston-accelerated protons increase the plasma density ahead of the piston and cause a secondary magnetic field compression there, which starts to reflect upstream protons (population Ⅱ) when its amplitude is sufficiently large (). The existence of the shock-reflected protons (population Ⅱ) is a very important signature in experiments to confirm whether a supercritical shock is formed, but care must be taken first to distinguish them from population Ⅰ. A significant difference between populations Ⅰ and Ⅱ can be observed from the phase space: for those protons that are reflected by the piston (population Ⅰ), their velocities in the direction never become negative during the reflection, as they are reflected directly by the electric field in the direction. However, for those protons that are reflected later (population Ⅱ), they are accelerated in the direction by the motional electric field first, and then gyrate in the direction. When the protons of population Ⅱ start to gyrate back downstream (), two isolated peaks appear for both the magnetic field and electron density, which is considered as the beginning of the separation between piston and shock as shown in figure 1. A shock reformation process can be observed after the separation: at the piston-accelerated protons (population Ⅰ) enter the downstream region completely and are separated from the shock. A ring distribution (vortex) in phase space is formed at the shock ramp by the previous shock-reflected protons (population Ⅱ), while a new bunch of protons start to get reflected at the edge of this distribution (population Ⅲ). A new shock front is formed by population Ⅲ later (), which is considered as a well-formed shock completely free from the piston effects in figure 1. The shock evolves quasi-periodically for about every after that. The situations for ions are quite similar to those for protons and therefore are not shown. A substantial part of the upstream ions is reflected when the shock is formed, which may result in electrostatic ion–ion instability and Weibel instability in the shock foot region. This process has been found in previous experiments, but happens on a much longer timescale than our simulation ($t \sim 2000\omega _{{\rm{pi}}}^{ - 1} \sim 40{\rm{\Omega }}_{{\rm{cp0}}}^{ - 1}$) [20]. Therefore, the ion dynamics are hardly affected by the Weibel instability on the timescale of our simulation.

Figure  2.  Formation process at ${\mathrm{\Omega }}_{\mathrm{cp}0}t=$ 0.4, 0.9, 1.4, 2.0, 2.6. First row: ambient protons in $v_x-x$ space, where the velocity is relative to the upstream Alfvén speed $V_{\mathrm{A0}}.$ Second row: ambient protons in $v_z-x$ space. Third row: piston protons in $v_x-x$ space. Fourth row: magnetic field and electron density relative to their upstream values. Note that the colorbar in this figure represents the density distribution in phase space, which has no physical meaning.

DownLoad: Full-Size Img PowerPoint

The structure of a well-formed perpendicular shock is illustrated at ${\mathrm{\Omega }}_{\mathrm{cp}0}t=3.0$ (figure 3). Typical features of perpendicular shocks including overshoot, ramp, and foot are identified from the magnetic field. There is an electric field in the $+x$ direction in the foot and ramp region (which is known as the shock potential) due to the charge separation and a motional electric field in the $z$ direction produced by $E_z=-V_x{B}_y.$ The number densities of ${\mathrm{H}}^{+}$ and ${\mathrm{C}}^{5+}$ are almost identical in front of the piston and vary consistently with the strength of the magnetic field. Significant transverse heating of ions happens in the foot region due to their reflection, while electrons are heated in a relatively downstream region. Since this is a one-dimensional simulation, particles are almost not subjected to force in the $y$ direction; there is no parallel heating for both ions and electrons.

Figure  3.  Profile of the shock at ${\mathrm{\Omega }}_{\mathrm{cp}0}t=3.0.$ (a) $B_y$ relative to its upstream value $B_0,$ $E_x$ and $E_z$ relative to $E_0=V_{\mathrm{A0}}{B}_0.$ (b) Number densities of ${\mathrm{H}}^{+}$ and ${\mathrm{C}}^{5+}$ relative to their upstream values. (c) Parallel and vertical components of electron and proton temperature.

DownLoad: Full-Size Img PowerPoint

The shock evolution in the simulation shows similar characteristics to our experimental results [31], including the shock Alfvén number ($M_{\mathrm{A}}=7–11$ in the experiment), the separation time between the shock and piston (${\mathrm{\Omega }}_{\mathrm{cp0}}t~1$ ), and the shock structure in figure 3 (obvious 'overshoot', 'ramp', 'foot', etc). Moreover, the energy spectrum of shock-accelerated ions in our simulation is also consistent with our experimental result, which shows a quasi-monoenergetic distribution at about twice the shock velocity [31]. The validity of our simulation is thus proven and we can discuss the details of ion dynamics in this kind of shock below.

3.2   Ion dynamics

We traced the ambient protons and divided them into three species according to their trajectories: directly transmitted protons, once-reflected protons, and twice-reflected protons, which are shown in figure 4. The same three species are also found for ${\mathrm{C}}^{5+}$ ions.

Figure  4.  Three typical trajectories of the ambient protons. Dynamics of the once-reflected and twice-reflected protons in this figure are described in detail in figures 5 and 8 below.

DownLoad: Full-Size Img PowerPoint

The reflection process is analyzed in detail by showing the velocity components in three directions and the acting force, which is divided into Lorentz force and electric force. The first interaction between the chosen once-reflected proton and the shock happens at ${\mathrm{\Omega }}_{\mathrm{cp}0}t\approx 3.48–4.02$ (shadow area Ⅰ in figure 5). This proton enters the foot region first, getting accelerated by the electric field caused by charge separation in the $+x$ direction and the motional electric field in the $-z$ direction. As the proton goes further into the shock, the electric field in the $x$ direction becomes negative. However, the proton is still subjected to a growing Lorentz force in the $+x$ direction and finally gets reflected by this force. This reflected proton gyrates back in the upstream region and interacts with the shock again at ${\mathrm{\Omega }}_{\mathrm{cp}0}t\approx 5.12–5.4$ (shadow area Ⅱ in figure 5). Since its velocity in the $z$ direction is positive this time, the Lorentz force reverses in the $x$ direction and it penetrates the shock directly.

Figure  5.  Dynamics of the once-reflected proton: (a) the trajectory of the proton, (b) the velocities in three directions, (c) and (d) the forces in the $x$ and $z$ directions. The forces are divided into electric force ($F_{\mathrm{e}x}$ and $F_{\mathrm{e}z}$ ) and Lorentz force ($F_{\mathrm{m}x}$ and $F_{\mathrm{m}z}$ ).

DownLoad: Full-Size Img PowerPoint

SDA and SSA can be distinguished by judging which force dominates the reflection [19]. The proton shown in figure 5 obviously experienced SDA, as the Lorentz force contributed more than the electric force in the $+x$ direction during its reflection. The SSA process in which electric force dominates the reflection is also found in our simulation (figure 6). The proton selected in figure 6 is reflected at ${\mathrm{\Omega }}_{\mathrm{cp}0}t\approx 3.05–3.65$ (shadow area in figure 6), where the electric force is generally larger than the Lorentz force in the $+x$ direction. It is of interest to compare the incidence of these two different acceleration mechanisms in laser-produced shocks. Figure 7 shows a 2D histogram of ambient protons reflected by the shock at ${\mathrm{\Omega }}_{\mathrm{cp}0}t=3.0–6.0$ in the $W_E-W_{V\times B}$ space, where $W_E$ and $W_{V\times B}$ are the work done by the electric force and Lorentz force in the $+x$ direction during the reflection of each proton, respectively. The orange line $W_{V\times B}=W_E$ divides the reflected ions into two parts: the SSA part where $W_{V\times B}<W_E$ and the SDA part where $W_{V\times B}>W_E.$ SDA takes place more often in our simulation, with the ratio of these two species being $N_{\mathrm{SDA}}/N_{\mathrm{SSA}}\approx 2.72.$ These two kinds of protons reach a similar speed after the reflection, as most of them are reflected once. The overall distribution can be fitted as $W_{V\times B}+W_E=85.1m_{\mathrm{p}}{V}_{\mathrm{A}0}^2,$ corresponding to a velocity $V_x\approx 13V_{\mathrm{A}0}$ after the reflection (about 1.6 times the shock speed). For the ${\mathrm{C}}^{5+}$ situation, the ratio is even larger: $N_{\mathrm{SDA}}/N_{\mathrm{SSA}}\approx 3.3.$ However, the velocity they get after the reflection is still similar to that of protons.

Figure  6.  Dynamics of an SSA proton. The same physical quantities in figure 5 are presented in four columns. The shadow area (${\mathrm{\Omega }}_{\mathrm{cp}0}t\approx 3.06–3.68$ ) marks out the reflection process, where electric force is clearly the dominant force in the $+x$ direction. This means that the proton is shock-surfing-accelerated.

DownLoad: Full-Size Img PowerPoint

Figure  7.  2D histogram of reflected upstream protons in $W_E-W_{V\times B}$ space. $W_E$ and $W_{V\times B}$ are the work done by the electric force and Lorentz force in the $+x$ direction during the reflection of each proton, respectively. For each proton, $W_E$ and $W_{V\times B}$ are calculated from the time it starts to be accelerated to the time it reaches its highest $v_x$ in the upstream region. In order to avoid the piston effects, only protons reflected at ${\mathrm{\Omega }}_{\mathrm{cp}0}t=3.0–6.0$ are presented in this figure.

DownLoad: Full-Size Img PowerPoint

Very few protons (< 0.1% ) can be reflected more than once in our simulation. A twice-reflected proton is traced in figure 8. When it first interacts with the shock at ${\mathrm{\Omega }}_{\mathrm{cp}0}t\approx 2.66–3.14$ (shadow area Ⅰ in figure 8), the proton penetrates the shock potential and finally gets reflected by the Lorentz force, similar to the proton in figure 5. However, as the reflected proton gyrates back from the upstream region and interacts with the shock for the second time at ${\mathrm{\Omega }}_{\mathrm{cp}0}t\approx 4.04–4.18$ (shadow area Ⅱ in figure 8), it is reflected again by the shock potential. The second reflection is obviously different from the previous one, as the Lorentz force remains negative in the $x$ direction throughout the process, and the electric force dominates the reflection. In other words, this proton experiences both SDA and SSA in sequence. The difference between the two reflections can be clearly seen from their trajectories in the $x‐z$ plane in the shock rest frame (figure 9). The displacement of the proton in the z direction is obtained by integrating its velocity over time. The proton gyrates around the shock overshoot for the first reflection, which is dominated by the Lorentz force (SDA). For the second time, there is no gyration and the proton is reflected at a relatively upstream position as it is reflected by the electric field caused by charge separation in the shock foot region (SSA). This twice-reflected proton is accelerated to very high energy and has a gyration diameter even larger than the distance between the piston and the shock when it enters the downstream region. However, this proton is still constrained in the downstream region by the electric field associated with the piston (shadow area Ⅲ in figure 8).

Figure  8.  Dynamics of the twice-reflected proton. The same physical quantities in figure 5 are presented in four columns.

DownLoad: Full-Size Img PowerPoint

Figure  9.  Trajectories of the twice-reflected proton in the $x‐z$ plane in the shock rest frame. The orange dashed line marks out the position of the shock overshoot, and the two reflections are also pointed out. The displacement of the proton in the $z$ direction is obtained by the integration of velocity.

DownLoad: Full-Size Img PowerPoint

A comparison is made in figure 10 between the two protons selected in figures 5 and 8. For the first reflection, they are accelerated while gyrating. They reach a similar velocity when they start to interact with the shock for the second time ($v_x\approx 7.85V_{\mathrm{A}0},v_z\approx 13.4V_{\mathrm{A0}}$ ). The twice-reflected ion then gets accelerated in the $+x$ direction abruptly by the shock potential and reaches a velocity about 2.5 times the shock speed when it enters the downstream region. It is quite obvious that the particle can reach higher energy by experiencing more reflections. Unfortunately, ions that are reflected more than twice are very difficult to find in the simulation, since those twice-reflected ions are already very rare (we only found about 80 of them in 300 000 tracked macroparticles in the upstream region). However, in experiments, the real particle number is much larger than that used in the simulation. It is possible to diagnose some energetic ions being reflected by the shocks many times in experiments.

Figure  10.  Comparison between two protons in the $v_z‐v_x$ space. The orange line represents the once-reflected proton selected in figure 5. The blue line represents the twice-reflected proton selected in figure 8.

DownLoad: Full-Size Img PowerPoint

4.   Conclusions and discussion

In conclusion, we have performed a 1D simulation of supercritical perpendicular shock using achievable parameters for laser experiments. A shock with Alfvén Mach number $M_{\mathrm{A}}~6.85-9.35$ is produced by a laser-driven piston. In the beginning, a small part of the ambient ions is accelerated by the electric field associated with the piston, causing the formation of a shock precursor ahead of the piston. This shock precursor finally separates from the piston and becomes a well-formed shock that can reflect part of the upstream ions. Shock-reflected ions are distinguished from piston-accelerated ions in the phase space, which can be applied in experiments to exclude piston effects from real shock physics.

Shock-reflected ions are accelerated by two different mechanisms. Our statistics show that about 73% of the reflected protons are SDA, while the others are SSA. Most of the protons are reflected only once, which can reach an average velocity of about 1.6 times the shock speed after the reflection. Very few protons are reflected twice, which can be accelerated to a velocity larger than twice the shock velocity. Our results predict that this high-energy population of ions may also be observed in future laser experiments.

Non-thermal proton spectra have been observed in previous experiments, which are believed to be produced by SSA alone [22]. However, by carefully analyzing the contribution of the electric force and the Lorentz force during the reflection of ions, we find that SDA can also take place in laser-produced shocks. It should be noted that the incidence of these two mechanisms depends largely on the shock parameters such as the Mach number, the plasma β and so on. Therefore, a scan over a range of laser-produced shock parameters is worth doing in the future.