# Beam pointing of relativistic high-order harmonics generated on a nonuniform pre-plasma

Chaoneng Wu, Yiming Xu, Andre Kalouguine, Jaismeen Kaur, Antoine Cavagna, Zuoye Liu, Rodrigo Lopez-Martens, Cangtao Zhou, Philippe Zeitoun, Stefan Haessler, et al.

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

To cite this version:

Chaoneng Wu, Yiming Xu, Andre Kalouguine, Jaismeen Kaur, Antoine Cavagna, et al. Beam pointing of relativistic high-order harmonics generated on a nonuniform pre-plasma. Plasma Physics and Controlled Fusion, 2025, 67 (4), pp.045009. 10.1088/1361-6587/adb1c hal-04986897

HAL Id: hal-04986897

https://hal.science/hal-04986897v1

Submitted on 11 Mar 2025

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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Beam Pointing of Relativistic High-order Harmonics Generated on a Nonuniform Pre-plasma

CHAONENG WU<sup>1,2</sup>, YIMING XU<sup>1,3</sup>, ANDRE KALOUGUINE<sup>4</sup>, JAISMEEN KAUR<sup>4</sup>, ANTOINE CAVAGNA<sup>4</sup>, ZUOYE LIU<sup>3</sup>, RODRIGO LOPEZ-MARTENS<sup>4</sup>, CANGTAO ZHOU<sup>1,2</sup>, PHILIPPE ZEITOUN<sup>4</sup>, STEFAN HAESSLER<sup>4,*</sup>, LU LI<sup>1,*</sup>

The use of tunable pre-pulse is a common technique to enhance the high-order harmonic generation from surface plasma. The shape and dynamic of the electron density, the degree of ionization and its rate, and the plasma heating are influenced by the pre-pulse properties. A non-uniform pre-pulse can generate a spatially varying pre-plasma density, which offers the surface for spectral up-conversion and reflection. The geometrical characteristics and plasma dynamics under intense laser fields subsequently affect the properties of the emitted harmonics. In this study, the variation in harmonic beam pointing due to the electron density shape was investigated. Particle-in-cell simulations demonstrated that both plasma hydrodynamics and geometrical optical effect induce the deviation of harmonic beam from specular reflection. This research contributes to the understanding of the surface plasma dynamics during high harmonic generation process.

## 1. Introduction

The advancement of laser technology has facilitated the generation of relativistic electromagnetic fields capable of generating extreme ultraviolet (EUV) radiation via laser-driven plasma surfaces. Surface high-order harmonic generation (SHHG) has the advantage of having no upper driving-intensity limit [1] and offering high conversion efficiency [2, 3], thus potentially allowing to deliver powerful EUV attosecond pulses. Extensive research has been dedicated to elucidating the mechanisms of SHHG [3-8], including coherent wake emission (CWE), relativistic oscillating mirror (ROM), and coherent synchrotron emission (CSE). In the past fifteen years, SHHG source development has evolved from physical understanding to characterization

and field manipulation. Spectral, spatial, and temporal characterization and modulation via laser-plasma interaction control are crucial for future applications. Ptychography [5] and Spectral Phase Interferometry for Direct Electric-Field Reconstruction (SPIDER) [9] methods have been proposed to accurately characterize both the spatial and spectral phases of SHHG sources. Concerning harmonic source modulation, the conversion efficiency has been boosted through the temporal shaping of multi-color driving-field-cycles [10-12]. The SHHG process is particularly sensitive to the plasma density gradient at the plasma-vacuum boundary: it is a key parameter for the conversion rate [13, 14], the temporal properties on the femtosecond scale [6, 7] as well as the spatial properties of SHHG [15] through the dynamic denting of the plasma surface [15, 16]. Methods were proposed to compensate for the spatial degradations including plasma truncation with an extra pre-pulse [17], moving the target out of focus for wave front compensation [18] and using a convex target to achieve a flat surface despite plasma denting-induced concavity [19]. Compensation of the CWE intrinsic phase via the pre-plasma distribution was also demonstrated [20, 21].

