# (Sub-)picosecond surface correlations of femtosecond laser excited Al-coated
multilayers observed by grazing-incidence x-ray scattering

L. Randolph<EMAIL_ADDRESS>Department Physik, University Siegen, 57072
Siegen, Germany  European XFEL, 22869 Schenefeld, Germany  M. Banjafar
European XFEL, 22869 Schenefeld, Germany  Technical University Dresden, 01069
Dresden, Germany  T. Yabuuchi  Japan Synchrotron Radiation Research Institute
(JASRI)  RIKEN SPring-8 Center, Sayo, Hyogo 679-5148, Japan  C. Baehtz
Helmholtz-Zentrum Dresden-Rossendorf, 01328 Dresden, Germany  M. Bussmann
Helmholtz-Zentrum Dresden-Rossendorf, 01328 Dresden, Germany  Center for
Advanced Systems Understanding (CASUS), 02826 Görlitz, Germany  N. P. Dover
The John Adams Institute for Accelerator Science, Imperial College London,
London SW7 2BW, United Kingdom  L. Huang  Helmholtz-Zentrum Dresden-
Rossendorf, 01328 Dresden, Germany  Y. Inubushi  Japan Synchrotron Radiation
Research Institute (JASRI)  RIKEN SPring-8 Center, Sayo, Hyogo 679-5148, Japan
G. Jakob  Institute of Physics, Johannes Gutenberg-University, 55099 Mainz,
Germany  M. Kläui  Institute of Physics, Johannes Gutenberg-University, 55099
Mainz, Germany  D. Ksenzov  Department Physik, University Siegen, 57072
Siegen, Germany  M. Makita  European XFEL, 22869 Schenefeld, Germany  K.
Miyanishi  RIKEN SPring-8 Center, Sayo, Hyogo 679-5148, Japan  M. Nishiushi
Kansai Photon Science Institute, National Institutes for Quantum Science and
Technology, Kyoto 619-0215, Japan  Ö. Öztürk  Department Physik, University
Siegen, 57072 Siegen, Germany  M. Paulus  Fakultät Physik / DELTA, TU
Dortmund, 44221 Dortmund, Germany  A. Pelka  Helmholtz-Zentrum Dresden-
Rossendorf, 01328 Dresden, Germany  T. R. Preston  European XFEL, 22869
Schenefeld, Germany  J.-P. Schwinkendorf  European XFEL, 22869 Schenefeld,
Germany  Helmholtz-Zentrum Dresden-Rossendorf, 01328 Dresden, Germany  K.
Sueda  RIKEN SPring-8 Center, Sayo, Hyogo 679-5148, Japan  T. Togashi  Japan
Synchrotron Radiation Research Institute (JASRI)  RIKEN SPring-8 Center, Sayo,
Hyogo 679-5148, Japan  T. E. Cowan  Helmholtz-Zentrum Dresden-Rossendorf,
01328 Dresden, Germany  T. Kluge  Helmholtz-Zentrum Dresden-Rossendorf, 01328
Dresden, Germany  C. Gutt  Department Physik, University Siegen, 57072 Siegen,
Germany  M. Nakatsutsumi<EMAIL_ADDRESS>European XFEL, 22869
Schenefeld, Germany

###### Abstract

Femtosecond high-intensity laser pulses at intensities surpassing
${10}^{14}\text{\,}\mathrm{W}\mathrm{/}\mathrm{c}\mathrm{m}^{2}$ can generate
a diverse range of functional surface nanostructures. Achieving precise
control over the production of these functional structures necessitates a
thorough understanding of the surface morphology dynamics with nanometer-scale
spatial resolution and picosecond-scale temporal resolution. In this study, we
show that individual XFEL pulses can elucidate structural changes on surfaces
induced by laser-generated plasmas, employing grazing-incidence small-angle
x-ray scattering (GISAXS). Using aluminum-coated multilayer samples we can
differentiate between ultrafast surface morphology dynamics and subsequent
subsurface density dynamics, achieving nanometer-depth sensitivity and
subpicosecond temporal resolution. The observed subsurface density dynamics
serve to validate advanced simulation models depicting matter under extreme
conditions. Our findings promise to unveil novel avenues for laser material
nanoprocessing and high-energy-density science.

Keywords: grazing-incidence x-ray scattering; ultrafast surface dynamics;
laser processing; XFEL

##  I Introduction

Intense, ultrashort laser-solid interactions at intensities ranging from
${10}^{13}\text{\,}$ to
${10}^{16}\text{\,}\mathrm{W}\mathrm{/}\mathrm{c}\mathrm{m}^{2}$ are of
paramount importance for laser nanoprocessing to achieve functional surfaces.
Among these interactions, one of the most notable examples is the generation
of Laser-Induced Periodic Surface Structures (LIPSS) Bonse and Gräf (2020),
which find applications, e.g., in antibacterial coatings, optical devices,
chemical sensing, and tribology Lutey  _et al._ (2018); *San-Blas_LIPSS_20;
*Li_LIPSS_20; *Bonse_Tribology18. Understanding the dynamics of both surface
and sub-surface phenomena at nanometer or even atomic scales, within the
requisite temporal frame, is pivotal for comprehending the underlying physics
responsible for creating desired surface structures in a controlled manner.
Despite notable advancements in theoretical frameworks and models in recent
years, aimed at elucidating intricate mechanisms involved in the self-
organization of nanostructures under ultrashort laser irradiation within
relevant temporal scales Rudenko  _et al._ (2020), direct experimental
visualization of this process remains scarce. Thus far, the majority of
experimental findings regarding surface manufacturing have relied on post-
mortem analyses, which lack temporal dynamics. Time-resolved experiments
employing optical lasers, such as optical reflectometry, interferometry, and
spectroscopy, suffer from limited spatial resolution and a lack of bulk
sensitivity. Laser processing involves a complex chain of various physical
processes occurring across different temporal and spatial scales. Initially,
ultrashort laser interactions with metals excite electrons within the surface
skin layer, typically spanning tens of nanometers. These excited electrons
then propagate into the bulk at Fermi velocity before rapid thermalization
occurs within sub-ps timescales via collisions. Depending on the excitation
strength and material properties, this process triggers electron-lattice/ion
thermalization, coherent phonon oscillation, and lattice heating in picosecond
timescales. This can lead to subsequent thermal or non-thermal melting,
followed by ablation through spallation in sub-ns timescales. These phenomena
collectively contribute to later crystallization or amorphization, leading to
specific surface nanostructures. Therefore, understanding above-mentioned
early time phenomena in (sub-)picosecond dynamics is crucial. To address this
need, here, we demonstrate experimental capabilities of visualizing physical
processes at both surface and subsurface levels with picosecond and nanometer
resolutions using an X-ray Free Electron Laser (XFEL) operating in a grazing-
incidence geometry. Grazing-incidence X-ray small-angle scattering (GISAXS) is
a well-established technique for probing lateral structures and correlations
along the surface normal. It has been extensively utilized over the past
decades at synchrotron x-ray facilities to characterize material structures
across length scales ranging from sub-nanometers to microns Roth (2016); Su
_et al._ (2017), albeit limited to the millisecond timescale due to photon
accumulation requirements.
Recently, we demonstrated the applicability of GISAXS to XFEL to track
nanometric multilayer dynamics at unprecedented picosecond scales, which is
six orders of magnitude faster than previously achievable Randolph  _et al._
(2022). This breakthrough is enabled by the XFEL’s ability to deliver an
immense number of photons, equivalent to those typically accumulated over a
second at third-generation synchrotron x-rays, within a single pulse lasting a
few tens of femtoseconds. Despite inherent time smearing due to grazing-
incidence, providing different parts of the x-ray arriving at the surface at
different times, achieving approximately picosecond time resolution is
feasible.
Utilizing the in-plane scattering signal, characterized by the wavevector
transfer along the sample depth ($Q_{z}$), which correlates closely with
specular reflectivity Holý  _et al._ (1993), our experiments revealed the
compressed, heated, intermixed, and ablated multilayer dynamics after the
laser irradiation. On the other hand, out-of-plane diffuse scattering along
the $Q_{y}$ direction provides insights into ultrafast changes in the lateral
distribution of surface roughnesses. This can be described with the help of
the height-height correlation function

