# Particle energization in colliding subcritical collisionless shocks investigated in the laboratory

A. Fazzini, W. Yao, K. Burdonov, J. Béard, S. N. Chen, A. Ciardi, E.

d'Humières, R. Diab, E. D. Filippov, S. Kisyov, et al.

► To cite this version:

A. Fazzini, W. Yao, K. Burdonov, J. Béard, S. N. Chen, et al. Particle energization in colliding subcritical collisionless shocks investigated in the laboratory. Astronomy and Astrophysics - A&A, 2022, 665, 10.1051/0004-6361/202243277. insu-03839650

HAL Id: insu-03839650

https://insu.hal.science/insu-03839650

Submitted on 4 Nov 2022

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
# Particle energization in colliding subcritical collisionless shocks investigated in the laboratory

A. Fazzini<fnref n="1" />, W. Yao<fnref n="1" />,<fnref n="2" />, K. Burdonov<fnref n="1" />,<fnref n="2" />,<fnref n="3" />, J. Béard<fnref n="4" />, S. N. Chen<fnref n="5" />, A. Ciardi<fnref n="2" />, E. d'Humières<fnref n="6" />, R. Diab<fnref n="1" />, E. D. Filippov<fnref n="3" />,<fnref n="7" />, S. Kisyov<fnref n="5" />,*, V. Lelasseux<fnref n="1" />, M. Miceli<fnref n="8" />,<fnref n="9" />, Q. Moreno<fnref n="6" />,<fnref n="10" />, S. Orlando<fnref n="9" />, S. Pikuz<fnref n="3" />,<fnref n="11" />, X. Ribeyre<fnref n="6" />, M. Starodubtsev<fnref n="7" />, R. Zemskov<fnref n="7" />, and J. Fuchs<fnref n="1" />

Received 5 February 2022 / Accepted 18 April 2022

## ABSTRACT

**Context.** Colliding collisionless shocks appear across a broad variety of astrophysical phenomena and are thought to be possible sources of particle acceleration in the Universe.

**Aims.** The main goal of our experimental and computational work is to understand the effect of the interpenetration between two subcritical collisionless shocks on particle energization.

**Methods.** To investigate the detailed dynamics of this phenomenon, we performed a dedicated laboratory experiment. We generated two counter-streaming subcritical collisionless magnetized shocks by irradiating two Teflon (C₂F₄) targets with 100 J, 1 ns laser beams on the LULI2000 laser facility. The interaction region between the plasma flows was pre-filled with a low-density background hydrogen plasma and initialized with an externally applied homogeneous magnetic field perpendicular to the shocks. We also modeled the macroscopic evolution of the system via hydrodynamic simulations and the microphysics at play during the interaction via particle-in-cell (PIC) simulations.

**Results.** Here, we report our measurements of the plasma density and temperature during the formation of the supercritical shocks, their transition to subcritical, and their final interpenetration. We find that in the presence of two shocks, the ambient ions reach energies around 1.5 times of those obtained with single shocks. Both the presence of the downstream zone of the second shock and of the downstream zone common for the two shocks play a role in the different energization: the characteristics of the perpendicular electric fields in the two areas indeed allow for certain particles to continue being accelerated or, at least, to avoid being decelerated.

**Conclusions.** The findings of our laboratory investigation are relevant for our understanding of the energy distribution of high-energy particles that populate the interplanetary space in our solar system and the very local interstellar medium around the heliopause, where observations have indicated evidence of subcritical collisionless shocks that may eventually go on to collide with one another.

**Key words.** shock waves – acceleration of particles – interplanetary medium

## 1. Introduction

Subcritical collisionless shocks are a class of shocks that are able to satisfy the Rankine-Hugoniot jump conditions using only dispersive and resistive dissipation mechanisms (Balogh & Treumann 2013). In subcritical shocks, the downstream flow velocity exceeds the sound speed behind the shock, but it is lower than the magnetosonic speed (defined as $c_{ms} = \sqrt{c_s^2 + v_A^2}$, where, $c_s$ and $v_A$ are the ion sound velocity and Alfvénic velocity, respectively). The limiting pre-shock magneto-

sonic Mach number $M_{ms} = v_s/c_{ms}$ for these conditions to be satisfied depends on the $\beta$ of the plasma and on the shock obliquity (i.e., the angle between the shock velocity and the upstream magnetic field). Subcritical shocks thus need to satisfy $M_{ms} \lesssim M_{ms}^{cr}$, with $M_{ms}^{cr}$ ranging between $M_{ms}^{cr} = 1$ for quasi-parallel shocks, up to $M_{ms}^{cr} = 2.76$ for perpendicular shocks, in the limit $\beta \to 0$ (Edmiston & Kennel 1984). In astrophysics, we can find these so-called subcritical shocks in a variety of scenarios. When a high-Mach number flow meets a dense medium, it becomes heavily “mass-loaded” and slows down to a velocity that allows for the formation of a subcritical shock. This is expected to happen when the solar wind interacts with the interstellar medium and forms the termination shock (Treumann 2009). Moreover,

<footnote n="1">¹ LULI - CNRS, CEA, UPMC Univ Paris 06 : Sorbonne Université, Ecole Polytechnique, Institut Polytechnique de Paris, 91128 Palaiseau Cedex, France
e-mail: julien.fuchs@polytechnique.edu</footnote>

<footnote n="2">² Sorbonne Université, Observatoire de Paris, Université PSL, CNRS, LERMA, 75005 Paris, France</footnote>

<footnote n="3">³ JIHT, Russian Academy of Sciences, 125412 Moscow, Russia</footnote>

<footnote n="4">⁴ LNCMI-T, CNRS, 31400 Toulouse, France</footnote>

<footnote n="5">⁵ "Horia Hulubei" National Institute for Physics and Nuclear Engineering, 077125 Bucharest-Magurele, Romania</footnote>

<footnote n="6">⁶ University of Bordeaux, Centre Lasers Intenses et Applications, CNRS, CEA, UMR 5107, 33405 Talence, France</footnote>

<footnote n="7">⁷ IAP, Russian Academy of Sciences, 603950 Nizhny Novgorod, Russia</footnote>

<footnote n="8">⁸ Università degli Studi di Palermo, Dipartimento di Fisica e Chimica E. Segrè, Via Archirafi 36, 90123 Palermo, Italy</footnote>

<footnote n="9">⁹ INAF, Osservatorio Astronomico di Palermo, Piazza del Parlamento 1, 90134 Palermo, Italy</footnote>

<footnote n="10">¹⁰ ELI-Beamlines, Institute of Physics, Czech Academy of Sciences, 5 Kvetna 835, 25241 Dolni Brezany, Czech Republic</footnote>

<footnote n="11">¹¹ NRNU MEPhI, 115409 Moscow, Russia</footnote>

<footnote>* Present address: Lawrence Livermore National Laboratory, Livermore, California 94550, USA.</footnote>
some astrophysical supercritical shocks evolve into subcritical ones in the course of their interaction with the upstream medium and consequent loss of energy, as it happens for solar coronal mass ejections (CMEs; Bemporad & Mancuso 2011). The collision of subcritical shocks is expected to occur between forward and reverse shocks in the solar wind, and also between solar wind shocks and planetary bow shocks (Whang & Burlaga 1985).

As in the case of their supercritical counterpart, namely, shocks with $M_{ms} > M_{ms}^{cr}$, where ions can be accelerated through a variety of mechanisms (Balogh & Treumann 2013; Marcowith 2016), subcritical shocks can also accelerate ions and induce thermal heating, although particle acceleration does not play a significant dissipative role in subcritical shocks. Both ion acceleration and heating have been observed in satellite crossings (Mellott 1984). The underlying ion acceleration mechanism(s) is still up for debate, but it has been suggested to include $v \times B$ heating (Ohsawa & Sakai 1985), ion reflection to a small degree from the shock front (Lee et al. 1987), and from other wave-particle interactions (Balikhin & Wilkinson 1996). As for particle acceleration from the collision of two subcritical shocks, no significant ion acceleration, with respect to the energies reached by particles accelerated by supercritical shocks, was observed in the simulations (Cargill et al. 1986). Accelerated ions with energy in the tens of MeV have been measured from the collision of two subcritical shocks at a small angle (Dudkin et al. 2000), however, their numbers are extremely small and there is still an ongoing effort to determine the acceleration mechanism.

