# PIC simulation of a nonoscillatory perturbation on a subcritical fast magnetosonic shock wave

M. E. Dieckmann, C. Huete, F. Cobos, A. Bret, D. Folini, B. Eliasson, R.

Walder

► To cite this version:

M. E. Dieckmann, C. Huete, F. Cobos, A. Bret, D. Folini, et al., PIC simulation of a nonoscillatory perturbation on a subcritical fast magnetosonic shock wave. Physica Scripta, 2024, 99, 10.1088/1402-4896/ad87c2. insu-04838838

HAL Id: insu-04838838

https://insu.hal.science/insu-04838838v1

Submitted on 15 Dec 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

<figure><img src="image_3.png" /></figure>

PAPER • OPEN ACCESS

# PIC simulation of a nonoscillatory perturbation on a subcritical fast magnetosonic shock wave

To cite this article: M E Dieckmann *et al* 2024 *Phys. Scr.* **99** 115606

View the [article online](#) for updates and enhancements.

## You may also like

- Application of fractional modified taylor wavelets in the dynamical analysis of fractional electrical circuits under generalized caputo fractional derivative
Ashish Rayal, Monika Anand and V K Srivastava

- Image encryption algorithm based on matrix projective combination-combination synchronization of an 11-dimensional time delayed hyperchaotic system
Jyotsna Kumari Bharti, P Balasubramaniam and K Murugesan

- Laser modulation to plateau region and intensity of high-order harmonic generation in monolayer $MoTe_2$
Zhiqiang Ji, Chenglong Wu, Xiaoyu Liu et al.
<figure><img src="image_4.png" /></figure>

CrossMark

OPEN ACCESS

RECEIVED

24 May 2024

REVISED

3 October 2024

ACCEPTED FOR PUBLICATION

16 October 2024

PUBLISHED

28 October 2024

Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence.

Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

<figure><img src="image_5.png" /></figure>

PAPER

# PIC simulation of a nonoscillatory perturbation on a subcritical fast magnetosonic shock wave

M E Dieckmann¹, C Huete², F Cobos³, A Bret³,⁴, D Folini⁵, B Eliasson⁶ and R Walder⁵

¹ Dept. of Science and Technology (ITN), Linköping University, Campus Norrköping, SE-60174 Norrköping, Sweden

² University Carlos III Madrid, Grp Mecan Fluidos, Leganes 28911, Spain

³ ETSI Industriales, Universidad de Castilla-La Mancha, 13071, Ciudad Real, Spain

⁴ Instituto de Investigaciones Energéticas y Aplicaciones Industriales, Campus Universitario de Ciudad Real, 13071, Ciudad Real, Spain

⁵ University of Lyon, ENS de Lyon, Univ Lyon 1, CNRS, Centre de Recherche Astrophysique de Lyon UMR5574, F-69230, Saint-Genis-Laval, France

⁶ Physics Department, University of Strathclyde, SUPA, Glasgow G4 0NG, Scotland, United Kingdom

E-mail: mark.e.dieckmann@liu.se

Keywords: collisionless shock, magnetized plasma, PIC simulation, shock deformation

Supplementary material for this article is available online

## Abstract

We use a two-dimensional particle-in-cell (PIC) simulation to study the propagation of subcritical fast magnetosonic shocks in electron-nitrogen plasma and their stability against an initial deformation. A slab of dense plasma launches two planar blast waves into a surrounding ambient plasma, which is permeated by a magnetic field that points out of the simulation box and is spatially uniform at the start of the simulation. One shock propagates into a spatially uniform ambient plasma. This reference shock has a Mach number of 1.75, and the heating of ions only along the shock normal compresses the ions that cross the shock to twice the upstream density. Drift instabilities lead to rapidly growing electron-cyclotron harmonic waves ahead of the location where the shock's density overshoot peaks, and to slowly growing lower-hybrid waves with a longer wavelength behind it. The second shock wave enters a perturbation layer that deforms it into a sine shape. Once the shock leaves the perturbation layer, the deformation is weakly damped and non-oscillatory, and the shock remains stable. Even without an external perturbation, and for the plasma parameters considered here, drift instabilities will cause ripples in the shock wave. These instabilities lead to a spatially and temporally varying compression of the plasma that crosses the shock.

## 1. Introduction

Shocks form if plasma collides at a speed that exceeds that of the ion density mode. In the absence of binary collisions between particles, the shock is mediated by the electromagnetic fields that are induced due to the collective motion of plasma particles [1]. Such a plasma, along with its associated structures and waves, is referred to as collisionless. Collisionless shocks are widely examined in the laboratory [2–5]. The Earth’s bow shock [6] is an example of a collisionless shock. It slows down the Solar wind plasma to a subsonic speed before it enters the shock’s downstream region, known as the magnetosheath. The shock’s structure depends on how the direction of the magnetic field, which is transported with the solar wind, is oriented relative to the normal of the shock plane. Here, we consider a perpendicular orientation.

Collisionless perpendicular shocks [7] are characterized by a compact transition layer between the freely flowing upstream plasma and the downstream plasma, which has been heated and compressed by the shock crossing. They are mediated by waves on the branch that become fast magnetosonic waves when the wavelength is long and the frequency is low. At higher frequencies and shorter wavelengths, these waves become partially electrostatic, compressive, and dispersive lower-hybrid waves [8, 9]. Lower-hybrid waves can sustain the density jump across the shock, provided that the shock speed is not too high. Subcritical shocks are mediated by the
electric cross-shock potential associated with the density jump between the upstream and downstream plasmas. Their supercritical counterparts require ion rotation in the magnetic field and a dense shock-reflected ion beam that dissipates some of the directed flow energy of the upstream plasma. Typically, the Earth's bow shock is supercritical [7] though it can turn subcritical if the solar wind is slow [10].

Observations by the magnetospheric multiscale mission (MMS) have revealed oscillations in the magnetic overshoot region of the bow shock that propagated along the magnetic field [11]. They have also been found in hybrid simulations that approximated the ions by a kinetic model and the electrons by a massless fluid [12]. Such oscillations are not limited to collisionless plasma shocks but have also been derived from models of shocks in (magneto)hydrodynamic fluids [13–19]. For an ideal gas, the equivalent adiabatic index above which self-sustained oscillations are possible is $\gamma > 1 + \sqrt{2}$ [20]. Shock oscillations require a restoring force, which can counteract the deforming one. In the hybrid simulations, the perturbations were caused by density clumps that resulted from interactions between the shock-reflected ion beam and the upstream plasma. As these clumps were convected into the perpendicular shock, they deformed its magnetic field that was oriented in the two-dimensional simulation plane. Magnetic tension provided the restoring force. A Fourier analysis of the shock's magnetic field direction revealed oscillations involving a wide range of wavenumbers. Their dispersion relation was that of Alfvén waves.

Motivated by the possibilities offered by laboratory astrophysics, we studied in [21] using particle-in-cell (PIC) simulations shock oscillations in plasma with conditions that are representative of some laser-plasma experiments [3]. Since PIC codes also approximate electrons by a kinetic model, they resolve additional wave modes and instabilities not captured by hybrid simulations, albeit at a higher computational cost. Our first series of simulations resolved the x-y plane. The magnetic field pointed along y while the subcritical fast magnetosonic shock expanded along x. We deformed the shock by letting it propagate across a perturbation layer with a limited extent along x, where the density of mobile ions varied sinusoidally with y. Once the shock left the perturbation layer and entered the spatially uniform upstream plasma, it performed damped oscillations around its equilibrium distribution. The frequency of the oscillations was in the range where the fast magnetosonic mode becomes quasi-electrostatic. The damping rate of the oscillations decreased with increasing box size along y. Here we use the same initial conditions as in the largest simulation performed in [21], but we rotate the magnetic field into the z-direction. This rotation removes magnetic tension as a means to stabilize the shock. Keeping the other plasma conditions unchanged, we can compare the shock we discuss here with that in our previous work to better understand the effects caused by the magnetic field.

We obtained the following results. Two fast magnetosonic shocks were launched by a slab of dense plasma and propagated in opposite directions into ambient plasma with identical conditions. One shock propagated into a spatially uniform plasma and served as the reference shock. It formed much faster than the inverse ion gyrofrequency and propagated at a speed of about 1.75 times the fast magnetosonic speed in the upstream frame of reference. The shock kept this subcritical speed until the simulation's end. The thin beam of shock-reflected ions did not drive strong waves ahead of the shock [22], and the shock transition layer remained narrow.

The magnetic ramp ahead of the shock was wider than an electron thermal gyroradius. Ambient electrons could be trapped magnetically by the ramp, while ions, having a much larger gyroradius, were unaffected by it. Ahead of the shock, the gradient of its density overshoot and the transport of trapped electrons across the ambient plasma resulted in an electric field that pointed normal to the shock into the upstream direction. Electrons trapped by the magnetic ramp underwent an ExB guiding center drift [23] along the shock front. Their electric current drove electrostatic waves with a wave vector that was aligned with their drift direction. Their rapid growth and short wavelength implied that they were electron-cyclotron harmonic waves. In time, electrostatic waves emerged in the density ramp between the shock's density overshoot and the downstream plasma. Their large wavelength suggests that they were lower-hybrid waves driven by electrons drifting in the ambipolar electric field of the density ramp. The electric field of the waves deflected the upstream ions that crossed the shock, giving rise to a spatially nonuniform ion density distribution downstream.

The second shock was deformed by its passage through a perturbation layer. After leaving this layer, the amplitude of the sinusoidal deformation was similar to that reported in [21] and small compared to its wavelength. The perturbation was nonoscillatory and weakly damped, and the shock was thus stable. The deformation of the shock front partially suppressed the drift waves behind it and led to an accumulation of ions near the interval of the shock boundary that lagged behind most. The absence of shock boundary oscillations in our simulation is evidence for the involvement of magnetic tension in the oscillations observed in [21].

Our manuscript is structured as follows. Section 2 discusses the numerical scheme of the PIC code, our initial conditions, and some aspects related to the shock. Section 3 presents the simulation results. Section 4 discusses our results and their potential relevance for shock observations by the MMS mission.
## 2. Background

### 2.1. Numerical scheme and initial conditions of the simulation

Particle-in-cell codes approximate each plasma species *i* by computational particles (CPs) with the charge $q_i$ and mass $m_i$ and a value $q_j/m_j$, which equals the charge-to-mass ratio $q_i/m_i$ of the represented species. The ensemble of all CPs that represent species *i* gives the phase space density distribution $f_i(\mathbf{x}, \mathbf{v})$. Each CP has a position $\mathbf{x}_j$ and velocity $\mathbf{v}_j$, from which we compute its contribution $\propto q_j\mathbf{v}_j$ to the macroscopic current density $\mathbf{J}$. The electric field $\mathbf{E}$ is updated with $\mathbf{J}$ via Ampère's law

$$ \nabla \times \mathbf{B} = \mu_0 \mathbf{J} + \mu_0 \epsilon_0 \frac{\partial \mathbf{E}}{\partial t}, \qquad (1) $$

($\epsilon_0$, $\mu_0$: vacuum permittivity and permeability) and subsequently the magnetic field $\mathbf{B}$ through Faraday's law

$$ \nabla \times \mathbf{E} = -\frac{\partial \mathbf{B}}{\partial t}. \qquad (2) $$

The equations for $\mathbf{E}$, $\mathbf{B}$, and $\mathbf{J}$ are solved on a numerical grid. Once $\mathbf{E}$ and $\mathbf{B}$ have been updated, they are interpolated to the position $\mathbf{x}_j$ of each CP and its velocity is advanced in time using the relativistic Lorentz force equation. The EPOCH code [24] we use is based on Esirkepov's scheme [25], which fulfills Gauss's law and the magnetic divergence law $\nabla \cdot \mathbf{B} = 0$ to round-off precision. In our simulation setup, CPs interact only via collective electromagnetic fields, and the absence of correlations between pairs of particles implies that velocity becomes an independent variable like the position. The plasma is thus collisionless. In our simulation, we initialize the particles of each species with a Maxwellian velocity distribution.

We study the evolution of shocks, which emerge if the velocity changes in plasma exceed the speed of the relevant density wave. The plasma must thus be set in motion and the magnetohydrodynamic momentum equation

$$ \rho \left( \frac{\partial}{\partial t} + \mathbf{v} \cdot \nabla \right) \mathbf{v} = \mathbf{J} \times \mathbf{B} - \nabla p = \frac{(\mathbf{B} \cdot \nabla) \mathbf{B}}{\mu_0} - \nabla \left( \frac{B^2}{2\mu_0} \right) - \nabla p \qquad (3) $$

illustrates how we can accomplish that. The mean velocity of the conducting fluid is $\mathbf{v}$, and $p$ is its thermal pressure. Let us consider the case where $\mathbf{v} = 0$ at the time $t = 0$. The fluid can be accelerated by any of the three terms on the right-hand side. The first two are the magnetic tension and the gradient of the magnetic pressure $P_B = B_0^2/2\mu_0$. In the case we consider here, the thermal pressure $p$ has contributions from the electrons and ions, but the latter is small due to the lower temperature and number density of the fully ionized nitrogen. Since the plasma is collisionless, the pressure is not mediated by binary collisions but through electric fields that arise from different mobilities and degrees of magnetization of electrons and ions.

Our ambient plasma consists of electrons with the number density $n_{e0} = 10^{21} \text{m}^{-3}$, temperature $T_{e0} = 1$ keV, and thermal pressure $P_{th} = n_{e0}k_B T_{e0}$ ($k_B$: Boltzmann constant). They have the electron plasma frequency $\omega_{pe} = (e^2n_{e0}/\epsilon_0 m_e)^{1/2} \approx 1.8 \times 10^{12} \text{s}^{-1}$, the thermal speed $v_{th,e} = (k_B T_{e0}/m_e)^{1/2} \approx 1.3 \times 10^7 \text{m/s}$ ($e, m_e$: elementary charge, and electron mass), and the Debye length $\lambda_D = v_{th,e}/\omega_{pe} \approx 7.4 \text{ µm}$. Nitrogen, which is used in experimental settings [3] because it is easier to handle than hydrogen, is the carrier of positive charge. It has the ionization state $Z = 7$, number density $n_{i0} = n_{e0}/Z$, the temperature $T_{i0} = T_{e0}/5$, and mass $m_i \approx 2.6 \times 10^4 m_e$. The ion plasma frequency is $\omega_{pi} = (Z^2e^2n_{i0}/\epsilon_0 m_i)^{1/2} \approx 3 \times 10^{10} \text{s}^{-1}$. A spatially uniform magnetic field $\mathbf{B}_0 = (0, 0, B_0)$ with $B_0 = 0.85 \text{ T}$ gives the electron gyrofrequency $\omega_{ce} = eB_0/m_e$ with $\omega_{ce}/\omega_{pe} = 0.084$ and the plasma $\beta = P_{th}/P_B = 0.56$. We normalize space by the electron skin depth $\lambda_e = c/\omega_{pe} \approx 170 \text{ µm}$ (c: speed of light) of the ambient plasma. We have only one positively charged particle species. Solar wind consists mostly of protons, and we may expect similar behavior in both plasmas. However, several characteristic frequencies and velocities, such as the ion's thermal velocity, depend on $m_i$. Our results are relevant for laser-plasma experiments but not necessarily for Solar wind plasma.

Our simulation box with periodic boundary conditions is oriented in the x-y plane as illustrated in figure 1.

We resolve the length $L_x = 180$ by 9000 grid cells and $L_y = 36$ by 1800 grid cells. The box is filled with spatially uniform ambient plasma everywhere except in two intervals. A dense plasma is placed in the center of the box along $x$ and is uniform along $y$. It has the same composition and ion temperature as the ambient plasma, 60 times its density, and the electron temperature $1.5 T_{e0}$. The direction of the magnetic field $\mathbf{B}_0$ defines the unresolved positive z-direction in the right-handed coordinate systems. The electron thermal pressure of $90 P_{th}$ in the dense plasma accelerates ions in both x-directions, as marked by the red and blue solid arrows. Blast waves formed by the accelerated ions drive shocks into the ambient plasma.

We define one right-handed coordinate system for each half-space and use the center of the dense plasma to define $x = 0$. The largest $x$ value is 90, and the initial boundary between the dense and ambient plasma is located at $x = 3$. The y-coordinate runs from 0 to $L_y = 36$. The shock in the blue half-space propagates through a spatially uniform ambient plasma and serves as the reference shock, while the shock in the red half-space
<figure><img src="image_6.png" /></figure>

propagates through a perturbation layer. The number density of mobile ions in the perturbation layer, $8.9 \le x \le 20.8$, is $n_{i,mob}(y) = 0.7 + 0.3 \sin(2\pi y/L_y)$. Since we set $E=0$ and keep the electron density constant in the perturbation layer at the time $t=0$, the net charge $Zn_{i,mob}(y) - n_{e0} \le 0$ is compensated by an immobile positive charge that serves as a grating for the shock.

We represent ions and electrons by 25 CPs per cell each and use the correct mass ratio between both species. The simulation evolves the plasma for a time of $t_{sim} = 1.2 \times 10^{-8}$ s or $\omega_{ci}t_{sim} = 0.48$, with the ion gyrofrequency $\omega_{ci} = ZeB_0/m_i$.

## 2.2. Relevant plasma processes and prior work

A rarefaction wave propagates into the dense plasma and accelerates the ions in the opposite direction, creating a blast wave. If the blast wave expands freely, the mean velocity of its ions increases, and their density decreases with an increasing distance from the rarefaction wave front [26]. The ambient plasma reacts to the blast wave by forming a charge density wave, which can steepen into a shock. The type of wave that grows depends on the plasma conditions. In unmagnetized plasma, and for low frequencies $\omega$, the dispersion relation of ion-acoustic waves $\omega = c_s k$ connects $\omega$ to their wavenumber $k$. For the initial conditions in the ambient plasma and the adiabatic constants $\gamma_e = 5/3$ and $\gamma_i = 3$, the ion-acoustic speed is $c_s = (k_B(Z\gamma_e T_{e0} + \gamma_i T_{i0})/m_i)^{1/2} \approx 2.8 \times 10^5$ m/s. The phase speed $\omega/k$ of ion-acoustic waves decreases as $\omega$ approaches $\omega_{pi}$. The presence of the perpendicular magnetic field $B_0$ modifies this dispersion relation, giving $\omega = v_{fms}k$ with the fast magnetosonic (FMS) speed $v_{fms} = (c_s^2 + v_A^2)^{1/2} \approx 5 \times 10^5$ m/s and the Alfvén speed $v_A = B_0/(\mu_0 n_{i0} m_i)^{1/2} \approx 4.1 \times 10^5$ m/s. The phase speed of FMS waves equals $v_{fms}$ for low frequencies and decreases as $\omega$ approaches the lower-hybrid frequency $\omega_{lh} = ((\omega_{ce}\omega_{ci})^{-1} + \omega_{pi}^{-2})^{-1/2}$ [8, 9].

Wave dispersion stabilizes the shock. Large amplitude waves steepen over time [27, 28], which shifts their wave spectrum to higher $k$. Once the wave spectrum of the steepening wave reaches a range of $k$ where the phase speed starts to decrease, short waves can no longer keep up with the shock, thereby halting further steepening (see [21, 29] for related case studies). Steepened ion-acoustic waves or FMS waves with frequencies close to $\omega_{lh}$ can sustain a shock, provided that the shock speed relative to the upstream plasma does not exceed $c_s$ or $v_{fms}$ by 2-3 times [30–34]. In our case, such a subcritical shock can evolve on time scales $\omega_{lh}^{-1} \approx \omega_{ci}^{-1}/60$. Supercritical FMS shocks evolve on a time scale $\omega_{ci}^{-1}$ [35–37].

One aspect of collisionless shocks, which sets them apart from hydrodynamic ones, is that they cannot fully thermalize the ions that traverse them. Consider a subcritical shock that is sustained by a jump in the electric potential across its boundary. The electric field, which slows down the inflowing upstream ions, points along the shock normal. If the upstream ions are warm, their velocity components along the shock normal have a thermal spread. How much their velocity changes during the shock crossing depends on how their kinetic energy in the rest frame of the shock compares to the shock’s electric potential jump. Fast ions lose a smaller fraction of their
kinetic energy than slow ones when they cross the shock, increasing the ions' thermal velocity spread. Thus, ions are heated by the shock along its normal. However, the electric field of the shock does not alter the ion velocity along the shock plane. Heating ions only in one direction reduces the compression ratio to 2 [38], and the thermal anisotropy of the ions drives collisionless instabilities downstream of the shock [39, 40]. Ions can also be heated stochastically [41].

Simulations of supercritical FMS shocks in [12, 37] demonstrated that the orientation of $B_0$ relative to the simulation plane affects their structure. Several factors are at play. Firstly, the magnetic tension in equation (3) is only important if magnetic field lines can be deformed, which requires them to be oriented in the simulation plane. Secondly, field lines perpendicular to the simulation plane can be pushed apart allowing plasma to enter the demagnetized region. Thirdly, the distribution of the thermally anisotropic ions in the shock transition layer, and the instabilities they drive, depend on the magnetic field direction relative to the resolved wavevectors.

Another aspect relevant to our PIC simulation is the evolution of the electric current $J$ near the shock front. It is caused by a drift of the electrons relative to the ions. We can understand this drift with the following idealized model. In the absence of an electric field, an electron with the speed $v \perp B$ moves along a closed circular trajectory perpendicular to a magnetic field $B$, which is constant in space and time. Let a unidirectional electric field $E$, which is constant in space and time and points orthogonally to $B$, act on this electron. In the rest frame of $B$, the electron is accelerated if $v \cdot E < 0$ and decelerated otherwise. The electron's gyroradius becomes large after acceleration and small after deceleration, causing a net drift of the guiding center of the electron trajectory perpendicularly to $E$ and $B$. The electron's trajectory perpendicular to $B$ is a cycloid. The speed $|E \times B|/B^2$ of the guiding center is known as the ExB drift speed, which is discussed in [23] along with other drift mechanisms. In principle, ions can also drift in the same direction. However, this is only true if the electromagnetic field patches are large compared to an ion's gyroradius, which is not the case in our simulations. The net drift of the electrons drives drift instabilities with wavevectors aligned with the drift velocity vector.

If the magnetic field is oriented in the simulation plane, the current is due to electron drift along the numerically unresolved direction perpendicular to the simulation plane. When the spatially varying magnetic field points out of the simulation plane, electrons stream relative to ions in the simulation plane at the speed $-v_D$ (see figure 1), and the unstable waves are resolved in the simulation. In the rest frame of the electrons, with an ion drift speed of $v_D$, the instability is driven by the interaction of the moving ions with the plasma eigenmodes. Several eigenmodes exist, each with characteristic frequencies. We discuss the two drift instabilities most relevant to our work, analyzed in an idealized one-dimensional geometry.

The electron-cyclotron drift instability leads to the growth of electron-cyclotron harmonic waves. Typically, the fastest-growing electron-cyclotron harmonic waves have frequencies close to the upper-hybrid frequency $\omega_{uh} = (\omega_{pe}^2 + \omega_{ce}^2)^{1/2}$. This instability [42–44] requires that the electron's guiding center drift speed $|E \times B|/B^2 > v_{th,e}$. It is quenched if the wavelength of the electron-cyclotron waves approaches the Debye length. This instability yields electrostatic waves with short wavelengths and a growth rate that is several percent of $\omega_{ce}$ for a plasma with $\omega_{pe} = 10\omega_{ce}$ and $T_{e0} > T_{i0}$. It has been observed ahead of planar supercritical shocks [45]. Although electron-cyclotron harmonic waves are caused by electron oscillations, their low phase speed in the ion rest frame modulates the ion density if their amplitude is large [22]. Ion density oscillations propagate in the direction of the drifting electrons.

The lower-hybrid drift instability [46–50] is also driven by a drift between electrons and ions perpendicular to the magnetic field and the waves it drives have been observed at the supercritical bow shock [51]. In the rest frame of the ions, the ion density waves propagate in the direction of the drifting electrons. Ion density waves can grow for drift speeds as low as $v_D > [k_B(T_e + T_i)/m_i]^{1/2} (T_e, T_i$: electron and ion temperatures) [47], which is below $10^{-2} v_{th,e}$ in our ambient plasma and thus much lower than the drift speed that triggers the electron-cyclotron drift instability. The work in [46] provides a criterion for instability if the electrons are cold. The frequency of each ion-cyclotron harmonic wave branch was approximated by $\omega_l(k) \approx l\omega_{ci}$ (l: integer). For a plasma density $n$ we get the condition for instability

$$ \frac{n}{|\nabla| n} \frac{\omega_{ci}}{v_{th,i}} \le (1/2l)(m_i/2m_e)^{1/2}. \quad (4) $$

The exponential growth rate of the wave is $\gamma_l/\omega_{ci} \approx l(m_e/m_i)^{1/4}$.

## 3. Results

We first discuss the evolution of the reference shock in the blue half-space in figure 1, followed by that of the perturbed shock in the red half-space. In what follows, we express the ion density $n_i$ in units of $n_{i0}$, the electric field $E_x$ in units of $cB_0$, and the magnetic field $B_z$ in units of $B_0$. The normalized spatial distribution of the ion density is $n_i(x, y)$, that of the electric field is $E_x(x, y)$, and that of the magnetic field is $B_z(x, y)$. Their averages over $y$ are denoted by $\langle n_i \rangle_y$, $\langle E_x \rangle_y$, and $\langle B_z \rangle_y$.
<figure><img src="image_7.png" /></figure>

### 3.1. The reference shock in the blue half-space

Figure 2 shows the relevant density and field distributions at the time $\omega_{ll}t = 10$. The ion density $n_i(x, y)$ in figure 2(a) increases from its upstream value to the peak value of the shock's density overshoot at $x \approx 25$. The density distribution near this structure is fragmented, and striation patterns extend to the lowest displayed values of $x$. Its average, $\langle n_i \rangle_y$, is plotted in figure 2(b). The density of the blast wave, $\langle n_i \rangle = 8$ at $x=8$, decreases to $\langle n_i \rangle = 2$ at $x = 18$. At larger values of $x$, the density remains close to 2 until $x = 23$, then it increases to 2.5 at $x \approx 25$ and subsequently decreases, eventually reaching 1. We identify the interval $18 \le x \le 23$ as the downstream region of the shock that separates it from the unperturbed upstream plasma. The density change across the shock is sustained by the electric field pulse in figure 2(c). The pulse is fragmented, and electric field oscillations extend upstream of the shock. Figure 2(d), which plots $\langle E_x \rangle_y$, confirms that $E_x(x, y)$ peaks where the density changes most and that it extends up to $x \approx 28$. A weaker negative electric field is visible in figure 2(d) for $24 \le x \le 25$, where the density grows from its downstream value to that at the shock's density overshoot. Figure 2(e) shows that the shock crossing amplifies the upstream magnetic field to twice its value. According to figure 2(f), this amplification starts at $x \approx 28$ and ends behind the shock's location. It decreases slowly in the downstream region $18 \le x \le 23$ of the shock, and $\langle B_z \rangle_y$ changes rapidly for $x \le 18$; the expanding blast wave pushes the magnetic field out.

The width of the magnetic field ramp $\approx 3$ in the interval $25 \le x \le 28$, also referred to as the shock foot, exceeds the electron's thermal gyroradius $v_{th,e}/(\omega_{ce}\lambda_e) \approx 0.5$. The foot's magnetic field, which moves in the rest frame of the ambient plasma, traps electrons upstream of the shock and pushes them to increasing $x$. Ions are unmagnetized and can thus not compensate for the electric current of the trapped electrons. This current induces the electric field ahead of the shock, which extends up to $x \approx 28$ in figure 2(d). The ambipolar electric field due to the density gradient at the shock's density overshoot adds to this field, creating the large spike near $x = 25$ in figure 2(d). The magnetically trapped electrons start drifting in the electric field. The condition $|\mathbf{E} \times \mathbf{B}|/B^2 \ge v_{th,e}$, which must be met for the growth of electron-cyclotron harmonic waves, reduces in our geometry to $E_x/B_z \ge v_{th,e}$. At $x = 26$, the values $B_z \approx 1.2B_0$ and $E_x \approx 0.05cB_0$ give $E_x/B_z \approx v_{th,e}$, which can lead to the growth of electron-cyclotron harmonic waves.

Figure 3 shows the density and field distributions at $\omega_{ll}t = 25$. The shock front has progressed to $x \approx 57$, which amounts to an average speed $1.75v_{fms}$.

The rapidly oscillating ion density striations ahead of the shock, which are compressed before entering the downstream plasma, remain present. The distribution of $n_i(x, y)$ now also reveals patches with a density $\approx 3$ that extend far downstream of the shock front, reaching widths of several electron skin depths $\lambda_e$. These large structures leave imprints on $E_x(x, y)$ and $B_z(x, y)$, as shown in figures 3(c), (e). The distributions of $\langle n_i \rangle_y$, $\langle E_x \rangle_y$, and $\langle B_z \rangle_y$ in figure 3 resemble those in figure 2. The shock transition layer, characterized by an increase in these quantities from upstream to downstream, retains a thickness of a few $\lambda_e$. Movie 1, which animates the distributions of $n_i(x, y)$, $|E_x(x, y)|$, and $B_z(x, y)$ for $10 \le \omega_{ll}t \le 25$, demonstrates the stability of the shock and the convection of the large density striations with it.
<figure><img src="image_8.png" /><figcaption>**Figure 3.** The reference shock at $\omega_{lh}t = 25$: Panel (a) shows the ion density $n_i(x, y)$ and (b) $\langle n_i \rangle_y$. The dashed red line marks the value 2. Panels (c, d) show $|E_x(x, y)|$ and $\langle E_x \rangle_y$. Panel (e) shows $B_x(x, y)$ and (f) $\langle B_x \rangle_y$.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 4.** The time-evolution of the y-averaged quantities that characterize the reference shock: Panel (a) shows $\langle n_i \rangle_y^{1/2}$. The color scale is clamped to the value of 5. The dashed red line indicates a speed of $1.75v_{fms}$. Panels (b, c) show $|\langle E_x \rangle_y|$ and $\langle B_z \rangle_y$.</figcaption></figure>

Figure 4 illustrates the time evolution of $\langle n_i \rangle_y$, $\langle E_x \rangle_y$, and $\langle B_z \rangle_y$. Figure 4(a) demonstrates how the thermal expansion of the dense plasma at $x \le 3$ reduces its peak density and drives the blast wave into the ambient plasma. The shock front forms over a time interval $\approx 2\pi\omega_{lh}^{-1}$ and expands at a constant speed, close to the previous estimate of $1.75v_{fms}$ based on the comparison between figures 2 and 3. The strongest pulse in $\langle E_x \rangle_y$ in figure 4(b) marks the position of the shock. Two less pronounced electric field pulses correspond to the rarefaction wave front, which propagates into the dense plasma, and the pulse visible for $\omega_{lh}t \le 7$, marking the location where the blast wave's density changed most rapidly. Figure 4(c) shows $\langle B_z \rangle_y$, illustrating that the expanding blast wave expels the magnetic field from the region to the left and piles it up ahead of it. The magnetic field amplitude doubles upon crossing the shock. By assuming the compressed magnetic field is stationary in the downstream plasma, the speed of the downstream plasma can be estimated by comparing the speeds of the left and right boundaries of the compressed magnetic field. At $\omega_{lh}t = 20$, the magnetic field with $\langle B_z \rangle_y \ge 1.2$ spans the interval $24 \le x \le 40.5$, and by $\omega_{lh}t = 25$, it has expanded to $32 \le x \le 58.5$. The shock front moves 2.25 times faster than the downstream plasma, suggesting that the downstream plasma speed is approximately $0.8v_{fms}$, while the shock speed is about $0.95v_{fms}$ faster.

Having determined the shock speed, we can now examine the growth and evolution of the striations in the shock's rest frame. Figure 5 shows the distributions of $n_i(x, y)$ near the left-most patch in figure 3 at times $\omega_{lh}t = 22, 25$, and 28. A moving window was applied, traveling with the shock speed of $1.75v_{fms}$ and keeping the overshoot region centered. Movie 2 depicts the evolution of $n_i(x, y)$ within this window for $10 \le \omega_{lh}t \le 25$.
<figure><img src="image_10.png" /><figcaption>**Figure 5.** The distributions of $n_i$ at the times $\omega_{lh}t = 22$ (a), 25 (b), and 28 (c). The moving window moves with the speed $1.75v_{fms}$ and keeps the shock overshoot centered. The diagonal black lines, which cover the intervals 3 along $y$ and 2.5 along $x$, are placed at the same locations in the moving window frame. The horizontal lines in (b) and (c) cover the intervals 1 and 2 along $y$, which shows that the structure propagated at the speed $\approx 3.7v_{th,i}$ along the shock front.</figcaption></figure>

In figure 5(a), the thin ion density striations ahead of the shock have wavelengths $\lambda \le 0.5$ along $y$, which is an order of magnitude larger than the Debye length $\lambda_D = 0.044$ of the ambient plasma. The ion density striations move leftward as they approach the shock front, as observed at $0 \le y \le 2$ and $x \approx 64$ in figure 5(c) and Movie 2. The brief time $\approx \lambda_e \omega_{ci} / v_{fms} = 0.014$ for ions to cross the shock's foot rules out ion rotation due to $B_0$ as a cause of the drift. Since $B_0$ points into the image plane, the striations move in the direction of electron drift. This, combined with their short wavelength, suggests that these are ion density waves driven by strong electron-cyclotron harmonic waves. Given that the growth rate of the electron-cyclotron drift instability is a few percent of $\omega_{ce}$ and $\omega_{ce} \approx 60\omega_{lh}$, this instability likely drives the ion density waves ahead of the shock in figure 5 and at $\omega_{lh}t = 10$ in figure 2. The ions in the thin striations are compressed and do not visibly expand after crossing the shock in figure 2(a), as $v_{th,i} / (\omega_{lh} \lambda_e) \approx 0.1$.

Figure 5 shows that the large density patch grows behind the shock overshoot, spanning an interval of width 3. It propagates rightward by a distance of 1 over an interval of 3, giving a speed of $3.7v_{th,i}$ relative to the ions' rest frame. This large ion density structure drifts in the opposite direction of the ion density striations ahead of the shock. It is confined to $54 \le x \le 57$ at $\omega_{lh}t = 25$, where $\langle n_i \rangle_y$ increases with $x$ in figure 3(b). The ambipolar electric field caused by this density change will point opposite to that ahead of the shock and be of lower magnitude. Even with the stronger magnetic field behind the shock, the ExB-drift speed will be slower than ahead of the shock, quenching the electron-cyclotron drift instability. Lower-hybrid drift waves, with a wavelength $\lambda \approx 2\pi$, are observed in figure 5.

We can estimate with the help of equation (4) and figure 3(b) whether these waves are driven by the lower-hybrid drift instability in the standard form, where the plasma is stationary to begin with. The ion density behind the overshoot is $n_i \approx 2n_{i0}$ and it changes by $0.5n_{i0}$ over $2\lambda_e$, giving $n_i / |\nabla n_i| \approx 8\lambda_e$. The ion's thermal gyroradius is $5.4\lambda_e$. The lower-hybrid frequency $\omega_{lh} \approx 60\omega_{ci}$ falls into the ion-cyclotron harmonic wave branch with $l=60$. Introducing these values into equation (4) gives $1.5 \le 1$. The inequality predicts that the lower-hybrid wave should not grow. However, the values on both sides are close, and equation (4) was derived for cold electrons and did not take into account the ion motion perpendicular to the electron drift direction near the shock. If we replace $v_{th,i}$ with the actual ion speed $\approx v_{fms}$ in the rest frame of the shock, equation (4) becomes $0.11 \le 0.95$. The growth rate estimate by [46] would give $\gamma / \omega_{ci} \approx 17$ or $\gamma \approx 0.3\omega_{lh}$, which would explain why this wave can grow and saturate during the simulation time.

## 3.2. Perturbed shock

Figure 6 presents the density and field distributions of the shock after it left the perturbation layer and entered the uniform ambient plasma. Its front has a sinusoidal deformation with an amplitude $\approx 1$. The value of $n_i(x, y)$ is largest at $x \approx 25$ and $y \approx 9$, where the density of mobile ions in the perturbation layer peaked. Figure 6(b) shows values $\langle n_i \rangle_y < 2$ in the perturbation layer because we averaged only over the mobile ions. The density behind the shock and outside the perturbation layer is, however, close to $\langle n_i \rangle_y = 2$, as in the case of the reference shock. Averaging smoothed out the overshoot that is visible in $n_i(x, y)$. Figure 6(c) shows a strong peak in $E_x(x, y)$ at the shock, which gives rise to the broad peak of $\langle E_x \rangle_y$. Two more peaks of $\langle E_x \rangle_y$ are located at $x \approx 21$ and 22.
<figure><img src="image_11.png" /><figcaption>**Figure 6.** The perturbed shock at $\omega_{\parallel}t = 10$: Panel (a) shows the ion density $n_i(x, y)$ and panel (b) $\langle n_i \rangle_y$. The dashed red line marks the value $\langle n_i \rangle_y = 2$. Panels (c, d) show $E_x(x, y)$ and $\langle E_x \rangle_y$. Panel (e) shows $B_z(x, y)$ and (f) $\langle B_z \rangle_y$.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 7.** Ion phase space densities $f_i(x, v_x)$ at the times $\omega_{\parallel}t = 5$ (left column) and $\omega_{\parallel}t = 10$ (right column). Panels (a, d) show the distributions at $y = 9$, panels (b, e) those at $y = 18$, and panels (c, f) those at $y = 27$. The vertical red lines mark the boundaries of the perturbation layer $8.9 \le x \le 20.8$. All distributions have been integrated over a spatial interval along $y$ with width 0.9 centered on the given $y$ value, normalized to the peak value in the ambient plasma at $t = 0$, and displayed on the same 10-logarithmic color scale.</figcaption></figure>

They are equal in magnitude and enclose a negatively charged potential, which suggests the presence of an ion phase space vortex. As in the case of the unperturbed shock, $B_z(x, y)$ is compressed to twice its value by the shock crossing.

Figure 7 examines the ion phase space density distribution at the times when the shock is about to enter the perturbation layer and at the time corresponding to figure 6. We selected the locations along $y$ where the density of the mobile ions in the perturbation layer has a maximum of 1, an average of 0.7, and a minimum of 0.4. Movie 3 animates these distributions for $0 \le \omega_{\parallel}t \le 10$.

Figures 7(a)–(c) correspond to $\omega_{\parallel}t = 5$, and the ion distributions are practically identical up to $x = 8.9$. Ions near the center of the dense plasma at $x = 0$ have a mean velocity of 0. The rarefaction wave, which propagates at a speed $\approx c_s$ into the dense plasma, has crossed a distance $\approx 4$ during $\omega_{\parallel}t \approx 5$, which is comparable to the thickness of the dense plasma. The density near $x = 0$ is still close to the initial value, but it starts to decrease after this time. Ions of the blast wave gain speed with increasing $x$, and the ambipolar electric field tied to the density change has expelled the ambient ions close to the initial boundary of the dense plasma. The ambient ions and the cooler, denser, and faster blast wave ions start to overlap at $x \approx 5$. The density contribution of the expelled
<figure><img src="image_13.png" /><figcaption>**Figure 8.** The function $f_x(y, t) = X(t) + A(t)\sin(2\pi y/L_y)$ that fits the shock best. Panels (a, b) plot $X(t)$ and $A(t)$. The vertical red lines mark the times $\omega_{\parallel}t = 12.5, 15.8, 19.2$, and 22.5. Panels (c-f) compare $|E_x(x, y)|$ and $f_x(y, t)$ (red curves) using the values from panels (a, b) at the marked times.</figcaption></figure>

ambient plasma reduces the density gradient and, thus, the ion acceleration by the ambipolar electric field in the interval $5 \le x \le 6$. At $x \approx 7$, we observe a hybrid structure composed of an electrostatic shock in the ambient plasma and a double layer for blast wave ions that reach the hybrid structure [52]. The crest of a steepening lower-hybrid wave is located at $x \approx 12$ and $v_x \le v_{fms}$. Ions between the hybrid structure and this crest are accelerated, and this acceleration increases with decreasing density of mobile ions, which we may understand as follows. The expanding blast wave ahead of the hybrid structure pushes the magnetically trapped electrons to larger $x$. Their density in the transition layer does not change with $y$. If the magnetic field in the foot of the hybrid structure moves at the same speed in figures 7(a)–(c), mobile ions must be accelerated to a higher speed in intervals with a lower mobile ion density to balance the electric current of the drifting electrons. The blast wave ions extend to $x \approx 18$ and speeds of $3.2v_{fms}$. Their phase space density and velocity spread are both an order of magnitude less than those of the ambient plasma. For $x \ge 15$, the blast wave density is only about $\sim 10^{-2}$, and the instabilities this ion beam can drive will be weak and grow slowly.

At the time $\omega_{\parallel}t = 10$, the steepening wave has changed into a shock that has left the perturbation layer. The shock in figure 7(d) has progressed farthest and is located at $x = 25$. At this time, the ion phase space distribution in the interval $23 \le x \le 25$ consists mostly of ambient ions that crossed the shock and were compressed by it; this structure is no longer a hybrid structure but an electrostatic shock. The mean speed of the ambient ions downstream of the shock is $0.8v_{fms}$, which matches the speed estimate obtained from figure 4. A population of blast wave ions with a low density is found at higher speeds; the blast wave ions and the shocked ambient ions have not yet been spatially separated. This is accomplished by a tangential discontinuity that separates the strongly magnetized shocked ambient plasma from the dense unmagnetized blast wave plasma [53]. The distribution of fast ions ahead of the shock in $25 \le x \le 45$ consists of a beam of blast wave ions with $x > 30$ that crossed the shock at early times and a thin beam of shock-reflected ambient ions to its left. Both beams can be distinguished by the level of granularity. Computational ions, which represented the dense plasma at $t = 0$, have a larger numerical weight than those in the ambient plasma. Figures 7(e), (f) show shocks that lag behind that in figure 7(d). Ions speed up in figures 7(e), (f) as they leave the perturbation layer at $x = 20.8$ and slow down at $x = 22$, which matches the positions and directions of the electric field spikes in figure 6(d).

Figure 8 tracks the evolution in time of the shock front's deformation. We determine the position along $x$ where $E_x$ reaches its maximum value as a function of $y$ and for all times, giving $s_x(y, t)$. We Fourier transform the function over space, giving $S_x(k_y, t)$, extract $S_x(0, t)$, and the imaginary component of $S_x(k_1, t)$ with $k_1 = 2\pi/L_y$. The first gives $X(t)$ and the second $A(t)$ of the function $f_x(y, t) = X(t) + A(t)\sin(k_1 y)$ that approximates the shock front best. Figure 8(a) shows $X(t)$. Its speed of $1.8v_{fms}$ is practically constant in time. Figure 8(b) shows that $A(t) = 0.85$ just after the shock left the perturbation layer, it reaches its maximum $0.97$ at $\omega_{\parallel}t = 12.3$ and decreases thereafter. Figures 8(c)–(f) compare $|E_x|$ with $f_x(y, t)$ at 4 selected times. The shape of the shock front changes in time but it does not oscillate. Figure 8(f) demonstrates a straightening of the shock front between its extrema, which involves spatial harmonics $mk_1$ (integer $m > 1$) of the sinusoidal deformation with period $k_1$.

Figure 9 shows the distributions of $n_i(x, y)$, $|E_x(x, y)|$, $B_z(x, y)$, and their y-averages at $\omega_{\parallel}t = 15.8$. Movie 4 shows the time evolution of $n_i(x, y)$, $|E_x(x, y)|$, and $B_z(x, y)$ for $10 \le \omega_{\parallel}t \le 25$. The density overshoot across the
<figure><img src="image_14.png" /><figcaption>**Figure 9.** The perturbed shock at $\omega_{\text{th}}t = 15.8$: Panel (a) shows the ion density $n_i(x, y)$ and panel (b) $\langle n_i \rangle_y$. The dashed red line marks the value $\langle n_i \rangle_y = 2$. Panels (c, d) show $|E_x(x, y)|$ and $\langle E_x \rangle_y$. Panel (e) shows $B_x(x, y)$ and (f) $\langle B_x \rangle_y$.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Figure 10.** The perturbed shock at $\omega_{\text{th}}t = 25$: Panel (a) shows the ion density $n_i(x, y)$ with the red box centered on $y = 27$ and covering $56 \le x \le 58$. Panel (b) plots $\langle n_i \rangle_y$. The dashed red line marks the value $\langle n_i \rangle_y = 2$. Panels (c, d) show $|E_x(x, y)|$ and $\langle E_x \rangle_y$. Panel (e) shows $B_x(x, y)$ and (f) $\langle B_x \rangle_y$.</figcaption></figure>

shock in figure 9(a) is located at $x \approx 37$. It is trailed by structures at $y \approx 9$ and $y \approx 27$, which are similar in size and density to those behind the unperturbed shock. Movie 4 reveals that the elongated structures with high ion density well behind the shock are caused by the compression of the inflowing upstream plasma by the non-stationary, non-planar shock front. The tangential discontinuity at $x \approx 21$, which separates the dense unmagnetized blast wave plasma from the shock-compressed magnetized ambient plasma, is not planar. It was probably deformed by drift instabilities. Figure 9(b) demonstrates that the density of the ambient plasma that crosses the shock is still doubled. Figure 9(c) shows a strong electric field peak at the shock and oblique waves ahead of it, as in the case of the unperturbed shock. Density gradients of the ion density structures downstream of the shock yield the electric field near $y = 9$ and $x = 36$. Figures 9(e), (f) show that the magnetic compression is the same as for the reference shock.

The structure of the shock remains qualitatively unchanged at the time $\omega_{\text{th}}t = 25$ shown in figure 10. We display here a much larger range along $x$ to visualize the change from the blast shell plasma to the shock-compressed ambient plasma. A high-density structure is located just behind the red mark in figure 10(a). It is the
<figure><img src="image_16.png" /><figcaption>**Figure 11.** Ion phase space densities $f_i(x, v_x)$ at the times $\omega_{li}t = 15.8$ (left column) and $\omega_{li}t = 25$ (right column). Panels (a, d) show the distributions at $y = 9$, panels (b, e) those at $y = 18$, and panels (c, f) those at $y = 27$. All distributions have been integrated over a spatial interval along $y$ with width 0.9 centered on the given $y$ value, normalized to the peak value in the ambient plasma at $t = 0$, and displayed on the same 10-logarithmic color scale. The vertical red line is the border of the perturbation layer $x = 20.8$.</figcaption></figure>

source of a thin stripe of dense plasma that extends far downstream of the shock. The tangential discontinuity, which is characterized by the front of the dense plasma, is now located at $x \approx 30$. Its propagation speed for the time $15.8 \le \omega_{li}t \le 25$ is thus about $0.8v_{fms}$. Small structures with high density are found downstream of the shock in figure 10(a). Movie 4 shows that they formed behind the shock and ahead of the discontinuity. The post-shock density $\langle n_i \rangle_y$ is now close to 1.9, suggesting a weakening of the shock. Its front is still marked by a narrow electric field pulse in figure 10(c), but with an amplitude below that at $\omega_{li}t = 15.8$. The rear end of the shock-compressed magnetic field near $x \ge 30$ in figure 10(e) does not follow the shape of the density structure, which separates the blast wave from the shocked ambient plasma. The nearby plasma is not in thermal equilibrium, which is also evident from the nonuniform density distribution.

Figure 11 displays the ion phase space density distributions at the times $\omega_{li}t = 15.8$ and 25 along the three slices $y = 9, 18,$ and 27. Movie 5 shows the time evolution of the distributions for $10 \le \omega_{li}t \le 15.8$, and Movie 6 animates the distributions for $15.8 \le \omega_{li}t \le 25$, thus between both columns in figure 11. The dense blast wave plasma with a mean speed $\approx v_{fms}$ has expanded from $x \approx 22$ in figure 11(a) to about $x \approx 31$ in figure 11(d), which is in line with what figures 9(a) and 10(a) showed. Structures resembling shocks are present to the right of the perturbation layer at low velocities in figures 11(b), (c). Movie 3 shows that they grow just to the right of $x = 20.8$ when the shock front enters the unperturbed ambient plasma and are convected with the flow to increasing $x$. At $\omega_{li}t = 15.8$, the blast wave ions overlap with the slower shocked ambient ions in the interval $25 \le x \le 33$. In this interval, the density of the blast wave ions decreases and that of the slower shocked ambient ions increases with $x$. For $x \ge 34$, we find almost exclusively ambient ions. The shock front in all three slices is trailed by an ion phase space vortex, which is seen best in figure 11(a) near $x = 36$, where the shocked ambient ions form a hemicircle. Ions in this vortex are accelerated by the electric field of the wave train that sustains the shock. This oscillatory wave train propagates in the downstream plasma with a speed just below $v_{fms}$. Its amplitude peaks at the shock's position and decreases with $x$, having a second maximum of the positive field near $x \approx 34$, evidenced by ions that gain speed at this location. At $\omega_{li}t = 25$, the ion phase space vortex has fully formed in the slice $y = 9$ near $x = 55$. For $x > 45$, the shocked ambient plasma is an order of magnitude denser than the blast wave plasma, and the latter no longer affects the shock dynamics. Figure 11(f) shows a single dense population of ions at $x \approx 54$ and $v_x \approx v_{fms}$ just behind the shock and in the interval behind the red box in figure 10(a), where the dense ion structure is located.

## 4. Discussion

We examined the expansion of a pair of subcritical fast magnetosonic shocks into an ambient plasma permeated by a spatially uniform magnetic field perpendicular to the simulation box. Our initial conditions produced two blast waves that propagated from a dense plasma slab into the surrounding ambient plasma. The ambipolar electric field associated with the blast wave's density gradient accelerated the ambient ions. Initially, hybrid structures formed, consisting of a combination of a shock in the ambient ions and a double layer for the blast
wave ions [52]. The ambient ions, reflected by the electrostatic shock, and the blast wave ions that crossed it pushed the magnetic field and magnetically trapped electrons ahead of the shock. The net electric current from the moving electrons induced an electric field, which drove a lower-hybrid wave in the ambient ions through a process similar to that described in [54]. Its saturation led to the formation of an electrostatic shock, well ahead of the initial hybrid structure and largely decoupled from the blast wave ions. One blast wave drove an electrostatic shock that propagated through a spatially uniform ambient plasma. This reference shock remained stable, with a Mach number of 1.75, until the simulation's end. A second electrostatic shock propagated through a perturbation layer, where the mobile ion density varied sinusoidally along the shock boundary, deforming it. We tracked the amplitude of the sinusoidal deformation and the shock speed over time after it left the perturbation layer. The shock propagated at the same speed as the reference shock. The boundary perturbation was non-oscillatory for the duration of the simulation, and the amplitude of the sinusoidal deformation decreased by about 20% by the simulation's end.

This differs from our previous case study [21], where the background magnetic field pointed within the simulation plane. In that case, the shock boundary's perturbation exhibited damped oscillations around its equilibrium, with a frequency just below the lower-hybrid frequency. This suggests that the oscillations were caused by magnetic tension due to deformation of the magnetic field with a wavevector along its direction. In the present case, where the magnetic field points out of the simulation box, the wavevector of lower-hybrid waves is always perpendicular to the magnetic field, making them linearly undamped. This was not the case in [21], where wave damping of lower-hybrid waves not strictly perpendicular to the magnetic field may have contributed to the different shock behavior.

Assuming magnetic tension is responsible for the shock oscillations, we can combine the results from this study and [21]. In a three-dimensional simulation that resolves the magnetic field direction and allows for magnetic tension, shock perturbations would give rise to waves propagating along but not orthogonal to the magnetic field, causing the shock to wrinkle in the latter direction.

We observed charge density waves ahead of and behind the shock boundary, which were not observed in [21], where the instability mechanism was geometrically suppressed. Ahead of the shock, the ExB drift of electrons resulted in the growth of electron-cyclotron harmonic waves through a mechanism similar to the electron-cyclotron drift instability [42]. This instability has a high growth rate if the drift speed is high and the electrons are hotter than the ions. A lower electron drift speed behind the shock suppressed this growth, and instead, lower-hybrid waves grew through an instability similar to the lower-hybrid drift instability. The criterion for this instability, derived under idealized assumptions in [46], which differ from those behind the shock in our simulation, was nearly fulfilled. Replacing the thermal speed of the ions in equation (4) with their mean speed in the shock frame suggests that this instability would grow. Electric field oscillations due to the lower-hybrid waves behind the shock front caused angular deflections of ions crossing the shock. Given that the ion density distribution downstream of the shock was less uniform than in [21], we conclude that these waves are responsible for the nonuniformities.

In both simulations, the ion density downstream of the shock was twice as high as upstream, and the lower-hybrid waves did not change the average shock compression. This low compression is caused by ion heating only along the shock normal [38]. The relatively low number of shock-reflected ions and the short simulation time, much shorter than the inverse ion gyrofrequency, implied that strong waves could not be driven ahead of the shock [45, 55].

Lower-hybrid waves behind the shock may contribute to the ripples observed by the MMS satellites at a subcritical fast magnetosonic shock [56]. However, verifying the relevance of the mechanism discussed here for the shock observed by the MMS mission is beyond the scope of our work.

Boundary oscillations that propagated along the magnetic field direction of Earth's bow shock were observed by the MMS mission [11]. This shock had a higher Mach number than the one studied here, and the oscillation frequency was lower, but the underlying mechanism might be the same. It is important to note that, apart from the supercritical speed of the shock discussed in [11], the shock dynamics in nitrogen plasma may differ from those in Solar wind plasma, and the spatial and temporal scales resolved by our simulation are orders of magnitude smaller than those needed to capture Alfvén waves.

Future work should examine subcritical shocks in proton plasma under conditions similar to those found in the Solar wind to determine when lower-hybrid waves occur behind the shock front. Larger simulation boxes could provide more detailed information about the spectrum of electron-cyclotron harmonic waves and lower-hybrid waves, which could then be compared to existing analytical work. Another interesting aspect would be to test the validity of the analytical estimate for the stability of the lower-hybrid drift instability, discussed in section 2, under Solar wind plasma conditions. Additionally, investigating shock boundary oscillations in laser-plasma experiments would be worthwhile. The wavelength at which shock oscillations become undamped is of the order of a few tens of electron skin depths, which may be within reach for laser-plasma experiments.
Acknowledgments

The simulations were performed on resources provided by the National Academic Infrastructure for Supercomputing in Sweden (NAISS) at the National Supercomputer Centre partially funded by the Swedish Research Council through grant agreement no. 2022-06725 and on the centers of the Grand Equipment National de Calcul Intensif (GENCI) under grant number A0090406960. The first author also acknowledges financial support from a visiting fellowship of the Centre de Recherche Astrophysique de Lyon. A.B. acknowledges support from the Ministerio de Economía y Competitividad of Spain (Grant No. PID2021-125550OBI00).

Data availability statement

All data that supports the findings of this study are included within the article (and any supplementary files). The
data that support the findings of this study are available upon reasonable request from the authors.

Conflict of interest

The authors declare that they have no conflict of interest.

ORCID iDs

M E Dieckmann https://orcid.org/0000-0003-4055-0552

C Huete https://orcid.org/0000-0002-3227-8520

F Cobos https://orcid.org/0000-0001-5953-4080

B Eliasson https://orcid.org/0000-0001-6039-1574

References

[1] Marcowith *A* et al 2016 The microphysics of collisionless shock waves *Rep. Prog. Phys.* 79 046901

[2] Romagnani L *et al* 2008 Observation of collisionless shocks in laser-plasma experiments *Phys. Rev. Lett.* 101 025004

[3] Ahmed H *et al* 2013 Time-resolved characterization of the formation of a collisionless shock *Phys. Rev. Lett.* 110 205001

[4] Schaeffer D B, Fox W, Habersberger D, Fiksel G, Bhattacharjee A, Barnak D H, Hu S X and Germaschewski K 2017 Generation and evolution of high-mach-number laser-driven magnetized collisionless shocks in the laboratory *Phys. Rev. Lett.* 119 025001

[5] Fazzini *A* et al 2022 Particle energization in colliding subcritical collisionless shocks investigated in the laboratory Astron A&A 665 A87

[6] Eastwood J P, Hietala H, Toth G, Phan T D and Fujimoto M 2015 What controls the structure and dynamics of earths magnetosphere? *Space Sci. Rev.* 188 251–86

[7] Burgess D, Mbius E and Scholer M 2012 Ion acceleration at the earths bow shock space *Sci. Rev.* 173 5–47

[8] Verdon A L, Cairns I H, Melrose D B and Robinson P A 2009 Properties of lower hybrid waves *Proc. Int. Astron. Union* 257 569–73

[9] Verdon A L, Cairns I H, Melrose D B and Robinson P A 2009 Warm electromagnetic lower hybrid wave dispersion relation *Phys. Plasmas* 16 052105

[10] Mellott M M 1984 The physical mechanisms of subcritical collisionless shock-wave formation *Adv. Space Res.* 4 245–53

[11] Johlander *A* et al 2016 Rippled quasiperpendicular shock observed by the magnetospheric multiscale spacecraft *Phys. Rev. Lett.* 117 165101

[12] Lowe R E and Burgess D 2003 The properties and causes of rippling in quasi-perpendicular collisionless shock front *Ann. Geophys.* 21 671

[13] Freeman N C 1955 A theory of the stability of plane shock waves *Proc. R. Soc. Lond. A. Math. Phys. Sci.* 228 341

[14] Zaidel P M 1960 Shock wave from a slightly curved piston *J. Appl. Math. Mech.* 24 316

[15] Dyakov S P 1954 Shock wave stability *Zh. Eksp. Teor. Fiz.* 27 288–304

[16] Kontorovich V M 1957 On the shock waves stability *Zh. Eksp. Teor. Fiz.* 33 1525-6

[17] Gardner C and Kruskal M 1964 Stability of plane magnetohydrodynamic shocks *Phys. Fluids* 7 700

[18] Wetta N, Pain J-C and Heuzé O 2018 Dyakov-Kontorovitch instability of shock waves in hot plasmas *Phys. Rev. E* 98 033205

[19] Huete C, Cobos-Campos F, Abdikamalov E and Bouquet S 2020 Acoustic stability of nonadiabatic high-energy-density shocks *Phys. Rev. Fluids* 5 113403

[20] Calvo-Rivera, Huete C, Garcia-Rubio F, Velikovich A L, Betti R and Tzeferacos P 2023 Stability of perpendicular magnetohydrodynamic shocks in materials with ideal and nonideal equations of state *Phys. Rev. E* 108 035203

[21] Dieckmann M E, Huete C, Cobos F, Bret A, Folini D, Eliasson B and Walder R 2023 PIC simulations of stable surface waves on a subcritical fast magnetosonic shock front *Phys. Scripta* 98 095603

[22] Dieckmann M E, McClements K G, Chapman S C, Dendy R O and Drury L O C 2000 Electron acceleration due to high frequency instabilities at supernova remnant shocks *Astron. Astrophys.* 356 377–88

[23] Cary J R and Brizard A J 2009 Hamiltonian Theory of guiding-center motion *Rev. Mod. Phys.* 81 693–738

[24] Arber T D *et al* 2015 Contemporary particle-in-cell approach to laser-plasma modelling *Plasma Phys. Control. Fusion* 57 113001

[25] Esirkepov T Z 2001 Exact charge conservation scheme for Particle-in-Cell simulation with an arbitrary form-factor *Comput. Phys. Commun.* 135 144–53

[26] Grismayer T and Mora P 2006 Influence of a finite initial ion density gradient on plasma expansion into a vacuum *Phys. Plasmas* 13 032103
[27] Dawson J M 1959 Nonlinear electron oscillations in a cold plasma *Phys. Rev. Lett.* **113** 383

[28] Shukla P K, Eliasson B, Marklund M and Bingham R 2004 Nonlinear model for magnetosonic shocklets in plasmas *Phys. Plasmas* **11** 2311–3

[29] Dieckmann M E, Folini D, Walder R, Romagnani L, d'Humieres E, Bret A, Karlsson T and Ynnerman A 2017 Emergence of MHD structures in a collisionless PIC simulation plasma *Phys. Plasmas* **24** 094502

[30] Marshall W 1955 The structure of magneto-hydrodynamic shock waves *Proc. R. Soc. A* **233** 367–76

[31] Bardotti G and Segre S E 1970 Laminar electrostatic shock waves in a plasma *Plasma Phys.* **12** 247–52

[32] Edmiston J P and Kennel C F 1984 A parametric survey of the 1st critical Mach number for a fast MHD shock *J. Plasma Phys.* **32** 429–41

[33] Dieckmann M E, Ahmed H, Sarri G, Doria D, Kourakis I, Romagnani L, Pohl M and Borghesi M 2013 Parametric study of non-relativistic electrostatic shocks and the structure of their transition layer *Phys. Plasmas* **20** 042111

[34] Ly M N, Sano T, Sakawa Y and Sentoku Y 2023 Conditions of structural transition for collisionless electrostatic shock *Phys. Rev. E* **108** 025208

[35] Winske D and Quest K B 1988 Magnetic field and density fluctuations at perpendicular supercritical collisionless shocks *J. Geophys. Res.* **93** 9681–93

[36] Chapman S C, Lee R E and Dendy R O 2005 Perpendicular shock reformation and ion acceleration *Space Sci. Rev.* **121** 5–19

[37] Burgess D and Scholer M 2007 Shock front instability associated with reflected ions at the perpendicular shock *Phys. Plasmas* **14** 012108

[38] Livadiotis G 2015 Shock strength in space and astrophysical plasmas *Astrophys. J.* **809** 111

[39] McKean M E, Omidi N and Krauss Varban D 1996 Magnetosheath dynamics downstream of low Mach number shocks *J. Geophys. Res.* **101** 20013–22

[40] Maruca B A *et al* 2018 MMS observations of beta-dependent constraints on ion temperature anisotropy in earths magnetosheath *Astrophys. J.* **866** 25

[41] Stasiewicz K and Eliasson B 2020 Quasi-adiabatic and Stochastic Heating and Particle Acceleration at Quasi-perpendicular Shocks *Astrophys. J.* **903** 57

[42] Forslund D W, Morse R L and Nielson C W 1970 Electron cyclotron drift instability *Phys. Rev. Lett.* **25** 1266–70

[43] Forslund D, Morse R, Nielson C and Fu J 1972 Electron cyclotron drift instability and turbulence *Phys. Fluids* **15** 1303–18

[44] Muschietti L and Lembege B 2013 Microturbulence in the electron cyclotron frequency range at perpendicular supercritical shocks *J. Geophys. Res.* **118** 2267–85

[45] Muschietti L and Lembege B 2017 Two-stream instabilities from the lower-hybrid frequency to the electron cyclotron frequency: application to the front of quasi-perpendicular shocks *Ann. Geophys.* **35** 1093–112

[46] Mikhailovskii A B and Timofeev A V 1963 Theory of cyclotron instability in a non-uniform plasma *JETP* **17** 919–21 http://jetp.ras.ru/cgi-bin/e/index/e/17/3/p626?a=list

[47] Krall A N and Liewer P C 1971 Low-frequency instabilities in magnetic pulses *Phys. Rev. A* **4** 2094–103

[48] Davidson R C, Gladd N T, Wu C S and Huba J D 1977 Effects of finite plasma beta on the lowerhybriddrift instability *Phys. Fluids* **20** 301–10

[49] Huba J D, Gladd N T and Papadopoulos K 1978 Lower-hybrid-drift wave turbulence in the distant magnetotail *J. Geophys. Res.* **83** 5217–26

[50] Drake J F, Huba J D and Gladd N T 1983 Stabilization of the lowerhybriddrift instability in finite plasmas *Phys. Fluids* **26** 2247–9

[51] Stasiewicz K 2020 Stochastic ion and electron heating on drift instabilities at the bow shock *Mon. Not. R. Astron. Soc.* **496** L133–7

[52] Hershkowitz N 1981 Double-layers and electrostatic shocks *J. Geophys. Rev.* **86** 3307–10

[53] Roth M, DeKeyser J and Kuznetsova M M 1996 Vlasov theory of the equilibrium structure of tangential discontinuities in space plasmas *Space Sci. Rev.* **76** 251–317

[54] Gueroult R, Ohsawa Y and Fisch N J 2017 Role of magnetosonic solitons in perpendicular collisionless shock reformation *Phys. Rev. Lett.* **118** 125101

[55] Kim S, Ha J-H, Ryu D and Kang H 2021 Microinstabilities in the transition region of weak quasi-perpendicular intracluster shocks *Astrophys. J.* **913** 35

[56] Hanson E L M, Agapitov O V, Mozer F S, Krasnoselskikh V, Bale S D, Avanov L, Khotyaintsev and Giles B 2019 Cross-shock potential in rippled versus planar quasi-perpendicular shocks observed by MMS *Geophys. Res. Lett.* **46** 2381–9