Efficient generation of turbulent collisionless shocks in laser-ablated counter-streaming plasmas

Laser-ablated high-energy-density (HED) plasmas offer a promising route to study astrophysically relevant processes underlying collisionless shock formation, magnetic field amplification, and particle acceleration in the laboratory. Using large-scale, multi-dimensional particle-in-cell simulations, we explore the interpenetration of laser-ablated counter-streaming plasmas for realistic experimental flow profiles. We find that the shock formation and its structure are substantially different from those of more idealized and commonly considered uniform flows: shock formation can be up to 10 times faster due to the transition from small-angle scattering to magnetic reflection and the shock front develops strong corrugations at the ion gyroradius scale. These findings have important consequences for current experimental programs and open exciting prospects for studying the microphysics of turbulent collisionless shocks with currently available high-energy laser systems.


I. INTRODUCTION
Collisionless shocks are ubiquitous across a wide range of scales and astrophysical environments, from galaxy clusters [1] and supernova remnants [2] to the Earth's bow shock [3,4]. They are mediated by plasma instabilities, which dissipate the flow energy by heating the plasma, amplifying magnetic fields, and accelerating particles. The microphysics governing the nonlinear shock dynamics is complex; it cannot be directly resolved by astrophysical observations and is still not fully understood. In particular, how the structure of the shock front and the amplification of turbulent magnetic fields affect particle injection into diffusive shock acceleration remains a very important open question in high-Mach number (M 10) astrophysical shocks [5][6][7][8][9][10][11][12]. The advent of high-power lasers has opened a new avenue to produce and study high-Mach number collisionless shocks in the laboratory. The main goal is to produce shocks with self-generated magnetic turbulence, where the mechanisms of magnetic field amplification and particle injection can be directly studied in a controlled environment. In particular, there is a strong interest in studies of the Weibel instability [13,14], as both simulations [6,9,11,12,15,16] and spacecraft observations [17] suggest this instability is dominant in driving magnetic turbulence at shocks in a wide range of scenarios, from planetary shocks to young supernova remnants to gamma-ray bursts. Experimental studies can greatly complement spacecraft observations and play a critical role in benchmarking numerical simulations and theoretical models.
This prospect has led to a large experimental effort in high-energy-density (HED) laser facilities in the last decade. In these experiments, high-Mach number counter-streaming flows are produced from the ablation of a solid-density target by kilojoule, nanosecond laser pulses. Important progress has been achieved in the characterization of the plasma flow properties and interaction [18][19][20][21], leading to the demonstration of the magnetic field amplification by the Weibel instability [22][23][24][25], and to a better understanding of the interplay between collisional and collisionless effects [26,27].
The generation of Weibel-mediated shocks requires a large flow interpenetration distance that poses a significant experimental challenge. Theoretical models and kinetic simulations have been largely limited to uniform plasma flows and suggest that shock formation requires system sizes exceeding 10 3 ion skin depths [15,28,29], which are beyond those attained at current experimental facilities. However, very recently, experiments at the National Ignition Facility (NIF) have demonstrated the formation of high-Mach number shocks mediated by electromagnetic instabilities [30]. The shock formation time observed is significantly faster than previous models predict. It is thus crucial to understand the physics of shock formation in the laboratory experiments and in which conditions can current systems allow the study of the relevant microphysics of astrophysical turbulent shocks.
In this paper, we study the formation of collisionless shocks in laser-ablated counter-streaming plasmas using large-scale two-(2D) and three-dimensional (3D) particle-in-cell (PIC) simulations that include the intrinsic density and velocity inhomogeneity of experimental plasma flows. We find that realistic plasma profiles impact the interaction in two important ways, which have not been previously recognized: (i) shock formation is significantly faster than previously found for homogeneous plasmas because the coherence length of the amplified magnetic field becomes comparable to the ion gyroradius, leading to a transition from small-angle scattering to magnetic reflection that efficiently slows down the flows, and (ii) enhanced magnetic field advection towards the shock leads to ion-gyroradius-scale corrugations of the shock front. An analytical model for the shock formation time and minimum system size required to study collisionless shocks with laser-ablated plasma flows is introduced and shown to be in good agreement with both PIC simulations and recent NIF experiments. These results indicate that state-of-the-art HED facilities, such as the NIF and Laser Megajoule (LMJ), can produce collisionless shocks with strong turbulence driven by the Weibel instability and where the shock front uniformity can be controlled, opening the possibility to probe their influence on particle injection and test relevant astrophysical shock models.