Recently, high-power lasers and externally controlled magnetic field generation have opened the door to investigations of astrophysically relevant collisionless shock studies on particle acceleration (Li et al. 2019; Fiuza et al. 2020; Yao et al. 2021, 2022).

In this work, we detail our laboratory study of subcritical perpendicular collisionless shocks, that is, shocks within an external magnetic field perpendicular to the shock propagation direction. We characterize in detail their global spatio-temporal dynamics using multiple diagnostics. Moreover, we investigate the head-on encounter of two such shocks, in order to determine if and how this could modify the conditions under which ions can be energized in such a configuration.

The shocks were characterized in the laboratory by interferometry and Thomson scattering (TS) measurements, performed at different times, which provided the electron density map and local electron density, as well as the local ion and electron temperature. We then performed three-dimensional (3D) hydrodynamic simulations with the magnetohydrodynamic (MHD) FLASH code (Fryxell et al. 2000), which reproduced the global dynamics of both the expanding plasmas driving the shock and the shock itself. Next, we studied the event using the 1D3V fully kinetic particle-in-cell (PIC) code SMILEI (Derouillat et al. 2018), where we again used the experimentally obtained parameters as the initialization values. In these kinetic simulations, we observed particle acceleration of the ambient particles. The acceleration is initially due to the electrostatic field associated with the shock front $E_x$, then to the inductive electric field $E_y \sim v_x B_z$, where $v_x$ is the flow velocity and $B_z$ the perpendicular magnetic field. During the interaction between the two subcritical shocks, we note that both the presence of the downstream zone of the second shock and the creation of a downstream zone common for the two shocks play a role in the higher energization of the ambient ions: the characteristics of the perpendicular electric fields of these two areas allow, indeed, for certain protons to be accelerated or to avoid being decelerated.

<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Setup of the experiment, conducted at LULI2000 by having two high-power lasers (1 ns, 100 J at $1\omega$, $1.6 \times 10^{13}$ W cm$^{-2}$ on target) irradiate two solid (Teflon, $C_2F_4$) targets to investigate the interpenetration of two magnetized shocks. An auxiliary beam of 15 J was used to perform Thomson scattering (TS) and an additional low-energy beam (not shown) probed the plasma along a line titled 9° upwards with respect to the z-axis in order to measure the integrated plasma electron density.</figcaption></figure>

As a result, ambient ions were energized to 1.5 times the energy of the single shock case. This is consistent with space measurements performed in-situ of ions accelerated ahead of outward propagating interplanetary shocks (Gosling et al. 1984), or in the interaction of an interplanetary shock with the bow shock of the Earth (Hietala et al. 2011).

## 2. Experiment

### 2.1. Setup and diagnostics

The setup employed in our experiment is shown in Fig. 1. We started by irradiating two Teflon ($C_2F_4$) targets with two high-power laser pulses (1053 nm wavelength, 1 ns, 100 J, $1.6 \times 10^{13}$ W cm$^{-2}$ each). The targets were tilted in a way that allowed the laser beams to reach them and such that the two plasma flows would encounter each other (as detailed in Fig. 2). The two targets were separated by a 9 mm distance. The region in between the targets was pre-filled with hydrogen at low density ($n_0 \sim 10^{18}$ cm$^{-3}$), injected by a gas nozzle and magnetized using an externally applied magnetic field of 20 T provided by a pulsed coil (Albertazzi et al. 2013), directed along the z-axis.

A focusing spectrometer with spatial resolution (FSSR; Faenov et al. 1994) was utilized to register x-ray ion emission of the plasma with and without the ambient medium. It allowed us to characterize both the plasma initiated by the laser interaction at the surface of each target and the heating of the ambient medium induced by the expanding plasmas. The spectrometer was equipped with a spherically bent mica crystal with a lattice spacing $2d = 19.9149$ Å and curvature radius of $R = 150$ mm. It was able to measure He-like (transitions 3p–1s, 4p–1s, 5p–1s etc.) and H-like (transition 2p–1s and its satellites) lines of Fluorine in the range of wavelengths between 13 and 16 Å with a spatial resolution about 0.1 mm along the axis which joins the centers of both targets. The presence of sulfur impurities in the targets also allowed us to register a corresponding He-like doublet (2p–1s transition) in the third order of reflection with a Li-like satellite structure being sensitive in our range of plasma parameters. The spectral resolution was achieved better than $\lambda/d\lambda = 1000$. The spectra were recorded using Fuji-film Imaging Plates of type TR, which were placed in a cassette
<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Side (left) and top (right) view of the targets.</figcaption></figure>

holder protected from the visible optical radiation. The signal is time-integrated. The analysis of x-ray spectra was done by comparison of the experimental line ratios with simulated ones using the collisional-radiative code PrismSPECT (MacFarlane et al. 2003) and by comparison of emissivity profiles in different conditions.

A high-energy auxiliary beam (527.5 nm wavelength, 1 ns, 15 J, focused over ~40 µm along the y-axis and propagated throughout the plasma, see Fig. 1) was used to perform Thomson scattering (TS) measurements off the electron and ion waves in the plasma. It was used in a mode where the plasma was sampled in a collective mode (Froula et al. 2011). The collection of the scattered light was performed at 90° (along the z-axis) from the incident direction of the laser probe (the y-axis). The light scattered off the ion (TSi) and electron (TSe) waves in the plasma was analyzed by means of two different spectrometers, set to different dispersions (3.1 mm nm⁻¹ for TSi and 7.5 × 10⁻² mm nm⁻¹ for TSe), which were coupled to two streak-cameras (Hamamatsu for TSe, and TitanLabs for TSi, both equipped with S-20 photocathode to be sensitive in the visible part of the spectrum, and both with typical 30 ps temporal resolution), allowing us to analyze the evolution of the TS emission in time. The scattering volumes sampled by the instruments were: 120 µm along the x- and y-axes, 40 µm along the z-axis for TSi; 100 µm along the x- and y-axes, 40 µm along the z-axis for TSe. The analysis of the Thomson scattered light was performed by way of comparison between the experimental images (recorded by the streak cameras) and the theoretical curves of the scattered spectrum for coherent TS in non-collisional plasmas, with the instrumental function width of 5.9 nm for the electron spectrometer and 0.12 nm for the ion spectrometer taken into account. We point out that the TS laser probe induces some heating in the hydrogen ambient gas (details can be found in Yao et al. 2021). With TS, we can get a spatially and temporally resolved measurement of the plasma density and temperature. In addition, another optical probe beam (λ = 530 nm) passed with a 9° angle with respect to the B-field lines through the interaction zone, allowing for a measurement of the integrated electron density through interferometry.

## 2.2. Experimental results

The electron temperature on the target surface was measured via FSSR by recording the emission of sulfur lines and by simulating

<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Experimental x-ray spectrum (black line) measured by the FSSR spectrometer as emitted from a C₂F₄ target. What is recorded is the spectrum of sulfur impurities in the third order of reflection. Overlaid are simulations performed using the PrismSPECT code (red, olive, and blue curves) for the target surface region, using the group of satellites sensitive to the plasma parameters. For all temperatures shown in the figure, the electron density was $N_e = 7 \times 10^{20}$ cm⁻³. All curves are normalized to the S Heα line. The best fitting corresponds to the red curve. The inset shown in the top left corner demonstrates the detailed fitting of the satellites of the experimental spectrum. The arrows point to the lines with the best fit.</figcaption></figure>

this emission in a steady-state approach using the code PrismSPECT (MacFarlane et al. 2003). This is shown in Fig. 3, yielding for the surface plasma a temperature $T_e = 550$ eV at almost critical density $N_e = 7 \times 10^{20}$ cm⁻³. We point out that this measurement is relative to the laser-target interaction, that is, to the collisional part of the system.

After the plasmas are generated at the surface of each target, they expand into the ambient medium. This expansion is monitored by optical probing. This is displayed in Fig. 4, which shows the measurements, at successive times, of the integrated (along the line-of-sight of the probe beam) electron density of the plasmas expanding from both targets. We point out that these images were obtained on different shots. As seen in our previous experiment, where we created one single magnetized shock (Yao et al. 2021), two structures develop out of each target: a piston front and a shock front characterized by two separated bumps of higher electron density (identified by arrows in Figs. 4b and e). Each piston is the result of the expansion of the plasma ablated from the solid target by each laser. The plasma flows expand in the low-density ambient hydrogen, which is quickly ionized by the x-rays produced by the irradiated targets, and shocks are generated as a result of the combined action of the supersonic piston expansion and of the externally applied B-field. In fact, the strong external magnetic field of 20 T is critical in providing additional pressure so that a magnetized shock can form in the hydrogen plasma (Yao et al. 2021), as in its absence we would get a shock only in the presence of a denser background plasma. As a result, for early times, that is, before ~12 ns, we observe two well-developed shocks propagating against each other.

