PIC simulations of stable surface waves on a subcritical fast magnetosonic shock front

PIC codes approximate each plasma species i by resolving its velocity distribution function with a cloud of computational particles (CPs). Each CP j of species i is characterized by a charge Q_j and mass M_j with a value Q_j /M_j , which matches the charge-to-mass ratio Q_i /M_i of a particle of species i. Every CP has a position x_j and velocity v_j and, hence, a spatially localized current density ∝Q_j v_j . The summation of the current density contributions of all CPs yields the macroscopic current density J(x, t). The electric field E(x, t), magnetic field B(x, t) and current density J(x, t) are defined on a numerical grid and the latter updates E(x, t) and B(x, t) according to discretized forms of Ampère's law

$$\begin{eqnarray}&&{\rm{\nabla }}\times {\bf{B}}-\displaystyle \frac{1}{{c}^{2}}\displaystyle \frac{\partial {\bf{E}}}{\partial t}={\mu }_{0}{\bf{J}},\end{eqnarray} \tag{ 1 }$$

where μ₀, c are the vacuum permeability and speed of light, and Faraday's law

$$\begin{eqnarray}&&{\rm{\nabla }}\times {\bf{E}}+\displaystyle \frac{\partial {\bf{B}}}{\partial t}=0.\end{eqnarray} \tag{ 2 }$$

The EPOCH code fulfills ∇ · B = 0 and Gauss' law ∇ · E = ρ/₀ (ρ, ₀: charge density and vacuum electric permittivity) to round-off precision [41]. Once the electromagnetic fields have been updated, they are interpolated to the position of each CP and update its momentum according to the relativistic Lorentz equation. All particle velocity components and, hence, all components of E, B, and J are updated also in simulations that resolve fewer than 3 spatial dimensions. This computational cycle is repeated for as many time steps Δ_t as necessary to cover the time scale of interest. More details of the numerical scheme of the EPOCH code are given elsewhere [42].

All simulations cover the interval −L_x /2 ≤ x ≤ L_x /2. The 2D simulations also resolve an interval 0 ≤ y ≤ L_y . Values of L_x and L_y vary between simulations. Boundary conditions are periodic in all directions. We model fully ionized nitrogen at the correct ion-to-electron mass ratio because it is widely used in laser-plasma experiments, for example in [3]. Figure 1 sketches the initial number density distribution of the ions along x. In the ambient plasma, we set the electron number density to n_e0 and that of the ions to n_i0 = n_e0/7. The electron plasma frequency ${\omega }_{{pe}}={({n}_{e0}{e}^{2}/{\epsilon }_{0}{m}_{e})}^{1/2}$ (e, m_e , c: elementary charge, electron mass, and light speed) sets the electron skin depth λ_e = c/ω_pe , with which we normalize space. At the time t = 0, the electric and magnetic fields are set to E = (0, 0, 0) and B = (0, B₀, 0) with eB₀/m_e ω_pe = 0.084. All ions have a temperature 200 eV. Electrons outside (inside) the dense cloud have a temperature 1000 eV (1500eV). Thermal diffusion lets more electrons stream from the dense into the ambient plasma than in the opposite direction, which yields an ambipolar electric field that points from the dense to the diluted plasma. This electric field lets the ions of the dense plasma expand into the ambient plasma. It also accelerates ambient electrons into the dense cloud. The accelerated ambient electrons form a beam that interacts with the electrons of the dense cloud. We give the latter a higher temperature in order to reduce the effects of a two-stream instability between both electron populations.

Zoom In Zoom Out Reset image size

