Time-resolved optical shadowgraphy of solid hydrogen jets as a testbed to benchmark particle-in-cell simulations

Particle-in-cell (PIC) simulations are a superior tool to model kinetics-dominated plasmas in relativistic and ultrarelativistic laser-solid interactions (dimensionless vectorpotential $a_0>1$). The transition from relativistic to subrelativistic laser intensities ($a_0 \lesssim 1$), where correlated and collisional plasma physics become relevant, is reaching the limits of available modeling capabilities. This calls for theoretical and experimental benchmarks and the establishment of standardized testbeds. In this work, we develop such a suitable testbed to experimentally benchmark PIC simulations using a laser-irradiated micron-sized cryogenic hydrogen-jet target. Time-resolved optical shadowgraphy of the expanding plasma density, complemented by hydrodynamics and ray-tracing simulations, is used to determine the bulk-electron temperature evolution after laser irradiation. As a showcase, a study of isochoric heating of solid hydrogen induced by laser pulses with a dimensionless vectorpotential of $a_0 \approx 1$ is presented. The comparison of the bulk-electron temperature of the experiment with systematic scans of PIC simulations demostrates that, due to an interplay of vacuum heating and resonance heating of electrons, the initial surface-density gradient of the target is decisive to reach quantitative agreement at \SI{1}{\ps} after the interaction. The showcase demostrates the readiness of the testbed for controlled parameter scans at all laser intensities of $a_0 \lesssim 1$.


