The importance of temperature-dependent collision frequency in PIC simulation on nanometric density evolution of highly-collisional strongly-coupled dense plasmas
Mohammadreza Banjafar Lisa Randolph Lingen Huang S. V. Rahul Thomas R. Preston Toshinori Yabuuchi Mikako Makita Nicholas P. Dover Sebastian Göde Akira Kon James K. Koga Mamiko Nishiuchi Michael Paulus Christian Rödel Michael Bussmann Thomas E. Cowan Christian Gutt Adrian P. Mancuso Thomas Kluge Motoaki Nakatsutsumi








mohammadreza.banjafar@xfel.eu
European XFEL, Holzkoppel 4, 22869, Schenefeld, Germany


Department Physik, Universität Siegen, Walter-Flex-Str. 3,  57072, Siegen, Germany


Helmholtz-Zentrum Dresden-Rossendorf, Bautzner Landstraße 400, 01328, Dresden, Germany


European XFEL, Holzkoppel 4, 22869, Schenefeld, Germany


European XFEL, Holzkoppel 4, 22869, Schenefeld, Germany


Japan Synchrotron Radiation Research Institute (JASRI), Sayo, Hyogo, 679-5198, Japan
RIKEN SPring-8 Center, Sayo, Hyogo 679-5148, Japan


European XFEL, Holzkoppel 4, 22869, Schenefeld, Germany


QST-Kansai, KPSI,8-1-7 Umemi-dai, Kizugawa-city, Kyoto, 619-0215, Japan 


European XFEL, Holzkoppel 4, 22869, Schenefeld, Germany


QST-Kansai, KPSI,8-1-7 Umemi-dai, Kizugawa-city, Kyoto, 619-0215, Japan 


QST-Kansai, KPSI,8-1-7 Umemi-dai, Kizugawa-city, Kyoto, 619-0215, Japan 


QST-Kansai, KPSI,8-1-7 Umemi-dai, Kizugawa-city, Kyoto, 619-0215, Japan 


Technical University Dortmund, Otto-Hahn Straße 4, 44227 Dortmund, Germany


Technical University Darmstadt, Schlossgartenstraße 9, 64289 Darmstadt, Germany


Helmholtz-Zentrum Dresden-Rossendorf, Bautzner Landstraße 400, 01328, Dresden, Germany


Technical University Dresden, 01069 Dresden, Germany
Helmholtz-Zentrum Dresden-Rossendorf, Bautzner Landstraße 400, 01328, Dresden, Germany


Department Physik, Universität Siegen, Walter-Flex-Str. 3,  57072, Siegen, Germany


European XFEL, Holzkoppel 4, 22869, Schenefeld, Germany
Department of Chemistry and Physics, La Trobe Institute for Molecular Science, La Trobe University, Melbourne, Victoria 3086, Australia


Helmholtz-Zentrum Dresden-Rossendorf, Bautzner Landstraße 400, 01328, Dresden, Germany


European XFEL, Holzkoppel 4, 22869, Schenefeld, Germany
Open and Transdisciplinary Research Institute, Osaka University, Suita, Osaka 565-0087, Japan




    Particle-in-Cell (PIC) method is a powerful plasma simulation tool for investigating high-intensity femtosecond laser-matter interaction. However, its simulation capability at high-density plasmas around the Fermi temperature is considered to be inadequate due, among others, to the necessity of implementing atomic-scale collisions. Here, we performed a one-dimensional with three-velocity space (1D3V) PIC simulation that features the realistic collision frequency around the Fermi temperature and atomic-scale cell size. The results are compared with state-of-the-art experimental results as well as with hydrodynamic simulation. We found that the PIC simulation is capable of simulating the nanoscale dynamics of solid-density plasmas around the Fermi temperature up to ∼2 ps driven by a laser pulse at the moderate intensity of 10^14-15 W/cm^2, by comparing with the state-of-the-art experimental results. The reliability of the simulation can be further improved in the future by implementing multi-dimensional kinetics and radiation transport.
    


The importance of temperature-dependent collision frequency in PIC simulation on nanometric density evolution of highly-collisional strongly-coupled dense plasmas
    Motoaki Nakatsutsumi
    April 29, 2024
==================================================================================================================================================================





