Detailed characterization of laboratory magnetized super-critical collisionless shock and of the associated proton energization

Collisionless shocks are ubiquitous in the Universe and are held responsible for the production of non-thermal particles and high-energy radiation. In the absence of particle collisions in the system, theoretical works show that the interaction of an expanding plasma with a pre-existing electromagnetic structure (as in our case) is able to induce energy dissipation and allow for shock formation. Shock formation can alternatively take place when two plasmas interact, through microscopic instabilities inducing electromagnetic fields which are able in turn to mediate energy dissipation and shock formation. Using our platform where we couple a fast-expanding plasma induced by high-power lasers (JLF/Titan at LLNL and LULI2000) with high-strength magnetic fields, we have investigated the generation of magnetized collisionless shock and the associated particle energization. We have characterized the shock to be collisionless and super-critical. We report here on measurements of the plasma density, temperature, the electromagnetic field structures, and particle energization in the experiments, under various conditions of ambient plasma and B-field. We have also modeled the formation of the shocks using macroscopic hydrodynamic simulations and the associated particle acceleration using kinetic particle-in-cell simulations. As a companion paper of \citet{yao2020laboratory}, here we show additional results of the experiments and simulations, providing more information to reproduce them and demonstrating the robustness of our interpreted proton energization mechanism to be shock surfing acceleration.