The temperature T_e = 1000 eV of the electrons in the ambient plasma sets their thermal speed ${v}_{{te}}={({k}_{B}{T}_{e}/{m}_{e})}^{1/2}$ (k_B : Boltzmann constant). The ion temperature T_i = 200 eV and T_e set the ion acoustic speed c_s in the ambient plasma. On the time scales of ion-acoustic oscillations, electrons have three degrees of freedom and ions have one. This gives the adiabatic constants γ_e = 5/3 for electrons and γ_i = 3 for ions. The ion-acoustic speed becomes ${c}_{s}={({k}_{B}(7{\gamma }_{e}{T}_{e}+{\gamma }_{i}{T}_{i})/{m}_{i})}^{1/2}$. The Alfvén speed ${v}_{A}={B}_{0}/{\left({\mu }_{0}{n}_{i0}{m}_{i}\right)}^{1/2}$ and c_s define the FMS speed ${v}_{{fms}}={({c}_{s}^{2}+{v}_{A}^{2})}^{1/2}$. Our plasma has the electron plasma beta $\beta ={n}_{e0}{k}_{B}{T}_{e}/({B}_{0}^{2}/2{\mu }_{0})=0.55$. Relevant plasma parameters in the uniform ambient plasma and their values are listed in table 1. Although we selected parameters in our simulation setup, which are representative of some laser-plasma experiments where physical scales matter, we use normalized units scaled to the electron skin depth λ_e , the electron plasma frequency ω_pe , and the speed of light c.

Table 1. The relevant normalized frequencies, speeds, and spatial scales for the physical values n_e0 = 10²¹m⁻³ and B₀ = 0.85 T we used in our simulation setup and for ${t}_{{sim}}=14500$ (in units of ${\omega }_{{pe}}^{-1}$).

Electron gyrofrequency ωce = eB0/me ωpe	0.084	${\omega }_{{ce}}{t}_{{sim}}=1200$
Ion gyrofrequency ωci = 7eB0/mi ωpe	2.3 × 10−5	${\omega }_{{ci}}{t}_{{sim}}=0.33$
Ion plasma frequency ${\omega }_{{pi}}={(49{e}^{2}{n}_{i0}/{\epsilon }_{0}{m}_{i})}^{1/2}/{\omega }_{{pe}}$	0.0165	${\omega }_{{pi}}{t}_{{sim}}=240$
Lower-hybrid frequency ${\omega }_{{lh}}={\left({\left({\omega }_{{ce}}{\omega }_{{ci}}\right)}^{-1}+{\omega }_{{pi}}^{-1}\right)}^{-1/2}/{\omega }_{{pe}}$	0.0014	${\omega }_{{lh}}{t}_{{sim}}=20$
Electron thermal speed ${v}_{{te}}={({k}_{B}{T}_{e}/{m}_{e})}^{1/2}/c$:	0.044	${v}_{{te}}{t}_{{sim}}/{\lambda }_{e}=640$
Ion-acoustic speed ${c}_{s}={({k}_{B}(7{\gamma }_{e}{T}_{e}+{\gamma }_{i}{T}_{i})/{m}_{i}{c}^{2})}^{1/2}$	9.4 × 10−4	${c}_{s}{t}_{{sim}}/{\lambda }_{e}=13.7$
Alfvén speed ${v}_{A}={B}_{0}/{\left({\mu }_{0}{c}^{2}{n}_{i0}{m}_{i}\right)}^{1/2}$	1.4 × 10−3	${v}_{A}{t}_{{sim}}/{\lambda }_{e}=20.0$
FMS speed ${v}_{{fms}}={({c}_{s}^{2}+{v}_{A}^{2})}^{1/2}/c$:	1.7 × 10−3	${v}_{{fms}}{t}_{{sim}}/{\lambda }_{e}=24.3$
Electron thermal gyroradius vte /ωce λe	0.53
Electron Debye length vte /ωpe λe	0.044

We consider here FMS shocks that propagate through the ambient plasma. Wave dispersive properties of FMS shocks are set by the plasma conditions in the downstream plasma. However, the density and magnetic field amplitude do not vary much between the upstream (ambient) plasma and downstream plasma of a subcritical FMS shock. We discuss wave properties in the ambient plasma and assume that they are also representative of the downstream plasma.