§ INTRODUCTION



    The interaction of sub-picosecond, high-intensity short-pulse lasers with solids at moderate laser intensities (10^14-15 W/cm^2) generates strongly-coupled dense plasmas with a few eV electron temperatures, i.e., the electron thermal energy is comparable to the Fermi energy <cit.>. Under such conditions, atomic-scale collisions and quantum effects such as electron degeneracy and ion correlations play an essential role in electron dynamics <cit.>. Such plasmas can be found inside planets <cit.>, inertial confinement fusion (ICF) plasmas  <cit.>, femtosecond laser ablation for micro-processing <cit.> and early stage of ultra-high intensity laser-solid interaction at election relativistic regime. 
    
    The electromagnetic Particle-in-Cell (PIC) simulation method is widely used to model high-intensity laser-matter interactions <cit.>. PIC solves Maxwell's equations on a grid system and applies relativistic equations to all particles, each representing an average of many real particles <cit.>. This allows the evaluation of particle transport and kinetic effects in plasmas. To account for binary interactions between particles, a Monte Carlo model is used in the PIC algorithm <cit.>. Recently, an experimental platform to benchmark the collision model in the PIC simulations at the regime of sub-relativistic laser intensity of ∼10^18 W/cm^2 has been developed  <cit.>. However, the reliability of PIC simulations in highly collisional degenerate regimes was always in doubt. Challenges include uncertainties in the electron collision frequency models, unreliability of the classical approximation for the momentum transport cross-section, i.e., Coulomb logarithm <cit.>, and the necessity to consider quantum effects and large-angle scattering and macro-fields <cit.>. Among these challenges, the uncertainty of the electron collision frequency (ν) is at most important, which governs the laser energy absorption <cit.>, electron-to-ion energy transfer  <cit.>, heat diffusion <cit.>, and subsequent surface expansion and ablation. 

    The electron collision frequency around the Fermi temperature (T_F) and its effects on short-pulse laser absorption and electron thermalization has been studied experimentally <cit.> and theoretically <cit.>. Fisher et al. <cit.> demonstrated that it is essential to account for both electron-electron and electron-phonon collisions in the determination of the electron momentum relaxation rate in the interaction of femtosecond laser pulses with an intensity of I⩽10^15 W/cm^2 in metallic solids <cit.>. Mueller and Rethfeld <cit.> derived Boltzmann collision integrals to study electron excitation, thermalization, and electron-phonon coupling under non-equilibrium electron distribution. Petrov et al. <cit.> calculated the thermal conductivity and electron-ion heat transfer coefficients for the interaction of short-pulse lasers and solids in the degenerate plasma regime <cit.>. Meyer-ter-Vehn and Ramis derived a global analytical expression for the effective collision frequency ν(T_e, E_ph, n_e) <cit.>. Here, T_e, n_e, and E_ph are electron temperature, electron density, and absorbed photon energy, respectively. Furthermore, Pineau et al. <cit.> modeled the effective electron collision frequency over the entire temperature range of solid-to-plasma transitions in ICF experiments.
    
    In this article, we implemented an improved temperature-dependent binary Coulomb collision model around the Fermi temperature (T_F) in one-dimensional with three-dimensional velocity space (1D3V) PICLS code <cit.> to study nano-scale density dynamics and heatwave propagation inside femtosecond (fs) laser-excited solid-density plasmas. Angstrom-scale cell size is also implemented to resolve atomic-scale electron collisions. The results are then quantitatively compared with the recent experimental results. We also discuss the comparison with 1D MULTI-fs hydrodynamic code <cit.>. 

    