I. INTRODUCTION
The acceleration of energetic charged particles by collisionless magnetized shock is a ubiquitous phenomenon in astrophysical environments, among which the most energetic particles are the ultra-high-energy cosmic rays (UHECRs) accelerated in the interstellar medium (ISM) 2,3 . In this case, the source of collisionless dissipation is self-generated electromagnetic fields, resulting from kinetic instabilities such as the Weibel one. Besides, particles are also accelerated in our solar system due to collisionless magnetized shocks developed by the interaction of the solar wind with planetary magnetospheres 4,5 and, at larger distances, with the ISM 6 . In that case, the source of collisionless dissipation is the preexisting global electromagnetic structure. This will be the a) Electronic mail: yao.weipeng@polytechnique.edu case for the experiment detailed here, where we apply a global strong magnetic field onto a laser-ablated fast plasma. Since these shocks usually have their Magnetosonic Mach number M ms = v sh /v ms 2.7 (where v sh is the shock velocity, v ms = C 2 s + v 2 A is the Magnetosonic velocity, C s and v A are the ion sound velocity and Alfvénic velocity, respectively), they belong to the so-called super-critical regime 7,8 , which means the shock is not maintained by classical dissipation means alone. In order to help maintain a shock, the additional channel to expel energy is achieved by reflecting particles back upstream 9 .
A variety of acceleration mechanisms have been evoked as a way to transfer the energy from the shock waves to the particles, including shock surfing acceleration (SSA), shock drift acceleration (SDA), and diffusive shock acceleration (DSA). DSA requires high initial energy before further acceleration 10 , thus raising the so-called "injection problem" 11 ; while SSA and SDA are believed to be responsible for generating the pre-accelerated seed particles, i.e. for the initial accelerating process from thermal energies. Although it is still under debate whether SSA or SDA dominates the pre-acceleration process in various collisionless shock environments [12][13][14] , we can distinguish them by the following two aspects: On the one hand, in SSA, charged particles first get reflected at the shock front (due to the cross-shock potential electric field), then they surf along the shock front against the convective electric field (E = −v × B), and thus they gain energy. While in SDA, charged particles drift (due to the magnetic field gradient at the shock front) along the convective electric field and then gain energy 15 . On the other hand, SSA requires a thin shock width, compared to the Larmor radius of the charged particles, while SDA needs the opposite (so that the charged particles can gyrate and drift within the shock layer) 10,16 .
However, because of the immense spatial scales involved with collisionless phenomenon (e.g. the mean-free-path is λ m f p ∼ 1 AU in the Solar system), only a very small sampling of the shock formation and dissipation mechanisms can be realized. As a result, we still do not have a full understanding of the formation and evolution of collisionless shocks, and the question of the effectiveness and relative importance of SDA and SSA is still largely debated in the literature 17 . To further our understanding, laboratory experiments (and their simulations) have been proven to be an effective tool, providing highly-resolved, reproducible and controllable multi-dimensional datasets that can complement astrophysical observations 18,19 . Below, we will now briefly review the investigation of collisionless shocks via laboratory experiments.
The route that has been up to now most explored in the laboratory is to produce a shock (mediated by the Weibel filamentation instability) by colliding two ablative, unmagnetized flows driven by high-energy nanosecond lasers. This setup has yielded promising results at the Omega Laser Facility 20-22 and the National Ignition Facility (NIF) 23,24 , as well as at many other laser facilities all over the world [25][26][27] . Recently, experiments on collisionless shocks in plasma flows in which there was significant self-generated magnetic field showed, for the first time, the formation of magnetized collisionless shock, with the generation of Weibel instability and observation of electron acceleration in the turbulent structure 28 . Most recently, the dynamics of the ion Weibel instability has been characterized by local, quantitative measurements of ion current filamentation and magnetic field amplification in interpenetrating plasmas via optical Thomson scattering (TS) 29 . What's more, the generation of sub-relativistic shocks, together with relativistic electron acceleration, has been demonstrated to be within the reach of larger-scale, NIF-class laser systems 30 .
Another setup relies on a plasma expanding into a preformed ambient magnetized secondary plasma. Thanks to the magnetisation, the target ions create a collisionless magnetic piston that accelerates the ambient plasma to super-Alfvénic velocity, thus creating a high-Mach number shock with velocity of the order of 1000 km/s [31][32][33][34] . Recently, Schaeffer et al. have been able to make significant progress in characterizing the formation of collisionless shocks in terms of ion and elec-tron density and temperature, as well as electric and magnetic field strengths as a function of time at OMEGA 35 .
Besides, at the LULI laser facility at École Polytechnique (France), collisionless shock waves and ion-acoustic solitons have been investigated by proton radiography 36 . Moreover, significant electron pre-heating via lower-hybrid waves was also achieved in laboratory laser-produced shock experiments with strong magnetic field, providing a potential mechanism for the famous "injection" problem 37 . Additionally, at the VULCAN laser facility at the Rutherford Appleton Laboratory, the temporally and spatially resolved detection of the forming of a collisionless shock was achieved 38 .
In contrast to the above schemes, novel setups have been used with ultra-high-intensity lasers. For example, at the XingGuang III laser facility at the Laser Fusion Research Center in China, using a short (2 ps) intense (10 17 W/cm 2 ) laser pulse, an electrostatic (ES) collisionless shock, together with the filaments induced by ion-ion acoustic instability, could be observed via proton radiography 39 .
In our experimental campaigns at JLF/Titan and LULI2000 1 , we investigated shock formation combining laser-produced plasmas, a background medium and a strong ambient magnetic field (as detailed below). We chose to have an expanding plasma to drive a shock into an ambient gas in the presence of a strong external magnetic field. Contrary to Schaeffer et al. 33 , in our setup, the expanding plasma and the magnetic field were decoupled as the higher Z piston evacuates the magnetic field and was thus unmagnetized. This also allowed us to simultaneously have a highly magnetized ambient plasma (with homogeneous and steady magnetic field) and a high-β piston (β ≡ P thermal /P mag is the ratio of the thermal pressure to the magnetic one). Moreover, since our magnetic field strength was more than two times higher 40 , reaching 20 T comparing to the 8 T in Schaeffer et al. 33 , we were able to decouple more strongly the electrons from the ions 41 , and the shock was able to fully separate from the piston, which is crucial for its characterization 42 . As a result, we have been able to characterize the plasma density, temperature, as well as the E-field developed at the shock front, and more importantly, observe strong non-thermal accelerated proton populations for the first time.
In this paper, we will first show that laboratory experiments can be performed to generate and characterize globally mildly super-critical, quasi-perpendicular magnetized collisionless shocks in Section II, and detail their characteristics. Then, we will detail in Section III three-dimensional (3D) magneto-hydrodynamic (MHD) simulations reproducing the laser-driven piston generation and the following shock formation process. In Section IV, with the parameters characterized in the experiment, we will report the results of kinetic particlein-cell (PIC) simulations, which pinpoint that shock surfing acceleration (SSA) can be effective in energizing protons from the background plasma to hundred keV-level energies.