II. PLASMA FLOW PROFILES AND SIMULATION SETUP
Plasma flows produced by laser-ablation of soliddensity targets have inhomogeneous density and velocity profiles given by a well-established self-similar solution [31], consistent with experimental measurements, such as those in Refs. [20,26]. In a typical configuration, two counter-facing solid targets separated by a distance 2L 0 are irradiated by kJ-class lasers to produce counterstreaming plasma flows that interact at the mid-plane between the two targets, defined at x = 0. The velocity and density of one of the flows is v(x, t) = c s +(x+L 0 )/(t+τ 0 ) and n(x, t) =ñ exp [−(x + L 0 )/(c s (t + τ 0 )) − 1], with τ 0 the time the flow takes to travel the distance L 0 to the mid-plane where the two flows first meet at t = 0, c s = (γZk B T e /m i ) 1/2 the sound speed, γ the adiabatic index, m i and Z the ion mass and charge number, T e the electron temperature, andñ the target surface density. The opposite, counter-propagating flow has the symmetric profile with respect to x = 0.
In order to study the importance of more realistic plasma profiles on shock formation we have performed 2D and 3D fully kinetic PIC simulations with the code OSIRIS 3.0 [32,33] for both homogeneous and inhomogeneous counter-streaming non-relativistic plasma flows. In the homogeneous case (identified by the index H), plasma flows are initialized with uniform velocity v H = v 0 = 0.11 c, a sonic Mach number M = v 0 /c s = 21, and density n H = n 0 . In the inhomogeneous case (identified by the index I), the flow velocity profile follows the selfsimilar theory, with maximum velocity v 0 = 0.11 c near the interaction region (mid-plane of the simulation) at t = 0 (Fig. 1a). The density profile also closely resembles the self-similar solution (see Fig. 1a and Appendix A), but has been modified such that n I ( This guarantees that the same plasma energy density is delivered to the interaction region in both cases and allows for a direct comparison of the shock formation efficiency between both. The baseline simulations have a longitudinal size of L x = 135 c/ω pi and a transverse size of L y = 50 c/ω pi , with ω pi = [4π(2n 0 )Z 2 e 2 /m i ] 1/2 the ion plasma frequency corresponding to the total plasma density. The simulations resolve the electron skin depth at the interaction region (c/ω pe ) with 8 cells, use a time step of 0.087 ω −1 pe , and use 36 particles per cell per species. The boundary conditions are periodic along the transverse y-direction and open along the x-axis for both fields and particles. A third-order interpolation scheme is used for improved numerical accuracy. The baseline ion mass to charge ratio used is m i /(m e Z) = 128. The high flow velocity (∼ 0.1 c) and reduced mass ratio used are typical choices in the PIC modeling of shocks [15,22,29,34] that capture the non-relativistic nature of the flows and a large separation between electron and ion physics, while balancing computational expense. We note that the physics of pure electromagnetic instabilities, such as the Weibel instability discussed in this work, can be rigorously scaled between non-relativistic systems with different flow velocities [35].
The evolution of the magnetic field and plasma density for both the homogeneous and inhomogeneous cases is shown in Fig. 1. As expected, the early time dynamics is dominated by the development of the ion Weibel instability and is similar between both cases, which is consistent with previous experiments and simulations [22,23]. The instability gives rise to filamentary currents with transverse wavelength comparable to the ion skin depth c/ω pi and to the exponential amplification of the magnetic field with growth rate Γ 0.07 ω pi for both cases, in good agreement with linear theory [13,14]. The instability saturates at τ sat 10 Γ −1 140 ω −1 pi , with B sat 0.014 m i ω pi c/e, which is in good agreement with predictions based on magnetic trapping B sat = v 0 m i ω pi /(2πZe) [36].
After saturation of the Weibel instability, the transverse coherence length of the magnetic field increases in both cases due to filament merging [34,[37][38][39]. However, the longitudinal extent of the filaments starts to be significantly more limited for the inhomogeneous flows, with the stronger magnetic fields being confined to the central region near the mid-plane (Fig. 1e,f). This is a consequence of different effects associated with the inhomogeneity of the flow profiles. The growth rate of the Weibel instability is reduced away from the mid-plane due to the density and velocity asymmetry and due to the ion heating of the flow that has crossed the interaction region (detailed calculations in Appendix B). Moreover, the asymmetric flow profiles lead to enhanced magnetic field advection towards the mid-plane region (Fig. 1i,j). The advection velocity is v adv = (n + v + +n − v − )/(n + +n − ) [40], where the index +(−) refers to the flow moving in the positive (negative) x-direction. (Note that v adv = 0 for homogeneous symmetric flows.) The weaker magnetic fields produced in the upstream region at a distance L from the mid-plane are thus advected to the central region and compressed on a time scale τ adv = L/v adv .
At later times, we observe dramatic differences between the two cases, both in terms of fields (Fig. 1g,h) and density (Fig. 1k,l). In the case of laser-ablated inhomogeneous plasmas, the flows are strongly compressed in the interaction region, two collisionless shocks are formed and propagate against the incoming plasma flows. Each shock reaches the density compression dictated by the Rankine-Hugoniot jump conditions [41]  285 ω −1 pi . k,l) Density compression ratio n2/n1 between the transversely averaged ion density of interacting flows n2 and single flow n1.
(for a 2D system) at a shock formation time τ sf 1030 ω −1 pi (Fig. 2a), where n 2 is the compressed (downstream) density and n 1 is the single flow (upstream) density. We have verified that in 3D the density compression reaches n 2 /n 1 4 (see Appendix D). In the homogeneous plasma case, the density compression observed is significantly weaker and slower. We have run a larger simulation with L x = 1200 c/ω pi up to 2500 ω −1 pi and the maximum compression obtained was just n 2 /n 1 2.5 (Fig. 2a).
The faster shock formation observed for laser-ablated inhomogenous plasma flows is due to the transition from small-angle scattering to magnetic reflection. In the case of homogeneous flows, the magnetic fields produced by the Weibel instability have a wavelength that is much smaller than the ion gyroradius of the flows, λ B r i FIG. 2. a) Temporal evolution of the density compression ratio n2/n1 at the simulation mid-plane (x = 0) for the inhomogeneous (blue dots) and homogeneous (orange dots) flow profiles. The black curve corresponds to the density evolution predicted from Eq. 2 for t > τc. b) Temporal evolution of the ion gyroradius ri (crosses) and dominant magnetic wavelength λB (dots) for the inhomogenous (blue) and homogeneous (orange) flow profiles. The black curve corresponds to the filament merging model for λB(t) from Ref. [34].
( Fig. 2b). Ions slow down via many small-angle scatterings with the magnetic fields [15,28] and kink-like instabilities of the current filaments [29], which operate on long time scales compared to the ion cyclotron frequency. A different dynamics is observed in the case of inhomogeneous flows where ions are efficiently confined near the mid-plane. This is a consequence of the temporal decrease of the flow velocity (and associated ion gyroradius) arriving at the interaction region, intrinsic to laserablated plasmas, . As the ion gyroradius becomes comparable to the dominant Weibel-mode wavelength, r i ∼ λ B (see Fig. 2b), the flows are effectively slowed down and heated by largeangle magnetic reflections.
As the shock is formed in the inhomogeneous flow interaction, strong corrugations of the magnetic field (and density) are observed at the shock front on the scale of the ion gyroradius. The plasma flows primarily in alternating filaments near the mid-plane region. Weaker magnetic fields from the upstream region are advected towards the center through these filaments and slowed down on a scale comparable to the ion gyroradius. This leads to the compression of the magnetic fields near the mid-plane and to strong anti-symmetric modulations at the shock front, as observed in Fig. 1g. Since the magnetic field is frozen-in the electron fluid, the ratio of electron thermal to fluid velocity, v th /v 0 , controls the shock front corrugation level. Figure 3 shows that indeed the corrugations are significantly less pronounced when v th v 0 . Experimentally, v th can be controlled by varying the laser intensity and Z (through radiative effects [42]) to probe the impact of shock front corrugations on magnetic field amplification and particle injection [4,7,10].