I. INTRODUCTION
High-impact applications of high-intensity laser-solid interactions such as fast ignition in inertial-confinement fusion [1][2][3], laser-driven ion acceleration [4,5], the investigation of warm-dense matter [6][7][8][9][10] or secondary-source development [11][12][13][14][15][16] have developed into independent research fields over the last years.Gaining insight into the microscopic interaction picture is the domain of numerical modeling through simulations.In practice, the simulations are most often used retrospectively to guide the analysis and interpretation of exploratory experiments.
Prediction making and the development of long-term strategies remains challenging.In translational research in radiation oncology [17][18][19][20][21], as one prominent example, the extrapolation of realized proof-of-concept studies towards future high-energy laser-driven proton accelerators [22,23] has been delayed [24,25].The technical realization of optimal laser-target interaction conditions and the further development of available simulation tools to correctly capture all occurring physical processes during the transition of the target from an initially solid state to an ultrarelativisticplasma state represent the key aspects.State-of-the-art modeling of ultrarelativistic laser-solid interactions with Petawatt-class lasers (peak intensity ≳ 10 21 W/cm 2 ) is currently based on a staged approach of numerical simulations [26][27][28] that includes target preexpansion by the leading edge of the high-power laser [29].For intensities in the vicinity of the laser peak, i.e., dimensionless vectorpotentials a 0 ≳ 1, particle-in-cell (PIC) simulations are optimized to compute the kinetic regime of particle motion and thermal nonequilibrium.Adequate initial conditions of the PIC simulations are derived by radiation-hydrodynamics simulations, which compute the interaction in the subrelativistic intensity regime (a 0 ≪ 1) during the leading edge.However, the transition from relativistic to subrelativistic laser intensities, i.e., dimensionless vectorpotentials a 0 ≲ 1 (10 16 to 10 18 W/cm 2 for 800 nm light), is currently reaching the limits of available modeling capabilities.
The two most obvious approaches to cover this laser-intensity regime are being pursued; the extension of hydrodynamics-simulation tools [30] and the inclusion of correlated and collisional plasma physics into PIC-simulation frameworks [25,31].These developments call for standardized theoretical and experimental benchmarks under unified interaction scenarios [32].As microscopic parameters are usually not directly accessible from laser-plasma interactions, such benchmarks would ideally require macroscopic observables that allow for an unambiguous allocation of specific interaction conditions.
In this work, we present a testbed to experimentally benchmark PIC simulations based on a laser-irradiated micron-sized cryogenic hydrogen-jet target [33][34][35].The corresponding PIC simulations of the laser-target interaction benefit from the comparably low target density [36][37][38][39][40], the single-species composition, the negligible amount of Bremsstrahlung radiation (atomic number Z = 1) and simple ionization dynamics [41].Furthermore, the plasma composition of only protons and electrons enhances comparability to analytic calculations.The temporal evolution of the laser-heated plasma density is visualized by time-resolved off-harmonic optical probing [42] via two spectrally separated, ultrashort laser beams.A fitting approach by hydrodynamics and ray-tracing simulations completes the testbed.It enables the determination of the bulk-electron temperature evolution after the laser energy was absorbed by the target.The capabilities of our testbed are showcased in an experiment studying isochoric heating of solid hydrogen when irradiated by laser pulses with 37 fs duration and a 0 ≈ 1.The focal spot of the laser is kept larger than the spatial dimensions of the target to guarantee a homogeneous interaction.Isochoric heating of solids by high-intensity lasers is a widely-used approach to generate warm-dense and hot-dense matter [43][44][45].Therefore, the process is particularly well suited to benchmark PIC-simulation results.During isochoric heating, the ultra-short laser pulse accelerates electrons on the laser-plasma interface up to relativistic velocities [46][47][48][49][50]. Subsequently, the electrons traverse the target and generate a thermalized bulk-electron and bulk-ion temperature, mainly by resistive heating, drag heating or diffusion heating [51][52][53][54].Finally, the plasma undergoes an adiabatic expansion into vacuum and the electron and ion temperatures equilibrate [55].For lasers with a dimensionless vectorpotential a 0 ≲ 1, a transition from vacuum heating to resonance heating of electrons exists, which depends on the surface-density scalelength of the target [56][57][58][59][60].Furthermore, for a ≲ 1, the collisionality of particles increases in relevance.The bulk-electron temperature after isochoric heating represents a specifically feasible endpoint, i.e., a macroscopic observable, to benchmark the entire modeling chain of nonthermal equilibrium that can be computed by PIC simulations, including laser acceleration of electrons and their thermalization to a single temperature within the target bulk.Experimentally, the determination of the achieved bulk-electron temperature most often relies on spatially and temporally integrated x-ray spectroscopy [51,[61][62][63][64][65].Time-resolved optical shadowgraphy of the induced expansion provides a measure of the bulk-electron temperature subsequent to isochoric heating [66,67].
After introducing the testbed and the method for determining the bulk-electron temperature from the experimental data, the results are compared to systematic scans of PIC simulations.We find that, due to the interplay of vacuum heating and resonance heating of electrons, the bulk-electron temperature is highly sensitive to the surface-density gradient of the target.By demonstrating the capability of the testbed to quantitatively benchmark PIC simulations, it establishes the basis towards controlled parameter scans at laser intensities well below a 0 = 1.This transition regime from solid state to ultrarelativistic plasmas is highly relevant to all high-intensity laser-solid interaction and their applications.

II. TESTBED PLATFORM
This section presents the testbed to experimentally benchmark PIC simulations of high-intensity laser-solid interactions via the endpoint of the bulk-electron temperature after isochoric heating.The endpoint includes a chain of physical processes, i.e., target ionization, the laser-acceleration of electrons as well as their thermalization within the initially cold target bulk to a single Maxwellian electron temperature.The bulk-electron temperature well above the Fermi energy induces an adiabatic plasma expansion into vacuum, which is visualized by time-resolved optical diagnostics.Section II A presents the time-delay scan of two-color shadowgraphy in the tens-of-picosecond timeframe together with the utilized experimental setup.The determination of the corresponding electron-temperature evolution is explained in section II B. Section II B 1 assumes a homogeneous electron temperature throughout the central plane of the target and calculates the electron-density evolution via a hydrodynamics simulation.Section II B 2 simulates the corresponding shadow diameters versus delay by raytracing simulations and section II B 3 presents a χ 2 fit of the simulated to the measured shadow diameters.By this we derive the best-matching hydrodynamics simulation and the corresponding electron-temperature evolution that corresponds to the bulk-electron-temperature evolution after isochoric heating and thermalization.The combination of hydrodynamics simulations, ray-tracing simulations and the χ 2 fit is referred to as HD-RT fit in the following.

A. Experiment
The experiment is conducted with the Draco-150 TW laser at the Helmholtz-Zentrum Dresden-Rossendorf [68].The top view of the setup is shown in figure 1  (a).By placing a 16 mm-wide circular aperture in the collimated beam path, the pump-laser energy is reduced to 160 mJ.Each laser pulse features a duration of 37 fs full width at half maximum (FWHM) and a central wavelength of λ laser = 800 nm.It is focused by a f/16 off-axis parabola (OAP) and generates an Airy-pattern focus with a spatial FWHM of 14 µm of the central disk and a peak intensity of I laser = 1.6 × 10 18 Wcm −2 .A continuous, self-replenishing jet of cryogenic hydrogen is used as a target [33,34].The cross section of the target at the source is defined by a circular aperture with a diameter of (5 ± 1) µm.At the position of the laser-target interaction, the mean diameter is 4.4 µm with a standard deviation of 0.2 µm (details in Appendix B).The electron density of the fully ionized target is 5.2 × 10 22 cm −3 , which corresponds to 30 times the critical plasma density of 800 nm light n 800 nm c .The laser-target interaction is investigated by timeresolved off-harmonic optical shadowgraphy [42] at 90 • angle to the pump-laser propagation direction.Two copropagating backlighter pulses with the wavelengths λ of 515 nm and 258 nm are generated from a synchronized stand-alone laser system [69].The 258 nm probe and the 515 nm probe feature a pulse energy of approximately 0.3 µJ and 5 µJ and a pulse duration of 260 fs and 160 fs, respectively.The pump-probe delay is variable, the temporal resolution of a time-delay scan is 175 fs and the captured data is blurred by the pulse duration [42].Dispersion effects in the beam path [70] cause a fixed delay of 4.6 ps between both probe-laser pulses.All delays are given with respect to the arrival of the pump-laser peak on target at 0 ps delay.To reduce the influence of parasitic plasma emission, both probe beams are focused by a f/1 OAP.A long-working-distance infinite-conjugate microscope objective (designed for laser wavelengths of 266 nm and 532 nm) is used to image the shadow of the target simultaneously for both wavelength, i.e., at two different ranges of plasma density.The two images are separated by a dichroic mirror and recorded by different cameras, each equipped with the corresponding color filters (more details in Appendix A).A pump-probe time-delay scan in steps of 5 ps is performed.The inherent spatial jitter of the target is controlled by a secondary optical-probe axis (not shown here) [37].Only data from within ±2 µm of the object plane is presented.Exemplary shadowgrams at 0 ps, 10 ps and 30 ps delay are shown in figure 1  both increase with pump-probe delay until volumetric transparency sets in.For the 258 nm probe and the 515 nm probe, this occurs between 20 and 25 ps and between 45 and 50 ps, respectively.For all delays greater than zero, D 515 nm E is on average larger than D 258 nm E .This is expected, as the critical plasma density n c drops with increasing wavelength λ according to Assuming an exponentially decreasing plasma-density surface, the difference between D 258 nm E and D 515 nm E enables the calculation of the corresponding scalelength L 0 .For 0 ps delay, the surface scalelength calulates to values between 0 and about 0.42 µm.Figure 2 (a) shows shot-to-shot fluctuations of the shadow diameter.They are caused by inherent shot-to-shot variations of the target and the pump laser.Target variations include local changes of diameter and aspect ratio on the submicron scale as well as bow-like deformations along the jet axis on the micron scale (details in Appendix B).Furthermore, the rapid freezing of the jet introduces compositional and structural variations of the solid hydrogen [71].The main source of fluctuation of the pump-laser intensity is the peak power with a measured standard deviation of 12 %.The fluctuation of the pump-laser energy is about 1 % and the pointing jitter is negligible.

Hydrodynamics simulation -HD
We assume that the density evolution of the hydrogen plasma is driven by a two-temperature adiabatic expansion, which can be modeled via hydrodynamics simulations.In the following, the simulation tool FLASH [72,73] is used to perform two-temperature hydrodynamics simulations in two-dimensional radial symmetry with the hydrogen equation of state FPEOS [74] (details in Appendix D).Initially, there are three free parameters of the target model: the ion temperature T i0 , the electron temperature T e0 and the initial plasma diameter D 0 .Because the evolution of the shadow diameter is sensitive to the initial plasma diameter D 0 , which is subject to experimental fluctuations, we handle D 0 as a free parameter within the range of the experimental uncertainties.Compared to the effect of T e0 , the initial ion temperature T i0 ≪ T e0 has little effect on the expansion process.
As derived by PIC simulations in section III, T i0 can be approximated by 10 eV.Hydrodynamics simulations with different initial electron temperatures T e0 and different initial plasma diameters D 0 are conducted.Each simulation box has a length of 50 µm and a vacuum density of 1 × 10 −8 times the target density.An exemplary evolution of the electron density is shown in figure 3 (a) and the corresponding evolution of ion and electron temperature is shown in figure 3 (b).

Ray-tracing simulation -RT
To compare the evolution of the electron density in the hydrodynamics simulation with the experimentally measured shadowgraphy data, the shadow formation of each probe beam needs to be modeled.For optical shadowgraphy of solid-density plasmas, shadow formation is governed by refraction on the density gradients of the under-critical regions of the plasma [42].Each simulated density profile is transformed into a two-dimensional distribution of refractive indices ñ from the local electron density n e via the dispersion relation [75] The critical density n c depends on the wavelength and ñ258 nm and ñ515 nm are calculated separately.
The experimental imaging setup and the beam path of each probe are reproduced in a virtual optical setup with the software Zemax [76].The calculated spatial distribution of refractive indices ñ258 nm or ñ515 nm is inserted into the object plane of the corresponding setup, which is exemplary shown in figure 3 (c).The purple lines illustrate the refraction of the 258 nm-probe rays in the gradients of the refractive-index distribution.The light distribution in the image plane is calculated by Zemax and presented in figure 3 (d).The graph shows the formation of a shadow with sharp edges, similar to the experiment.Refraction leads to an increased intensity level at the rim of the shadow edges.The shadow diameter D 258 nm S is derived at half of the unperturbed background intensity (0.5 arb.units).By utilizing the virtual setup of the 515 nm probe, D 515 nm S is calculated accordingly.The simulated evolution of D 258 nm S (cyan) and D 515 nm S (blue) of a hydrodynamics simulation with the initial setting T e0 = 250 eV and D 0 = 3.5 µm (case α) is shown as solid lines in figure 2 (a).Like in the experiment, D S is set to zero at delays for which no sharp shadow edge is derived and volumetric transparency is observed instead.For the 258 nm probe and the 515 nm probe of case α, this occurs at 22.5 ps and 47.5 ps delay, respectively.

χ 2 fit
To find the best-matching hydrodynamics simulation to the experimental data, the initial electron tempera- ture T e0 is varied from 100 to 350 eV in steps of 50 eV.Furthermore, the initial plasma diameter D 0 is scanned from 3 to 5 µm in steps of 0.5 µm.The variance χ 2 of the difference between the experimental data D E (t, λ) and the simulation data Here, t is the pump-probe delay and λ is the wavelength of the probe.The resulting χ 2 map is shown in figure 2 (b).Minimum values of χ 2 are reached for case α (T e0 = 250 eV, D 0 = 3.5 µm) and case β (T e0 = 300 eV, D 0 = 4 µm).Case β has better agreement to the 258 nm-probing data and case α has better agreement to the 515 nm-probing data.The two cases give a lower and upper limit of the best-fitting T e0 and D 0 .A discussion of the method of the HD-RT fit is given in Appendix C.
In summary, the HD-RT fit allows to fit a heuristic electron-temperature evolution subsequent to isochoric heating and thermalization of the bulk electrons, which constitutes the endpoint for the comparison to PIC simulations in the presented testbed.For the here presented experimental data, the HD-RT fit yields a heuristic initial bulk-electron temperature T e0 between 250 and 300 eV at 0 ps delay.

III. PIC SIMULATION
This section exemplary demonstrates a comparison of the results of the testbed with a PIC simulation of the high-intensity laser-solid interaction.The comparable endpoint is the temporal evolution of the bulk-electron temperature.The presentation of the results of the PIC simulation details the processes of isochoric heating and thermalization of the bulk electrons (not captured by the hydrodynamics simulation of the HD-RT fit).Both finally lead to the evolution of the bulk-electron temperature after thermalization, which can be compared to the hydrodynamics simulation of the HD-RT fit.The simulation is performed with the 2D3V-PICsimulation tool PIConGPU [77].The target is modeled by a fully ionized spherical hydrogen column.The density distribution follows the model with n 0 = 30 n 800 nm c , a surface scalelength L 0 = 0.25 µm (mean of the experimental uncertainty), an unperturbed target diameter of D 0 = 4.4 µm (see Appendix B) and the reduced radius to fulfill mass conservation.The initialized total length of the surface plasma corresponds to ten times L 0 .The laser is incident along the y direction and polarized in x direction.It features a Gaussian shape in space and time, a pulse duration of 37 fs FWHM and a laser-spot size of 14 µm FWHM.The simulation box is 32 × 32 µm 2 with absorbing boundary conditions and the simulation uses a cell length of 0.8 µm/96 and 240 macro particles per cell.The simulation includes a relativistic-binarycollision mechanism. .Comparing the electron-density profiles at 0.2 ps, 1.0 ps and 2.1 ps, two characteristics are identified.The density region denoted by "1" shows a temporally increasing exponential scalelength of the plasma density.It results from adiabatic expansion, which is driven by the electron temperature T e inside the target bulk (refer to "Bulk plasma" in fig. 4 (b)).For the electron-density profile at 1.0 ps delay, however, the exponential scalelength is overlayed by a transiently occurring step of the density profile that is denoted by "2".It is caused by the thermal pressure of the coronal plasma that features a much higher temperature than the bulk plasma.The spatial distribution of the electron temperature T e at 1.0 ps delay is displayed in figure 4 (b).The electrons inside the bulk plasma feature a spatially constant temperature of 450 eV while the coronal plasma features temperatures above 10 keV.For comparison of figures 4 (a) and (b), the vertical dashed line indicates the radius of the strong increase of T e , which coincides with the step of the electron-density profile denoted by "2".As the hotter coronal plasma expands faster than the bulk plasma, the influence of the coronal plasma on the plasma densities in the range above 0.1 n 800 nm c decreases with delay and adiabatic plasma expansion most likely dominates for longer delays (compare the evolution of regions "1" and "2" of the density profiles at 1.0 ps and 2.1 ps delay).Note that the hydrodynamics simulations of the HD-RT fit considers adiabatic plasma expansion only and agreement between the density profiles of the PIC simulation and the HD-RT fit to the experiment is expected only at delays for which region "1" dominates the expansion process.It follows that, in the discussed case, the evolution of the electron temperature T e in the bulk plasma is a better suited endpoint of the comparison than a direct comparison of the density profiles.The temporal evolution of the bulk-electron temperature T e and the equivalent temperatures to the average kinetic energy T x , T y and T z are displayed in figures 4 (d) and (c).T n with n = x, y or z is calculated from the average kinetic energy of the electrons within the bulk plasma (circle with a radius of 2 µm).As the laser is polarized into the x direction and propagates along the y direction, T x and T y include the contribution of non-thermal laser-accelerated electrons.In contrast, T z is generated by collision of particles only.For 2D3V-PIC simulations, the bulk-electron temperature T e is close to the temperature component T z .Throughout this work, the quantity T e of the PIC simulations is derived by a Maxwellian fit to the electronvelocity distribution of the bulk plasma.The maximum of T x and T y of 20 keV is reached at 67 fs.Together with the approximately constant electrondensity profile between −0.1 and 0.2 ps, this indicates that the heating of the plasma bulk happens isochorically.Subsequently, as the plasma starts to expand, the temperatures T x and T y decrease while T z increases.The term thermalization refers to the circumstance that all electrons get the same Maxwellian temperature distribution into all spatial dimensions via collisions.As figure 4 (c) shows, T x , T y and T z equal each other after about 1.5 ps, i.e., the plasma thermalized within about 1.4 ps after the termination of heating.From the analytic electron-electron collision rate of hot electrons and assuming a plasma with an electron density of 30 n 800 nm c and a temperature of 20 keV we find a thermalization time [78] which is in close agreement to the PIC simulation.We emphasize that the testbed utilizes a pure hydrogen target, which allows for the comparison to analytical calculations without further approximations.Figure 4 (d) compares the evolution of the bulk-electron temperature T e of the PIC simulation (red line) and the electron-temperature evolution of the HD-RT fit to the experiment (black line).The lower and upper limit of the black line correspond the cases α and β.Up to 1.5 ps the PIC-simulation results show the process of isochoric heating and thermalization.Subsequently, T e declines because of adiabatic plasma expansion.The HD-RT fit, however, shows adiabatic plasma expansion only, which artificially starts with a heuristic initial electron temperature at 0 ps delay.Both approaches are comparable only after thermalization, i.e., at delays later than 1.5 ps.Although the trend of adiabatic cooling by plasma expansion is present for both approaches, the PIC simulation overestimates the bulk-electron temperature compared to the electron temperature of the HD-RT fit.To demonstrate the feasibility of the testbed to benchmark PIC simulations, the next section discusses the disagreement by presenting systematic scans of PIC simulations.

IV. DISCUSSION
As exemplified in the previous section, the results of the testbed are feasible to be used for quantitative comparisons to PIC simulations.In this section we discuss the effect of different parameters of the PIC simulations with respect to the chosen endpoint of the bulk-electron temperature after thermalization.References [56][57][58][59][60] show that resonance absorption dominates the heating of electrons for laser irradiances I laser λ 2 laser ≲ 1 × 10 18 Wµm 2 /cm 2 and targets with a surface-density scalelengths L 0 ≳ 0.1 λ laser = 80 nm.Vacuum heating of electrons, however, dominates for targets with a surface scalelengths L 0 ≲ 0.1 λ laser = 80 nm.Furthermore, Ref. [79] shows that the existence of preplasma changes the laser-absorption efficiency.The here-utilized peak irradiance is 1.0 × 10 18 Wµm 2 /cm 2 .Consequently, we expect a dependence of the heating of electron to the surface-density scalelength of the target.According to the experimental uncertainties, PIC simulations with different initial surface scalelengths L 0 between 0 and 500 nm are conducted.Figure 5 (b) displays the derived bulk-electron temperatures T e at 1.8 ps delay.The two limiting cases are T e = 1087 eV for L 0 = 0 and T e = 353 eV for L 0 = 500 nm, which corresponds to a variation of almost a factor of 3.An exponential fit to the data (dashed line) suggests the saturation of the temperature decrease with increasing initial scalelength.As explained in the following, the reduction of the bulk-electron temperature for increased surface scalelengths L 0 is caused by a change of the laser-absorption mechanism from vacuum heating at small L 0 to resonance absorption at higher L 0 .Figure 5 (a) shows the velocity distribution of electrons in laser-propagation direction for the cases of L 0 = 0 and L 0 = 250 nm at 20 fs delay.Both distributions feature a forward-moving electron current with velocities between ∼ 0.1 c and 0.5 c.However, in the case of L 0 = 250 nm, the overall number of forward-moving electrons is decreased, which is a signature of the transition from vacuum heating to resonance absorption.A second signature of the transition to resonance absorption at higher L 0 is given by an increased temperature of the coronal plasma (red-shaded area in fig. 4 (b)) for higher L 0 .The total temporal peak of the forward-moving electron current decreases from 32.4 kA/µm 2 at L 0 = 0 to 10 kA/µm 2 at L 0 = 500 nm.A temporally resolved analytic calculation of the electron-temperature increase by resistive returncurrent heating based on the simulated forward-moving electron current demonstrates good agreement to the increase of the bulk-electron temperature of the PIC simulations up to 70 fs delay.It follows that the decrease of the forward-moving electron current contributes to a reduction of the bulk-electron temperature as the laser-absorption mechanism changes from vacuum heating at low L 0 to resonance absorption at high L 0 .
The Appendix of this work summarizes other parameter scans and assumptions of the PIC simulations that influence the bulk-electron temperature, each on a sub-20 % level.Appendices E and F show that the assumption of an initially fully ionized target and the negligence of radiative cooling are reasonable approximations with insignificant effect on the final evolution of the bulk-electron temperature.Appendix G investigates the difference of 2D3V-PIC and 3D-PIC simulations with a high number of particles per cell.The results reveal that, at 0.4 ps delay, the bulk-electron temperature of the 3D-PIC simulation is 14 % lower than the bulk-electron temperature of the 2D3V-PIC simulation.To furthermore check the influence of the pressure gradient of hot electrons along the z axis in the experimental parameter range, Appendix J presents 3D-PIC simulations with a larger box size and lower resolution.The results demonstrate that the influence of three-dimensional effects of the hot-electron density distribution are negligible within the FWHM of the pump-laser focal spot.PIConGPU features a Coulomb-binary-collision model, which is given by the Spitzer-resistivity model of return-current Joule heating with a constant cutoff of the collision frequency at low temperatures.However, as shown in Ref. [80], even if a constant cutoff is applied in the low-temperature regime, the Spitzer resistivity strongly deviates from the Lee-More, Redmer, and Monte-Carlo results.
Therefore, a low-temperature collision-frequency correction is needed for electron temperatures lower than the Fermi temperature [30].A calculation of the relevance for the here discussed parameter range is given in Appendix H.During the leading edge of the laser pulse, at bulk-electron temperatures below 100 eV, slight deviations between the calculation and the results of the PIConGPU simulation are observed.However, as the collision rate remains unaffected at the hundreds-of-eV temperature level, the low-temperature collision correction has negligible influence on the bulk-electron temperature after thermalization.To furthermore test the influence of the analytic uncertainty of the Coulomb logarithm ln Λ of O(ln Λ −1 ), PIC simulations with a fixed Coulomb logarithm are performed in Appendix I.A decrease of the Coulomb logarithm from 15 to 3.5 decreases the bulk-electron temperature at 1 ps delay by about 20 %, which is small compared to the dependence on the inital surface-density scalelength.
In summary, systematic scans of PIC simulations show that mainly the dependence on the initial surface-density scalelength contributes to the overestimation of the bulkelectron temperature by the PIC simulation in section III.The observed transition over different heating mechanisms of electrons confirms previous work at similar laser-intensity levels [56][57][58][59][60]. Furthermore, the discussion underlines the importance of the exact initial distribution of the target density in PIC simulations that try to model experimental scenarios, which is known, for example, from laser-driven proton acceleration [26,81].

V. FUTURE PROSPECTS
The here-presented showcase of the testbed utilizes isochoric heating of solid hydrogen by an ultrashort laser pulse with a dimensionless vectorpotential a 0 ≈ 1.A simple reduction of the pump-laser energy directly leads to a reduction of a 0 .With that, the showcase demonstrates the readiness of the testbed for controlled parameter scans in experiment and simulation at all laser intensities of a 0 ≲ 1 and varied laser-pulse duration.Furthermore, the testbed is able to investigate the transition from the thermal-driven regime of plasma expansion (a 0 < 1) to the hot-electron sheath-driven regime of plasma expansion (a 0 > 1) by increasing the pump-laser energy.We expect similarities to the plasma physics of laser-heated nanoparticles that show aspects of hydrodynamic expansion together with Coulomb-explosion [82][83][84][85][86].At intensities approaching 10 22 W/cm 2 (a 0 ≫ 1) similar investigations like in the presented showcase suggest that relativistically induced transparency becomes relevant [42].In the present implementation, the testbed features probe beams at 258 nm and 515 nm wavelength.The probe beams are sensitive to plasma-density gradients around an electron density of about 0.1 n 258 nm c ≈ 1 n 800 nm c and 0.1 n 515 nm c ≈ 0.2 n 800 nm c [42].A comparison of both densities to the density profiles of the PIC simulation in figure 4 (a) shows that the shadowgraphy diagnostic does not image the density of the hot coronal plasma, which features electron densities between 0.1 n 800 nm c and 10 −3 n 800 nm c (decreasing with delay).A modification of the experimental setup to probing wavelengths in the near-infrared spectral range (e.g.2.5 µm wavelength) would allow to measure the respective densities and, by this, enable a complementary benchmark of PIC simulations with respect to non-thermal effects beyond the bulk-electron population.Cryogenic jets of different material and composition are readily available and frequently used [87].In the future, the effect of multi-species mixtures of hydrogen and deuterium will be studied [88] and cryogenic Argon-jets will allow to benchmark ionization and recombination dynamics as well as plasma opacitiy, e.g., by probing with extreme-ultraviolet backlighters [89][90][91].Finally, we would like to emphasize that the testbed is ready to be used in combination with other laser-driven secondary sources that induce isochoric heating, for example laser-accelerated ion beams [92].

VI. SUMMARY AND CONCLUSION
In summary, we introduce a testbed to experimentally benchmark PIC simulations based on laser-irradiated micron-sized cryogenic hydrogen-jet targets.
Timeresolved optical shadowgraphy by two spectrally seperated laser beams measures the temporal evolution of the plasma density.A fitting approach by hydrodynamics and ray-tracing simulations enables the determination of the bulk-electron temperature evolution after the laser energy was absorbed by the target (HD-RT fit).A showcase of the testbed studies isochoric heating of solid hydrogen by laser pulses of 37 fs duration and a dimensionless vectorpotential a 0 ≈ 1.The HD-RT fit yields a bulk-electron temperature between 250 and 300 eV after absorption of the laser energy, which is supported by systematic scans of PIC simulations.The results confirm that, due to the interplay of vacuum heating and resonance heating of electrons, an exact determination of the surface-density gradient of the target is key to achieve quantitative agreement of experiments of high-intensity laser-solid interactions with PIC simulations in the regime of a 0 ≈ 1.The showcase demonstrates the readiness of the testbed for controlled parameter scans at all laser intensities of a 0 ≲ 1, which is particularly relevant as PIC-simulation tools develop towards the inclusion of physics models at subrelativistic laser-intensity levels.By extending the platform with additional multi-color probes and by including diverse atomic target species in the future, the platform establishes a path towards a sophisticated, versatile testbed to systematically explore plasma opacity as well as ionization and recombination dynamics in laser-heated plasmas.
dynamics simulations to calculate the plasma-expansion process, which is discussed in this section.We first calculate the Knudsen number Kn. Taking the temperature T e = 300 eV, the electron density n e = 30 n 800nm With the diameter of the target L = 4.4 µm as a characteristic spatial scale of the system, the Knudsen number of the bulk plasma is which is in the applicable range of hydrodynamics equations [93].
The Debye length λ D = 1 nm gives the plasma parameter which shows that the plasma is weakly coupled.The PIC simulation in section III yields magnetic fields B between 1 T to 100 T in the single-picosecond timeframe and the hydrodynamics simulations yield ion temperatures between 10 eV and 100 eV in the tens-of-picosecond timeframe.The characteristic Lamor radius shows that the investigated plasma is weakly magnetized during all times.e is the elementary charge.We conclude that hydrodynamics simulations are feasible to calculate the plasma-expansion process.
The simulation tool FLASH is commonly used to model high-energy-density physics [94,95].Uncertainties mainly arise from the hydrogen equation of state and the assumption of two-dimensional radial symmetry.Radial symmetry is supported by the experimental observation of similar expansion in a secondary optical-probing axis antiparallel to the pump-laser axis.Lateral heat diffusion by transverse temperature gradients along the z axis is, however, not considered by two-dimensional radial symmetry.Two-dimensional cylinder-symmetric simulations in Appendix J shows that the influence is negligible in the investigated parameter range.Different equations of state are considered in Appendix D.
As we did not use any laser-plasma-interaction model for the HD-RT fit, the approach promises to be robust and versatile.The method does not depend on specific laser and target parameters and can be readily applied to other laser-target systems.To estimate the influence of radiative cooling, we use the results of the PIC simulation (section III) and calculate the corresponding energy losses from the bulkelectron-temperature evolution.There are mainly three kinds of radiation.Bound-bound radiation P Ir accounts for the radiation that is emitted from atomic ionization, free-bound radiation accounts for the radiation by electron recombination (F1) Z corresponds to the ionization state.The free-free radiation P Br accounts for Bremsstrahlung radiation by electrons.For a hydrogen-like plasma, B Br is given by The total radiation loss is given by the sum of all processes P r = P Ir + P Br + P Rr . (F3) Comparing equations F1 and F2, P Rr is is the ionization energy of hydrogen.As the relevant bulk-electron temperatures are in the multihundred-eV range, P Rr is in the percent range of P Br and is neglected in the following.For fully ionized hydrogen atoms, P Ir is zero.It follows that the total emitted power is approximately equivalent to the power of Bremsstrahlung P r ≈ P Br . (F4) In the following, we calculate the total radiation loss E r as a function of time t 0 from the temperature evolution of the PICLS simulation in figure 4 (d) by With N i = N e = 30 n 800 nm c × 1 cm 3 and Z = 1, the total radiation within 2 ps is 1.76 × 10 3 J, which is 0.14 eV per electron.Compared to the electron temperature of the PIConGPU simulation at 2 ps, the radiation loss amounts to 0.04 % only.The temporal evolution is shown by the black line in figure 7.As radiation loss scales with Z 2 N e N i , the result is in agreement with Ref. [51].The calculation is supported by a calculation with the non-local-thermal-equilibrium tool FLYCHK [100], for which we use the temperature evolution of the PIC simulation.The result is shown in figure 7. The bound-bound, bound-free, free-free and total radiation energy are given by the blue, orange, green and red line.Bound-bound and bound-free radiation contribute much less than freefree radiation, which confirms the approximation of equation F4.The total radiation loss per electron as calculated by FLYCHK is 0.13 eV.In summary, the main mechanism of radiation loss is Bremsstrahlung (free-free) radiation.Because of the low atomic number of Z = 1, the radiation loss of an electron is lower than 0.1 % of the relevant thermal electron energies and radiative cooling can be neglected.

Appendix G: 2D3V-PIC simulation versus 3D-PIC Simulation
To identify possible differences between 2D3V-PIC and 3D-PIC simulations, PIConGPU simulations with a high spatial resolution are compared in the following.To account for limited computing resources, a hydrogen target with 1 µm diameter is used.The 2D3V simulation uses a box size of 32 × 32 µm 2 (x, y) and the 3D simulation uses a box size of 8 × 8 × 12 µm 3 (x, y, z).The total length of the hydrogen target is 4 µm into the z direction.To resolve high electron densities in all dimensions, the cell size is 0.8 µm/96 and the number of macro particles per cell is 240 for both simulations.The 3D simulation is stopped at about 400 fs.For both approaches, a Maxwell-Boltzmann distribution is fitted to the electron-velocity distribution to derive the thermal bulk-electron temperature.This eliminates the The comparison of the bulk-electron temperature between the 2D3V-PIConGPU and the 3D-PIConGPU simulation is shown in figure 8. Before 100 fs, there is a high agreement between the two simulations.After 100 fs, the temperature in the 2D3V case is increased to slightly higher absolute values.Between 200 fs and 400 fs the difference between the two cases amounts to about 14 %.

Appendix H: Low-temperature-collision correction
The collision frequency of the PIC simulations does not include electron-phonon scattering, which is important at electron temperatures below the Fermi energy.However, we can use the temperature evolution as simulated by the PIC simulation, calculate the corresponding current of hot electrons, and reversibly derive the expected temperature evolution by taking the low-temperature-collision correction into account.By comparing the original temperature evolution from the PIC simulation with the artificially calculated temperature evolution, an estimate of the influence of the low-temperature-collision correction can be given.According to the hot-electron scaling in Ref. [101], the hot-electron temperature is  The three terms on the right-hand side are diffusion heating (dif f ), resistive return-current heating (res), and drag heating (dra).n e is the electron density, j h is the current density of hot electrons, n h is the density of hot electrons, T h is the average kinetic energy of hot electrons, σ Te is the electron conductivity, τ e is the collision time of bulk electrons and K Te is the thermal conductivity of bulk electrons Here, m e is the electron rest mass, e is the elementary charge, k B is Boltzmann's constant and ln Λ is the Coulomb logarithm.
In the PICLS simulation (section III), a total amount of energy of η absorb = 25 % is absorbed from the laser.The total energy of laser pulse of E laser = 160 mJ and T h = 90.5 keV gives the number of hot electrons N h that is generated by the laser-target interaction We assume that all the absorbed energy is converted into hot-electron energy.We derive a total number of hot electrons of N h ≈ 2.8 × 10 12 .By assuming a uniform distribution of hot electrons in the plasma column and in the laser spot, we calculate the hot electron density from the corresponding volume V by The current density of hot electrons is From Ref. [52] we derive the fraction of temperatures that are generated by resistive versus drag heating T res/dra and resistive heating versus diffusion heating T res/dif f : and with α = 12.9, n c,23 = 0.51 and L T = 4.4 µm.The equations show that resistive return-current heating is dominant for the heating of bulk electrons.Consequently, equation H2 is simplified to The electron conductivity is given by The electron-ion collision frequency ν ei is given by Spitzer's formula Z av e 4 m e n e (m e k B T e ) 3/2 ln Λ .(H11) Z av = 1 in the here considered case of a hydrogen plasma.
As the mean free path length of electrons f cannot be smaller than the average ion distance r 0 with the ion density n i , a cut-off at low collision frequencies is introduced v T h and v T F are the electron thermal velocity and the Fermi-temperature velocity.The electron-ion collision frequency including the cutoff correction is The low-temperature-collision correction applies to plasma temperatures lower than the Fermi temperature.Here, Spitzer's formula of electron-ion collisions (equation H11) is invalid, because the plasma is in a degenerate state.The collision frequency depends on the scattering of electrons with phonons.The corresponding collision frequency is given by [30] k s is a constant value that is estimated from experiments, T i is the ion temperature, ℏ is the reduced Plank's constant.The electron-ion collision frequency ν eic of a plasma with a temperature smaller than the Fermi temperature is The overall electron-ion collision frequency based on equations H11 to H16 is shown in figure 9 (a).For this graph, the Coulomb logarithm is set to 5 and the electron density is set to 30 n 800 nm c .The trend of the collision frequency is reversing for temperatures below the Fermi temperature.
From equation H9 we derive a formula of the returncurrent density with the temporal iteration step i of the following calculation.For each timestep, the electron temperature T e is derived from the PIC simulation (ln Λ = 5 in figure 10).The upper limit of t is set to 50 fs.The retrieved current density of hot electrons is presented in figure 9 (b).From the evolution of the hot-electron current j h , the influence of the low-temperature-collision correction on the bulk-electron temperature is calculated and compared to the PIC simulation in figure.9(c).The artificially derived electron temperature evolution with the temperature evolution of the PIC simulation.Minor differences are observed between −30 fs and −10 fs, i.e., before the laser peak arrives on target.After that and due to the rapid Joule heating, both approaches show the same results.It follows that the low-temperaturecollision correction has negligible effect on the final bulkelectron temperature.The presented 2D3V-PIC simulations simulate the temperature evolution in the x-y plane only.They ignore the dynamics caused by the hot-electron-pressure gradient in the z direction (along the jet axis).Full 3D-PIConGPU simulations are conducted to investigate the effect of the density gradient of hot electrons to the bulkelectron-temperature distribution within the laser-spot region.The size of the simulation box is 28 µm in z direction and 20 µm in x and y direction.The initialized target is a cylindrical hydrogen rod of 20 µm total length (z direction) and a radius of 2.2 µm.The grid size is 0.8 µm/24 in all directions and the number of macro particles per cell is 240.The initialized hydrogen plasma is fully ionized and features a density of 30 n 800 nm c without surface scalelength.The laser parameters are equivalent to section III.The supplementary figure S 3 shows sectional planes through the center of the target at 426 fs after the laser peak.The bulk electrons are already thermalized at this time.Temperature gradients are observed in all directions.The temperature is decreased from 1.7 keV in the laser-spot region to 1 keV above and below the laser spot (z-axis).Along the laser-propagation direction (y-axis), i.e., into the target bulk, the temperature decreases from 1.7 keV on the front to 1.6 keV on the rear side of the target.It follows that the temperature distribution in the x-y plane (supplementary figure S 3 (c)) is homogeneous within 20 %.This is consistent with the 2D3V-PIC simulation in section III, for which a uniform temperature distribution is formed after about 500 fs.
To estimate the contribution of the hot-electronspressure gradient, we calculate the dependency of the local-bulk-electron temperature on the local laser intensity and subsequently compare the calculation to the 3D-PIC-simulation result.In Appendix H we show that Joule heating by hot electrons is dominant for the heating of bulk electrons.The laser intensity on target is The average bulk-electron temperature variation along the z-axis of the 3D-PIConGPU simulation is shown in figure 11 by the blue line.It is calculated from the supplementary figure S 3 (b) by averaging along the x axis.The green and the orange curve in figure 11 show the intensity distribution of the laser I(r) 0.4 and I(r), each scaled to the maximum of the blue line.Within the FWHM of the laser (14 µm), T e of the PIC simulation coincides with the scaling I(r) 0.4 .Beyond this region, the scaling of temperature is slightly different from I(r) 0.4 .This shows that the pressure gradient of hot electrons is relevant only outside the focal-spot FWHM.As the 2D3V-PIC simulations refer to the central x-y plane of the interaction at z = 0, the influence three-dimensional effects of the hot-electron-density distribution can be neglected.

Lateral heat transfer by diffusion
Temperature gradients along the jet axis (z axis) can influence the hydrodynamic expansion of the plasma in the central x-y plane at z = 0 by lateral heat diffusion.To verify the assumption of two-dimensional radial symmetry of the hydrodynamics simulations of our experimental scenario, we conduct three-dimensional hydrodynamics simulations and compare the plasma density on the tens-of-picosecond timescale for different initial temperature gradients.
The simulations utilize FLASH and the simulation box is 20 µm in height (z-axis) and 20 µm in radius.The initial diameter of the plasma column is 4.4 µm.Following the previous subsection, three different temperature gradients are initialized and shown in the supplementary figure S 4 (a).The blue line shows a uniform temperature distribution of 300 eV, the green line shows a temperature distribution proportional to the intensity distribution of the laser I(r) with a peak temperature of 300 eV and the orange line shows a temperature distribution proportional to I(r) 0.4 , similar to the result of the 3D-PIC simulation of the previous subsection.The results of the simulations are presented in the supplementary figures (b), (c), and (d) for 10 ps, 20 ps, and 30 ps delay and at the position of the laser peak at z = 0.For comparison to the hydrodynamics simulations in section II B 1, the red-dashed line shows the results of a two-dimensional radial-symmetric simulation with T e0 = 300 eV and D 0 = 4.4 µm.All simulations show high agreement to each other.Small deviations of the density profiles occur at densities below 0.1 n 800 nm c only and amount to 5 % at maximum for the case of an initial temperature distribution proportional to I(r).The overall agreement of all three-dimensional cylindersymmetric simulations to the two-dimensional radialsymmetric simulation shows the negligible influence of lateral heat diffusion in our scenario and justifies the utilization of two-dimensional radial-symmetric hydrodynamics simulations by the HD-RT fit.

FIG. 1 .
FIG. 1. a) Experimental setup: The 37 fs-pump laser with a peak intensity of 1.6×10 18 W/cm 2 interacts with a cylindrical hydrogen-jet target.Time-resolved high-resolution shadowgraphy utilizes a single microscope and two simultaneous backlighters with a pulse duration of 260 fs (258 nm probe) and 160 fs (515 nm probe).b) Exemplary shadowgrams of a timedelay scan as recorded by the 258 nm probe.The pump-probe delay is given in the upper-right corner of each shadowgram and the shadow diameter D 258 nm E is measured between the cyan arrows.For 30 ps delay, volumetric transparency is observed.The spatial extent and fluence of the parasitic plasma emission varies from shot to shot.The colorscales and spatial scales are consistent between the images.