In our situation, the shocks are perpendicular, that is, the angle between the magnetic field and the shock propagation direction is $\theta_{Bn} \approx 90^\circ$ and the shocks are characterized by a
<figure><img src="image_5.png" /></figure>

**Fig. 4.** Temporal sequence of integrated (along the z-axis, with a 9° tilt) electron density measurements showing the evolution of the interpenetration of two magnetized shocks, at 6 ns (a), 11 ns (b), 13 ns (c), and 16 ns (d), after the main laser pulses hit the targets. The aperture of the magnetic coil structure restrained the passage of the optical probe, diminishing slightly the field of view (FoV) in *a–d*. (*e*) shows the line-out of the integrated electron density at 11 ns (blue solid line) and at 16 ns (red dashed line), along the lines shown on the relative map – (b) and (d), respectively; the location of the targets at 9 mm distance from each other is also shown, while the gray dashed areas represents the zones out of the FoV. Before the interpenetration of the two plasma flows, a shock front and a piston develop out of each target, as indicated by the two arrows in (b) and (e): the dotted blue arrow points at the left-drifting piston and the solid green arrow at the left-drifting shock front. When the two shocks collide, the density has a spike at the meeting point at around 3.5 mm in (e).

$\beta = P_{\text{therm}}/P_{\text{mag}} \approx 0.1$; hence, the critical Mach number has a value of $M_{\text{ms}}^{\text{cr}} \sim 2.6$ (Edmiston & Kennel 1984). The shocks obtained in our experiment are supercritical up to 3–4 ns after the laser beams hit the targets and turn into subcritical for later times (Yao et al. 2021). Indeed, they propagate with an initial velocity

of $v_s \approx 1500 \text{ km s}^{-1}$, which corresponds to $M_{\text{ms}} \approx 3.3 > M_{\text{ms}}^{\text{cr}}$, and when they eventually interact, they have a velocity of a few hundreds of km/s, which gives, for $v_s \approx 500 \text{ km s}^{-1}$, $M_{\text{ms}} \approx 1.1 < M_{\text{ms}}^{\text{cr}}$. We point out that the measurements of shock velocity are obtained from the interferograms by measuring the positions at different times, which correspond also to different shots. As we can observe, the structures developing from the two targets are of different sizes: they start forming at the same time, but their different distance from the gas nozzle exhaust makes them propagate in a medium of slightly different density which has a visible impact on their propagation velocities. Moreover, the fact that the two shocks do not propagate directly against each other, but perpendicularly to the targets, has also been taken into account while calculating the velocity. Indeed, the interferometry view corresponds to the side view (apart for a ~9°) of Fig. 2, which is a projection of the displacements along the z-axis. Hence, the distances extracted from the interferometry figures are multiplied by a factor of $1/\cos(60°) = 2$.

As for the collisionality, we find that the mean-free-path of the drifting ions with respect to the ambient ones is $\lambda_{mfp}^{i-i(d-a)} \approx 33 \text{ mm}$ (calculated according to Braginskii 1965), which is much larger than the characteristic length over which the interaction takes place (~hundreds of µm), thus making the shock collisionless.

Moreover, we measured the plasma Thomson scattering of the plasma thermal waves to assess the plasma characteristics. Figure 5 shows two examples of TS spectra from electron plasma waves (a, c) and from ion acoustic waves (b, d), corresponding to 15 ns and the period from 13 ns to 16 ns, respectively. The temporal evolution of the electron density, the electron temperature, and ion temperature is shown in Fig. 6. We observe that after around 13 ns the TSi signal suddenly broadens, which corresponds to the time of the collision between the two shocks. This broadening is attributed to the heating of the ions in the plasma due to the energy released when the two plasma bubbles collide. As is shown in the time evolution in Fig. 6, after the collision, the electron density slightly increases, while the electron temperature remains initially unperturbed (around 80 eV). The interpenetration of the two plasma shocks heats the ions up to temperatures $\approx 135 \text{ eV}$, according to an adiabatic gas compression. Electrons are then heated at a slower rate by ion-electron collisions. Further increase of ion and electron temperature as well as electron density is observed when the pistons collide at a later time (~14.5 ns). Here, we will only focus on the shock-shock collision, before the encounter of the pistons.

Complementary to the TS measurements, Fig. 7 shows the x-ray emissivity profiles of the plasma located in between the two targets. We compare three cases: (1) when the applied magnetic field and the ambient medium are present and the plasmas expand from the two targets (thin red line); (2) still with two plasmas, but in the absence of the ambient medium (dashed black line); and (3), when only one plasma is flowing from either the right or the left target, but in the same magnetic field and ambient medium conditions as in case (1) – thin gray lines, with areas filled by patterns. What we observe is that the collision between the two plasma flows results in an emission enhancement in the zone between the two targets (compare the red curve to the filled areas). Here, the left target has a lower intensity due to the positioning of the corresponding part of the spectrometer in the “shadow zone” of the right target. The electron temperature was measured in this region as $T_e = 240 \text{ eV}$ (a lower-limit estimate) at an electron density $n_e = 10^{18} \text{ cm}^{-3}$ using the ratio between the resonance lines $Ly_\alpha$ and $He_\beta$ by the method described in Khiar et al. (2019), in reasonable agreement
<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Thomson scattering measurements of the plasma density and temperatures in the region of shock collision. Spectra of Thomson scattering off electron plasma waves (a, c) and ion acoustic waves (b, d). Then, (a) and (b) show the spectra profiles, corresponding to 15 ns, while (c) and (d) show the temporal evolution of the scattering spectra over a time period from 13 ns to 16 ns for the electron plasma and ion acoustic waves, correspondingly. Black solid lines (in a and b) are for experimental data profiles, while red solid lines are for theoretical spectra, composed of a superposition of narrow (black dotted lines) component relative to the ambient medium, having density $1.5 \times 10^{18}$ cm$^3$, electron temperature 100 eV, ion temperature 200 eV, and broad (black dashed lines) component relative to the piston plasma, having density $6 \times 10^{18}$ cm$^3$, electron temperature 300 eV, ion temperature 100 eV. The ratio between the magnitudes of the narrow and broad components is 3.5. We note that the deep central dip in the experimental spectra (a, b) and the white vertical region in the streak-camera images (c, d) is related to a filter (a black aluminum stripe) which is positioned right before the entrance of the two streak cameras (recording respectively the light scattered off the electron and ion waves). This filter is used to block the very intense and unshifted laser wavelength (the Rayleigh-scattered light), which otherwise would saturate the cameras. Thus, no signal is recorded in this zone, which is materialized by the gray dashed box.</figcaption></figure>

with the TS measurements, knowing that the X-ray diagnostics is time-integrated. In addition, one can note that the emissivity drops significantly faster between the target and the middle zone in the case when the ambient medium is applied, which is most probably related to a faster recombination rate as well as to a higher confinement of particles close to the target (Filippov et al. 2021).

The main plasma parameters extracted from the experimental measurements are summarized in Table 1. These values are used to initialize our simulations detailed below, which we use to further investigate the particle acceleration during the shock collision.

## 3. Numerical simulations

Our simulation effort is two-fold: the first step was to undertake MHD simulations of the laser-driven plasma expansion and interaction with the ambient gas and magnetic field, leading to the experimentally observed piston. In the second step, we used the results of the experimental diagnostics as a starting point for kinetic simulations that allow us to investigate in detail the microphysics of the shocks colliding and the underlying particle acceleration mechanisms.

<figure><img src="image_7.png" /><figcaption>**Fig. 6.** TS measurements of the temporal evolution of electron density (a), electron temperature (b), and ion temperature (c). The narrow and broad configurations are respectively related to the ambient and piston plasmas (for more details, see Fig. 5 caption).</figcaption></figure>

### 3.1. MHD simulations