A. Experimental setup
The experiments were performed at the JLF/Titan (LLNL, USA) and LULI2000 (France) laser facilities with similar laser conditions but using complementary diagnostics, which was mostly linked with the availability of different auxiliary laser beams at each facility.
In the experiment at JLF/Titan, as is shown in Fig. 1, the collisionless shock was generated by sending a plasma, generated by having a high-power laser (1 µm wavelength, 1 ns duration, 70 J energy, and 1.6 × 10 13 W/cm 2 on-target intensity) irradiating a solid target (Teflon, CF 2 ), into a low-density (10 18 cm −3 ) H 2 ambient gas pulsed from a nozzle prior to the shot, and in the presence of a 20 T magnetic field that is homogeneous and steady-state at the time scale of the experiment. As shown in Fig. 1, the magnetic field created by a Helmholtz coil system 40,43 is oriented along the y-or z-axis.

B. Density characterization through optical interferometry
Using an interferometry setup 44 , the plasma electron density is recorded by optically probing the plasma (with a mJ, 1 ps auxiliary laser pulse). In Fig. 2, we present the overall electron density recorded in three different cases.
For the case with both ambient gas and B-field shown in Fig. 2 45,46 in the xy-and xz-plane, respectively. (e) Cases with only ambient gas but without B-field in the xy-plane (the xz-plane will be the same). Each image corresponds to a different laser shot, while the color scale shown at the top applies to all images. (f) The lineouts along the thin dark lines shown in each image. The laser comes from the right side and the piston source target is located at the left (at x = 0). Yellow arrows indicate the piston front, while green arrows indicate the shock front. of a hot plasma (the piston) that propagates along the x-axis and the collisionless shock is formed as a consequence of the plasma piston propagating in the magnetized ambient gas 33 . We can clearly see both the piston front and the shock front (indicated by the orange and green arrows, respectively), and indeed they are well detached from each other, enabling us to characterize them separately.
A lineout of the plasma density is shown in Fig. 2 (f), where the piston and shock fronts are also well identified by the abrupt density changes. The piston front is steepened by the compression of the magnetic field (see also below). Besides, we can clearly see a "foot" structure ahead of the shock front in the upstream (US) region for the case with both ambient gas and B-field, indicating the formation of the magnetized shock 47 .
In contrast, for the case with only B-field but without ambient gas 45 shown in Fig. 2 (c) and (d), due to the lack of ambient gas, no collisionless shock is formed ahead of the piston. For the case with only ambient gas but without B-field in Fig. 2 (e), no shock is formed as well in the ambient gas. From the corresponding lineout in (f), it is clear that only a smooth plasma expansion into the ambient (the green dashed line) can then be seen.

C. Piston compression characterization through X-ray spectroscopy
To further characterize the piston, the x-ray ion emission of Fluorine compressed within the expanding piston was measured by a Focusing Spectrometer with high Spatial Resolution (FSSR) 48 at both laser facilities. It was based on a spherically-bent mica (2d = 19.9376 Å) crystal with a curvature radius of R = 150 mm. Spatial resolution of 100 µm per pixel was achieved along the plasma expansion. Image Plate (Fujifilm TR BAS) was used as a fluorescent detector. The implemented scheme resulted in 13-16 Å spectral range with a high resolution (λ /dλ is higher than 1000). It covers spectral lines of Fluorine: resonance H-like (2p-1s transition) and Helike (3p-1s, 4p-1s, 5p-1s etc.) transitions as well as dielectronic satellites to Ly α . The diagnostic allowed us to measure electron density and temperature profiles of the piston expansion using a quasi-stationary approach 49 . The method is based on analysing the relative intensities of spectral lines of the same charge state and also takes into account the recombining plasma with a "frozen" ion charge. Figure 3 (a) shows that obviously the piston encounters stronger hindrance in the case with both ambient gas (H 2 ) and B-field (B z ) (see the green diamonds), comparing with other cases (i.e. the case with only B z in red dots and the case with only H 2 in blue triangles). We also see in Fig. 3 (b) that the electron temperature in the case of B z + H 2 becomes the highest at the piston front (between 4 and 7 mm), comparing with other cases. In addition, at the position of 4.5 mm, the evaluated electron density for the case of B z + H 2 is around 2 − 3 × 10 18 cm −3 and the electron temperature is about 65 eV, which are well-reproduced by our FLASH simulations, see Fig. 8 (a) and (b).