FIG. 2 .
FIG. 2. a) Measured shadow diameter from two-color optical shadowgraphy and HD-RT fit to the experimental data (case α, see section II B 3). b) The total χ 2 distribution between the experimental data and all HD-RT simulation data.
(b).The images are recorded by the 258 nm probe and show a sideview of the target.The pump laser impinges from the left and the spatial FWHM of the focal spot is indicated by the dashed lines on the left-hand side.Plasma emission occurs mostly from within the FWHM of the focal spot and from the position of the unperturbed target-front surface.The shadow of the target features sharp edges.The shadow diameter D 258 nm E measures the width of the shadow at the position of the pump-laser peak, as illus-trated by the cyan arrows in the shadowgrams at 0 ps and 10 ps delay.At 30 ps, the edges of the shadow are diffuse and prevent the measurement of D 258 nm E .Consequently, D 258 nm E is set to zero for this delay.The shadowgram shows penetration of probe light through central parts of the target.In the following we refer to this observation as volumetric transparency.The whole scan of shadow diameter versus pump-probe delay is shown by the markers in figure 2 (a).D 258 nm E and D 515 nmE

FIG. 3 .
FIG. 3. Results of a hydrodynamics simulation and ray-tracing simulations: a) Radial electron density of the hydrodynamics simulation with an initial setting of Te0 = 300 eV and D0 = 4 µm (case β in figure 2 (b)).The corresponding shadow diameters D 258 nm S are derived by ray-tracing simulations and are given by the dashed arrows on the top.b) Electron and proton temperature of the target bulk at zero radius versus time (case β in figure 2(b)).c) Top view of the distribution of refractive indices ñ258 nm (gray colorscale), which is calculated from the radial electron-density profile at 15 ps delay.The distribution is placed in the object plane of the ray-tracing simulation and refraction of the 258 nm-probe rays is visualized by the purple lines.d) Light distribution in the image plane of the ray-tracing simulation.The shadow diameter D 258 nm S is measured at half of the unperturbed background intensity.