$\displaystyle C(R)=<h(0)h(R)>=\sigma^{2}\exp{-(R/\xi_{\parallel})^{2H}}$ (1)

where $R$ is the spatial separation of two points, $\sigma$ is the rms
roughness, $\xi_{\parallel}$ is the lateral correlation length and $H$ the
hurst parameter Sinha  _et al._ (1988). The correlation function along the
sample depth can be written as

$\displaystyle<h_{j}(0)h_{k}(R)>=$
$\displaystyle\frac{1}{2}\left[\frac{\sigma_{k}}{\sigma_{j}}C_{j}(R)+\frac{\sigma_{j}}{\sigma_{k}}C_{k}(R)\right]$
$\displaystyle\times\exp{\left(-\frac{z_{j}-z_{k}}{\xi_{\perp}}\right)}.$ (2)

Here, $C(R)$ is the height-height correlation function, $z$ denotes the depth
inside the sample and $\xi_{\perp}$ the cross-correlation length. The indices
$j$ and $k$ denote different interfaces of the ML sample Schlomka  _et al._
(1995).

Additionally, by incorporating a low-Z capping layer on top of a high-Z
multilayer (ML), GISAXS facilitates measurements of energy transport several
hundred nanometers within the material. Upon ultrashort laser irradiation,
strongly confined temperature gradients induce compression (shock) waves,
where the progression of temperature gradients depends on the mean-free-path
of excited electrons which is a function of both temperature and density. Our
result show that increasing laser excitation strength transitions the dominant
energy transfer mode from compression wave dominance to the electron-thermal
wave dominant regime.

##  II Experimental methods and results

Figure 1:  Schematics of the experimental setup to investigate the surface and
subsurface solid-density plasma dynamics with GISAXS using single femtosecond
x-ray FEL pulses. (a) A ML sample consists of 5 repetitions of tantalum (Ta)
and copper nitride (Cu3N) of $4.5$ and $8.5\text{\,}\mathrm{nm}$ thickness
covered by a $200\text{\,}\mathrm{nm}$ aluminium capping layer. The samples
were irradiated by an optical laser with a central wavelength of
$800\text{\,}\mathrm{nm}$, intensity of $(8\pm 0.2)\times${10}^{14}\text{\,}$$
and $(8\pm
0.2)\times${10}^{15}\text{\,}\mathrm{W}\mathrm{/}\mathrm{c}\mathrm{m}^{2}$$
with $40\text{\,}\mathrm{fs}$ pulse duration. The laser is irradiated by an
incident angle of $17\text{\,}\mathrm{\SIUnitSymbolDegree}$ from the surface
normal in $p$ polarization. The x-ray pulses with $8.81\text{\,}\mathrm{keV}$
photon energy, $7\text{\,}\mathrm{fs}$ FWHM duration are irradiated on the
sample at the grazing-incidence angle of
$\alpha_{i}=$0.75\text{\,}\mathrm{\SIUnitSymbolDegree}$$, i.e., slightly above
the critical angle of external total reflection of the layer materials. The
laser beam is defocused to obtain a
$\sim$500\text{\,}\mathrm{\SIUnitSymbolMicro m}$$ spot diameter to cover the
x-ray footprint of $\sim$300\text{\,}\mathrm{\SIUnitSymbolMicro m}$$ at the
surface. Scattered photons are recorded by 2 modules of the MPCCD area
detector placed around the specular direction. The strong specular peak at
$Q_{z}=$1.16\text{\,}\mathrm{n}\mathrm{m}^{-1}$$ is blocked. (b) In-plane
signal along $Q_{z}$ (at $Q_{y}=0$) for different time delays between $0$ and
$14\text{\,}\mathrm{p}\mathrm{s}$ after the laser intensity peak.