III. MODEL OF SHOCK FORMATION WITH LASER-ABLATED FLOWS
A model for shock formation in laser-ablated counterstreaming plasma flows can be derived by estimating the critical time τ c at which the magnetic wavelength of the Weibel mode becomes comparable to the ion gyroradius at the interaction region, r i (τ c ) ∼ λ B (τ c ). After saturation of the Weibel instability, the evolution of the dominant magnetic wavelength due to filament merging follows [34] pi is the typical merging time, and v 0 = L 0 /τ 0 is the initial flow velocity. The ion gyroradius is r i (t) = m i cv(t)/(ZeB sat ) = 2π(c/ω pi )L 0 /(L 0 + v 0 t). The critical time τ c is the solution of For the conditions of our simulation, this yields τ c 565 ω −1 pi , in good agreement with our results (Fig. 2b). We note that filament merging ceases at late times for both homogeneous and inhomogeneous flows (Fig. 2b) due to the development of the drift-kink instability [29]. However, in practice, this occurs at times comparable to or larger than τ c and thus the use of the filament merging model to estimate τ c is appropriate (see Appendix C).
The critical time τ c marks the transition from smallangle scattering to magnetic reflection. At this stage, the ions reaching the center are slowed down and confined over a region of extension ∼ r i . The density compression in the mid-plane region between τ c and the shock formation time τ sf can then be described as where we have assumed that at t = τ c the two flows overlap without yet significant compression, i.e. n 2 (τ c ) = 2n 1 (τ c ). The shock formation time τ sf is obtained from the density jump condition n 2 (τ sf ) = (2+δ)n 1 (τ sf ), where δ = 1 in 2D and δ = 2 in 3D. The shock formation times predicted for homogeneous flows τ sf−iso (black line) [28] and τ sf−kink (blue line) [29] largely exceed that observed for laser-ablated flows. The shaded areas indicate the range of parameters achievable at current HED laser facilities.
We have confirmed that Eq. 2 describes well the density compression observed in the simulations of laser-ablated plasma flows after τ c , as shown in Fig. 2a. The predicted shock formation time τ sf 1060 ω −1 pi is also in good agreement with the simulation results. In the limit where the initial single flow density varies weakly over an ion gyroperiod, the shock formation time can be approximated by with Ω c = v(t)/r i (t) the ion cyclotron frequency. For typical experimental flow velocities [19-24, 26, 27] v 0 ∼ 1000 − 2000 km/s c, the critical time τ c largely exceeds the time of flow compression by magnetic reflection, and dominates the shock formation time τ sf ∼ τ c . The minimum flow interpenetration distance required for shock formation is L sf ∼ τ sf v 0 ∼ τ c v 0 and thus corresponds to a minimum target separation 2L 0 ∼ 2τ c v 0 . By setting τ c = τ 0 and taking the appropriate limits τ c τ m and τ c τ sat we can calculate the minimum shock formation time as In practice, laser-ablated fully ionized plasmas have A/Z 2, for which the required target separation to reach shock formation is 2L 0 250 c/ω pi and is independent of the flow velocity v 0 . This analysis considered the longitudinal inhomogeneity of laser-ablated planar plasma flows. On time scales longer than L 0 /c s , 3D divergence effects can modify the plasma density profiles. It is thus important to guarantee that the shock formation time τ sf = L sf /v 0 < L 0 /c s or, equivalently, L 0 > L sf /M . This condition is always comfortably met for system sizes large enough to study high-Mach number shocks (L 0 ≥ L sf and M 1) and thus the profiles considered are reasonably justified. The transverse inhomogeneity of the plasma and obliquity of the flows away from their axes (due to divergence) can also affect plasma interpenetration. These effects should not impact the physics of shock formation and corrugation (both on the ion gyroradius scale) if L 0 r i ∼ 0.064/ n 0 [10 19 cm −3 ] cm. For the typical system size, L 0 ∼ 0.4 − 1.25 cm, and plasma density, n 0 ∼ 5 × 10 18 − 5 × 10 19 cm −3 , of current experiments [20,24,26,30], this condition is always met. We have performed additional 2D simulations with spherically divergent flows and 3D simulations with modified density profiles and confirmed that the shock formation time and structure predicted by our model are robustly observed in all cases (see Appendix D).