The experiment was first modeled with the 3D MHD code FLASH (Fryxell et al. 2000). We model the formation and the propagation of pistons and shocks generated by the laser interaction with two Teflon targets having the same arrangement as shown in Fig. 2. However, to reduce the computational cost, the separation between the targets is here limited to 6.5 mm, instead of 9 mm as in the experiment. As in the experiment, the targets are embedded inside an ambient hydrogen gas-jet within an
<figure><img src="image_8.png" /><figcaption>**Fig. 7.** X-ray emissivity profiles measured by the FSSR spectrometer in different cases. The curves obtained when both plasmas expand out of the targets are represented by the dashed black line and the full red line, respectively for the cases of absence or presence of the ambient medium. The gray curves with a pattern correspond to a single target case as a reference. In all cases, the magnetic field is present. All curves are normalized to the right target emissivity. The emissivity for the left target is the inverted right one with a multiplier taking into account the signal reduction due to the location of a part of the spectrometer in the shadow zone.</figcaption></figure>

<table><caption>**Table 1.** Parameters extracted from our measurements at ~11 ns, i.e., right before the interpenetration of the two shock structures.</caption><thead><tr><th colspan="2">Characterized ambient plasma conditions</th></tr></thead><tbody><tr><td>Upstream elec. number density n<sub>e</sub> [cm<sup>-3</sup>]</td><td>1.0 &times; 10<sup>18</sup></td></tr><tr><td>Upstream elec. temperature T<sub>e</sub> [eV]</td><td>80</td></tr><tr><td>Upstream ion temperature T<sub>i</sub> [eV]</td><td>20</td></tr><tr><td>Downstream elec. temperature T<sub>e</sub> [eV]</td><td>130</td></tr><tr><td>Downstream ion temperature T<sub>i</sub> [eV]</td><td>200</td></tr><tr><td>Shock velocity at meeting point v<sub>s</sub> [km s<sup>-1</sup>]</td><td>~500</td></tr><tr><td>Upstream magnetic field strength B<sub>z</sub> [T]</td><td>20</td></tr><tr><th colspan="2">Calculated parameters</th></tr><tr><td>Ion collisional mean-free-path &lambda;<sup>i-i(d-a)</sup><sub>mfp</sub> [mm]</td><td>33</td></tr><tr><td>Flow ion Larmor radius r<sub>L,i,f,l</sub> [mm]</td><td>0.26</td></tr><tr><td>Upstream plasma thermal beta &beta;<sub>ther</sub></td><td>0.10</td></tr><tr><td>Mach Number M</td><td>4.42</td></tr><tr><td>Alfvénic mach number M<sub>A</sub></td><td>1.15</td></tr><tr><td>Magnetosonic mach number M<sub>ms</sub></td><td>1.12</td></tr></tbody></table>

**Notes.** $\lambda_{mfp}^{i-i(d-a)}$ is the collisional mean free path between drifting and ambient ions.

external magnetic field. The laser intensity, the hydrogen gas-jet density, and the external magnetic field strength are the same as in the experiment.

Figure 8 shows the electron density time evolution in the case with the external magnetic field (20 Tesla), that is, before the shocks collision at $t = 5.8$ ns, at collision time $t = 7.2$ ns and after the collision at $t = 7.8$ ns, respectively. Due to the reduced distance between the targets used in the simulation, to scale it with the experiment, the collision time should be scaled by a factor 1.4, resulting in a scaled collision time of 10 ns, which is quite close to the experimentally observed one (~12 ns). When

<figure><img src="image_9.png" /><figcaption>**Fig. 8.** Simulation, using the 3D MHD code FLASH, of the volumetric electron density plotted at $t = 5.8$ ns (a), $t = 7.2$ ns (b), and $t = 7.8$ ns (c), along the laser beam direction, with external magnetic field $B = 20$ T.</figcaption></figure>

the two shocks collide, the electron density increases only of 20%. The pistons expand more slowly due to the increase of the magnetic pressure behind the shock.

These hydrodynamic simulations performed with FLASH are capable of describing the overall dynamic of the system, but they were not able to quantitatively reproduce the temperatures measured in the experiment (Yao et al. 2021), likely due to the fact that the kinetic effects associated to our collisionless system cannot be taken into account. That is why we also performed PIC simulations, so that we may take them into account.

## 3.2. Particle-in-cell simulations

The interaction between the two subcritical shocks has been modeled via the fully kinetic particle-in-cell (PIC) code, SMILEI (Derouillat et al. 2018), for which we used profiles of the plasma density, temperature, and magnetic field extracted from the experimental data as initial conditions (see Table 1).

We simulated such a system in a 1D3V geometry, as the scale of the shock front interaction with the ambient medium is much smaller across the shock (a few hundreds of microns) than along the shock (a few mm). We point out that our PIC simulations are dedicated to capture only the kinetic effects of the shock colliding process. The laser-target ablation and piston formation are well-reproduced by the FLASH simulations and the shock formation and transition from supercritical to subcritical is detailed in our previous papers (Yao et al. 2021, 2022).

In order to understand the effects of the collision of two shocks, we simulated both a single drifting shock configuration and a double counter-streaming shocks scenario. The initial configurations are shown in Fig. 9. The box has a length of $L_x = 2048$ $d_e \approx 11$ mm and the spatial resolution is $d_x = 0.2$ $d_e \approx 1.1$ µm, where $d_e = c/\omega_{pe} \approx 5.3$ µm is the electron inertial length, and $\omega_{pe} = \sqrt{n_0 q_e^2/(m_e \epsilon_0)} \approx 5.6 \times 10^{13}$ rad s$^{-1}$ is the electron plasma angular frequency. Here, c is the velocity of light, $n_0 = 1.0 \times 10^{18}$ cm$^{-3}$ is the electron (and proton) number density of the ambient plasma, and $m_e$, $q_e$ and $\epsilon_0$ are the electron mass, the elementary charge, and the permittivity of free space, respectively. We point out that the x-axis we are talking about here in the case of PIC simulation does not correspond to the one relative to the experimental setup (used in Fig. 1, 2, 4) and to the FLASH simulations (Fig. 8). In our PIC simulation, the x-axis is the axis along which the two shocks propagate. Each cell has 1024 particles plus 1 tracked particle for each species. Moreover, an external uniform magnetic field $B_{z0} = 20$ T is set in the z-direction perpendicular to the plasma
<figure><img src="image_10.png" /><figcaption>**Fig. 9.** 1D PIC simulations initialization setups. (a) Single-shock case: a hydrogen plasma ($n_1 = 2n_0 = 2 \times 10^{18} \text{ cm}^{-3}$, $T_{e1} = 130$ eV and $T_{i1} = 200$ eV) drifts through a background hydrogen plasma ($n_1 = n_0 = 1 \times 10^{18} \text{ cm}^{-3}$, $T_{e2} = 80$ eV, and $T_{i2} = 20$ eV) with a drifting velocity of $v_d = 350 \text{ km s}^{-1}$ (see Table 2). (b) Double-shock case: a background plasma ($n_2 = n_0 = 1 \times 10^{18} \text{ cm}^{-3}$, $T_{e2} = 80$ eV, and $T_{i2} = 20$ eV) is set at rest between two counter-streaming denser plasmas ($n_1 = n_3 = 2n_0 = 2 \times 10^{18} \text{ cm}^{-3}$, $T_{e1} = T_{e3} = 130$ eV and $T_{i1} = T_{i3} = 200$ eV). The drifting velocity $v_d = 350 \text{ km s}^{-1}$ imposed to the protons in both configurations leads to a shock velocity of $v_s \approx 640 \text{ km s}^{-1}$. The simulation is initialized ~11 ns after the lasers start ablating the targets, with a distance of 1.8 mm in between the two shock fronts, while the box has a total length $L_x = 11$ mm.</figcaption></figure>

velocity ($\omega_{ce}/\omega_{pe} = 0.06$, where $\omega_{ce} = q_e B_{0z}/m_e$). The simulation lasts for $1.5 \times 10^5 \omega_{pe}^{-1} \approx 2.5$ ns, with an initial time that corresponds to ~11 ns in the experiment, i.e., ~1 ns before the shocks collide. In short, the simulation covers time between 11 and 13.5 ns of the experiment. These times will be used below for better comparison with the experiment.