The experiment was performed at the BL2 EH6 station at the SACLA XFEL facility
in Japan Yabuuchi  _et al._ (2019). The basic setup is similar to Randolph
_et al._ (2022). A metallic ML sample, consisting of $5$ repetitions of
tantalum (Ta, $4.5\text{\,}\mathrm{nm}$) and copper nitride (Cu3N,
$8.5\text{\,}\mathrm{nm}$) capped with a $200\text{\,}\mathrm{nm}$ thick
Aluminum top layer was irradiated by an optical laser with a central
wavelength of $800\text{\,}\mathrm{nm}$, intensities of
$8\text{\times}{10}^{14}\text{\,}$ and
$8\text{\times}{10}^{15}\text{\,}\mathrm{W}\mathrm{/}\mathrm{c}\mathrm{m}^{2}$,
pulse duration of $40\text{\,}\mathrm{fs}$ and a focus size of
$\sim$500\text{\,}\mathrm{\SIUnitSymbolMicro m}$$ (Fig. 1a). After a variable
delay time, the surface state of the sample was probed via surface sensitive
x-ray scattering employing ultrafast XFEL pulses of x-ray photons of energy
$8.81\text{\,}\mathrm{keV}$. The x-ray beam with a pulse energy of
$\sim$0.1\text{\,}\mathrm{m}\mathrm{J}\mathrm{/}\mathrm{p}\mathrm{u}\mathrm{l}\mathrm{s}\mathrm{e}$$
was incident under a grazing incidence angle of $0.75\text{\,}\deg$ and the
diffusely scattered intensity was recorded by a multi-port CCD (MPCCD) area
detector Kameshima  _et al._ (2014) with the intense specular peak (incident
angle equals exit angle,
Q${}_{\text{specular}}=$1.16\text{\,}\mathrm{n}\mathrm{m}^{-1}$$) being
blocked by a beam stop. Fig. 1 a) displays a typical single pulse GISAXS
pattern. The shallow angle-of-incidence, in combination with the beam size of
$4\text{\,}\mathrm{\SIUnitSymbolMicro m}$ leads to a large footprint
($\sim$300\text{\,}\mathrm{\SIUnitSymbolMicro m}$$) of the beam on the sample
surface, limiting the temporal resolution of the detected x-ray signal to
$\sim$1\text{\,}\mathrm{ps}$$ integration time.

Figure 2:  In-plane scattering signal along $Q_{z}$ around the Bragg-like peak
at $2.1\text{\,}\mathrm{n}\mathrm{m}^{-1}$. (a) Circular dots represent the
experimental data for various time delays, while solid lines correspond to a
Gaussian model fit. (b) Peak intensity derived from the Gaussian model plotted
against time delay. Circular dots denote the maximum intensity values from
(a), while the dashed line represents a refinement using an error function.
Blue color denotes the lower laser intensity of
$8\text{\times}{10}^{14}\text{\,}\mathrm{W}\mathrm{/}\mathrm{c}\mathrm{m}^{2}$
while red indicates the higher laser intensity of
$8\text{\times}{10}^{15}\text{\,}\mathrm{W}\mathrm{/}\mathrm{c}\mathrm{m}^{2}$.
(c) $Q_{z}$-position of the Bragg-like peak. The dashed line is a guide to the
eyes.

Fig. 1b) illustrates the in-plane scattering signal ($Q_{y}=0$) as a function
of $Q_{z}$ for various delay times after laser irradiation. Because of the
peculiar geometry of GISAXS, the in-plane scattering signal contains momentum
transfers with components in both the normal ($Q_{z}$) and surface-parallel
($Q_{x}$) directions. As a result, the scattering signal reflects vertical
density correlations along the z-direction ($Q_{z}$) of the ML structure,
roughness correlations along the surface plane ($Q_{x}$) direction, and cross-
correlations between different interfaces. For clarity, only the $Q_{z}$
values are plotted here, as $Q_{x}\ll Q_{z}$. The in-plane signal at
$Q_{z}>Q_{\text{specular}}$ is closely associated with the specular
reflectivity curve Holý  _et al._ (1993), which we characterized ex-situ. The
spacing of the intense Bragg-like peaks at $Q_{z}=$1.6\text{\,}$$,
$2.1\text{\,}$ and $2.6\text{\,}\mathrm{n}\mathrm{m}^{-1}$ corresponds to the
typical length scale, i.e. $13\text{\,}\mathrm{nm}$ thickness of each Ta/Cu3N
double layer. The Kiessig fringes Kiessig (1931), represented as smaller peaks
between the Bragg peaks, are a fingerprint of the number of double-layer
repeats in the sample. Upon laser irradiation, substantial changes in the in-
plane scattering signal are observed during the first
$12\text{\,}\mathrm{ps}$. The number and intensity of Kiessig fringes remain
relatively stable during the initial $\sim$9\text{\,}\mathrm{ps}$$ after laser
irradiation, thereafter exhibiting a gradual decrease, concomitant with a
simultaneous broadening and reduction of the intense Bragg-like peaks. After
$12\text{\,}\mathrm{ps}$, the Kiessig fringes are nearly gone, leaving behind
only a broad residual of the Bragg-like peaks.

The most notable structural changes inside the ML sample are manifested by the
change of intensity and position of the intense peaks at $Q_{z}=1.6$, $2.1$
and $2.6\text{\,}\mathrm{n}\mathrm{m}^{-1}$. Fig. 2a) illustrates the temporal
evolution of the peak at $Q_{z}=$2.1\text{\,}\mathrm{n}\mathrm{m}^{-1}$$
following laser irradiation with an intensity of
$8\text{\times}{10}^{14}\text{\,}\mathrm{W}\mathrm{/}\mathrm{c}\mathrm{m}^{2}$.
While only the peak at $Q_{z}=$2.1\text{\,}\mathrm{n}\mathrm{m}^{-1}$$ is
displayed here, qualitatively similar behaviour is observed for the peaks at
$Q_{z}=$1.6\text{\,}$$ and $2.6\text{\,}\mathrm{n}\mathrm{m}^{-1}$. The solid
lines represent refinements using a Gaussian function. The peak amplitude as a
function of time delay, extracted from the Gaussian refinement, is depicted in
Fig. 2b), with blue dots, which shows an ultrafast increase within the first
picosecond. Subsequently, the intensity remains constant, followed by a
gradual decrease at around $\sim$10\text{\,}\mathrm{ps}$$ associated with a
loss of structural correlation between the layers. This decay is modelled
using an error function

$\text{Intensity}=A\cdot\text{erf}\left(-\frac{(x-t_{0})}{\Delta t}\right)-B$
(3)