§ TEMPERATURE-DEPENDENT ELECTRON COLLISION FREQUENCY AT LOW TEMPERATURES

    Binary Coulomb collision between weighted particles has been implemented in PICLS <cit.>, based on the Takizuka and Abe model <cit.>. Initially, particles are grouped in each cell, and then the code randomly determines which pairs of particles are likely to collide. In the next step, the alterations in the velocity of particles caused by the binary collision during the time interval Δ t are computed and applied to advance the particle velocities. In an ideal plasma (high temperature and low density), the collision frequency of the small-angle scattering between two sample particles α and β can be obtained using the Spitzer formula as  <cit.>;
  
    
    ν_plasma = 4 π (q_αq_β)^2 n_l lnΛ/p_rel^2 v_rel,

    
    where q_α and q_β are the charge of the particles, p_rel = m_αβ v_rel is the relative momentum between two particles, m_αβ is the reduced mass, v_rel is the relative velocity, n_l is the lower density among n_α and n_β, and lnΛ = ln(λ_De/λ_dB) is the Coulomb logarithm. 
    Here, λ_D_e = (k_B T_e/4π n_e e^2)^1/2 is the Debye length and λ_dB=h/2m_e v_e is de-Broglie's wavelength, where k_B is the Boltzmann constant, e and v_e are the electron charge and velocity, respectively, and h is the Planck's constant.
    Considering the thermal velocity of electrons v_th = √(3k_BT_e/2m_e), where m_e is the electron mass, the denominator in Eq. <ref> is proportional to T_e^3/2. 
    Therefore, the Spitzer collision frequency diverges at low T_e. To avoid this, the Fermi temperature T_F was considered as the threshold temperature for the degenerate plasma regime, below which the collision frequency was set constant as ν_const = 4m_eZe^4lnΛ/3πħ^3 <cit.> (see the orange line in Fig <ref>), where Z is the ionization state and ħ is the reduced Planck constant.  

    
    To improve the model at low-temperature and high-density plasmas around and below T_F, we implemented the temperature-dependent collision frequency at T ≤ T_F into PICLS based on the model using an ad hoc interpolation as  Ref. <cit.>:
    
    
    ν_interp^-1 = ν_plasma^-1 + ν_metal^-1.

    
    Eq. <ref> smoothly bridges the collision frequency at high T ideal plasmas (ν_plasma) and cold metallic solids ν_metal (T_e << T_F). The ν_metal is the summation of the electron-electron ν_ee and electron-phonon ν_el-ph contribution:
    
    
    ν_ee_metal(T_e) = E_F/ħ(T_e/T_F)^2,
       ν_el-ph(T_i) = K_se^2 k_B T_i/ħ v_F,

    
    where E_F = k_B T_F = m_e v_F^2 / 2 is the Fermi energy, v_F = ħ (3 π^2 n_e)^1/3 / m_e is the Fermi velocity, T_i is the ionic temperature, and K_s is an empirical parameter which is adjusted using experimental observations <cit.>.
    At T_e >> T_F, the Spitzer formula (Eq. <ref>) includes all types of collisions; electron-electron (ν_ee_plasma), electron-ion (ν_ei_plasma), and ion-ion (ν_ii_plasma). The Coulomb logarithm in Eq. <ref> is also modified to be <cit.>
    
    
    lnΛ = 1/2ln[1 + max(1,b_max/b_min)]

    
    where b_max and b_min are upper and lower cutoffs on the impact parameter for Coulomb scattering <cit.>. Eq. <ref> ensures the avoidance of negative values for the Coulomb logarithm at low temperatures. The b_max is set by the screening effect or Debye-Hückel <cit.>. Considering the effect of electron degeneracy, the Debye-Hückel length of an electron distribution can be approximated by

    
    1/b_max^2 = 1/λ_DH^2    = 1/λ^2_De + 1/λ^2_Di,
        = 4π n_e e^2/k_B(T_e^2 + T_F^2)^1/2 + 4π n_i Z^2 e^2/k_B T_i,


    where λ_Di is the ionic Debye shielding length, and n_i is the ion density <cit.>. In dense plasmas, the Eq. <ref> breaks down as the screening length (λ_DH) becomes less than inter-atomic distances R_0 = (3/4π n_i)^1/3 <cit.>. To avoid this nonphysical phenomenon, R_0 was set as the minimum cutoff for the electron mean free path λ_e <cit.>. Accordingly, the cutoff collision frequency is given as:
    
     
    ν_cutoff = electron  total  velocity/ion  sphere  radius = √(v_th^2+v_F^2)/R_0.

     
     The lower cutoff b_min can be set by λ_dB at high energy, according to the Heisenberg uncertainty principle, or by the classical distance of closest approach Ze^2/3k_BT_e at low energy <cit.>. Taking into account the electron degeneracy and a partial ionization in high-density degenerate plasmas, b_min can be reasonably approximated by <cit.>
    
    
    b_min = max[ h/2√(3m_ek_B √(T_e^2+T_F^2)), Ze^2/3k_B√(T_e^2+T_F^2)].

     
         

     The blue and red solid lines in Fig. <ref> a) illustrate the temperature-dependent interpolated collision frequency ν_interp for electron-ion (ei) and electron-electron (ee) collisions implemented in 1D PICLS code for the pre-ionized aluminium (Al^5+). T_F was considered as the threshold temperature for the degenerate plasma regime, and ν_el-ph was calculated assuming the local thermodynamic equilibrium (LTE) condition T_e = T_i. 
     
     The evolution of the electron and ion temperature of Al^5+ plasmas using our new collision frequency model in PIC simulation is summarized in Fig. <ref> b) with solid lines. The initial electron and ion temperatures were set to 58 and 0.5 eV, respectively, and the ionization was turned off. The system is confined between periodic boundaries, and no temperature gradient or heat diffusion is assumed. The thermalization dynamics are compared with the simple analytical calculation (dashed lines), where the relaxation rate is determined by the difference in electron and ion temperatures  <cit.>:

     
    dT_e/dt = -n_i/n_eT_i-T_e/τ_ei,
       dT_i/dt = T_e-T_i/τ_ei,


    where τ_ei is the characteristic time of ion heating <cit.>;
     
     
    τ_ei = m_i/2m_e ν_ei,

     
     with m_i denoting the ion mass. The electron-ion collision frequency was chosen as ν_ei = 1.8 × 10^16 s^-1, to be consistent with its value in PIC simulation at 58 eV. Regarding the ion heating, a notable agreement between the PIC simulation and the analytical predictions is observed. However, the electron temperature exhibits a more rapid decay in the PIC simulation due, likely, to the radiation loss in the PIC simulation, a topic to be discussed in greater detail in the subsequent section.


     It is worth emphasizing that although the collision model correction is demonstrated in a one-dimensional PIC code, the momentum transfer accounting that results from collisions is performed in a comprehensive three-dimensional velocity space.
        