The FMS mode is not dispersive for low wavenumbers k. Since shocks form by a steepening of the wavefront, the wavenumbers of the waves that form the shock increase in time. Eventually, they reach values where the FMS mode becomes dispersive. Simplified expressions for the dispersion relation of FMS waves in collisionless plasma can be obtained by either neglecting space charge, which gives ω_EM (k) = v_fms k, or by taking the electrostatic limit where the electric field is tied to oscillations of the charge density. The latter gives the lower-hybrid mode ${\omega }_{{ES}}(k)={\left(3{v}_{{ti}}^{2}{k}^{2}+{\omega }_{{pi}}^{2}({\omega }_{{ce}}^{2}+{v}_{{te}}^{2}{k}^{2})/({\omega }_{{pe}}^{2}+{\omega }_{{ce}}^{2}+{v}_{{te}}^{2}{k}^{2})\right)}^{1/2}$ with the ion thermal speed ${v}_{{ti}}={({k}_{B}{T}_{i}/{m}_{i})}^{1/2}$.

We use the noise distribution of the PIC code to track the full wave branch and show how it goes over into limits ω_EM (k) and ω_ES (k). Charge- and current-density fluctuations due to the moving CPs yield electric and magnetic field fluctuations. These fluctuations will have values of k and ω, which are connected to the particle motion. If CPs create fluctuations with such values of k, ω, they can often also absorb them. Hence, strong fluctuations of the electric and magnetic fields tend to reveal locations in k, ω-space where waves resonate with particles. They are strongest close to eigenmodes of the plasma (See [43] for a related discussion). The FMS mode compresses the background magnetic field and we can use fluctuations of B_y to track FMS modes in k, ω-space. In our 1D geometry, fluctuations in E_x are always tied to electrostatic charge density fluctuations and we use E_x to identify the lower-hybrid mode in k, ω-space.

Figure 2 shows the power spectra of the fluctuations in the magnetic B_y and electric E_x components. We computed them by running a 1D PIC simulation that resolved x with the plasma parameters of the ambient plasma discussed above. We Fourier-transformed B_y (x, t) and E_x (x, t) over space and time and multiplied the result with its complex conjugate giving ∣B_y (k, ω)∣² and ∣E_x (k, ω)∣².

Zoom In Zoom Out Reset image size

The noise at low wavenumbers k is magnetic and follows ω_EM . At large k, the electric field noise maps out lower-hybrid waves. The waves that connect ω_EM with ω_ES at k ≈ 1 have significant electric and magnetic components. They are thus not represented correctly by either limit.

Shocks form in our simulations just before the expanding ions reach the perturbation layer. We change the number density of mobile ions in this interval in figure 1 and keep that of mobile electrons unchanged. Since we set E = 0 at t = 0, Gauss' law ∇ · E = ρ/₀ implies initially a zero total plasma charge density ρ = 0 everywhere, which places an immobile positive charge density in our perturbation layer. Its electric field cancels out the one caused by the jump in the number density of mobile ions.

The 1D simulations and the 2D simulations 1 and 2 resolve L_x =160 with 8000 grid cells and end at the time ${t}_{{sim}}=14500$. The 2D simulation 3 uses 9000 grid cells to resolve L_x,3 = 180 and follows the shock for the time ${t}_{{sim},3}=1.25{t}_{{sim}}$. The 2D simulations 1-3 resolve L_y,1 = 12 by 600 cells, L_y,2 = 24 by 1200 cells, and L_y,3 = 36 by 1800 cells. The data is averaged over 4 cells in the 1D simulations and over patches of 2 × 2 cells in the 2D simulations. In the 2D (1D) simulations, ions and electrons are resolved by 25 CPs (375 CPs) per cell each.