Experimental control of steep pre-plasma profiles is typically achieved by diverting a small fraction of the laser energy from the main beam using a holey or pick-up mirror prior to focusing [14, 22]. The small pre-pulse beam and the larger main beam are then recombined and focusing by the same optic, ensuring the small main-pulse focus to be accurately positioned with respect to the larger pre-pulse focus. The SHHG signal is collected in the specular reflection direction, hundreds of millimeters distant from the laser-plasma interaction point, within the far-field region of the harmonic beam, which typically exhibits a divergence angle ranging from tens to 100 mrad [5, 14, 15, 22]. For metrology and applications, controlling the beam pointing is crucial. The harmonic beam characteristics, including directionality, have been the subject of theoretical and experimental studies over the past few decades. However, a common limitation of current theories is the

<footnote n="1"><sup>1</sup>Shenzhen Key Laboratory of Ultraintense Laser and Advanced Material Technology, Center for Intense Laser Application Technology, and
College of Engineering Physics, Shenzhen Technology University, Shenzhen 518118, China</footnote>

<footnote n="2"><sup>2</sup>College of Applied Sciences, Shenzhen University, Shenzhen 518060, China</footnote>

<footnote n="3"><sup>3</sup>Frontiers Science Center for Rare Isotopes and School of Nuclear Science and Technology, Lanzhou University, 730000 Lanzhou, China</footnote>

<footnote n="4"><sup>4</sup>Laboratoire d'Optique Appliquée, Institut Polytechnique de Paris, ENSTA Paris, Ecole Polytechnique, CNRS, Palaiseau, France</footnote>

<footnote>*stefan.haessler@ensta-paris.fr</footnote>

<footnote>*lilu@sztu.edu.cn</footnote>
assumption of a 1D pre-plasma density distribution, typically modelled as a transversely uniform exponential decay profile. Practically, the pre-pulse induced density can exhibit significant transverse spatial variation, leading to an initial reflecting surface with a tilt. When the main pulse interacts with this tilted surface, the harmonic beam deviates from the specular direction.

## 2. Optical Tilt on A Gaussian Pre-plasma Profile

To demonstrate this effect, a transversely (i.e. parallel to the target surface) Gaussian pre-plasma density profile is considered, with an exponential decay along the target normal direction [23]. This is realistic for a Gaussian pre-pulse profile since the plasma expansion velocity is known to be proportional to the pre-pulse fluence [22, 23]. The target is placed on the x-z plane, while the y axis is perpendicular to the target. We solve the problem in the x-y laser incidence plane. Note that we can simply extend the following theory to obtain comprehensive 3D results. In figure 1(a), the pre-pulse and main-pulse are illustrated by the lighter and darker red cones, respectively. The blue highlighted region is where the plasma density exceeds the critical density $n_c \cos^2 \theta$, where $\theta$ signifies the angle of incidence and $n_c = (\epsilon_0 m_e / e^2) \omega^2$, with $\omega$ being the laser angular frequency, e the electron charge, $m_e$ the electron mass and $\epsilon_0$ the vacuum permittivity. We model the electron density profile $n_e(x, y)$ as:

$$n_e(x, y) = n_0 \exp \left[ -\frac{y}{L(x)} - 1 \right], L(x) = L_0 \exp \left[ -\frac{x^2}{w_{\text{pre}}^2} \right], \quad (1)$$

where $n_0$ is the peak density reached in the target bulk, x is the transverse distance to the pre-pulse center, and y is the distance along the target normal from the initial solid surface before thermal expansion.

<figure><img src="image_3.png" /><figcaption>Fig. 1. (a) Schematic for the origin of pre-plasma tilt: the interaction between a Gaussian distributed pre-plasma and the main pulse. The lighter red cone, darker red cone and purples cone represents the pre-pulse, main pulse and HHG beam, respectively. For clarity in the figure, the main pulse has been exaggerated in size. The purple lines signify the reflection directions with and without pre-plasma. (b) angle deviation induced by the geometrical tilt at different focal positions for $w_{\text{pre}} = 10\lambda_L$ and $25\lambda_L$.</figcaption></figure>

To elucidate the impact of the Gaussian pre-pulse on the reflected beam, we initially assume the critical surface as a static and perfect reflector. Using the scenario without a pre-plasma ($L_0 = 0$) as a benchmark to measure beam deviation caused by geometric tilt, the EUV beam propagates along the specular direction, indicated by the purple dashed lines starting from $y = 0$. To address the beam deviation when a pre-plasma is present ($L_0 \neq 0$), the reflection point ($x_r, y_r$) is determined by calculating the intersection between the geometric critical surface $y_c(x)$ (where $n_e(x, y) =$

$n_c \cos^2 \theta$ for p-polarized laser [24]) and the light ray path $y_{\text{ray}}(x)$, which are represented as:

$$y_c(x) = L_0 \exp \left[ -\frac{x^2}{w_{\text{pre}}^2} \right] \ln \left( \frac{n_0 / e}{n_c \cos^2 \theta} \right), y_{\text{ray}}(x) = \frac{x - x_0}{\tan \theta}, \quad (2)$$

Here, e is Euler's number (according to Eq. 1, the density at the original solid-vacuum boundary, $y=0$, remains constant at $n_0 / e$ over the whole laser heating area), and $x_0$ is the transverse distance of the main pulse focal spot to the center of the pre-plasma distribution, as shown in Fig. 1(a). Note the domain of definition of $y_{\text{ray}}$ is from the far-field beam center on the focusing mirror to the intersection point, which serves as the secondary reflection surface on the plasma.

The inclination of the reflecting surface due to the pre-plasma expansion, as depicted in Fig. 1(a), causes a deviation of the reflected beam from the expected specular direction by the angle

$$\alpha = 2 \arctan \left. \frac{dy_c(x)}{dx} \right|_{x=x_r} \approx 2 \ln \left( \frac{n_0 / e}{n_c \cos^2 \theta} \right) \delta(x_r), \quad (3)$$

where the slope $\delta(x_r)$ is defined as $\delta(x) = dL/dx|_{x=x_r}$. It is a more universal indicator for the pre-plasma tilt than the specific L and $w_{\text{pre}}$ in Eq. 1, and for a Gaussian pre-plasma can be written as:

$$\delta = -\frac{2x_r}{w_{\text{pre}}^2} L(x_r), \quad (4)$$

Fig. 1(b) shows the deviation angle $\alpha$ at different main-pulse positions for two different pre-plasma conditions: a narrow profile with $w_{\text{pre}} = 10\lambda_L$ and a normal one with $w_{\text{pre}} = 25\lambda_L$, $\theta = \pi/4$ and $L = 0.1\lambda_L$ are the same in these two cases. Assuming a bulk plasma density such that $n_0 / e = 100n_c$, which is close to the case of a fully ionized SiO$_2$, the deviation angle $\alpha$ reaches tens of milliradians at $x_r = w_{\text{pre}} / \sqrt{2}$. This is very similar to typically observed SHHG beam divergences [5,14, 15, 22] and can thus be obviously observed in experiments. While the deviation angle at $x_0 = 0$ remains fairly small ($\alpha \approx 10$ mrad) even with the narrow pre-plasma profile, for longer scale lengths $L > 0.1\lambda_L$, significant beam deviations may occur. Around the main-pulse focus position $x_0$, L and $\delta$ are used as two independent indicators for describing the pre-plasma condition. The transverse scale-length variation $\delta$ at the flank of the pre-plasma, with a maximum value of $\delta \approx 8 \times 10^{-3}$ for the narrow pre-plasma profile considered above, is identified as the origin of an effective surface tilt, deviating the reflected beam.

## 3. Plasma Correction on the Beam Drift

The above model describes the geometrical angle deviation. However, upon the incidence of a relativistic, tightly focused short pulse, the reflection does not happen at a fixed surface but is caused by a dynamic compressed electron sheet driven around the critical-density region. This compression modifies the reflection spot relative to the original surface position, introducing a correction to the geometrical deviation. To assess the impact of this plasma effect on the directionality of the harmonic source, particle-in-cell (PIC) simulations were conducted with the WarpX code [25]. They show that the actual tilt angle observed on the reflected pulse is reduced compared to the geometrical tilt angle, suggesting a correction coefficient <1 should be applied to the right-hand side of Eq. (3). Once this correction is known, the easily measurable modulation of
the SHHG beam direction directly indicates the critical surface tilt, which in turn can be employed to characterize the plasma expansion process within a pump-probe like scheme, where the pump is the pre-pulse and the probe is the SHHG-driving pulse. This approach can provide valuable insights into the behaviour of the plasma under the influence of the driving laser field.

The results of the PIC simulations are illustrated in Fig. 2(a). A polarized laser is injected into the simulation box, with a beam waist of $w_0 = 1.5\lambda_L$ and a duration of $\tau = 2T_0$. Here, $\lambda_L = 800$ nm and $T_0 = \lambda_L/c$ represent the fundamental laser wavelength and period, respectively. The normalized vector potential is $a_0 = 2.5$, and