§ SIMULATION RESULTS

    We tested the constant ν_const and interpolated ν_interp collision frequency models in a 1D PIC simulation of the laser-excited Al thin foil. The simulation includes ionization dynamics and non-LTE conditions. This pure PIC simulation does not incorporate quantum effects or electron degeneracy within the algorithm of the PIC framework. The sample is a 200 nm thick aluminium (Al^3+) foil at ambient temperature initially (T_e=T_i=300 K). A p-polarized laser with 50 fs duration at full-width half maximum (FWHM), 800 nm wavelength, and 6 × 10^14  W/cm^2 peak intensity irradiates on the Al sample surface.  
    To resolve λ_e ∼ 3 Å around T_F, the cell size was set to be Δ x = 1.25 Å corresponding to the time-step of Δ t = 4.16 × 10^-19 s. Each cell initially contains 12 ions (Al^3+) and 36 free electrons. We also used the fourth-order particle shape to minimize the numerical noise <cit.>. The ion number density was set to the density of Al (n_Al = 34.6n_c), where n_c = m_e ω_L^2/(4 π e^2) = 1.742 × 10^21cm^-3 is the critical plasma density at the laser wavelength of λ_L = 800 nm. Here, ω_L is the laser angular frequency. The ionization dynamics due to the intense laser EM-field and the electron-ion collisions were modeled by the Landau-Lifshitz field and direct-impact (DI) ionization models, respectively.

    
    For the comparison, we also performed a similar simulation using the 1D MULTI-fs hydrodynamic simulation code, specifically designed for short (≤ ps) pulse high-intensity laser-solid interactions below < 10^17 W/cm^2 laser intensity <cit.>. MULTI-fs calculates an explicit solution to Maxwell's equations for the fs laser pulse interactions with steep plasma density gradients. It includes the temperature-dependent collision frequency as well as a thermal conductivity from metallic solids up to high-temperature ideal plasmas and a separate equation-of-state (EOS) for electrons and ions (two-temperature model). We used EOS and ionization tables for Al provided by the MULTI-fs package. 500 cells were used for simulating the 200 nm thick Al foil. Finer cells were used closer to the front and rear surfaces compared to the middle of the foil to reduce the numerical error. The electron collision frequency was modeled using the same interpolated collision frequency model as implemented in the PICLS code. The radiation transport module was switched off, and the heat flux inhibition parameter was set to f = 0.1, which is defined as the ratio between the actual heating flux and the free streaming flux. The latter refers to the heat flow carried by electrons when they travel freely at their thermal velocity along the temperature gradient <cit.>. The value of f equal to 0.1 has been chosen based on comparisons between hydrodynamic code calculations and experiments on heat wave propagation into solids, which have shown good agreement for this value of f <cit.>. 

    
    
     
    Fig. <ref> a) illustrates the amplitude of the transverse oscillating electric field of the laser (E_y) as a function of the position in 1D PIC simulations with ν_const and ν_interp models when the peak laser intensity reaches the Al foil front surface on the left side. The grey shaded area indicates the foil. For the ν_const case, there is considerable leakage of the laser field due to the overestimation of ν at low temperatures. This is caused by the suppression of reflection due to excessive collisions  (λ_e < 1 nm), which prevent electrons from coherent oscillations in the laser field. This increases the collisional skin depth as
    
    
    δ = i/√(ε)c/ω_L,

     
     where ε is the dielectric function of the plasma medium defined by <cit.>
     
     
    ε = 1-ω_p^2/ω_L^2(1+iν_ei/ω_L).

     
    Fig. <ref> b) illustrates the skin depth as a function of ν_ei, for the propagation of an EM-wave with a wavelength of 800 nm inside a solid-density Al^+3 plasma. The skin depth increases with ν_ei, resulting in a more extended penetration of the EM field into solids. The horizontal dashed line represents the collisionless skin depth (δ_collLess = c/ω_p). On the other hand, the laser EM-field is highly shielded by high-frequency plasma oscillations in the case with ν_interp. As a result, the laser is reflected more from the surface and is highly damped inside solids.

    
    
    
    Fig. <ref> summarized the retrieved electron (a-c) and ion (d-f) temperature profiles inside the Al foil as a function of time and depth in 1D PIC with ν_const (left), ν_interp (middle), as well as the MULTI-fs (right). The zero time and depth correspond to the peak intensity of the laser pulse and the initial front surface of the target, respectively. Due to the short duration of the laser pulse and the high density of the Al foil target, the surface expansion during laser irradiation is almost negligible. The laser's energy is deposited primarily at a solid density and instantly elevates the electron temperature at the surface, resulting in well-separated electron and ion temperatures. In this laser intensity regime, the main absorption mechanism in metallic solids is the inverse Bremsstrahlung, which directly depends on the collision frequency <cit.>. The laser absorption observed in the PIC simulation with ν_interp is around 29% and 35% in MULTI-fs, which are close to experimental measurements reported in <cit.>. To compare with the theory, the absorption coefficient of a solid surface A, with a reflection coefficient of R, irradiated by a linearly polarized monochromatic EM wave is given by <cit.> 
    
    
    A = 1 - R,      R = |1-√(ε)/1+√(ε)|^2.

     
    For ν_ei≈ 1. × 10^16 s^-1, the Eq. <ref> yields the absorption coefficients of A ≈ 0.3 for the ionization state in PIC with ν_interp (Al^6.5+) and A ≈ 0.35 for ionization state in MULTI-fs (Al^4.5+), which agree with values measured by simulations. In PIC with ν_const, 28% of the incident light was transmitted, and 35% was reflected from the front surface. Thus, one can estimate a 37% laser energy deposition into the target. The higher absorption for ν_const can be attributed to the broader skin depth area and absorption of the laser in deeper parts of the target. Both PIC simulations also show a faster surface expansion due to the higher temperature and, therefore, higher pressure (P=P_e+P_i=k_B(n_eT_e+n_iT_i)) inside the Al foil, compared to the MULTI-fs. 
    
    Comparing the electron temperature T_e in the PIC simulation with ν_interp and MULTI-fs (b and c), one can see a generally good agreement in the surface heating and subsequent heat transport into the target. During laser excitation, the surface temperature reaches around 40 eV in MULTI-fs, and 50 eV in PIC with a spatial gradient inside the target of size ∼ 18 nm in both simulations, which agrees with the collisional skin-depth c/ω_p√(ν_ei/2ω_L) <cit.>, for Al^3+ and ν_ei=10^16 s^-1. The higher T_e in the PIC simulation at around 0 ps can be attributed to kinetic effects (e.g., electron oscillations in the laser field), which are excluded in the MULTI-fs simulation. The temperature gradient is quickly smoothed out by heat transfer via thermal conduction. Bulk heating due to the ballistic transport of oscillating electrons does not play an essential role here due to high collisionality with sub-nm λ_e. Therefore, the dominant energy transport is driven by the heat diffusion expressed by the two-temperature energy conservation equations as
     
    
    C_e ∂ T_e/∂ t = -∇·q - γ(T_e - T_i) + Q(z,t),
        C_i ∂ T_i/∂ t = γ(T_e - T_i),

    
    where C_e, C_i are, respectively, the volumetric heat capacity of electrons and ions, q = κ∇ T_e describes the electron thermal heat flow with κ(ν) being the heat conductivity, and Q(z,t) = ∂ I_abs⁄∂ z is the power density deposited by the laser with I_abs being the absorbed laser flux. The energy transfer rate from electrons to ions is expressed by γ = C_i τ_ei^-1. 

    
    
    
    Both simulations also exhibit a similar depth of heatwave propagation. However, PIC with ν_interp shows a relatively faster electron cooling than MULTI-fs. In PIC, a considerable part of such fast electron cooling can be attributed to the collisional ionization process (roughly 28% of the absorbed laser energy is used to ionize bound electrons). We also observed that almost one-third of the electron kinetic energy was left from the simulation box due to radiation energy loss via Bremsstrahlung. The blue line in fig. <ref> illustrates the spectrum of radiation loss at the rear side of the target in PIC with ν_interp. The spectrum peak intensity is around 1.5×10^16 Hz, which coincides with the Planckian electromagnetic energy density ρ_ν at T_e = 22 eV, as shown by the red area. This temperature is close to that observed in PIC, up to 1 ps. In addition, one can see a lower radiation power loss at frequencies below 10^16 Hz compared to the Planckian distribution due to plasma shielding.
    
    In PIC simulation with ν_const, the overestimated collision frequency and the excessive leakage of the laser field lead to a higher temperature (4-5 eV) in deeper parts of the target during laser excitation, compared to two other simulations. The overestimation of collisional heating also results in a higher surface temperature than in other cases. However, the exceeding electron collision frequency leads to very limited heat conduction into the target after the laser excitation <cit.>. We also observed a higher radiation loss compared to the case with ν_interp, resulting in faster surface cooling.
    It is useful to mention that the highest level of agreement with MULTI-fs simulations was achieved for the empirical parameter K_s = 78.8 for electron-phonon collisions in the PIC simulation. We observed much-limited heat conduction for lower K_s values in PIC simulations.

    
    The evolution of the ion temperature T_i via electron-ion energy transfer is summarized in Figs. <ref> (d-f). To focus on the heating via collisions, only the transverse kinetic energy, i.e., perpendicular to laser propagation is shown. Here, one can also see that the overall profile of the ion temperatures agrees between MULTI-fs and 1D PICLS with ν_interp. The maximum ion temperature of around 10 eV is reached in both simulations. It is also visible that the PICLS shows a faster ion cooling, which is likely to be attributed to the overestimated radiative energy loss. On the other hand, PICLS with ν_const (Fig. <ref> d)) shows a faster ion heating and faster cooling. The ion temperature quickly reaches above 10 eV near the surface and 4 eV in deeper parts already during laser excitation due to the overestimation of the ν_ei. However, the temperature decays very quickly via excessive collisions and radiative loss.

    