We test with two 1D simulations how shocks react to changes in the number density of mobile ions in the perturbation layer. One 1D simulation represents 0.7n_i0 by mobile ions and the other 1D simulation represents 0.4n_i0 by mobile ions. According to figure 1, one shock will be launched at the density jump to the right of the dense cloud and a second one at the density jump to its left. Both will propagate away from the dense cloud and into the ambient plasma. Only the shock that moves to the right propagates through the perturbation layer. The left-moving shock in one of the simulations is taken as the reference shock. We invert the sign of the position x and velocity v_x in our plots so that we can compare the ion phase space density distribution of the left-moving unperturbed shock with that of the right-moving perturbed shock.

Figure 3 shows ion phase space density distributions. The supplementary movie 1 animates their evolution in time until ${t}_{{sim}}=14500$. At t₁ = 2000, a localized oscillation of the ion phase space density distribution is growing at x ≈ 8. We may interpret it as a steepening wave. It breaks shortly after this time and changes into a shock [19]. The magnetic field of the shock traps the electrons ahead of it and pushes them upstream. The current of the moving electrons and their space charge induces an electric field. It accelerates the practically unmagnetized ions with ${\omega }_{{ci}}{t}_{{sim}}\ll 2\pi$ ahead of the shock until their current balances the electronic one, giving rise to the shock foot. The change in the ion number density near the boundary x = 8.9 of the perturbation layer in figures 3(b), (c) is compensated by a faster motion of the ions in the perturbed interval.

Zoom In Zoom Out Reset image size

At t₂ = 5000 or ω_lh t₂ = 2.2π, qualitative differences can be observed between the three shocks. Most downstream ions in the interval 12 ≤ x ≤ 16 in figure 3(d) move at a speed below v_fms while the mean speed of those in figure 3(e) is about v_fms . The downstream ions in figure 3(f) are confined to a smaller interval 13 ≤ x ≤ 15 and most ions move faster than v_fms . Despite its faster downstream ions, the shock at x = 15 in figure 3(f) is trailing those in figures 3(d), (e). Fewer ions, which cross the shock, reduce the thermal pressure behind it and, hence, the shock speed in the downstream frame. Figure 3(f) also shows that the slowdown of the perturbed shock decreased the speed of the reflected ions. The ions at x ≈ 20, which were reflected by the shock when it just entered the perturbation layer, have a speed ≈3v_fms . Those at x ≈ 17, which were reflected after the slowdown of the shock, reach a speed of only 2.2v_fms .

In figure 3(d), the mean velocity of the ions is approximately constant for 9.5 ≤ x ≤ 16. The ion density change at x ≈ 11 is caused by a tangential discontinuity, which balances the high thermal pressure of the blast shell plasma against the combined thermal and magnetic pressure of the shocked ambient plasma. As the shocks in figures 3(e), (f) move into the perturbed interval, the pressure ahead of the tangential discontinuity decreases. Blast shell ions accelerate for 10 ≤ x ≤ 11 in figures 3(e), (f) in the ambipolar electric field of the blast shell's density gradient and decelerate at larger x, transferring their momentum to the downstream plasma.

Figures 3(d), (e), (f) show a structure in the shock-reflected ion beam at x ≈ 21 and v_x /v_fms ≈ 2.8. It develops when the shock-reflected ion beam, which has not yet developed at the snapshots for t₁ = 2000, catches up with the ions at the front of the rarefaction wave that is visible in figures 3(a), (b), (c) at high speeds for x > 8 (See supplementary movie 1).

Figure 4 shows the ion phase space density distributions at later times. We also plot the distributions of the ion densities that correspond to those of the phase space density. At t₃ = 7000 (ω_lh t₃ = 3.1π), the shocks are about to leave the perturbation layer. As before, the mean velocity of the downstream ions increases, and the shock speed in their rest frame decreases with a decreasing number density of the mobile ions ahead of the shock. A substantial lag is observed in particular for the shock in figure 4(c). The ion densities for this time are plotted in figure 4(d). In the overshoot of the unperturbed shock, the ion density reaches the value 3. It decreases to just below 2 downstream of the shock. The shocks, which move through the perturbation layer, have a lower density of their downstream plasma, and the ion density at their overshoot reaches about 3 times that of the mobile ions in the perturbation layer. At ${t}_{{sim}}=14500$ (${\omega }_{{lh}}{t}_{{sim}}=6.4\pi$), the shock speed and the mean velocity of the downstream ions behind the shocks are similar. The ion phase space vortices [44] in figures 4(f), (g) have been separated from the shock by the influx of shocked upstream plasma. The density distributions near the overshoots in figure 4(h) are similar for all shocks apart from a displacement along x.