We initialized the system with different portions of plasma, all composed of protons and electrons with $m_p/m_e = 1836$. For a single subcritical shock, we set, from 0 to $L_x/2$ (Region 1 in Fig. 9a), a hydrogen plasma drifting towards positive $x$ at a velocity of $v_x = v_d = 350 \text{ km s}^{-1}$ and density $n_1 = 2n_0 = 2 \times 10^{18} \text{ cm}^{-3}$, while between $L_x/2$ and $L_x$ (Region 2) we put a background hydrogen plasma with density $n_2 = n_0$ at rest. In the case of double shock, we set two counter-streaming hydrogen plasmas with densities $n = 2n_0$ moving in the x-direction at velocities $v_d$ and $-v_d$ between 0 and $5L_x/12$ (Region 1 in Fig. 9) and between $7L_x/12$ and $L_x$ (Region 3), respectively. Moreover, a background hydrogen plasma with density $n_0$ was set at rest in between (Region 2). As for the temperature, we used the results of the TS diagnostic in the experiment for both simulations, namely, $T_{e1} = T_{e3} = 130$ eV and $T_{i1} = T_{i3} = 200$ eV for the drifting plasmas, and $T_{e2} = 80$ eV and $T_{i2} = 20$ eV for the background plasma. The list of parameters used to initialize our simulations is summarized in Table 2.

Figure 10 shows the spatial profiles of ion density, magnetic field $B_z$, and electric fields $E_x$ and $E_y$, for both single shock and double shock simulations. We point out that the profiles relative to the single shock are shifted by $\Delta x = -L_x/12$ in order to simplify the comparison.

In Fig. 10a at 1 ns ($\approx 12$ ns for the experiment), we clearly see that the shock formation and propagation happens in the same way for both the single shock case and the double shock case, as the profiles relative to the single shock overlap the ones of the right-drifting shock in the double shock case. In other words, the presence of the left-drifting shock has no effect on the evolution of the right-drifting one yet and vice versa. In the downstream regions, we note the spontaneous formation of fast magnetosonic waves propagating away from the shock fronts (Moreno et al. 2019). We point out that in the double-shock case, the electric field components $E_x$ and $E_y$ of the right-drifting shock and the left-drifting one are going in opposite directions: the shock com-

<table><caption>**Table 2.** List of parameters we used to initialize our PIC simulations.</caption><thead><tr><th rowspan="2">PIC initialization parameters</th><th colspan="2">Regions</th></tr><tr><th>1 and 3</th><th>2</th></tr></thead><tbody><tr><td>$v_d$ [km s<sup>-1</sup>]</td><td>350</td><td>0</td></tr><tr><td>$v_s$ [km s<sup>-1</sup>]</td><td>640</td><td>0</td></tr><tr><td>$B_{z0}$ [T]</td><td>20</td><td>20</td></tr><tr><td>$T_e$ [eV]</td><td>130</td><td>80</td></tr><tr><td>$T_i$ [eV]</td><td>200</td><td>20</td></tr><tr><td>$n_i$ [10<sup>18</sup> cm<sup>-3</sup>]</td><td>2</td><td>1</td></tr><tr><td>$v_A$ [km s<sup>-1</sup>]</td><td>308</td><td>436</td></tr><tr><td>$c_s$ [km s<sup>-1</sup>]</td><td>144</td><td>102</td></tr><tr><td>$c_{ms}$ [km s<sup>-1</sup>]</td><td>340</td><td>448</td></tr><tr><td>$r_{Li}$ [mm]</td><td>0.26</td><td>0.39</td></tr><tr><td>$\lambda_{mfp}^{i-i(d-a)}$ [mm]</td><td colspan="2">33</td></tr><tr><td>$M_{ms}$</td><td colspan="2">1.43</td></tr></tbody></table>

**Notes.** Regions 1 and 3 correspond to the drifting hydrogen plasmas, while region 2 is relative to the background hydrogen plasma at rest. $\lambda_{mfp}^{i-i(d-a)}$ is the mean free path relative to the collisions between drifting and ambient ions, while $M_{ms}$ is the magnetosonic Mach number of the shock wave moving at speed $v_s = 640 \text{ km s}^{-1}$ in the ambient plasma characterized by a magnetosonic speed $c_{ms} = 448 \text{ km s}^{-1}$. We point out that the values relative to regions 1 and 3 refer only to the initial situation, before the shocks completely form, hence, they are not to be confused with the parameters of the downstream region at a certain time.

ing from the left is characterized by oscillations starting with a positive peak of $E_x$ at the shock front and by a downstream region with $E_y > 0$, while the shock coming from the right has oscillations starting with a peak $E_x < 0$ at the shock front and a field $E_y < 0$ in the downstream region.

In Fig. 10b, we can see that at 1.5 ns ($\approx 12.5$ ns for the experiment), the interaction of the two shocks has begun: the ion density and magnetic field start overlapping, and, at longer times not shown here, they keep piling up, reaching values of $n_i \approx 5 \times 10^{18} \text{ cm}^{-3} = 5n_0$ and $B_z \approx 50 T = 2.5 B_{z0}$ at $\sim 2.2$ ns ($\approx 13.2$ ns in the experiment). The x-component of the electric field $E_x$ fluctuates around 0 and has peaks of $\sim 100 \text{ MV m}^{-1}$ associated with the shocks. After the interpenetration of the two shocks, $E_x$ reaches values of $\sim 250 \text{ MV m}^{-1}$. The only contributions to the y-component $E_y$ is the inductive electric field $E_y = -(v \times B)_y = v_x B_z$ (Ilie et al. 2017) and it, too, has a fluctuating profile, centered on 0 in the upstream regions and on $\sim \pm 8 \text{ MV m}^{-1}$ in the downstream ones. Between these two zones, $E_y$ passes gradually from $8 \text{ MV m}^{-1}$ (or $-8 \text{ MV m}^{-1}$) to 0, as the plasma velocity distribution decreases (increases) and the magnetic field increases (decreases), even after the two shocks have met.

By tracking a set of representative protons, we were able to understand the energization mechanism undergone by the most energetic ones, namely, the ones reaching a kinetic energy of $E_k > 10$ keV. In Figs. 11a and b, we show the motion in the $v_x-v_y$ space of protons from the drifting and the ambient plasma, respectively. For clarity, we plotted, for each plasma, only one particle for each case, but we checked that these trajectories are well representative of all other tracked particles coming from the same populations. We observe that the protons of the drifting plasma, having initially a bulk velocity $v_x = 350 \text{ km s}^{-1}$, rotate in the $v_x-v_y$ space, without showing any special difference between the single and double shock cases.
<figure><img src="image_11.png" /><figcaption>**Fig. 10.** Profiles of ion density, magnetic field $B_z$, and electric fields $E_x$ and $E_y$ at times 1 ns (a) and 1.5 ns (b) after the beginning of the simulation. At 1 ns ($\approx$12 ns for the experiment), the two fast shocks propagate without perturbing each other, as their profiles correspond to the one of the single shock case. At 1.5 ns ($\approx$12.5 ns for the experiment), the interaction between the two shocks makes the ion density and the magnetic field $B_z$ increase and pile up in the middle; the electric field $E_x$ presents several spikes also due to the interaction between the downstream fast magnetosonic waves; the inductive $E_y$ does not drastically change structure or magnitude, but it presents opposite signs in the downstream regions depending on the propagation direction of the shock. We point out that the profiles relative to the single shock are shifted of $\Delta x = -L_x/12$ in order to simplify the comparison.</figcaption></figure>

The situation is definitely different for the protons of the ambient plasma: after starting at rest, the protons are accelerated by the shocks up to kinetic energies which, in the double shock case, are around 1.5 times higher than the single shock case. This difference is well presented by the proton spectra at the final simulation time shown in Fig. 11c and, zoomed, in Fig. 11d, where to higher energies are associated higher distribution values in the double-shock case.

By comparing the cases with one and two subcritical shocks, we could understand the reason of the higher energization in

<figure><img src="image_12.png" /><figcaption>**Fig. 11.** Comparison of the trajectories in the $v_x - v_y$ space of two typical energetic tracked drifting (a) and ambient (b) protons, for the double shock case in red and for the single shock case in blue. (We note that their initial velocities might be different from the species bulk speed ($v_x = 350$ km s$^{-1}$), since the protons are initiated with an initial temperature.) (c) Final energy spectra of the ions for the two configurations and (d) zoom on the range from 10 to 30 keV.</figcaption></figure>

the double shock case by analyzing the dynamics of such energetic tracked particles. Let us start considering the protons of the ambient plasma whose dynamics is reported in Fig. 11b. For the proton from the single shock case, we analyzed its motion in the $x-t$ space over a map of $E_x$ and a map of $E_y$ (Figs. 12a1 and b1, respectively), and in the $v_x - v_y$ space (Fig. 12c1), where the color of the proton trajectory follows a scale based on its kinetic energy $K$. Moreover, we compared the temporal evolution of the proton kinetic energy and the work done by $E_x$ and $E_y$ on it as shown in Fig. 12d1. On these four graphs (a1–d1), we distinguish the presence of three phases corresponding to different regimes experienced by the proton.