resulting in values for $t_{0}=11.4\pm$0.3\text{\,}\mathrm{ps}$$ and a width
of $\Delta t=2.1\pm$0.4\text{\,}\mathrm{ps}$$. Here, $A$ represents a
stretching factor, $t_{0}$ denotes the time at which the error function
decreased to half of its initial value, $\Delta t$ is the duration of the
decrease, and $B$ stands for a constant offset. In comparison, the same
analysis for the higher laser intensity case of
$8\text{\times}{10}^{15}\text{\,}\mathrm{W}\mathrm{/}\mathrm{c}\mathrm{m}^{2}$,
is depicted with red circular dots. Similar to the lower intensity case, the
ultrafast increase in scattering intensity within the first
$\sim$1\text{\,}\mathrm{ps}$$ is evident. However, for the time delays
$>\sim$1\text{\,}\mathrm{ps}$$ the amplitude quickly begins to decline again.
The refined error function yields a time constant of
$t_{0}=2.8\pm$0.3\text{\,}\mathrm{ps}$$ and a width of $\Delta
t=1.0\pm$0.5\text{\,}\mathrm{ps}$$. On top of the change in the Bragg
intensity, we observe a steady increase in the Q-position with time delay for
both laser intensities, as summarized in Fig. 2 c). This shift towards larger
Q-values implies compression of the double-layer structure as the typical
length scale decreases.

Figure 3:  In-plane signal along $Q_{z}$ around the dynamical diffraction area
between $0.6\text{\,}$ and $1.0\text{\,}\mathrm{n}\mathrm{m}^{-1}$. (a) In-
plane signal for different time delays at lower laser intensity case
($8\text{\times}{10}^{14}\text{\,}\mathrm{W}\mathrm{/}\mathrm{c}\mathrm{m}^{2}$).
(b) The circular dots display the integrated signal between
$Q_{z}=0.65-$0.79\text{\,}\mathrm{n}\mathrm{m}^{-1}$$ as a function of time
delay. The dashed line represents a guide to the eyes. Blue color stands for
the lower laser intensity of
$8\text{\times}{10}^{14}\text{\,}\mathrm{W}\mathrm{/}\mathrm{c}\mathrm{m}^{2}$
and red denote the higher laser intensity
$8\text{\times}{10}^{15}\text{\,}\mathrm{W}\mathrm{/}\mathrm{c}\mathrm{m}^{2}$.

The dynamic diffraction effects at $Q_{z}<Q_{\text{specular}}$ yield insight
into the dynamics of the top Al layer. In particular, the scattering signals
at exit angles equal to the critical angle for the total external reflection –
$Q_{z}$ = $0.75\text{\,}\mathrm{n}\mathrm{m}^{-1}$ for Al – are particularly
surface-sensitive, originating from an evanescent x-ray wave travelling
parallel to the surface. This so-called Yoneda peak Yoneda (1963), primarily
originated from interference within the topmost surface layer, serves as a
sensitive marker of its surface structure. A close-up of this Yoneda peak is
summarized in Fig. 3 a). The peak at
$Q_{z}=$0.75\text{\,}\mathrm{n}\mathrm{m}^{-1}$$ exhibits a progressive
decrease after laser irradiation, which persists for at least at
$t=$12\text{\,}\mathrm{ps}$$ before almost disappearing at
$t=$14\text{\,}\mathrm{ps}$$, indicative of the ongoing presence of the solid
density Al cover layer. Assuming the surface ablation speed approximates the
speed of sound,
$C_{s}=\sqrt{{(\gamma_{e}Z_{mean}k_{B}T_{e}+\gamma_{i}k_{B}T_{i})}/{M_{i}}}$,
the ablation of the $200\text{\,}\mathrm{nm}$ layer would require
$\sim$15\text{\,}\mathrm{ps}$$. Here, $T_{e}\sim 10\mathrm{eV}$ is the
electron temperature, $T_{i}\sim 1\mathrm{eV}$ is the ion temperature,
$Z_{mean}\sim 5$ is the mean ionization, $\gamma_{e}=1$ and $\gamma_{i}=3$ are
adiabatic index of electrons and ions, respectively and $M_{i}=26$ a.u. is the
ion mass for Al.
The integrated intensity of the Yoneda peak
($Q_{z}=0.65-$0.79\text{\,}\mathrm{n}\mathrm{m}^{-1}$$) is depicted in Fig. 3
b) for both low (blue) and high (red) intensity cases. In both cases, a rapid
initial decay in intensity is observed within the first picosecond, followed
by a period of quasi-constant signal intensity extended up to
$\sim$12\text{\,}\mathrm{ps}$$ and $\sim$5\text{\,}\mathrm{ps}$$ before
subsequent decrease, respectively. This temporal behaviour aligns
qualitatively with that observed at the Bragg-like peak at higher $Q_{z}$
shown in Fig. 2 b). In Fig. 3 a), two additional Yoneda peaks are also visible
at higher $Q_{z}$, corresponding to the Yoneda peaks for $\mathrm{Cu_{3}N}$
($Q_{z}=$0.86\text{\,}\mathrm{n}\mathrm{m}^{-1}$$) and Ta
($0.95\text{\,}\mathrm{n}\mathrm{m}^{-1}$), respectively. Given that these
materials are embedded within a $200\text{\,}\mathrm{nm}$ thick Al layer, the
time-dependency of their intensity is not as pronounced as that of the Al
Yoneda peak, except during the period between 12 and $14\text{\,}\mathrm{ps}$
when all peaks suddenly disappear. As we will discuss below, aided by plasma
simulations, this phenomenon can be attributed to the arrival of a compression
wave to the ML, which initiates significant modulation of the ML structure. An
additional noteworthy feature in the figure is the presence of high-frequency
fringes superimposed on the Ta and $\mathrm{Cu_{3}N}$ Yoneda regions. These
fringes arise from the interference between waves scattered from the Al
surface and those from the Al-ML interface. Following laser excitation and the
subsequent expansion of the top Al surface, these fringes rapidly vanish as a
result of blurring of the distinct interface between vacuum and Al.

Figure 4: GISAXS signals and corresponding retrieved real-space electron
density profile. (a) Lineout at $0\text{\,}\mathrm{ps}$ (white circular dots)
and $1\text{\,}\mathrm{ps}$ (grey circular dots) after the laser intensity
peak. Solid lines represent models using the program BornAgain. (b) Zoom into
the area of the Bragg-like peak at
$Q_{z}=$2.1\text{\,}\mathrm{n}\mathrm{m}^{-1}$$ and (c) zoom into the area of
dynamical diffraction around $Q_{z}=$0.75\text{\,}\mathrm{n}\mathrm{m}^{-1}$$.
(d) Retrieved real-space effective electron density profile as a function of
the depth (Z). The laser irradiates the sample from the left. (e) Generated
surface model for $H=0.6$, $\sigma=$2.3\text{\,}\mathrm{nm}$$ and (f)
$H=0.15$, $\sigma=$2.6\text{\,}\mathrm{nm}$$. (g) Lineout along $Q_{y}$ at the
Al Yoneda peak ($Q_{z}=$0.75\text{\,}\mathrm{n}\mathrm{m}^{-1}$$). Dashed
lines indicate a fit of the decay. (h) Lineout along $Q_{y}$ at
$Q_{z}=$2.1\text{\,}\mathrm{n}\mathrm{m}^{-1}$$.