<figure><img src="image_4.png" /><figcaption>Fig. 2 Simulation results for different scale length L and decay slope δ. (a) 2D spectrum of L = 0.04 and 0.14, (b) lineouts of far-field angular spectrum map driven by different scale lengths, the centroid of harmonic spectrum is shown in dashed lines, (c) the spectrum drift caused by increasing δ, the centroid of spectrum is shown in dashed lines. (d) the simulated angular drift by different decay slope, the fitting linear relationship is shown with error bar, which indicates for $d\alpha/d\delta = 6.1 \pm 0.9$.</figcaption></figure>

the incidence angle is $\theta = 45^\circ$. The target consists of a fully ionized carbon slab with a thickness of $2\lambda_L$, with a maximum density of $100n_c$. This maximum density remains constant at the initial slab surface and exponentially decays with rate L that varies along the slab surface. Since the plasma denting depth in the conditions considered here remains much smaller than the distance from the critical density surface to the $100n_c$-surface, this is equivalent to the situation considered for the geometrical model above, where Eq.1 fixed the plasma density at $100n_c$ at the initial vacuum-plasma interface. The simulation grid size is set to $\lambda_L/512$, and 100 particles per cell are used for resolving up to the 30th order harmonic signals. With a slope $\delta$, the scale length $L(x)$ of the pre-plasma varies as $L = L_c + \delta x$, where $L_c$ is the central scale length at $x = 0$. Typical results with $L_c/\lambda_L = 0.04, 0.07, 0.11, 0.14$ and $10^4\delta = 5, 10, 15, 20$ are shown in Fig. 2. The beam drift was quantified through the angular distribution of the simulated harmonics. A near-field-to-far-field transformation was applied to project the reflected harmonic beam from the simulation domain into the far-field. Subsequently, the drift angle $\alpha$ was determined by the displacement from the specular reflection direction of the

centroid of the high-order harmonic ($n > 10$). The simulated near-field spectrum at simulation box boundary for $L_c = 0.07, \delta = 0$ and $L_c = 0.14, \delta = 0$ are displayed in Fig. 2(a) as an example, where $k_{\text{spec}}$ is along the specular reflection direction. The corresponding far-field angular distributions for various $L_c$ values with $\delta = 0$ were calculated by integrating the spectral intensity along each far-field angular direction and are shown in Figure 2(b). The distributions' centroids are shown by dashed lines, which indicate the angular drift induced by the scale-length variation. The inherent angle $\alpha_0$ for the uniform target with $\delta = 0$ was found to be 13, -2, 21, 16 mrad for the scale lengths $L_c/\lambda_L = 0.04, 0.07, 0.11, 0.14$, respectively. In this scenario, few-mrad level angular drifts are observed, which are small compared to the overall angular divergence and are dominantly caused by a morphing beam shape, without a clear dependence on $\delta$. Note that a very tightly focused driving beam is adopted in this work, and when it comes to the typical experimental results with an angular divergence of ~50 mrad, the beam drift due to variations in $L_c$ can be expected to be less than 5 mrad.

While the scale length $L_c$ has a minimal impact on the beam directionality, the decay rate $\delta$ influences the angular drift. With $L_c$ set to $0.04\lambda_L$, the angular intensity distributions for various $\delta$ values are depicted in Figure 2(c). The angular profile shifts monotonously without changing its general shape, demonstrating a linear relationship between the angular drift $\alpha$ (with respect to $\delta = 0$) and the slope $\delta$, consistent with the prediction from Eq. (3). However, Figure 2(d) reveals, across different scale lengths $L_c$, a linear slope $k_{\text{PIC}} = d\alpha/d\delta = 6.3 \pm 0.7$, which is below the value derived from Eq. (3): $k_c = 10.6$. Note that the large error range is mainly caused by the morphing effect at different $\delta$. This discrepancy confirms our expectation of a correction factor <1 arising from the plasma effect that affects the geometrical angular deviation.

<figure><img src="image_5.png" /><figcaption>Fig. 3. Correction of decay slope on the geometrical profile by the plasma effect. (a) density maps driven by laser at different instants for $\delta = 0.002$ and (b) the corresponding lineouts at the equilibrium positions for $\delta = 0$ and 0.002, the dashed lines show the equilibrium surface dynamics along $x$ when laser swipes along target, blue represents for the flat target and green represents for the tilt one, (c) the reflecting surface motion around the end of compression phase for position 1, (d) denting depth difference between $\delta = 0$ and 0.002.</figcaption></figure>