In phase I, the proton is accelerated by the electrostatic field $E_x$ associated with the shock front, it gains a velocity $v_x > 0$ and starts rotating in the upstream plasma. While gyrating clockwise due to the applied magnetic field $B_z$, its velocity in the $x$-direction $v_x$ decreases and the proton meets again the shock front. After crossing it, in phase II, the proton is in the downstream region characterized by $E_y > 0$. Since it has a $v_y < 0$,
<figure><img src="image_13.png" /><figcaption>**Fig. 12.** Comparison of the trajectories of two typical energetic protons of the ambient plasma in the single- and double-shock cases in the *left column* (1) and *right column* (2), respectively. In their evolution between 0 and 2.5 ns we can distinguish different phases (I, II, II*a*, II*b*, III), that are detailed in the main text. In (a1–a2) and (b1–b2), the maps of the $E_x$ and of the $E_y$ fields in the *x–t* space, together with the trajectories of the chosen protons, are shown, respectively. Note the different scales for the electric fields $E_x$ and $E_y$ color maps. In (c1–c2), the proton trajectories in the $v_x – v_y$ space and the temporal points delimiting the different phases, are plotted. The color of the proton trajectories in (a1–a2), (b1–b2), and (c1–c2) follows the evolution of their kinetic energy $K$ on a color scale from 0 to 20 keV. In (d1–d2), are shown the temporal evolution of the kinetic energies of the two tracked protons (full green line), the work done by $E_x$ (dashed red line) and by $E_y$ (dashed blue line).</figcaption></figure>

the positive $E_y$ does a negative work on the proton and lowers its kinetic energy, as can be seen in Fig. 12d1. In phase III, the particle has again a velocity $v_y > 0$, hence the positive $E_y > 0$ does a positive work on the proton and keeps energizing it.

We conducted a parallel analysis on the proton from the ambient plasma in the double shock case and plotted the results on the right column of Fig. 12: (a2) and (b2) show its motion in the *x–t* space over a map of $E_x$ and $E_y$, respectively; (c2) shows its trajectory in the $v_x – v_y$ space; (d2) shows the temporal evolution of the particle kinetic energy and the work made on it by $E_x$ and $E_y$. The color of the proton trajectory in Figs. 12a2–c2

follows the same kinetic energy scale used for the single shock case.

We can now distinguish four different regimes for the tracked proton. In phase I, similarly to what happens in the single shock case, the proton is accelerated by the right shock front and rotates in the upstream region. Phase II*a* and II*b* are due to the presence of the second shock (coming from the right side) and hence are different to phase II in the single shock scenario. Specifically, in phase II*a*, the proton interacts with the shock front coming from the left, whose $E_x$ oscillates but does a net positive work on the proton, increasing its energy. In phase II*b*, the proton encounters the shock coming from the right and finds itself in the region downstream of both shocks, which is characterized by a low $E_y$ (due to the colliding of the two shocks) and by an oscillating $E_x$ (due to the interaction of the two magnetosonic waves), that do not change its energy in a net way. Phase III in the double shock case is similar to the one for the single shock: the proton is in the downstream zone of the shock coming from the right and has a velocity $v_y < 0$, hence the negative $E_y < 0$ makes a positive work and keeps energizing it.

Other energetic protons show a similar behavior: in general, the higher energization in the double shock case is mostly due to the fact that while it is gyrating because of the magnetic field, the proton can find areas with $E_y$ directed accordingly to its velocity $v_y$ or zones downstream (of both shocks) with low $E_y$ that does not decelerate it either. These two contributions have the effect of accelerating some of the ambient ions to higher energies in the case of two counter-propagating shocks than in the presence of only one.

Among the tracked background ambient protons, we can compare how many reach, for instance, at least 10 keV in the cases of single and double shock. In the presence of a single shock only 1 proton out of 5120 (≈0.02%) reaches 10 keV or higher energies, while with the second shock we have 20 out of 1706 (≈1.17%); hence, a much higher percentage becomes accelerated to higher energies. We note that the inequality of number of tracked protons in the two cases comes from the difference of ambient plasma size.

We point out that we are not able to compare the experimental proton spectrum with the simulated one: during the first 3–4 ns of evolution, the shocks are supercritical and the individual interaction of both of them with the background plasma already leads to some high particle energization, as shown in our previous work (Yao et al. 2021, 2022). As the shocks propagate, their velocity drops, leading to an interaction between two subcritical shocks. This interpenetration is indeed the only part that we are simulating here, hence, the resulting spectrum cannot take into account the protons previously accelerated by the supercritical shocks. The protons accelerated by supercritical shocks reach much higher energies, hence in the final spectrum obtained experimentally they “cover” the portion of protons accelerated to lower energies by the interaction between the two subcritical shocks. We also considered, as detailed in the appendix, the limit of validity of the 1D framework used here.

## 4. Astrophysical relevance

Although the vast majority of collisionless shocks in astrophysics are supercritical, the results of our experiment can be relevant for some phenomena observed in both interplanetary and astrophysical plasmas and involving subcritical shocks. In fact, observations of subcritical shocks are sparse and most of them are restricted to the interplanetary space, where colliding shocks can also be observed (e.g., Colburn & Sonett 1966).
<table><caption>**Table 3.** Dimensionless quantities relative to our laser-driven shock and the interstellar medium shock characterized by Voyager in *Burlaga et al.* (2013), as a weak subcritical resistive laminar shock.</caption><thead><tr><th>Parameters</th><th>Our shocks</th><th>Shocks as in (1)</th></tr></thead><tbody><tr><td><i>B</i><sub>DS</sub>/<i>B</i><sub>US</sub></td><td>1.5</td><td>1.4</td></tr><tr><td><i>β</i><sub>ther,US</sub></td><td>0.1</td><td>0.23</td></tr><tr><td>(<i>v</i><sub>A</sub>/<i>c</i><sub>s</sub>)<sup>2</sup><sub>US</sub></td><td>4.6</td><td>5</td></tr><tr><td><i>θ</i><sub>Bn</sub></td><td>90</td><td>85</td></tr><tr><td><i>M</i><sub>ms</sub></td><td>1.12<sup>(exp)</sup>−1.43<sup>(sim)</sup></td><td>≈1.9</td></tr></tbody></table>

**Notes.** *DS* and *US* refer to the downstream and upstream zones, respectively.

**References.** (1) *Burlaga et al.* (2013).

Examples are the interplanetary forward shocks convected with the solar wind that are expected to propagate with the wind mostly outward into the outer heliosphere, so that they can have a sufficiently low Mach number to be subcritical. Furthermore, subcritical shocks are expected in cometary environments, where the collision between the solar wind and the atmosphere of the comet can reduce the upstream flow velocity, resulting in low Mach number cometary bow shocks. CMEs can also produce subcritical shocks. The observations show that CMEs resulting from x-ray flares of the solar corona associated with type-II radio bursts (so-called radio-loud CMEs) are, in general, very fast and extended and, initially, they produce supercritical shocks that become subcritical at later times (e.g., Bemporad & Mancuso 2011, 2013); CMEs resulting from flares not producing a type-II burst (radio-quiet CMEs) produce subcritical shocks at all times (e.g., Bemporad & Mancuso 2011, 2013). All these subcritical shocks are present in the interplanetary space and can interact with each other or with planetary bow shocks, thus contributing to the acceleration of particles (electrons, protons, ions) up to near-relativistic energies. These shock-shock interactions are reproduced by our experiment, showing that ambient ions can be energized around 1.5 times more than in single shocks and that the particles are accelerated in different ways, depending on the areas of the shock-shock interaction region where the particles are located. The implication is that the population of high-energy particles in the interplanetary space may depend on the rate of occurrence of shock-shock interactions.