§ COMPARISON WITH THE EXPERIMENT AND DISCUSSION
  
    In this section, we discuss how the improved collision model in PIC simulation better describes the nanoscale solid-density plasma dynamics around the Fermi temperature by comparing it with state-of-the-art experimental results  <cit.>.

    
    
    The experiment was performed at the EH6 station at the SACLA X-ray free-electron laser (XFEL) facility in Japan <cit.>. 
    A schematic of the experimental setup is shown in Fig. <ref> a). 
    A p-polarized Ti:sapphire short-pulse laser with 800 nm central wavelength, 70 mJ pulse energy, 40 fs FWHM duration was irradiating with  17^∘ incident angle from the surface normal on multilayer (ML) samples. The ML samples consisted of 5 layers of tantalum (Ta) and copper nitride (Cu_3N) 4.3 / 11.5 nm thickness each. After a variable delay, an x-ray probe pulse with 7 fs duration and photon energy of 8.81 keV was focused to a 4 μm FWHM on the sample at a grazing incidence angle of 0.64^∘, which is slightly larger than the critical angle of total external reflection. Scattered photons were collected by a 2D area detector, multi-port charge-coupled device (MPCCD) <cit.>. The strong specular peak at the exit angle equal to the incidence angle was blocked. Due to the grazing incidence geometry, the x-ray footprint on the surface was elongated in one dimension to ∼ 360 μ m. Accordingly, the laser beam was defocussed to ∼ 500 μ m diameter to cover the x-ray footprint, resulting in a laser intensity of 3.6±0.2×10^14 W/cm^2. The x-ray scattering, as a consequence of the large probe area, represents a time-averaged signal of about 1.2 ps <cit.>. 

        
    Fig. <ref> b-e) summarized the lineout of the x-ray scattering intensity along the plane spanned by the incoming x-ray beam and the target surface normal, denoted as Q_z in Fig <ref> a), at different delays with respect to the laser intensity peak. 
    The circular dots represent the experimental data, while the solid lines are fits. This so-called in-plane scattering signal contains momentum transfers dominantly along the depth direction (z), which is closely associated with the specular reflectivity curve <cit.>. The Bragg-like peak at Q_z = 1.33  nm^-1 represents the typical length scale, i.e. 15.8 nm thickness of each Ta/Cu_3 N double-layer. The small peaks between Q_z = 1.0  nm^-1 and 1.33  nm^-1, are known as Kiessig fringes <cit.>, resembling the number of double-layer repeats in the ML sample. The scattering signal at the exit angle below the incident angle Q_z < Q_specular is particularly surface sensitive, where the evanescent x-ray wave travels parallel to the surface. The peaks at 0.80  nm^-1 and 0.87  nm^-1 (dashed vertical purple lines in Fig. <ref> b-e)) correspond to the solid densities of Cu_3 N and Ta, respectively, so-called Yoneda peaks <cit.>. Those signals originate primarily from the interference of the topmost surface layer and serve as a sensitive marker for the structure of the uppermost layers. 

    
    
    The change in the scattering signal is already visible at 0 ps (which is the integration over ± 0.60 ps around the laser intensity peak) compared to that of the cold sample in the Yoneda peaks region, but also in the number of Kiessig fringes and the Bragg peak intensity. Those modulations are more pronounced in later times; 1 and 2 ps. 
    By comparing these x-ray scattering signals with the reconstructed electron density profile, shown as grey-shaded areas with the dashed black lines in Fig. <ref>, this can be understood as the surface expansion (ablation) of the top Ta layer into vacuum and modulation of inner layers which intermix to each other. The details on the reconstruction from the x-ray scattering pattern to the real-space density profiles are discussed in detail in ref. <cit.>. Here, we focus on the comparison between the experimental data and simulations.  
    Note that in our experiment, x-ray scattering was originated mostly from the bound electrons due to the high x-ray photon energy (8.81 keV) and relatively low mean ionization Z_mean = 7<cit.>. The x-ray scatters from all electrons with binding energies below the photon energy, which corresponds to 27 electrons for Cu and 63 for Ta. This means we track here mostly the bound electrons or, equivalently, ion motion. We applied the same criteria for plotting the electron density (effective electron density: n_eff) for our simulations. 
      
    Fig. <ref> shows the comparison between the reconstructed density profiles from the experiment (the grey area with dashed black lines), PICLS using ν_const (red line) and ν_interp (blue line) and MULTI-fs (orange line). Each simulation result is averaged over 1.2 ps to be consistent with the time resolution of the experiment. 
    In Multi-fs, the simulation box included 1208 cells, with 833 cells for the ML and 375 cells for the Al substrate. Finer cells were used close to the layer interfaces to avoid the numerical noise. The EOS and ionization data for Ta and Cu were produced using FEOS code <cit.>. The radiation transfer module was switched off in the simulation. The free streaming limiting factor f of 0.6 was used here instead of f=0.1 used above, to be consistent with our previous publication <cit.>. We also verified that the selection of f between 0.1 and 0.6 led negligible impact on overall density dynamics. The incident angle of the laser was set to 17 from the target normal to be consistent with the experiment. For the substrate material, instead of silicon (Si) used in the experiment, aluminium (Al) was used in our simulations because of its well-known EOS. As the energy transport occurs from the ML to the substrate, no meaningful difference in ML dynamics between these two material choices is expected. Only the expansion of the last Ta layer into the substrate can be slightly different between these two materials due to the difference in their densities by about 20%. 

    In PICLS, each cell contains 15 virtual ion particles with initial charge states of 2, 1, and 3 for Ta, Cu, and Al, respectively. The fourth-order particle shape with different particle weightings for different ion species is used. The ion number densities are set to realistic densities of Ta (n_Ta=31.8n_c), Cu (n_Cu=48.7n_c), and Al (n_Al=34.6n_c). The ionization dynamics are modeled using the field and direct-impact ionization models. Here, the incident angle of the laser was set to normal to the surface, as the oblique incidence is not supported in 1D PICLS. We verified from Multi-fs that the difference in laser absorption between normal-incidence and 17 incidence was less than 1%, as has been reported previously <cit.>. 
    For both simulations, the laser with τ=50 fs duration, λ_L=800 nm wavelength, and 4×10^14W/cm^2 intensity impinges on the surface of the ML target from the left side. 

    It is important to note that the advantage of using the PIC here compared to the hydrodynamic simulation is, that despite its higher computational cost, the PIC allows simulating the particle exchanges between adjacent nm-scale layers, which leads to layer intermixing. This effect is excluded in Lagrangian hydrodynamic simulations as the boundary between cells needs to be preserved <cit.>. 

    
   
      
    The electron temperature (T_e) profiles in simulations at 0.1 and 1 ps after the laser-excitation are shown in figure <ref>a) and figure <ref>b), respectively. Here, different colors represent different simulations (same as fig. <ref>). PIC simulations show a higher surface temperature than the MULTI-fs due to higher laser absorption, as already discussed in the previous section. The laser absorption in PIC with ν_interp was 43%, which is in agreement with experimental measurements for Ta, reported in ref. <cit.>, and was 29% in the MULTI-fs. In PIC with ν_const, 37% of the laser energy was reflected from the surface, and around 13% was transmitted through the target, including the ML and 200 nm Al substrate. This means around 50% of the laser energy was deposited in the target. The overestimated laser propagation length into the target leads to a higher T_e in deeper parts (between 50 and 80 nm) of the ML target and the Al substrate. Furthermore, both PIC simulations here also show faster electron cooling than that observed in MULTI-fs due to excessive radiation energy loss.

    
    The highly localized surface heating instantaneously increases the pressure in the surface region. As a result, the top Ta layer expands rapidly into the vacuum within the first 1 ps, as can be seen in the experiment, PIC, and MULTI-fs simulations. The density modulation of the ML target starts from the surface and moves towards the substrate due to heat diffusion. The whole ML structure appears to be significantly modulated after ∼1 ps as indicated both by the experiment and simulations. 
    If we consider the position of 1/e of the surface temperature to be the position of the heat wavefront, the heatwave velocity of around 4-5×10^4 m/s is observed in both MULTI-fs and PICLS with ν_interp. This value was in good agreement with the ferocity of the density perturbation front observed in the experiment.
    
    In terms of the overall density dynamics, the PIC simulation with ν_interp shows a better agreement with the experiment compared to the PIC simulation with ν_const. In MULTI-fs, the initial ML structure persists due to the lack of intermixing between Lagrangian cells. 
    
    The ablation of the solid-density surface develops over time as the pressure wave travels into the target. If we define the ablation surface as the position of half of the initial solid Cu density, the PIC simulation with ν_interp predicts an ablation speed of ∼10 nm/ps up to 2 ps. This is in good agreement with the experiment and also with the sound speed at T_e=20 and T_i=5 eV, assuming an ideal gas EOS <cit.> 
    
    
    C_s = ( γ_e Z_mean k_B T_e + γ_i k_B T_i/m_i)^1/2

    
    where Z_mean is the mean ionization, γ_e and γ_i are adiabatic index of electrons and ions, respectively. M_i=181 a.u. and 63 a.u. is the ion mass for Ta and Cu, respectively, and γ = 1+2⁄n with n being the number of degrees of freedom. 
    In most cases, γ_e = 1 and γ_i = 3 can be used <cit.>. For (T_e,T_i)=(20,5)  eV and Z_mean = 5, Eq. <ref> yields C_s _ Ta∼ 8  nm⁄ps and C_s _ Cu∼ 13  nm⁄ps, respectively. 
    
    However, the PIC simulation with ν_const shows less ablation due to the faster surface cooling. 
     
    
     
    Furthermore, we investigated the impact of the selection of the grid size on PIC simulation under the strongly coupled plasma regime.
    Fig. <ref> shows the result obtained with the larger grid size of 6 Å, which is larger than the electron mean free path at T_e=10 eV; as the total electron velocity of v_e=√(v_th^2+v_F^2)=2.7×10^6 m/s and the collision frequency ν≈ 1. × 10^16 s^-1 leads λ_e = 2.7 Å. This resulted in lower collisional absorption of 32% compared to the case of 1.25 Å  grid size with 43% absorption. This led to lower T_e, thus slower subsequent density modulation. Consequently, the agreement with the experiment is lost, as shown in <ref> b-c). We found the maximum cell sizes should not exceed 1.5 times the electron mean free path. Also, the higher number of ions per cell results in more accurate and stable outputs.
     