Zoom In Zoom Out Reset image size

Figure 5 quantifies the impact of the changed density of mobile ions on the shock speed and position. We identify the position, where the shock is located, as the one with the largest curvature of the ion density and track it over time. We averaged the ion density over several grid cells to decrease statistical fluctuations of the ion density, which reduces the accuracy with which we can determine the spatial position of the shock's overshoot. The method, with which we determined the shock position, does also not always find the exact position of the shock in particular in the perturbation layer with the density of mobile ions 0.4n_io , which created several outliers in the data. Inaccuracies will be visible in particular in the velocity data that is obtained by differentiating the noisy position data. Nevertheless, the computed curves reveal trends that are confirmed by the supplementary movie 1. Figure 5(a) plots the separation of the reference shock, which moves through the unperturbed plasma, from the two shocks that move through the perturbed plasma. We start plotting the separation after t = 2500 when the shock fully reformed in figures 3(a)–(c). The shock moving through the perturbation layer with the mobile ion density 0.7n_i0 falls steadily behind the unperturbed shock until t ≈ 8000 when figure 4(b) shows that it crosses the upper boundary x = 20.8 of the perturbed interval. The shock, which moves through the plasma with the mobile ion density 0.4n_i0, leaves the perturbation layer last. Hence it is slowed down until t ≈ 10⁴. Figure 5(b) obtains the shock speed from the change of its position over a time interval ≈640. The speeds of the shocks remain above 1.3v_fms for all times. Given that the FMS wave with frequencies just below ω_lh is dispersive and that its phase speed is well below v_fms at wave numbers k ≈ 1 in figure 2, the Mach numbers are higher.

Zoom In Zoom Out Reset image size

In what follows, the term 'simulation j' with 1 ≤ j ≤ 3 refers to the 2D simulation with the length L_y,j along y. We express the number density of mobile ions n_i (x, y) in units of n_i0 and the in-plane magnetic field is $B{(x,y)=({B}_{x}^{2}+{B}_{y}^{2})}^{1/2}/{B}_{0}$. The magnetic B_z component remains at noise levels. The perturbation layer covers again 8.9 ≤ x ≤ 20.8. Within the perturbation layer of simulation j, the density of mobile ions varies as ${n}_{i}(x,y)=0.7\,+\,0.3\sin (2\pi y/{L}_{y,j})$ for 0 ≤ y ≤ L_y,j . We take the shock that moves in the direction of decreasing x < 0 as the unperturbed reference shock.

Figure 6 shows n_i (x, y) and B(x, y) at the time t = 7000, when the perturbed shocks in the three simulations have left the perturbation layer and entered the spatially uniform ambient plasma. Figures 6(a), (c), (e) show in each simulation j a localized peak of the ion density at x ≈ 22. All shocks lag behind the reference shock. The shape of the density peak varies between the simulations. With respect to the upstream region, the shock front is concave near y = 10 in simulation 3 and follows the isocontour of the magnetic amplitude. It is almost planar at y ≈ 5 in simulation 2 and convex at y ≈ 3 in simulation 1. We also observe a structure with a reduced ion density, which is a signature of an ion phase space vortex like the ones shown in figure 4, at x ≈ 16 and a non-planar front of the dense blast shell for x ≈ 14.

Zoom In Zoom Out Reset image size