Thanks to the Voyager 1 and 2 missions, as of 2012 it has been possible to probe the density of the very local interstellar medium with accurate in situ measurements. This has allowed for the presence of several shock waves in the interstellar plasma to be revealed, which are most likely interplanetary shocks originating from energetic solar events (e.g., CMEs) that traveled outward through the supersonic solar wind and, after colliding with the heliospheric termination shock, crossed through the heliopause into the interstellar medium (Burlaga et al. 2013; Gurnett et al. 2013; Liu et al. 2014). Burlaga et al. (2013) reported the first in situ measurement of a shock in interstellar plasma, whose characteristics are summarized and compared to the ones relative to our shocks in Table 3. Then, evidence of multiple shocks was reported in 2015 data collected with Voyager 1 (Ocker et al. 2021). According to the Voyager measurements these shocks are, in general, weak low beta and subcritical shocks (Mellott & Greenstadt 1984; Burlaga et al. 2013; Mostafavi & Zank 2018). Thus, we can argue that the results of our experiment may be applicable to the interactions between these subcritical shocks which populate the very local interstellar medium. Interestingly, these shocks of solar origin are char-

acterized by a precursor consisting of various disturbances in the intensity and anisotropy of galactic cosmic rays (Gurnett et al. 2015). Voyager missions revealed that these disturbances are typically preceded by bursts of high-energy (≈5–100 MeV) electrons, most likely due to the reflection and acceleration of cosmic-ray electrons by magnetic field jumps at the shock or due to interactions (or both) with upstream plasma waves or shocks (Gurnett et al. 2021). Our experiment shows that the interaction between these subcritical shocks has a direct effect on the way particles are accelerated by the shocks and on the maximum energization of particles.

Possible shock-shock interactions as those discussed above can also be present (and play a significant role in the acceleration of particles) in the environments around exoplanets. In addition to cases analogous to those we observe in our solar system, the cases of hot-Jupiters (gas giant exoplanets that should be similar to Jupiter but are in close proximity to their stars) are of particular interest, given the strong interaction with their host stars via the stellar wind, the magnetic field and the irradiation. Depending on the parameters of the star-planet system (distance between the two objects, masses of the star and the planet, wind velocity, stellar irradiation, etc.) complex flow structures can form from the colliding planetary and stellar winds, as bow shocks, cometary-type tails, and inspiraling accretion streams (e.g., Matsakos et al. 2015). In particular, the speed of the planetary wind is, in general, only marginally supersonic, so that the Mach number is low (e.g., Tremblin & Chiang 2013). Also the cometary-type tails around the exoplanets are advected by the stellar winds (so that possible shocks can have low Mach numbers to be subcritical) and can be highly perturbed (producing a highly variable complex pattern of shocks), depending on the parameters of the star-planet system. Under these conditions, interactions between subcritical shocks may develop and can be analogous to those produced in our experiment.

The results of our experiment may also be relevant to shed light on the origin of the so-called suprathermal halo which characterizes the velocity distribution of electrons in the solar wind (e.g., Pierrard et al. 2011). In fact, observational evidence suggests that this electron velocity distribution consists of two components: a dominant low-energy core of thermal (Maxwellian) origin and a high-energy non-thermal tail. The latter component corresponds to a suprathermal halo which is nearly isotropic in the velocity space and can be described by power-law kappa-like distribution functions (Pierrard et al. 2001, 2016; Štverák et al. 2009; Maksimovic et al. 2005; Kajdič et al. 2016; Berčič et al. 2019) and a strahl component which appears as a magnetic-field-aligned beam streaming away from the Sun and which is more prominent during solar energetic events (e.g., in the fast component of the solar wind or during CMEs; Pilipp et al. 1987; Anderson et al. 2012). The origin of the high-energy non-thermal tails is widely debated in the literature (e.g., Seough et al. 2015; Boldyrev & Horaites 2019; Horaites et al. 2019). In particular, as it concerns the halo, some authors have suggested that it may originate from the solar corona in the form of high-energy electrons that escape from the Sun with velocities much higher than those of the electrons that populate the core. The nearly isotropic distribution in the velocity space may be due to a combination of Coulomb collisions and magnetic defocusing (e.g., Boldyrev & Horaites 2019). An alternative mechanism that was suggested considers the halo generated from the strahl electrons through strong angular scattering due to local interaction with ambient plasma turbulence (e.g., Štverák et al. 2009; Boldyrev & Horaites 2019). Our work shows that a
single event of shock-shock collision can lead to enhanced heating and particle acceleration. An implication could be that such an event – and even more so a multiplicity of such events (the interplanetary space is expected to be populated by a multitude of shocklets) – could contribute to a turbulent high-energy spectrum of electrons, as observed in the solar wind, thus putting forth a possible explanation for the origin of the halo electrons. More in general, similar effects can be expected in turbulent stellar winds, in supernova remnants in the radiative phase, or in environments in which turbulence is a prominent component. In order to challenge this scenario, however, it is necessary to study the electron energy distribution resulting from a multiplicity of events similar to that investigated here. This would be feasible using advanced laser facilities, such as NIF (Moses & Meier 2008) or LMJ (Casner et al. 2015), which offer many individual laser beams. We consider this a promising scenario to investigate in a future work. Complementarily, we will also better characterize turbulence in shock-shock encounters using large-scale 2D PIC simulations, where both longitudinal and transverse kinetic instability (e.g., Bunemann instability and Weibel instability) can be induced, which may lead to the generation of turbulence (Yao et al. 2018).

## 5. Conclusions

In our experimental campaign presented in this paper, we investigated the interpenetration of two laser-driven collisionless subcritical shocks, relevant for astrophysical phenomena such as the interplanetary medium and the local interstellar medium. The data obtained from their characterization has been used to feed MHD and PIC simulations, respectively, run with the FLASH and the SMILEI codes. While the MHD simulations were used to describe the collisional part of the problem, that is, the piston, and the overall evolution of the system; PIC simulations have provided insights into the microphysics at play in such a scenario. We compared the cases of single and double shocks and observed an acceleration of the background ions which was up to 1.5 times higher in the case of double shocks. Such acceleration is initially due to the electrostatic field $E_x$ associated with the shock front (injection), followed by a contribution mostly due to the “surfing” effect on the inductive electric field ($E_y \sim v_x B_z$). The presence of a second shock benefits this second mechanism as it allows for the existence of zones with $E_y$ directed as the proton $v_y$, and thus enhance their acceleration. Unfortunately, we could not measure such a spectrum during our experimental campaign, as the previously formed supercritical shocks produced higher energy protons whose spectrum covered the one of the particles accelerated by the further stage of double subcritical shock interaction.

In spite of the fact that most astrophysical collisionless shocks are supercritical, these results can shed light on the less frequently investigated subcritical shocks which are still relevant in various space phenomena. In particular, we have shown that the interaction of two subcritical shocks could lead to a higher energization of the background protons. This is relevant information when determining the distribution of high-energy particles that populate the interplanetary space and the very local interstellar medium surrounding the heliopause where colliding subcritical shocks are present. Moreover, high-power laser-plasma experiments have been demonstrated to be an essential tool that facilitates a laboratory recreation of scaled astrophysical phenomena, whose characterization is crucial for initializing the relative numerical simulations.

**Acknowledgements.** The authors would like to thank the teams of the LULI2000 (France) laser facility for their expert support, as well as the Dresden High Magnetic Field Laboratory at Helmholtz-Zentrum-Dresden-Rossendorf for the development of the pulsed power generator. We thank the SMILEI development team for technical support; we thank P. Loiseau (CEA-France) for the Thomson scattering analysis code. This work was supported by funding from the European Research Council (ERC) under the European Unions Horizon 2020 research and innovation program (Grant Agreement No. 787539). S.O. and M.M. acknowledge financial contribution from the PRIN INAF 2019 grant “From massive stars to supernovae and supernova remnants: driving mass, energy and cosmic rays in our Galaxy” and the INAF mainstream program “Understanding particle acceleration in galactic sources in the CTA era”. The computational resources of this work were supported by GENCI. The FLASH software used in this work was developed in part by the DOE NNSA ASC – and DOE Office of Science ASCR-supported FLASH Center for Computational Science at the University of Chicago. Part of the experimental system is covered by a patent (1000183285, 2013, INPI-France). The work of JIHT RAS team was supported by The Ministry of Science and Higher Education of the Russian Federation (Agreement with Joint Institute for High Temperatures RAS No 075 15 2020 785).

## References

Albertazzi, B., Béard, J., Ciardi, A., et al. 2013, Rev. Sci. Instrum., 84, 043505

Anderson, B. R., Skoug, R. M., Steinberg, J. T., & McComas, D. J. 2012, J. Geophys. Res. (Space Phys.), 117, A04107

Balikhin, M. A., & Wilkinson, W. P. 1996, Geophys. Res. Lett., 23, 1063