§ CONCLUSION AND PERSPECTIVES

    We demonstrated that employing a correct collision frequency around the Fermi temperature with the Angstrom-scale cell size significantly improves the predictive capability of the current PIC simulations at the degenerated plasma regime. It is revealed that our new model reasonably describes the density dynamics of the laser-excited nanometric ML samples at a few eV temperatures up to ∼ 2 ps. However, the implemented collision model still suffers from several fundamental issues that may affect its predictive power. For instance, the binary collision model implemented in the PICLS code is suitable only for small-angle scattering, a reasonable assumption in high-temperature plasmas <cit.> but not at low temperatures and high densities. This limitation becomes particularly significant when dealing with low-temperature dense plasmas, where the assumptions of the binary collision model may not hold. Here, one needs careful attention when interpreting results and making conclusions about the electron-ion equilibrium rate and other essential plasma properties.

    Moreover, the relatively substantial radiation loss, approximately ∼30 % of the total energy, observed in the PIC simulation can result in a pronounced effect on the electron temperature and the subsequent density evolution. This is primarily attributed to Bremsstrahlung radiation — a consequence of electron collisional processes. To address this issue, incorporating a radiation transport module into the simulation is necessary. Such an enhancement would enable a more accurate characterization of the interplay between radiation loss and its impact on the dynamics of the laser-excited ML target.

    Our experimental results also indicate that the intermixing of the nanolayer within the ML target plays a crucial role in the density dynamics. This intermixing effect, however, is not fully captured in the current 1D PIC simulation employed in this paper. The 1D nature of the simulation may not be well-suited to capture the three-dimensional complexities of the system and its consequential impact on the density dynamics of laser-excited multilayer targets.
    
    Quantum mechanical effects and electron degeneracy, crucial factors in accurately describing the behavior of particles in low-temperature and high-density plasmas, are not inherently taken into account in PIC simulations. One may also encounter methodological challenges in exploring the parameter space studied here, which intersects with plasma parameters typically investigated using Time-Dependent Density Functional Theory (TDDFT) <cit.> and Molecular Dynamics (MD) <cit.> simulations. These challenges are especially pronounced under conditions featuring exchange-correlation effects in cold and dense environments. The adoption of PIC simulations in this realm is innovative, and the contributions outlined here allow for the study of this parameter space using a simplified surrogate model approach. This approach builds upon state-of-the-art collisional operators and ionization models. However, it introduces complexity when diverging from an LTE framework towards a non-LTE approach. Specifically, it implicitly assumes LTE by constructing dynamical means from the local phase space distribution, a process that overlooks the long- and short-range order within the material. At the same time, it provides access to study large-scale systems inaccessible by TDDFT/MD methods. While enabling simulations at this scale, it may neglect physical interactions that could impact the accuracy of our findings, which can best be evaluated in comparison to experiments.

    Including the above-mentioned physics in our PIC simulation would require prohibitive computational cost, especially for a multidimensional PIC simulation with angstrom scale cell size. The numerical effort for solving Maxwell's equations scales as (∝ 1/Δ x^4), making such an approach impractical <cit.>. This computational limitation can be effectively addressed by leveraging modern, stable, and accelerated PIC simulation frameworks like PIConGPU <cit.>. Additionally, innovative hybrid approaches, such as MicPIC <cit.>, combining MD and PIC techniques, offer a promising strategy to overcome this challenge. In MicPIC, MD handles short-range electrostatic interactions, while PIC is employed for long-range electromagnetic effects, providing a balanced and efficient solution.

    Experimental measurements using GISAXS could serve as reliable benchmarks for assessing the impact of implementing new physics in PIC simulations. Here, some advancements in measurement techniques can be explored to enhance the utility of these experimental benchmarks and improve their synergy with simulations. Specifically, adopting nanofocusing of XFEL beams in experiments holds the potential to enhance temporal resolution significantly. This improvement enables more accurate and precise capture of ultra-fast density dynamics in laser-excited targets. When incorporated into simulations, these improved experimental results can offer a more detailed understanding of crucial processes, such as laser absorption and subsequent heat propagation into the target. The intricacies of these processes, inherently dependent on collisions, become more discernible with better experimental precision.
    
    Moreover, the improved experimental time-resolution precision provides a unique opportunity to investigate the mechanisms that initiate the intermixing process between adjacent nanolayers of dense plasmas in the early stages after the laser excitation. It also enables a closer examination of how collisional effects influence the speed and dynamics of the intermixing process. Consequently, the refined experimental benchmarks not only validate simulation outcomes but also facilitate a more nuanced exploration of the underlying physics, contributing to a deeper understanding of laser-solid interactions and subsequent plasma dynamics.



§ ACKNOWLEDGMENTS

    M.B., C.G., and M.Na. acknowledge funding by DFG GU 535/6-1. We are particularly grateful for the invaluable contributions of our colleagues. Special thanks to Y. Sentoku (Osaka University) for offering essential insights into PIC simulations. We extend our gratitude to R. E. Brambrink (EuXFEL) for his expert advice on hydrodynamic simulations, and we appreciate the valuable discussions with J.-P. Schwinkendorf (HZDR), L. Wollenweber (EuXFEL), M. Habibi (EuXFEL), C. F. Grote (Max Planck Institut für Evolutionsbiologie), and C. Rödel (U. Schmalkalden).
    The success of the XFEL experiments at BL2 of SACLA, approved by the Japan Synchrotron Radiation Research Institute (JASRI) under Proposal No. 2018B8049, was a collaborative effort directed by C.G. and M.Na. The team, including L.R., M.Ba., T.R.P., M.M., N.P.D., S.G., T.M., M.Ni., C.G., and M.Na., was ably supported by T. Yabuuchi, Y. Inubushi, K. Sueda, and T. Togashi. We also acknowledge N.P.D. and J.K.K. for their assistance in generating EOS tables. The manuscript benefited from the comments and suggestions of all involved.  