FIG. 4 .
FIG. 4. a) Spatial distribution of the electron density in the PIC simulation.b) Spatial distribution of the electron temperature Te at 1 ps delay.c) Temporal evolution of the equivalent temperatures Tn (equation 6).Thermalization is reached at about 1.5 ps.d) Comparison of the PIC simulation and the HD-RT fit to the experiment.

Figure 4 (
a) shows the simulated electron-density profile for different delays.0 ps delay corresponds to the arrival of the pump-laser peak on the front side of the target.The comparison of the blue and the orange line shows that the electron density at 0.2 ps delay does not significantly differ from the initial state at −0.1 ps.The expansion of hot electrons into the surrounding vacuum occurs at densities below about 0.1 n 800 nm c

FIG. 5 .
FIG. 5.PIC simulation results: a) Electron-velocity distribution in the laser-propagation direction to the time of the maximum return current at 20 fs delay.b) Bulk-electron temperature Te at 1.8 ps delay versus the initial surface scalelength L0.The dashed line is an exponential fit to the data and serve as a guide for the eye.

2 e= 1 .
logarithm ln Λ = 3.55 , the electronelectron-collision rate (for thermalized temperatures T e ) calculates to νee = 2.91 × 10 −6 n e ln Λ T −3/04 × 10 14 s −1 .(C1) The thermal velocity of electrons is v the = k B T e m e = 7.26 × 10 6 m/s , (C2) with Boltzmann's constant k B and the electron rest mass m e .The mean free path length f of the electrons is