Balogh, A., & Treumann, R. A. 2013, Physics of Collisionless Shocks: space Plasma Shock Waves (New York: Springer)

Bemporad, A., & Mancuso, S. 2011, ApJ, 739, L64

Bemporad, A., & Mancuso, S. 2013, J. Adv. Res., 4, 287

Berčič, L., Maksimović, M., Landi, S., & Matteini, L. 2019, MNRAS, 486, 3404

Boldyrev, S., & Horaites, K. 2019, MNRAS, 489, 3412

Braginskii, S. I. 1965, Rev. Plasma Phys., 1, 205

Burgess, D., & Scholer, M. 2007, Phys. Plasmas, 14, 012108

Burlaga, L. F., Ness, N. F.,urn:https://doi.org/10.1007/s10043-016-0349-6

Cargill, P., Goodrich, C., & Papadopoulos, K. 1986, Phys. Rev. Lett., 56, 1988

Casner, A., Caillaud, T., Darbon, S., et al. 2015, High Energy Density Phys., 17, 2

Colburn, D. S., & Sonett, C. P. 1966, Space Sci. Rev., 5, 439

Derouillat, J., Beck, A., Pérez, F., et al. 2018, Comput. Phys. Commun., 222, 351

Dudkin, G. N., Nechaev, B. A., & Padalko, V. N. 2000, Plasma Phys. Rep., 26, 502

Edmiston, J. P., & Kennel, C. F. 1984, J. Plasma Phys., 32, 429

Faenov, A. Y., Pikuz, S. A., Erko, A. I., et al. 1994, Phys. Scr., 50, 333

Filippov, E. D., Makarov, S. S., Burdonov, K. F., et al. 2021, Sci. Rep., 11, 8180

Fuiza, F., Swadling, G. F., Grassi, A., et al. 2020, Nat. Phys., 16, 916

Froula, D. H., Luhmann, N. C., Jr, Sheffield, J., & Glenzer, S. H. 2011, Plasma Scattering of Electromagnetic Radiation: Theory and Measurement Techniques (Elsevier)

Fryxell, B., Olson, K., Ricker, P., et al. 2000, ApJS, 131, 273

Gosling, J. T., Bame, S. J., Feldman, W. C., et al. 1984, J. Geophys. Res., 89, 5409

Gurnett, D. A., Kurth, W. S., Burlaga, L. F., & Ness, N. F. 2013, Science, 341, 1489

Gurnett, D. A., Kurth, W. S., Stone, E. C., et al. 2015, ApJ, 809, 121

Gurnett, D. A., Kurth, W. S., Stone, E. C., et al. 2021, AJ, 161, 11

Hietala, H., Agueda, N., Andréeová, K., et al. 2011, J. Geophys. Res.: Space Phys., 116, A10

Horaites, K., Boldyrev, S., & Medvedev, M. V. 2019, MNRAS, 484, 2474

Ilie, R., Daldorff, L. K. S., Liemohn, M. W., Toth, G., & Chan, A. A. 2017, J. Geophys. Res.: Space Phys., 122, 5391

Kajdič, P., Alexandrova, O., Maksimovic, M., Lacombe, C., & Fazakerley, A. 2016, ApJ, 833, 172

Khiar, B., Revet, G., Ciardi, A., et al. 2019, Phys. Rev. Lett., 123, 205001

Lee, L. C., Mandt, M. E., & Wu, C. S. 1987, J. Geophys. Res., 92, 13438

Li, C. K., Tikhonchuk, V. T., Moreno, Q., et al. 2019, Phys. Rev. Lett., 123, 055002

Liu, Y. D., Richardson, J. D., Wang, C., & Luhmann, J. G. 2014, ApJ, 788, L28

MacFarlane, J. J., Golovkin, I. E., Woodruff, P. R., et al. 2003, Proceedings of Inertial Fusion and Science Applications (La Grange Park, IL: American Nuclear Society), 1

Maksimovic, M., Zouganelis, I., Chaufray, J. Y., et al. 2005, J. Geophys. Res. (Space Phys.), 110, A09104

Marcowith, A., & e. a., 2016, Rep. Prog. Phys., 79, 046901

Matsakos, T., Uribe, A., & Königl, A. 2015, A&A, 578, A6
Mellott, M. 1984, Adv. Space Res., 4, 245

Mellott, M. M., & Greenstadt, E. W. 1984, J. Geophys. Res., 89, 2151

Moreno, Q., Dieckmann, M. E., Ribeyre, X., & D'Humières, E. 2019, Plasma Res. Express, 1, 035001

Moses, E. I., & Meier, W. R. 2008, IEEE Trans. Plasma Sci., 36, 802

Mostafavi, P., & Zank, G. P. 2018, ApJ, 854, L15

Ocker, S. K., Cordes, J. M., Chatterjee, S., et al. 2021, Nat. Astron., 5, 761

Ohsawa, Y., & Sakai, J.-I. 1985, Geophys. Res. Lett., 12, 617

Pierrard, V., Maksimovic, M., & Lemaire, J. 2001, Astrophys. Space Sci., 277, 195

Pierrard, V., Lazar, M., & Schlickeiser, R. 2011, Sol. Phys., 269, 421

Pierrard, V., Lazar, M., Poedts, S., et al. 2016, Sol. Phys., 291, 2165

Pilipp, W. G., Miggenrieder, H., Montgomery, M. D., et al. 1987, J. Geophys. Res. Space Phys., 92, 1075

Seough, J., Nariyuki, Y., Yoon, P. H., & Saito, S. 2015, ApJ, 811, L7

Štverák, Š., Maksimovic, M., Trávníček, P. M., et al. 2009, J. Geophys. Res. Space Phys., 114, A5

Tremblin, P., & Chiang, E. 2013, MNRAS, 428, 2565

Treumann, R. A. 2009, A&ARv, 17, 409

Whang, Y., & Burlaga, L. 1985, J. Geophys. Res. Space Phys., 90, 221

Yang, Z., Lembège, B., & Lu, Q. 2012, J. Geophys. Res. Space Phys., 117, A7

Yao, W., Qiao, B., Xu, Z., et al. 2018, J. Phys., 20, 053060

Yao, W., Fazzini, A., Chen, S., et al. 2021, Nat. Phys., 17, 1177

Yao, W., Fazzini, A., Chen, S. N., et al. 2022, Matter. Radiat. Extremes, 7, 014402
## Appendix A: Limits of validity of the 1D simulations

<figure><img src="image_14.png" /><figcaption>**Fig. A.1.** Estimation of the maximum distance that a particle can travel in the y-direction before escaping the shock, considering the hemispherical geometry of the shock front. If we consider a radius $R \approx 2.5$ mm and a shock thickness $\delta R \approx 0.2$ mm, we have $\Delta y \approx 0.73$ mm.</figcaption></figure>

In our 1D particle-in-cell (PIC) simulations we have not taken the multidimensional effects into consideration, for instance, the shock front non-stationarity (Burgess & Scholer 2007), which might affect the proton dynamics due to the

ripping along the shock front (Yang et al. 2012). However, we verified in our previous work that the proton acceleration mechanism in the single shock case is not affected by the non-stationarity in the early few ns (Yao et al. 2022).

Moreover, the 1D geometry approximation has the limit of not considering the hemispherical profile of the shock. This can become a problem if the particles have gone too far in the y-direction and exit the shock. To quantify the maximum length that the protons can travel, we consider the radius and the thickness of the shock right before the collision: from the experimental characterization obtained via interferometry, we estimate a shock front radius $R \approx 2.5$ mm and a shock thickness $\delta R \approx 0.2$ mm.

As shown in Fig. A.1, we approximate the maximum length $\Delta y$ for which the particle can still be considered inside the shock as the distance that a particle moving only along the y-axis would travel between the middle of the shock and the external edge of it. This length would result in $\Delta y = \sqrt{R \delta R + \frac{3}{4} \delta R^2} \approx 0.73$ mm.

This leads to an imposition of such a limit on the background protons that have been accelerated by the double shock. Among the 20 protons with energies of 10 keV and higher, only 8 (= 40%) manage to remain confined in the shock all the time. This means that out of the total 1706 tracked ambient protons, only 8 (= 0.47%) reach at least 10 keV while remaining confined in the shock. Hence, even considering this limitation, the percentage of protons that are energized to 10 keV and more is considerably higher in the presence of a double shock structure than with a single one (1 out of 5120, i.e., ≈ 0.02%).