To elucidate the plasma correction, an analysis of the electron
dynamics in the pre-plasma was conducted. The SHHG process can be divided into three distinct phases based on the motion of the reflecting surface: compression phase, acceleration phase, and dephasing phase. According to the plasma denting theory [12], electrons in the pre-plasma reach an equilibrium position at the end of the compression phase. It can be considered as an indicator for the emission depth. In the case of short driving pulses, the effect of the positive ion motion can be neglected. The equilibrium depth that the reflecting surface ultimately reaches, relative to the critical surface position $y_c(x)$ (Eq. 2), is given as [12]:

$$y_e = -L \ln \left[ 1 + \frac{\lambda_L a_0 (1 + \sin \theta)}{\pi L} \frac{1}{\cos^2 \theta} \right], \quad (5)$$

Note that $y_c(x)$ is not limited to a Gaussian distribution. Every function has its own geometrical decay slope, and thus the geometrical linear factor can be calculated from Equation (3). The plasma denting effect can thus provide a correction to the reflection direction via $dy_e(x)/dx$. Since $y_e$ is approximately linearly related to $L$ within a main-pulse spot that is small compared to the pre-plasma, it results in a nearly constant correction to the angle deviation $\alpha$.

In physical terms, the plasma denting counteracts the tilting of the reflective plasma surface caused by the slope $\delta(x_r)$ because a larger scale length facilitates pushing the reflecting electron sheet further into the dense plasma region. When the laser swipes over the plasma surface in Fig. 3(a), the equilibrium depth difference between the tilt case with $\delta \neq 0$ and the uniform case with $\delta = 0$ can be regarded as the reason for the beam drift. The three panels show snapshots at different instants spaced by one laser cycle, where the center panel gives the density map when the electron sheet reaches $x = 0$. Vertical lineouts of the electron density at the three positions marked in Fig. 3(a) are plotted in Fig. 3(b). To exhibit how $\delta$ tilts the reflective surface, the lineouts at the same instants and positions are shown for $\delta = 0$ as a reference. The sheet in the flat case ($\delta = 0$) starts from a shallower position and ends at a deeper position, opposite to the tilt case ($\delta = 0.002$). The reflective-surface dynamics around the equilibrium position are depicted in Fig. 3(c). Using the equilibrium depth at $x = 0$ (position 2 in Fig. 3(a)) as a reference, if the equilibrium depth at position 1 was calculated by only considering the geometrical tilt of the critical surface, it should be located at the blue dotted line. If the correction by the plasma effect from $dy_e(x)/dx$ is considered, the equilibrium depth is predicted perfectly. That means that our combined theory with the geometrical tilt and the plasma effect can be used to calculate the equilibrium depth for different spatial positions. The relative difference between the equilibrium positions of a flat and a tilted surface at different transverse positions can be used to evaluate the correctness of Eq. 5, shown in Fig. 3(d). Good agreement is found between the combined theory and the PIC-simulation results. Quantitatively, the angle deviation induced by the plasma denting correction can be expressed here for our simulation parameters as:

$$k_e = \frac{d}{d\delta} \left( 2 \arctan \frac{dy_e}{dx} \right) \approx 4.8, \quad (6)$$

Note that $k_e$ varies with $a_0$, and here the result at $x = \sqrt{2}\lambda_L$ was used as an averaged value. The combined angular drift is thus $\alpha = (k_c - k_e)\delta = 5.8\delta$, which fits well the simulation results of $k_{PIC} = 6.1 \pm 0.9$. While the plasma correction is only verified here

with a tightly focused few-cycle pulse, more simulations were performed to verify its robustness. The dependence of the angular deviation on the laser parameters, including carrier envelope phase and intensity, remained consistent with our model and our conclusions can be applied to a wide range of SHHG driving conditions. Note that when a long driving pulse is used, a further correction should be considered to account for ion motion [15]. However, in our parameter range, from 2 to tens of femtoseconds, it has been verified in the simulations that the ion motion has little influence on the results above. The detailed correction from scale-length and ion motion is out of scope of this article, and it will be studied in our future work.

## 4. Conclusion