FIG. 6 .
FIG. 6.Comparison of the isotherms of SESAME EOS, FPEOS, FEOS and FLASH EOS of hydrogen for an electron temperature of a) 10 eV, b) 100 eV, and c) 1000 eV.Comparison of the electron-density profile at d) 10 ps, e) 20 ps, and f ) 30 ps as calculated by SESAME EOS and FPEOS.The initial conditions are Te0 = 300 eV and D0 = 4.4 µm.

FIG. 9 .
FIG. 9. a) Electron-ion collision frequency versus electron temperature with and without low temperature correction.b) Current density of hot electrons as calculated by equation H9 from the PIConGPU simulation with a Column logarithm of 5 in figure 10. c) Electron-temperature evolution from PIConGPU and as calculated from the current density in fig.(b) by taking into account the low-temperature-collision correction.

Appendix J: Lateral heat transfer 1 .
FIG. 11.Average bulk-electron-temperature distribution along the z axis (jet axis) in the 3D-PIConGPU simulation at 426 fs (blue line).The green and the orange line are the laser-intensity distribution I(r) and I(r) 0.4 , each scaled to the maximum of the blue line.

r 2 , 2 ( 2 0 2 − 1 ) m e c 2 ≈ a 2 0 4 m e c 2
(J1)    with the laser-spot radius r.The hot-electron temperature is denoted as T h (r).Assuming a constant absorption coefficient η, the density of hot electrons n h is approximated by equation H4,n h T h (r) dV = η I(r) τ 0 dS .(J2)Here, dV and ds are the differential volume and area and τ 0 is the laser-pulse duration.It follows thatn h ∝ I(r) T h (r) .(J3)For v h ≪ c, the current density of hot electrons isj h (r) = e n h v h (r) = c e n h 2 k B T h (r)m e c 2 H9, H10, and H11 we haveT e (r) ∝ j h (r) 2 σ Te ∝ j h (r) 2 T e (r) − 3 J6)andT e (r) ∝ j h (r) 0.8 .(J7)T e (r) is the bulk-electron temperature.For a 2 0 /2 ≪ 1, equation H1 givesT h (r) = ( 1 + a (J8)with a 0 ∝ I(r).We deriveT h (r) ∝ I(r) (J9)and together with equation J5j h (r) ∝ I(r) .(J10)Finally, from equation J7 we derive the proportionality of the bulk-electron temperature and the intensity distribution of the laser T e (r) ∝ I(r) 0.4 .(J11)