IV. COMPARISON OF SHOCK FORMATION MODELS WITH PIC SIMULATIONS AND EXPERIMENTS
Shock formation in counter-streaming laser-ablated plasmas is significantly faster than that predicted from previous models of Weibel-mediated shocks in homogeneous flows. The predicted shock formation time due to flow isotropization during filament merging is [28] 1867c/ω pi . Additionally, shock formation due to the disruption of the current filaments by kinktype instabilities predicts [29] τ sf−kink 5τ kink (5/3)(m i /Zm e ) 3/4 (c/v 0 )ω −1 pi , where τ kink is the growth time of the drift-kink instability, and corresponds to 2L 0 1572c/ω pi . The shock formation model proposed in this work is thus 6 − 8× faster than previous predictions. Figure 4 summarizes the comparison of our model predictions for the shock formation time and those of previous models [28,29] with the results of 2D and 3D PIC simulations and experimental measurements on NIF [30]. The PIC simulations performed over a range of flow velocities (v 0 = 0.03 − 0.1 c) and mass to charge ratios [m i /(m e Z) = 32 − 512] show good agreement with our model, confirming the scaling of τ sf ∝ 1/v 0 (shock formation distance independent of v 0 ) and the weak dependence on mass to charge ratio. Moreover, our model is in good agreement with the recent NIF experiments, which observed shock formation after 3 ns of interaction. For the measured plasma overlapping density 5×10 19 cm −3 and v 0 2000 km/s, the shock formation time predicted by our model is 2 × 10 4 ω −1 pi 3 ns, consistent with the experimental measurements.
The results presented here are important to guide the development and interpretation of experimental studies of turbulent collisionless shocks and in connecting them with astrophysical simulations and models. Figure 4 indicates the parameter space (shaded areas) that can be probed by the current HED facilities (see Appendix E for details). For kJ-class OMEGA experiments [19][20][21][22][23][24]26], with measured plasma densities n 0 ∼ 10 18 − 10 19 cm −3 , we predict τ sf 13 ns and 2L 0 2.5 cm, which largely exceed the interaction time of ∼ 5 ns and system size of ∼ 0.8 cm of previous experiments. However, our results indicate that shock formation can be comfortably reached at higher energy laser facilities (> 100 kJ), such as NIF and LMJ. The collisionless shocks will have a transition size dictated by the ion gyroradius and the magnetic turbulence in the shock foot is driven by the Weibel instability, which are critical properties of high-Mach number astrophysical shocks revealed by kinetic simulations [9] and spacecraft observations [17]. Moreover, by varying the electron temperature it will be possible to control and study the impact of shock front corrugations on magnetic field amplification and particle injection [4,7,10].