D. Electric field characterization through proton radiography
The single shock front was also probed with protons in order to measure the local electric field. The probing protons (accelerated by the Target Normal Sheath Acceleration process 51 from an auxiliary target and using the short-pulse arm of Titan) was sent parallel to the B-field, i.e. along the z-axis, as is shown in Fig. 1 (a).
As shown in Fig. 4 (a), we could clearly observe the same structures of the piston front and the shock front, consistent with those observed via optical probe, as shown in Fig. 2. By analysing the proton deflection structure, we could infer that we had a bipolar electric field at the shock front, with a total width of 0.4 mm and an amplitude of around 1 MV/m, as shown in Fig. 4 (b). To get this result, we imposed a certain 3D electric field map and simulated the proton dose that we would get on a detector. The electric field had been modulated in order to obtain a simulated dose (blue line) matching as much as possible the measured one (green line). Note that the proton deflection structure is accumulated along the z-direction. We will compare it with the particle-in-cell simulation results and discuss it in detail in Sec. IV.
Moreover, we compared the position of the shock structures seen in the electron density (via interferometry) with that in the electric field (via proton radiography) for the case with both external B-field and ambient gas. For the former, we have considered the point where the electron density had a sharp jump, as shown in Fig. 2f; as for the latter, we have taken into account the external edges of the proton dose accumulation. As is shown in Fig. 5, the evolution of the piston front and the shock front through both diagnostics are illustrated together (see legends for details), and it is clear that the results are quite consistent with each other. Note that when the target was not clearly visible in the radiography, i.e. for the series of points around 5 ns, we made use of the interferometry results to shift all the points of the right amount, while the distances between the piston and the shock fronts were kept constant. The original RCFs for the data points at various times are also shown.
FIG. 5. Piston and shock front position over time on the electron density (via interferometry) and on the electric field (via proton radiography). Images of proton radiography doses at different times are also shown, with dashed lines for piston front (in orange) and shock front (in green).

E. Temperature characterization through Thomson scattering
With a second high-energy auxiliary (0.5 µm wavelength, 1 ns, 15 J, focused over ∼ 40 µm along the z-axis and propagated throughout the plasma) available at LULI2000, we are able to characterize the plasma temperature by performing Thomson scattering (TS) off the electron and ion waves in the plasma (used in a collective mode 52 and analyzed by different spectrometers). Figure 6 shows the TS measurements in the region downstream (DS) compared to the shock front for cases with and without the external B-field. By comparing the experimental data profiles with the theoretical equation of the scattered spectrum for coherent TS in unmagnetized and non- collisional plasmas, with the instrumental function taken into an account, we are able to retrieve the local electron number density, as well as the electron and ion temperatures 53 . For the case without the B-field (i.e. with only ambient gas), both TSe and TSi give n e ∼ 1.5 × 10 18 cm −3 and T e ∼ 80 eV, and TSi also gives T i ∼ 40 eV in the DS region, as can be seen in Fig. 6 (a) and (c). However, for the case with B = 20 T, we see strong compression and heating in the DS region, indicated by the higher density and temperatures, i.e. n e ∼ 2.5 × 10 18 cm −3 , T e ∼ 115 eV, and T i ∼ 200 eV, as can be seen in Fig. 6 (b) and (d). With the characteristic feature of T i > T e , the effective formation of a shock can be inferred.

F. Evidence for proton energization
For the observation of the non-thermal proton spectrum, we use a standard magnetic spectrometer, with permanent magnets of 0.5 T strength. It was located close to the target (17.5 cm away) in order to maximize its collection efficiency, and it had its main axis along z, the main of the external B-field. That spectrometer has been calibrated precisely with a Hall probe and on many previous campaigns using filters to verify its energy dispersion. The protons are deflected by the magnetic field inside the spectrometer and landed after a short drift space onto Imaging plates (of TR type), the detector used here. These detectors are absolutely calibrated 54 .
The recorded proton spectrum is shown in Fig. 7   dots and blue error bars. Comparing it to the analytical thermal proton spectra (200 eV in red dash-dotted lines, as is observed in 1 through TS), it is clear that the proton energization is non-thermal. The cutoff energy reaches to about 80 keV. Also note that there is no signal recorded above the noise baseline for the case with only the B-field or the ambient gas.