Both Bragg-like peaks and the Yoneda peak indicate the presence of three
distinct time regimes. Immediately after laser irradiation
$\sim$1\text{\,}\mathrm{ps}$$, we observe an ultrafast increase in the
intensity of the Bragg-like peaks without any change in their
$Q_{z}$-positions, accompanied by the decay of the Yoneda peak. As we will
discuss further below, this phenomenon is attributed to a modification in the
surface structural properties. In the second time frame, spanning from $1$ to
$\sim$10\text{\,}\mathrm{ps}$$ or to $\sim$3\text{\,}\mathrm{ps}$$ for low and
high laser intensities, respectively, the intensities of the Bragg-like peaks
and Yoneda peak remain constant, while the $Q_{z}$ position of the Bragg-like
peaks steadily increases. This suggests that the embedded ML undergoes
compression, although ML structure remains constant. In the third time frame,
the intensities of the Bragg and Yoneda peaks are further reduced, accompanied
by a continued shift of the $Q_{z}$-position to even higher values. This
indicates further compression of the entire multilayer, as well as significant
modulation of the ML structure, resulting in the loss of correlation and
subsequent decrease in the x-ray scattering intensity.

##  III Discussion

To facilitate a more quantitative discussion of the ultrafast surface
dynamics, we employed the BornAgain Pospelov  _et al._ (2020) GISAXS analysis
program to model our experimental observations. The white circular dots in
Fig. 4a) represent the experimental data at $0\text{\,}\mathrm{ps}$ while the
grey dots depict the signal at a delay of $1\text{\,}\mathrm{ps}$. The solid
lines represent the BornAgain model corresponding to the effective electron
density profiles shown in Fig. 4d). Figures 4b) and c) provide enlarged views
of the intense peak at $2.1\text{\,}\mathrm{n}\mathrm{m}^{-1}$ and the
surface-sensitive Yoneda region, respectively. According to the model, the
observed variation between the blue ($0\text{\,}\mathrm{ps}$) and red
($1\text{\,}\mathrm{ps}$) lines in Figs 4a)-c) can be attributed to the
decrease in the Hurst parameter of the Al surface layer from $0.6$ to $0.15$,
indicating an increase in the spatial frequency of roughness, as illustrated
schematically in Figs. 4 e) and f). This observation is supported by analyzing
the out-of-plane signal along $Q_{y}$ at different $Q_{z}$ positions. Fig. 4g)
illustrates a lineout along $Q_{y}$ at the Al Yoneda peak. The decay starting
at $Q_{y}=$0.01\text{\,}\mathrm{n}\mathrm{m}^{-1}$$ is modeled via

$\text{Intensity}\sim Q_{y}^{-m}$ (4)

where $m$ defines the slope of the decay which is proportional to the Hurst
parameter $H$ Salditt  _et al._ (1994). Refining the lineouts for $0$ and $1$
ps indeed reveals a reduction of the Hurst parameter as $m$ has decreased.
Conversely, the lineout along $Q_{y}$ at larger
$Q_{z}=$2.1\text{\,}\mathrm{n}\mathrm{m}^{-1}$$ (Fig. 4h)) shows no change in
the Hurst parameter, indicating that the ultrafast change is localized to the
surface rather than at the ML interfaces. Furthermore, it appears that in
order to align the model with the experimental observation from $0$ to
$1\text{\,}\mathrm{ps}$, the vertical correlation length needs to be increased
from $100$ to $150\text{\,}\mathrm{nm}$, and the surface roughness RMS
$\sigma$ should be slightly increased from $2.3$ to $2.6\text{\,}\mathrm{nm}$.
We hypothesize that this ultrafast change in surface roughness properties to
be attributed to the presence of a thin aluminium oxide layer on the surface,
which typically forms within minutes when exposed to air Nguyen  _et al._
(2018). This oxide layer likely possesses a distinct surface morphology
compared to the intrinsic surface morphology parameters of bare aluminium,
leading to a reduction in spatial frequency of roughness as well as a decrease
in vertical correlation length. Given that this layer is typically extremely
thin, on the order of a nanometer, it evaporates instantaneously
($<$1\text{\,}\mathrm{ps}$$) after laser excitation revealing the intrinsic
surface properties of the initial Al layer.

Figure 5: MULTI-fs hydrodynamic simulations irradiated by a p-polarized laser
pulse with a duration of 50 fs, a wavelength of 800 nm, and a laser intensity
of
$8\text{\times}{10}^{14}\text{\,}\mathrm{W}\mathrm{/}\mathrm{c}\mathrm{m}^{2}$
at time delays of (a) 0.5 and (b) 7.75 ps, and with a laser intensity of
$8\text{\times}{10}^{15}\text{\,}\mathrm{W}\mathrm{/}\mathrm{c}\mathrm{m}^{2}$
at time delays of (c) 0.5 and (d) 2.7 ps. The ML structure matched that of the
experimental setup. (e and f) 1D Particle-in-cell PICLS simulation results at
time delays of (e) 1.5 and (f) 2.7 ps. In all figures, the yellow, grey, and
red regions represent the laser-induced pressure wave (in Mbar or Mbar$/10$),
mass density (in $\mathrm{g/cm^{3}}$), and electron temperature (in eV),
respectively. The dashed line denotes the initial target density. The red
arrow in (a) indicates the direction of the incoming laser.

To substantiate our experimental findings further, we also conducted two
advanced simulations: a 1D MULTI-fs hydrodynamic simulation and 1D PICLS
Particle-in-Cell (PIC) simulation. The 1D MULTI-fs hydrodynamic simulation
code is tailored for modeling short pulse high-intensity laser-solid
interactions which directly solves the Maxwell equations. The code
incorporates a temperature-dependent collision frequency, thermal conductivity
ranging from metallic solids to high-temperature ideal plasmas, and a distinct
equation-of-state (EOS) for electrons and ions (two-temperature model). The
details of these simulations’ setup is summarized in Methods section. The 1D
PICLS simulation incorporates the collisional 1d3v (one-dimensional in space
and three-dimensional in velocity) Sentoku and Kemp (2008). An interpolated
collision frequency akin to the MULTI-fs method was integrated into the code
to address Angstrom-scale collisions occurring at electron temperatures
($T_{e}$) around the Fermi temperature ($T_{F}$). The setup details for these
simulations are summarized in the Methods section.