V. CONCLUSIONS
In summary, we have shown for the first time that the inhomogeneity of laser-ablated plasma flow profiles has a significant impact in accelerating the formation of turbulent collisionless shocks and in controlling their dynamics. These findings have important implications as they enable the use of current HED experimental facilities to study the complex interplay between magnetic field amplification, shock front corrugation, and particle injection underlying high-Mach number astrophysical shocks.
The The plasma used in our inhomogeneous simulations has a velocity profile v(x, t) = c s + (x + L 0 )/(t + τ 0 ) that follows the self-similar solution for laser-ablated plasmas [31]. The density profile has been slightly modified from the self-similar solution in order to guarantee that the flow energy density, nv 2 , delivered at the midplane region (where the flows interact) is constant in time nv 2 = n 0 v 2 0 . This is convenient because it allows for a direct comparison of the shock formation efficiency between inhomogeneous (laser-ablated) and homogeneous flows. We compute the initial density profile n(x, t = 0) by integrating the continuity equation for the density, with boundary condition n(x = 0, t) and assuming ballistic propagation of a density element (i.e. dv/dt = 0). This gives the density profile showed in Fig. 5 (blue line), which is a reasonable approximation of the self-similar profile (red dashed line). We have performed simulations with different density profiles and checked that the main results of this work are relatively insensitive to the exact details of the density profile and determined primarily by the inhomogeneity of the flow velocity, which is responsible for the transition from small-angle scattering to magnetic reflection in the shock formation.