III. MHD SIMULATIONS WITH FLASH
We use the 3D MHD code FLASH 55 to study the dynamics of the plasma plume expansion and shock formation in the ambient gas with the strong magnetic field, using the same parameters as the JLF/Titan experiment. The simulations are initialized in 3D geometry, using three temperatures (two for the plasma, and one for the radiation) with the equation-ofstate of Kemp and Meyer-ter Vehn 56 and radiative transport, in the frame of ideal MHD and including the Biermann battery mechanism of magnetic field self-generation in plasmas 57 . Specifically, the laser beam is normal to a Teflon target foil and has an on-target intensity of 10 13 W/cm 2 ; the generated plasma plume expands in the hydrogen gas-jet having an uniform density of 10 18 cm −3 . Moreover, the plasma plume expands in the uniform external magnetic field of 20 Tesla (aligned along the z-axis, as in the experiment). Figure 8 shows the FLASH simulation results, i.e. the electron density, electron temperature and ion temperature from FLASH at t = 2 ns (after the laser irradiation), in two different cases (the upper row is for the case with only ambient gas but without B-field, while the lower row is for the case with both the ambient gas and the B-field). We can observe that the structures of both the hydrodynamic piston and the induced shock, which propagates inside the ambient, are qualitatively reproduced compared to the experiment. The Teflon expanding piston produces a forward shock in the ambient (around x = 1.4 mm), as well as a reverse shock inside the Teflon piston (around x = 0.8 mm). The electron density is ∼ 1.6 × 10 18 cm −3 in the forward shock in the gas and increases up to ∼ 5 × 10 19 cm −3 in the reverse shock. The electron temperatures are between 60 to 70 eV in the forward and reverse shocks. Both correspond quite well to what is measured in the experiment (see the FSSR measurements in Fig. 3 and the TS measurements in 1 ). The ion temperature is 15 eV in the forward shock and between 80 eV and 180 eV inside the reverse shock.
Concerning the electron temperature, the FLASH simulation results are two times lower compared to the TS measurements in the DS region shown in Fig. 6; while for the ion temperature, the situation is worse as it is ten times less in the forward shock compared to the TS measurements. Also note that we have not seen the foot structure ahead of the shock in the FLASH simulations. Such discrepancies between the MHD simulations and the experiments show the difficulties to reproduce the shock condition in our case. This points to the fact that the shock evolution is dominated by kinetic effects. This is why we have resorted to using PIC simulations, the initial conditions of which are taken from the experimental measurements. Nevertheless, we can still observe that the FLASH simulations reproduce well the dynamics of the piston that induces the shock.
Since FLASH has the ability to model magnetic field generation through the Biermann battery effect, it allows us to assess the importance of this effect in the present configuration. Biermann battery generation of magnetic field is typically important only close to the target surface (order of 1 mm), and it is localized over the steep temperature gradients generated by the laser beam and rapidly decays once the laser beam is off (see for example [58][59][60]. As the shock is induced by the piston in the ambient gas 1 mm away from the target surface after the laser is off (∼ 2 ns), as shown in Fig. 8, the Biermann battery effect is negligible, compared to that of the strong externally applied B-field.

IV. KINETIC SIMULATIONS WITH SMILEI
The proton energization via the collisionless shock is modelled with the kinetic PIC code SMILEI 61 . During the interaction between the shock front and the ambient plasma, as the scale across the shock (∼mm) is much larger than that along the shock (∼ µm), we can treat this quasi one-dimensional (1D) interaction via the 1D3V version of the code.
As is shown in Fig. 9, the ambient plasma lies in the right half of the simulation box, while the left half is for the shocked plasma, flowing towards the right with an initial velocity of v 1 = 1500 km/s. Both of them consist of electrons and protons, with the real mass ratio m p /m e = 1836. The simulation box size is L x = 2048d e = 11 mm, and the spatial resolution is d x = 0.2d e = 1.1 µm, in which d e = c/ω pe = 5.3 µm is the electron inertial length, and ω pe = (n e0 q 2 e /m e /ε 0 ) 1/2 = 5.6 × 10 13 s −1 is the electron plasma frequency. Here, c is the speed of light, n e0 = 1.0 × 10 18 cm −3 is the electron number density of the ambient plasma, and m e , q e and ε 0 are the electron mass, elementary charge, and the permittivity of free space, respectively. Note that the shock width is initialized to be equal to the ion inertial length d i = 200 µm. The magnetic field is homogeneously applied in the z-direction with B z = 20 T (ω ce /ω pe = 0.06, where ω ce = q e B/m e ). The simulation lasts for 1.5×10 5 ω −1 pe ∼ 2.5 ns. Inside each cell, we put 1024 particles for each species. From the perspective of the ion Larmor motion, the simulation size is more than 10 r Li , in which For the shocked plasma, the electron number density is n e1 = 2n e0 = 2.0 × 10 18 cm −3 , and the temperature is T e1 = 100 eV and T i1 = 200 eV, all inferred from the TS characterization 1 . The boundary conditions for both particles and fields are open, and enough room is left between the boundary and the shock, so that the boundary conditions do not affect the concerned physics. Given the initial low temperature of the ambient plasma in the simulation (T e0 = 50 eV), the Debye length is small compared to the grid resolution d x , i.e. λ De = (ε 0 kT e0 /n e0 q 2 e ) 1/2 ≈ 0.01d e = 0.05d x . However, we do run a series of simulations with different initial temperatures, showing that the energy conservation for those cases is limited around 0.05% and the physical results are almost the same. The mean-free-path of the presented case is λ m f p ≈ 1800d e , which is larger than the interaction scale, further confirming that the shock is collisionless.
We report in Fig. 10 the results of two PIC simulations, i.e. with and without B-fields. For the case with the applied B-field (on the left column), typical structures of a supercritical quasi-perpendicular collisionless shock can be seen 9 . For example, the overshoot in the DS region (on the left of the red dashed line), the ramps in the shock fronts (both the red dashed line and the cyan dotted line), and the foot in the upstream (US) region, as can be seen in Fig. 10 (a). This foot region is formed by the reflected protons at a distance within r L,i and modulated by the modified two-stream instability 62 . The proton density (n i ) in Fig. 10 (e) shows a compression ratio of n i,DS /n i,US ≈ 4, which agrees with the theoretical jump condition prediction 63 . This density profile, together with the transverse electric field E y (not shown here), also follows the distribution of the external applied B-field B z . The longitudinal electric field (E x ) in Fig. 10 (c) peaks right at the ramps, providing the electrostatic cross-shock potential to trap and reflect the protons, as can be seen in the phase-space distribution in Fig. 10 (g). Because the proton reflection is clearly due to the E x in our case, not the DS compressed B-field 10 , together with the fact that the ion Larmor radius (about 0.8 mm) is larger than the shock width (around 200 µm), the dominant particle acceleration mechanism is SSA, not SDA. Note that the cyan dotted line indicates one of the periodic shock reformation 9 . On the contrary, for the case without B-field (on the right column), the drifting plasma just penetrates through the ambient gas and no shock is formed, thus no proton energization can take place, which is in accordance with our experimental observation.
Note that in Fig. 10 (c), the PIC simulation gives a longitudinal electric field E x ∼ 5 × 10 8 V/m in the shock layer, which is two order-of-magnitude higher than the fitting of the proton radiography in Fig. 4 (b). This discrepancy may be due to two reasons: on the one hand, the bipolar electric field structure fitted in the proton radiography has a size of 0.4 mm, while the E x peaks in the PIC simulations are very sharp, with their width smaller than 0.02 mm. With a time-average of the PIC simulation over 0.2 ns, the E x profile around the shock front reaches a size of 0.4 mm, and its value drops down to 2 × 10 7 V/m. On the other hand, the PIC simulation represents the tip of the semi-sphere shaped expanding shock front at a single slice of z-direction, where the B-field is strictly perpendicular to the plasma flow and the shock is the strongest; however, the proton radiography covers the whole shock front with an integration along the z-direction. It includes all other plasma flow directions in the xz-plane, which are not perpendicular to the B-field and the corresponding shocks are weaker. Together with the above two aspects of reasons, it is understandable that the electric field fitted from the proton radiography shall be smaller than the PIC result.
Particle dynamics of a high-velocity shock (as well as the comparison with the low-velocity case) and of the subsequent shock surfing proton energization is detailed in our previous paper 1 , while here we focus on demonstrating the robustness of the SSA mechanism that is at play in our experiment via 2D simulations, taking the non-stationarity 64 into considera-tion. Due to the limitation of the computational resources, we reduce the 2D simulation scale to an acceptable level: the simulation box size are L x = 8 mm, L y = 0.8 mm, the simulation time t end = 0.7 ns, and the resolution is d x = 0.4d e . From Fig. 11 (a), we can clearly see that the transverse nonstationarity has already occurred, with 2D-stripes mostly positioned at/behind the shock layer; while for protons with kinetic energy above 30 keV, their trajectories show that they mainly appear at the shock front, travelling down the negative y-direction. Note that the convective electric field E = −v ×B is towards the positive y-direction, i.e. E y = v x B z ; and the drifting of the protons against the convective electric field serves as a distinctive feature that the dominate proton acceleration mechanism is SSA, not SDA 15 . Fig. 11 (b) shows the proton energy spectra at 0.7 ns of both the 1D and 2D cases, which are close to each other, and there is only a 2 eV difference in the highest energy cut, which can be caused by the numerical heating of the 2D case (with lower spatial resolution). Moreover, checking the energy evolution of the protons in the x − t diagram, overlaid on the transversely-averaged B-field map in the reference frame of the contact discontinuity (CD), it is clearly demonstrated that the accelerated proton is first reflected at (or, picked up by) the shock front in Fig. 11 (c), and then surfing along the shock front while keeping gaining energy in Fig. 11 (d). This is exactly the same picture as we have shown for the 1D simulations 1 , proving that the SSA is the dominating proton acceleration mechanism at play (even in the multi-dimensional case).
Nevertheless, the non-stationarity of the shock might further accelerate the proton at a later time, especially after the protons pass through the shock front and gyrate in the DS region. But unfortunately right now we do not have the computational resources to reveal that scenario. In short, our simulation shows that the non-stationarity does not prevent the protons being accelerated by SSA (reflecting and surfing), at least not at an early time.

V. CONCLUSIONS
In conclusion, we have shown that laboratory experiments can be performed to generate and characterize globally mildly super-critical, quasi-perpendicular magnetized collisionless shocks. More importantly, non-thermal proton spectra are observed for the first time, and the underlying acceleration mechanism is pinpointed to be SSA via kinetic simulations, which can remarkably reproduce the experimental proton spectra. Such laboratory studies for proton acceleration, as well as those for electrons reviewed above, can not only further our understanding of the shock formation and evolution by complementing spacecraft and remote sensing observations, but also help shed new light on solving the fundamental issue of injection for the UHECR production.
Our platform can be tuned in the future to perform a systematical study of collisionless shock with different B-field strength and orientation, enabling us to capture the transition of the magnetized collisionless shock from sub-critical regime to super-critical one, so that we can explore the triggering of  FIG. 11. 2D simulation results. (a) B-field maps at 0.7 ns, normalized to 20 T, with trajectories of protons (E k > 30 keV). Solid lines are for protons from the ambient plasma and dashed ones are for protons from the drifting plasma; blue squares are the starting position at 0.5 ns, while red dots are the ending position at 0.7 ns. (b) Energy spectra of both 1D and 2D simulation results at 0.7 ns. Red lines are for the 1D case (solid line for protons in the whole simulation box, dashed line for those which lie around the shock layer in the vicinity of 1.8 mm), while black lines are for the corresponding 2D one. (c) Trajectory of a proton reflected at the shock front in the x − t diagram, overlaid on the transversely-averaged B-field map in the reference frame of the contact discontinuity (the grey colorbar is for the B-field strength, while the colored one is for the proton kinetic energy). (d) Trajectories of two protons surfing along the shock front, also in the x − t diagram, overlaid on the transversely-averaged B-field map in the same reference frame.