Figures 6(b), (d), (f) show that the distributions of B(x, y) react to the density modulations of the shock fronts. The varying density of mobile ions in the perturbation layer affected the balance between the thermal- and magnetic pressures downstream of the shock and the ram pressure of the upstream ions. The magnetic field expanded in the upstream direction in intervals with a low density of mobile ions. We observe bent magnetic field lines at the shock and around the ion phase space vortex. The curves ${x}_{c,j}(y)=22.5-{A}_{j}\sin (2\pi y/{L}_{y,j})$ with A_j /L_y,j = 0.021 follow the modulation of the front of B(x, y). The supplementary movie 2 (simulation 1), movie 3 (simulation 2), and movie 4 (simulation 3) show the evolution of the shocks in the perturbation layer. It demonstrates that the density and magnetic field perturbations grow, saturate, and remain stationary until they leave it. This differs from what we observed in figure 5(a) where the distance between the reference shock and the perturbed shock steadily increased until the shock left the perturbation layer.

The density oscillation along y and around x ≈ 22 in figures 6(a), (c), (e) induces an electric field that points from the dense to the dilute plasma and accelerates ions. It acts as a standing surface wave with a wavenumber k_y,j = 2π/L_y,j . Electrons can move freely along the magnetic field and the perturbation could thus oscillate in the ion-acoustic mode with the frequency ω_cs,j = c_s k_y,j . If it does, it will oscillate with ω_cs,1 ≈ 0.36ω_lh , ω_cs,2 = 0.24ω_lh , and ω_cs,3 = 0.12ω_lh in simulations 1 to 3. The density perturbation along y is also coupled to the shock with its average normal along x. The upstream plasma that crosses the shock is compressed by a FMS mode with a frequency close to ω_lh . The large difference between the ion-acoustic frequency and ω_lh allows us to distinguish between both.

For every data time step, we determine the maximum value of the ion density along x as a function of y, which gives us a density distribution ${n}_{i,\max }(y,t)$. We also determine for each y the location where we reach B(x, y) = 1.67 first coming from the upstream region, which gives us the evolution in time of a curve x_B (y, t) that is similar to the solid curve in figures 6(b), (d), (f). The imaginary part of the fundamental frequency of the Fourier transform over space of x_B (y, t) and ${n}_{i,{\max }}(y,t)$ gives us the amplitudes x_B (k_y,j , t) and ${n}_{i,\max }({k}_{y,j},t)$. Figure 7 plots x_B (k_y,j , t)/L_y,j and ${n}_{i,{\max }}({k}_{y,j},t)/{L}_{y,j}$. In simulation 1 x_B (k_y,j , t) and ${n}_{i,{\max }}({k}_{y,j},t)$ are strongly damped. Both curves reach extrema at around t = 1.1 × 10⁴ and converge to a steady state after that. Weaker damping in simulation 2 yields pronounced extrema at t ≈ 10⁴ and weaker ones at t ≈ 1.3 × 10⁴. Both curves in simulation 3 oscillate in antiphase with the period ≈6000. The oscillation frequency is lowest in simulation 1 and it is similar in simulations 2 and 3, which is not what we expect from ion-acoustic oscillations with ω_cs,1 = 2ω_cs,2 and ω_cs,1 = 3ω_cs,3. Simulation 3 gives us the oscillation frequency 0.75ω_lh or 6.25ω_cs,3.

Zoom In Zoom Out Reset image size

Lower-hybrid waves have been analyzed in uniformly magnetized plasma [45, 46]. They trap electrons magnetically and the trapped electrons confine the ions electrically. Lower-hybrid waves require that the trapped electrons do not move far along the magnetic field during $2\pi {\omega }_{{lh}}^{-1}$. Their wavevector, which is aligned with the density gradient, therefore has to be almost perpendicular to the magnetic field. The perturbation in our 2D simulation rotates the magnetic field direction and the density gradient differently and both will not remain mutually orthogonal. This rotation sets an upper limit on the amplitude of the spatial oscillations of the magnetic field and, hence, on the perturbation amplitude in our 2D simulations.