Appendix B: Growth rate and saturation of the Weibel instability in inhomogeneous counter-streaming flows
The growth of the ion Weibel instability in inhomogeneous counter-streaming flows can be best described in a frame where the instability is purely growing (termed "Weibel frame"), defined as with n ± , v ± , T ± the local densities, velocities, and temperatures of the ion flows moving in the positive (+) and negative (-) direction. The corresponding dispersion relation, accounting for the ion dynamics in the nonrelativistic regime, reads pi is the ion plasma frequency at the local density, v th is the ion thermal velocity and the primed quantities are defined in the Weibel frame. The growth rate Γ = iω, obtained from the solution of Eq. S2, in the limit of cold flows and k ω * pi /c, can be expressed as FIG. 6. Comparison of the longitudinal spatial profile of the magnetic field measured in the inhomogenous flow simulation at t = 125 ω −1 pi (blue curve) with the prediction of the magnetic trapping mechanism (Eq. S4) for the growth rate obtained from Eq. S2 (black curve).
which is maximum at the center of interaction (x = 0), where the flows have equal density and opposite velocity, consistent with the simulation results. However, it is important to take into account temperature effects, since the heating of the flows as they cross the strong field region near the mid-plane will further decrease the growth rate away from it. This will also impact the saturation of the magnetic field, which, following the magnetic trapping mechanism [36], is given by We have compared the magnetic field profile obtained in the inhomogeneous simulation with that predicted by Eqs. S2 and S4, where we have solved Eq. S2 numerically for the ion temperature, density, and velocity extracted from the simulation. The results, shown in Fig. 6 near the saturation of the instability (t = 125 ω −1 pi ), indicate that the sharp decrease of the magnetic field amplitude away from the mid-plane in inhomogeneous flows can be reasonably well described by our model.