Figs. 5 a) and c) display the density profile simulated using MULTI-fs at a
delay of 0.5 ps. The localized heating induced by the laser within the skin
depth instantaneously elevates the surface electron temperature to a few tens
of eV and generates pressures exceeding $>$10\text{\,}\mathrm{Mbar}$$,
prompting immediate surface expansion. This observation is aligned with the
experiment demonstrating ultrafast changes in surface properties. The
resultant expansion causes an evaporation of the top few nanometers. This
value contributes to the roughness RMS parameter for the GISAXS pattern.

During the second time frame, spanning from $1$ to
$\sim$9\text{\,}\mathrm{ps}$$ in the case of lower laser intensity
($I_{L}=$8\text{\times}{10}^{14}\text{\,}\mathrm{W}\mathrm{/}\mathrm{c}\mathrm{m}^{2}$$),
the intensity of the Bragg-like peak remains constant (Fig.2 b)). As confirmed
by MULTI-fs simulations, it appears that between $1$ and
$8\text{\,}\mathrm{ps}$ the compression wave is only in the Al layer until it
arrives at the Al-ML boundary (Fig. 5 b), highlighted in yellow). Primarily,
only the Al layer undergoes modulation during this interval, while the ML
remains static with slight compression. Subsequently, after
$\sim$8\text{\,}\mathrm{ps}$$, the simulation reveals a strong pressure wave
compressing the ML. This aligns with the experimental findings (Fig. 2 c))
that demonstrated a further increase in the $Q_{z}$-position of the Bragg
peaks.

Figures. 5 c)-d) are the same simulation but performed at the higher laser
intensity
$I_{L}=$8\text{\times}{10}^{15}\text{\,}\mathrm{W}\mathrm{/}\mathrm{c}\mathrm{m}^{2}$$.
Here, the pressure peak arrives at the multilayer surface already after $\sim$
3 ps, which encounters strong modulation subsequently. On the contrary, the
experiment (Fig. 2 b)) shows the significant decrease in Bragg intensity
already started at $1-$2\text{\,}\mathrm{ps}$$. It implies that the
significant modulation of the ML occurs before the arrival of the strong
compression wave to the ML surface, which is not seen in the MULTI-fs
simulation.

Lagrangian hydrodynamic simulations like MULTI-fs are incapable of accurately
simulating the atomic mixing of adjacent layers Randolph  _et al._ (2022).
Therefore, we conducted kinetic 1D-PICLS simulation for the intensity case of
$I_{L}=$8\text{\times}{10}^{15}\text{\,}\mathrm{W}\mathrm{/}\mathrm{c}\mathrm{m}^{2}$$,
summarized in Fig. 5 e-f). The simulation shows that the embedded ML undergoes
significant layer intermixing from a thermal wave, even before the arrival of
the strong compression wave. This phenomenon begins already at 1.5 ps delay as
shown in Fig. 5 e). Even though the pressure peak is still in the middle of
the Al layer, the ML has already started to undergo modulation. This leads to
reductions in the Bragg peak in the x-ray scattering. By 2.7 ps delay, when
both MULTI-fs and PICLS indicate the arrival of the pressure peak (Fig. 5 f)),
PICLS demonstrates that the entire ML is intermixed, and the original ML
structure is entirely disrupted. Such a structure cannot produce any
pronounced Bragg-like peaks due to the absence of a statistically relevant
periodic structure.

In summary, for the low-intensity case
$I_{L}=$8\text{\times}{10}^{14}\text{\,}\mathrm{W}\mathrm{/}\mathrm{c}\mathrm{m}^{2}$$,
modifications to the ML structure are primarily driven by the compression
wave, as the electron temperature is low enough that the thermal wave cannot
significantly alter the ML before the arrival of the pressure wave.
Conversely, for the high-intensity case
$I_{L}=$8\text{\times}{10}^{15}\text{\,}\mathrm{W}\mathrm{/}\mathrm{c}\mathrm{m}^{2}$$,
significant heating occurs deep within the ML much before the arrival of the
pressure wave.

##  IV Conclusions

In conclusion, our investigation of surface and subsurface dynamics in high-
intensity laser-excited ML systems using ultrafast GISAXS has provided
important insights into the complex physical processes occurring in laser-
matter interactions on a picosecond time scale.
We can distinguish between ultrafast surface dynamics and subsequent
subsurface dynamics. The ultrafast modification in surface properties occurs
in sub-picosecond time scale induced by strong localized surface heating. The
slower subsurface dynamics manifest as changes in intensity and position of
the intense Bragg-like peak at
$Q_{z}=$2.1\text{\,}\mathrm{n}\mathrm{m}^{-1}$$. Comparison of the temporal
evolution of the density dynamics with 1D hydrodynamic and 1D particle-in-cell
simulations reveals qualitative agreement across different time regimes
observed in the experiment. Notably, the observed time scales align well with
experimental observations, indicating the arrival of the pressure wave at the
ML surface and subsequent intermixing of individual layers.
The presented GISAXS experiment was limited in time resolution to
$\sim$1.2\text{\,}\mathrm{ps}$$ due to the large footprint in grazing-
incidence. However, this can be improved by using a smaller x-ray focal spot,
e.g. $100\text{\,}\mathrm{nm}$, which is readily available at XFEL facilities.
This would imply only $\sim$30\text{\,}\mathrm{fs}$$ smearing at
$0.75\text{\,}\mathrm{\SIUnitSymbolDegree}$ grazing-incidence, allowing
experimental studies on surface and subsurface dynamics with nanometer spatial
and femtosecond temporal resolution to observe the ultrafast changes on the
surface or using laser intensities above
${10}^{16}\text{\,}\mathrm{W}\mathrm{/}\mathrm{c}\mathrm{m}^{2}$ to study e.g.
the effect of density oscillations Paschke-Bruehl  _et al._ (2023).

##  V Supplementary Materials

###  V.1 Multi-fs