The angle θ between the wavevector of an undamped lower-hybrid mode and the magnetic field is limited to the range $| 90^\circ -\theta | \lesssim {\theta }_{\max }$. In a plasma with equal temperatures of electrons and ions, the angular range can be estimated with ${\cos }^{2}{\theta }_{\max }={m}_{e}/{m}_{i}$, which yields ${\theta }_{\max }\approx 0.4^\circ$ for nitrogen ions.

The perturbed shock propagates into an upstream plasma, which is spatially homogenous and does therefore not favor a specific propagation direction of the boundary oscillation. Hence, the surface wave and the mode that compresses the upstream plasma constitute a standing wave, which is composed of modes that propagate obliquely to the shock normal and in opposite directions. We can estimate the propagation angle θ_j of these modes relative to the background magnetic field as follows. According to the ion velocity distribution in the interval 37 ≤ x ≤ 41 in figure 4(e), the wavelength of the lower-hybrid mode, which sustains the shock and forms the component of the wave that points along its normal, is k_x,s ≈ π. The wavelength of the surface wave that moves along the shock boundary is k_y,j = 2π/L_y,j . The propagation angles θ_j = k_x,s /k_y,j are thus θ₁ = 80.5°, θ₂ = 85.2°, and θ₃ = 86.8°.

Shock modes do not have to be undamped as they are continuously fed with energy from the inflowing upstream plasma. The wave vector of the shock mode may however be rotated into a direction with less damping. Another aspect is that lower-hybrid waves with frequencies ω ≈ ω_lh and with θ = 90° are dispersive [18]. If their phase velocity also changes with θ, we obtain an undamped surface wave only if the spread in wavenumbers along k_x and k_y is small. The higher resolution of k_y by simulation 3 allows the surface wave to involve wave modes with similar frequencies. Given that the shock fronts involve more than one wave mode, different spectral resolutions may also explain why their shapes differ in figures 6(a), (c), (e).

In addition to the requirement that lower-hybrid waves can only mediate a shock if they propagate almost perpendicularly to the background magnetic field, magnetic tension (B · ∇)B/μ₀ is likely to contribute to the saturation of the perturbation. In simulation 3 and at the time depicted in figure 6(f), its magnitude is comparable to that of the magnetic pressure gradient force density ∇B²/2μ₀ at x ≈ 22 and y ≈ 9 and about 20% of the thermal pressure gradient force density at the front of the shock in figure 6(e) at the same position (not shown). As we will see, the magnetic tension builds up over a much wider x-interval than the other forces and its impact on the saturation of the shock deformation in the perturbation layer is thus difficult to quantify.

We examine the distributions of the ion density and magnetic field amplitude at two times. Figure 8 shows them for simulations 1 and 2 at the time t = 10⁴, when the curves for simulation 2 went through their extrema in figure 7. The perturbation of the shock front in simulation 1, which is shown in figures 8(a), (b), has damped out. The front is planar and the density of the overshoot and the magnetic amplitude does not vary with y. The front is located one electron skin depth behind that of the reference shock and P_bx , which quantifies the deformation of the magnetic field in the simulation plane, is at noise levels. The perturbation of the shock front in simulation 2 has the opposite phase than the initial one in figures 6(c), (d); it is oscillating. In simulation 2, P_bx evidences a deformation of the magnetic field, which results in magnetic tension, in an interval with the width 3 centered on the shock.

Zoom In Zoom Out Reset image size

Figure 9 shows that at $t={t}_{{sim}}$, the shock perturbation has also damped out in simulation 2. The density distribution of simulation 1 in figure 9(a) shows another density maximum at x ≈ 40. Its separation from the leading density maximum at x ≈ 42 reveals that the wavelength of the lower-hybrid mode, which mediates the shock, is 2 and k_x,s = π as we observed in the 1D simulations. The magnetic field near the shocks in both simulations is aligned with y at this time. The lag between the perturbed and unperturbed shocks in figures 9(c), (f) is unchanged compared to that in figures 8(c), (f).