Appendix C: Interplay between filament merging and kink instability
Our PIC simulations (for both homogeneous and inhomogeneous flow profiles) show that right after the saturation of the Weibel instability the dynamics of the current filaments is primarily governed by filament merging, in good agreement with the model of Ref. [28] (see Fig. 2b). However, at late times (t 500 ω −1 pi in Fig. 2b), we observe that filament merging ceases and the dominant magnetic wavelength remains approximately constant. We have found that the stopping of the filament merging coincides with the development of longitudinal modulations of the filaments and the onset of magnetic turbulence, consistent with the development of kink-like instabilities described in Ref. [29]. We observe that this transition occurs at a time t 2τ kink (m i /Zm e ) 3/4 (c/v 0 ) ω −1 pi after the saturation of the Weibel instability (τ sat ), which is consistent across the range of simulations performed for different flow velocities and ion to electron mass ratios (Fig. 7). We have also verified that our estimate of the critical time τ c for the transition between small-angle scattering and magnetic reflection (Eq. 1) is in good agreement with the PIC simulation results (Fig. 7). We observe that for realistic ion to electron mass ratios the critical time occurs always before the saturation of filament merging (onset of the kink instability), thus validating the use of the merging model in the estimate of τ c .
Appendix D: Impact of flow divergence and 3D effects on shock formation Laser-produced plasma flows will have an associated divergence, i.e. they will have a radial velocity component away from the central axis (r = 0). The flow divergence at the interaction region can be controlled by the ratio of the laser focal spot size and the interaction distance L 0 . Earlier experiments used laser spot sizes significantly smaller than the system size, and the flow divergence was relatively large [20,23,26,27]. Very recently, experiments at the National Ignition Facility (NIF) [30] used a larger spot size comparable to the system size, and showed that the transverse size of the heated plasma region at the flow-flow interaction is comparable to the laser spot at the target, indicating that the flows are relatively well collimated, as considered in our simulations. In general, as discussed in the main text, provided that the system size is much larger that the ion gyroradius associated with the Weibel magnetic field, L 0 r i , the flow divergence is not expected to affect the shock formation physics and the use of collimated flows considered in the simulations is well justified.
In order to explore the impact of the plasma flow divergence, we performed a simulation with spherically divergent flows. The density and velocity profiles along the y = 0 axis are equal to the ones discussed in the manuscript and constant at a radial distance R from the target irradiated spot. The flow is transversely limited at ±30 • from the y = 0 axis to minimize the need of a very large transverse simulation box size. As can be seen in Fig. 8, the simulation with spherically expanding flows (right panel) shows a magnetic field amplitude and development of corrugations similar to the one presented in the main paper (left panel). This confirms that the dynamics is weakly affected by the flow divergence if the ion gyroradius is much smaller than the system size, as in the recent NIF experiments.
In addition, we have performed 3D simulations to verify that the shock physics and model prediction for the shock formation time is robust in 3D. We have used m i /(Zm e ) = 32 due to the very high computational cost. The velocity and density profiles are the same as in Fig.  1, with the difference that the density is kept constant after τ plateau = 700 ω −1 pi . This plateau in density profile is typically observed experimentally due to the flow divergence and/or finiteness of the laser duration [20,30]. In order to confirm that this does not impact significantly the shock formation physics, which is primarily determined by the flow velocity profile, we have simulated a case where τ sf > τ plateau . As can be seen in Fig. 9, shock formation is reached in 3D between t ∼ 900 − 1000 ω −1 pi with a density compression factor of 4, as expected from the shock jump conditions. The observed shock formation time is in very good agreement with our model prediction of τ sf 940 ω −1 pi .

Appendix E: Plasma conditions of current HED facilities
There has been a significant experimental effort in using high-energy-density laser facilities such as OMEGA and NIF to explore the physics of collisionless shocks in counter-streaming laser-ablated plasmas [19][20][21][22][23][24][25][26][27]30]. These experiments produce typical flow velocity v 0 ∼ 1000 − 2000 km/s, with the plasma density being primarily constrained by the laser energy and system size. Experiments at OMEGA [19][20][21][22][23][24]26] have been performed with up to 5kJ/target and have measured typical densities at the interaction region in the range 10 18 − 10 19 cm −3 . Recent experiments on the NIF used 250-450 kJ/target [27,30] and measured a density at the interaction in the range 10 19 − 10 20 cm −3 . While, to our knowledge, collisionless shock experiments on LMJ have not yet been reported, based on the results from NIF, at the currently available energy of 150 kJ/target it is reasonable to expect plasma densities in the range 5 × 10 18 − 5 × 10 19 cm −3 . The maximum plasma system sizes that can be produced at these conditions can be estimated from basic energy conservation arguments as where η is the energy conversion efficiency from the laser to the plasma flows and θ is the divergence angle of the flows. We have considered typical parameter ranges η ∼ 0.25 − 0.5 and θ ∼ 30 • − 45 • . The range of interaction times L 0 /v 0 obtained is shown in Fig. 3.
To study collisionless shocks, it is important to guarantee that Coulomb collisions will not affect the shock formation process. For collisions to have a negligible effect, one should guarantee that the ion-ion mean free path λ mfp = m 2 i v 4 r /(16πZ 4 e 4 n i log Λ), which is responsible for slowing down the flows, is much larger than the required shock formation distance L sf = 125 c/ω pi (as described in our model). Here n i = n 0 /Z is the ion density of each flow, v r = 2v 0 is the relative velocity between flows, and log Λ is the Coulomb logarithm. This condition can be written as where we have used A/Z = 2 and log Λ = 10, as typical of laser-ablated plasmas. We find that for the typical range of flow density and velocity values of experiments discussed above and Z = 6 (Carbon flows), we have λ mfp /L sf ∼ 6 − 280 and thus collisions are not expected to affect the shock formation process discussed in our work.