The 1D MULTI-fs hydrodynamic simulation code is tailored for modeling short
($\leq$ ps) pulse high-intensity laser-solid interactions at intensities below
$<10^{17}\,\mathrm{W/cm^{2}}$ Eidmann  _et al._ (2000); Ramis  _et al._
(2012). MULTI-fs directly solves the Maxwell equations, capturing ultrashort
laser plasma dynamics with sharp plasma density gradients. The code
incorporates a temperature-dependent collision frequency, thermal conductivity
ranging from metallic solids to high-temperature ideal plasmas, and a distinct
equation-of-state (EOS) for electrons and ions (two-temperature model). In the
MULTI-fs simulation setup, the simulation box comprises 1478 cells distributed
as follows: 400 cells for the 200 nm-thick Al top layer, 678 cells for the ML,
and 400 cells for a 200 nm-thick Al substrate. To mitigate numerical
artifacts, finer cell resolutions were employed near layer interfaces. We
employed the EOS and ionization tables for aluminium (Al) as provided by the
MULTI-fs package. For tantalum (Ta) and copper (Cu), the EOS and ionization
data were generated using the FEOS code Faik  _et al._ (2018). The radiation
transfer module was disabled in the simulation, and the free streaming
limiting factor was set to $f=0.6$. To maintain consistency with the
experimental setup, the laser angle of incidence was fixed at $17\degree$ from
the target normal. In our simulations, Al was chosen as the substrate material
instead of silicon (Si) used in the experiment. This choice was motivated by
aluminium’s well-established EOS. Given that the energy transfer occurs from
the ML to the substrate, we anticipate that the choice between these two
materials will have minimal impact on the dynamics of the ML.

###  V.2 1D PICLS

The 1D PICLS simulation incorporates the collisional 1d3v (one-dimensional in
space and three-dimensional in velocity) Sentoku and Kemp (2008). An
interpolated collision frequency akin to the MULTI-fs method was integrated
into the code to address Angstrom-scale collisions occurring at electron
temperatures ($T_{e}$) around the Fermi temperature ($T_{F}$). To accurately
simulate the microscopic particle collisions at the atomic scale, the spatial
resolution was configured with a cell size of $\Delta x=2$ Å corresponding to
a time–step of $\Delta t=6.67\times 10^{-19}$ s. Each computational cell
accommodated 30 virtual ion particles with initial charge states of 2, 1, and
3 for Ta, Cu, and Al, respectively. Additionally, a fourth-order particle
shape was employed, along with distinct particle weightings for different ion
species. The ion number densities are set to realistic values for Ta
($n_{\mathrm{Ta}}=31.8n_{c}$), Cu ($n_{\mathrm{Cu}}=48.7n_{c}$), and Al
($n_{\mathrm{Al}}=34.6n_{c}$), where $n_{c}=m_{e}\omega_{L}^{2}/(4\pi
e^{2})=1.742\times 10^{21}\mathrm{cm^{-3}}$ represents the critical plasma
density at the laser wavelength of $\lambda_{L}=800$ nm. Here, $m_{e}$ and $e$
denote the electron mass and charge, respectively, and $\omega_{L}$ is the
laser angular frequency. Ionization dynamics are modeled employing field and
direct-impact ionization models. The laser incident angle is set as normal to
the surface, as oblique incidence is not supported in 1D PICLS. Our test
simulations using the MULTI-fs have revealed a negligible difference,
approximately 1%, in laser absorption between normal incidence and a 17-degree
angle of incidence.

## Acknowledgements

The XFEL experiments were performed at the BL2 of SACLA with the approval of
the Japan Synchrotron Radiation Research Institute (JASRI) (Proposal No.
2019B8076). We thank E. Brambrink, S. Göde, S.V. Rahul, C. Rödel, A. Kon, J.
Koga and Y. Sentoku for various discussions and advices.

### Funding

C.G. and M.Na. acknowledge funding by DFG GU 535/6-1. G.J. and M.K.
acknowledge the support by the Deutsche Forschungsgemeinschaft (DFG, German
Research Foundation) Project No. 268565370 (SFB TRR173 Projects A01 and B02)
by TopDyn and the ForLab MagSens.

### Author contributions

C.G. and M.Na. conceived the study. L.R., M.Ba., T.R.P., M.M., N.P.D., S.G.,
T.M., M.Ni., A.P., J.-P.S., C.B., C.G., and M.Na. performed the experiment at
SACLA with support from T.Y., Y.I., K.S., and T.T. L.R. analyzed the GISAXS
data under the supervision of C.G. M.Ba. performed MULTI-fs and PICLS
simulations to interpret the data with extensive support of L.H. and T.K,
under the supervision of A.P.M., M. Bu., T.E.C., and M.Na. N.P.D. and J.K.K.
gave support for generating EOS tables. Multilayer samples were prepared by
G.J. and M.K. which were further characterized by L.R. and M.P. before the
experiment. J.K., M.M., and M.Na. organized the sample cutting. L.R, M.Ba,
C.G. and M.Na. wrote the initial draft which was extensively revised by
contribution from D.K. and Ö.Ö. All authors read commented on the manuscript.

### Data Availability Statement

The data are available upon reasonable request.

### Conflicts of Interest

The authors declare no conflicts of interest.