Zoom In Zoom Out Reset image size

Figure 10 shows the ion density and magnetic field distributions at the times t = 1.1 × 10⁴ and 1.66 × 10⁴, when the curves n_i (k_y,3, t) and x_B (k_y,3, t) go through extrema in figure 7.

Zoom In Zoom Out Reset image size

The oscillations of the shock fronts are in phase at the times t = 1.1 × 10⁴ and t = 1.66 × 10⁴ in simulation 3. Figure 10(a) reveals a density maximum at x ≈ 33 and y > 18, which is concave with respect to the upstream region. A weak density enhancement near y ≈ 8 bulges out into the upstream region. The position of the shock front is thus also a function of y and the variation is sinusoidal. The density oscillation is still present in figure 10(d). The front is almost flat, which indicates that some of the waves that sustained the shock have subsided. The in-plane magnetic field B(x, y) in figures 10(b), (e) is deformed near the shock front. It expands upstream in intervals, where the shock density is low. The supplementary movie 4 confirms that the shock front performs sinusoidal oscillations around its mean position along x, which are synchronized with those in the ion density along the shock front. The ion density oscillations are correlated with those of B(x, y).

The curves in figures 10(c), (f) show that the average ion density downstream of the shock is about 2; this shock does not compress the plasma much because the shock speed in the downstream frame is not small compared to the relative speed between the downstream and upstream plasma. According to P_bx , deformations of B(x, y) are strongest near the shock front. They extend upstream and downstream and decrease exponentially with distance from the maximum. Note that the curve has its maximum behind the front because the strong magnetic field does not fill out the simulation box for all values of y at larger x. The peaks in the ion density have a much smaller spread along x than their magnetic counterparts; charge density oscillations are shielded on Debye length scales and magnetic ones on skin depth scales. Figure 11 shows the evolution in time of the y-averaged density and field distributions in simulation 3. We have transformed these distributions into a reference frame that moves with the speed v_s = 1.6v_fms in the direction of the perturbed shock. Figure 11(a) shows that the shock is located at $\tilde{x}=x-{v}_{s}t\approx 3.5$. The position of its front oscillates in time. Its evolution is determined by a simultaneous oscillation of the shock front position, which varies with y and t, and of the density along the front. A strong in-plane electric field marks the front of the shock in figure 11(b) where the ion density gradient is largest. According to figure 11(c), the magnetic field is amplified by the shock crossing to more than twice its upstream value. The contribution of B_x to the magnetic pressure shown in figure 11(d) shows oscillations, which extend far into the upstream and downstream regions. The first and strongest maximum at early times is associated with the forced deformation of the shock as it propagated across the perturbation layer. After the shock left the perturbation layer, the shock performed free oscillations. They are coherent along the normal direction of the shock. Typical amplitudes of B_x are weak compared to B_y and we can thus interpret them as perturbations of an almost uniform guiding magnetic field that is aligned with y.

Zoom In Zoom Out Reset image size

The supplementary movie 5 animates the ion phase space density distribution in time for the interval x > 0 and figure 12 shows its final frame. The shock reforms at the time t ≈ 1500 and it remains stable after this time. It remains stationary in its comoving frame, which is typical for subcritical shocks. Some ions behind the shock are accelerated to high energies forming phase space vortices with an axis that is approximately aligned with the y-axis. One such vortex survives at x ≈ 36 in figure 12. A phase space vortex is sustained by the ambipolar electric field, which is associated with localized ion density depletion. The shock is located at x ≈ 53 at this time and it accelerates a small fraction of the upstream ions to a speed ≈3v_fms . The density of the shock-reflected ion beam is affected by the shock's density oscillation as can be seen from several patches with a reduced density. Since the upstream plasma outside of the perturbation layer has a uniform density, oscillations in the density of the shock-reflected ion beam must lead to a change in the number of ions that cross the shock and enter its downstream region.

Zoom In Zoom Out Reset image size