In conclusion, this study has thoroughly examined the effects of SHHG beam drifting for various pre-plasma distributions. For SHHG driven by a few-cycle laser, the deviation of the harmonic emission angle from the specular direction is primarily governed by the slope of the critical surface. The theory and simulations indicate that the drifting of the harmonic beam can serve as a probe for the critical surface condition at the main pulse focal spot. When a high-contrast laser system is utilized, the geometrical tilt dynamics of the pre-plasma can be continuously measured with varying pre-pulse delays. Furthermore, SHHG beam drifting as well as the far-field wave front information carries information on the relativistic plasma at micrometer and femtosecond to picosecond scales.

**Funding.** National Natural Science Foundation of China (Grant No. 12205203), LASERLAB-Europe (LOA0024038, 654148, H2020-EU.1.4.1.2), Agence Nationale de la Recherche (ANR-22-CE30-0005-01 BANDITO), Guangdong Province Key Construction Discipline Scientific Research Capacity Improvement Project (Grant No. 2021ZDJS107)

**Acknowledgments.** This research used the open-source particle-in-cell code WarpX https://github.com/ECP-WarpX/WarpX. Primary WarpX contributors are with LBNL, LLNL, CEA-LIDYL, SLAC, DESY, CERN, and TAE Technologies. We acknowledge all WarpX contributors.

**Disclosures.** The authors declare no conflicts of interest.

**Data availability.** Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.

## References

1. C. Thaury, F. Quéré, J.-P. Geindre et al., Nature Phys 3, 424 (2007).

2. T. Baeva, S. Gordienko, and A. Pukhov, Phys. Rev. E 74, 046404 (2006).

3. C. Thaury and F. Quéré, J. Phys. B: At. Mol. Opt. Phys. 43, 213001 (2010).

4. U. Teubner, and P. Gibbon, Rev. Mod. Phys. 81, 445-479 (2009).

5. A. Leblanc, S. Monchocé, C. Bourassin-Bouchet et al., Nature Phys 12, 301 (2016).

6. D. an der Brügge and A. Pukhov, Physics of Plasmas 17, 033110 (2010).

7. A. A. Gonoskov, A. V. Korzhimanov, A. V. Kim et al., Phys. Rev. E 84, 046403 (2011).

8. M. R. Edwards and J. M. Mikhailova, Sci Rep 10, 5154 (2020).

9. C. Wu, L. Li, M. Yeung et al., Opt. Express 30, 389 (2022).

10. M. Yeung, S. Rykovanov, J. Bierbach et al., Nature Photon 11, 32 (2017).

11. L. Chopineau, G. Blaclard, A. Denoeud et al., Phys. Rev. Research 4, L012030 (2022).

12. M. R. Edwards and J. M. Mikhailova, Phys. Rev. Lett. 117, 125001 (2016).
13. S. Kahaly, S. Monchocé, H. Vincenti et al., Phys. Rev. Lett. **110**, 175001 (2013).

14. S. Haessler, M. Ouillé, J. Kaur et al., Ultrafast Sci **2022**, 2022/9893418 (2022).

15. H. Vincenti, S. Monchocé, S. Kahaly, et al., Nat Commun **5**, 3403 (2014).

16. H. Vincenti, Phys. Rev. Lett. **123**, 105001 (2019).

17. B. Y. Li, F. Liu, M. Chen et al., Phys. Rev. Lett. **128**, 244801 (2022).

18. Y. Zhang, C. Zhong, S. Zhu, et al., New J. Phys. **24**, 033038 (2022).

19. J. Gao, B. Li, F. Liu, et al., Phys. Rev. E **101**, 033202 (2020).

20. E. Porat, H. Yehuda, I. Cohen et al., Phys. Rev. Research **3**, L032059 (2021).

21. H. Yehuda, E. Porat, I. Cohen et al., Opt. Lett. **46**, 4674 (2021).

22. J. Kaur, M. Ouillé, D. Levy et al., Review of Scientific Instruments **94**, 113002 (2023).

23. P. Mora, Phys. Rev. Lett. **90**, 185002 (2003).

24. P. Gibbon, Short pulse laser interactions with matter : an introduction, (Imperial College Press, 2005).

25. L. Fedeli et al., *SC22: International Conference for High Performance Computing, Networking, Storage and Analysis*, Dallas, TX, USA, pp. 1-12 (2022).