## References

  * Bonse and Gräf (2020) J. Bonse and S. Gräf, “Maxwell Meets Marangoni – A Review of Theories on Laser‐Induced Periodic Surface Structures,” Laser and Photonics Reviews 14, 2000215 (2020).
  * Lutey  _et al._ (2018) Adrian H.A. Lutey, Laura Gemini, Luca Romoli, Gianmarco Lazzini, Francesco Fuso, Marc Faucon, and Rainer Kling, “Towards laser-textured antibacterial surfaces,” Scientific Reports 8, 1–10 (2018).
  * San-Blas  _et al._ (2020) A. San-Blas, M. Martinez-Calderon, J. Buencuerpo, L. M. Sanchez-Brea, J. del Hoyo, M. Gómez-Aranzadi, A. Rodríguez, and S. M. Olaizola, “Femtosecond laser fabrication of LIPSS-based waveplates on metallic surfaces,” Applied Surface Science 520, 146328 (2020).
  * Li  _et al._ (2020) Chao Yu Li, Sai Duan, Jun Yi, Chen Wang, Petar M. Radjenovic, Zhong Qun Tian, and Jian Feng Li, “Real-time detection of single-molecule reaction by plasmon-enhanced spectroscopy,” Science Advances 6, 1–9 (2020).
  * Bonse  _et al._ (2018) Jörn Bonse, Sabrina V. Kirner, Michael Griepentrog, Dirk Spaltmann, and Jörg Krüger, “Femtosecond laser texturing of surfaces for tribological applications,” Materials 11, 1–19 (2018).
  * Rudenko  _et al._ (2020) A. Rudenko, A. Abou-Saleh, F. Pigeon, C. Mauclair, F. Garrelie, R. Stoian, and J.P. Colombier, “High-frequency periodic patterns driven by non-radiative fields coupled with marangoni convection instabilities on laser-excited metal surfaces,” Acta Materialia 194, 93–105 (2020).
  * Roth (2016) Stephan V Roth, “A deep look into the spray coating process in real-time—the crucial role of x-rays,” Journal of Physics: Condensed Matter 28, 403003 (2016).
  * Su  _et al._ (2017) Bo Su, Volker Körstgens, Yuan Yao, David Magerl, Lin Song, Ezzeldin Metwalli, Sigrid Bernstorff, and Peter Müller-Buschbaum, “Pore size control of block copolymer-templated sol–gel-synthesized titania films deposited via spray coating,” Journal of Sol-Gel Science and Technology 81, 346–354 (2017).
  * Randolph  _et al._ (2022) L. Randolph, M. Banjafar, T. R. Preston, T. Yabuuchi, M. Makita, N. P. Dover, C. Rödel, S. Göde, Y. Inubushi, G. Jakob, J. Kaa, A. Kon, J. K. Koga, D. Ksenzov, T. Matsuoka, M. Nishiuchi, M. Paulus, F. Schon, K. Sueda, Y. Sentoku, T. Togashi, M. Bussmann, T. E. Cowan, M. Kläui, C. Fortmann-Grote, L. Huang, A. P. Mancuso, T. Kluge, C. Gutt, and M. Nakatsutsumi, “Nanoscale subsurface dynamics of solids upon high-intensity femtosecond laser irradiation observed by grazing-incidence x-ray scattering,” Phys. Rev. Research 4, 033038 (2022).
  * Holý  _et al._ (1993) V. Holý, J. Kuběna, I. Ohlídal, K. Lischka, and W. Plotz, “X-ray reflection from rough layered systems,” Phys. Rev. B 47, 15896–15903 (1993).
  * Sinha  _et al._ (1988) S. K. Sinha, E. B. Sirota, S. Garoff, and H. B. Stanley, “X-ray and neutron scattering from rough surfaces,” Phys. Rev. B 38, 2297–2311 (1988).
  * Schlomka  _et al._ (1995) J.-P. Schlomka, M. Tolan, L. Schwalowsky, O. H. Seeck, J. Stettner, and W. Press, “X-ray diffraction from si/ge layers: Diffuse scattering in the region of total external reflection,” Phys. Rev. B 51, 2311–2321 (1995).
  * Yabuuchi  _et al._ (2019) T. Yabuuchi, A. Kon, Y. Inubushi, T. Togahi, K. Sueda, T. Itoga, K. Nakajima, H. Habara, R. Kodama, H. Tomizawa, and M. Yabashi, “An experimental platform using high-power, high-intensity optical lasers with the hard X-ray free-electron laser at SACLA,” J. Sync. Rad. 26, 585–594 (2019).
  * Kameshima  _et al._ (2014) T. Kameshima, S. Ono, T. Kudo, K. Ozaki, Y. Kirihara, K. Kobayashi, Y. Inubushi, M. Yabashi, T. Horigome, A. Holland, K. Holland, D. Burt, H. Murao, and T. Hatsui, “Development of an X-ray pixel detector with multi-port charge-coupled device for X-ray free-electron laser experiments,” Review of Scientific Instruments 85, 033110 (2014).
  * Kiessig (1931) H. Kiessig, “Untersuchungen zur Totalreflexion von Röntgenstrahlen,” Annalen der Physik 10, 769 (1931).
  * Yoneda (1963) Y. Yoneda, “Anomalous Surface Reflection of X Rays,” Phys. Rev. 131, 2010–2013 (1963).
  * Pospelov  _et al._ (2020) Gennady Pospelov, Walter Van Herck, Jan Burle, Juan M. Carmona Loaiza, Céline Durniak, Jonathan M. Fisher, Marina Ganeva, Dmitry Yurov, and Joachim Wuttke, “BornAgain: software for simulating and fitting grazing-incidence small-angle scattering,” Journal of Applied Crystallography 53, 262–276 (2020).
  * Salditt  _et al._ (1994) T. Salditt, T. H. Metzger, and J. Peisl, “Kinetic roughness of amorphous multilayers studied by diffuse x-ray scattering,” Phys. Rev. Lett. 73, 2228–2231 (1994).
  * Nguyen  _et al._ (2018) Lan Nguyen, Teruo Hashimoto, Dmitri N. Zakharov, Eric A. Stach, Aidan P. Rooney, Benjamin Berkels, George E. Thompson, Sarah J. Haigh, and Tim L. Burnett, “Atomic-scale insights into the oxidation of aluminum,” ACS Applied Materials & Interfaces 10, 2230–2235 (2018), pMID: 29319290, https://doi.org/10.1021/acsami.7b17224 .
  * Sentoku and Kemp (2008) Y. Sentoku and A.J. Kemp, “Numerical methods for particle simulations at extreme densities and temperatures: Weighted particles, relativistic collisions and reduced currents,” Journal of Computational Physics 227, 6846 – 6861 (2008).
  * Paschke-Bruehl  _et al._ (2023) Franziska-Luise Paschke-Bruehl, Mohammadreza Banjafar, Marco Garten, Lingen Huang, Brian Edward Marré, Motoaki Nakatsutsumi, Lisa Randolph, Thomas E. Cowan, Ulrich Schramm, and Thomas Kluge, “Heating in multi-layer targets at ultra-high intensity laser irradiation and the impact of density oscillation,” New J. Phys. 25, 043024 (2023).
  * Eidmann  _et al._ (2000) K. Eidmann, J. Meyer-ter Vehn, T. Schlegel, and S. Hüller, “Hydrodynamic simulation of subpicosecond laser interaction with solid-density matter,” Phys. Rev. E 62, 1202–1214 (2000).
  * Ramis  _et al._ (2012) R Ramis, K Eidmann, J Meyer-ter Vehn, and S Hüller, “Multi-fs–a computer code for laser–plasma interaction in the femtosecond regime,” Computer Physics Communications 183, 637–655 (2012).
  * Faik  _et al._ (2018) S. Faik, A. Tauschwitz, and I. Iosilevskiy, “The equation of state package FEOS for high energy density matter,” Comput. Phys. Comm. 227, 117 – 125 (2018).