Resolving transient temperature and density during ultrafast laser ablation of aluminum

To understand the dynamics of ultrashort-pulse laser ablation, the interpretation of ultrafast time-resolved optical experiments is of utmost importance. To this end, spatiotemporally resolved pump-probe ellipsometry may be utilized to examine the transiently changing dielectric function of a material, particularly when compared to two-temperature model simulations. In this work, we introduce a consistent description of electronic transport as well dielectric function for bulk aluminum, which enables unambiguous quantitative predictions of transient temperature and density variations close to the surface after laser excitation. Potential contributions of these temperature and density fluctuations to the proposed optical model are investigated. We infer that after the thermal equilibrium of electrons and lattice within a few picoseconds, the real part of the dielectric function mostly follows a density decrease, accompanied by an early mechanical motion due to stress confinement. In contrast, the imaginary part is susceptible to a complicated interaction between time-varying collision frequency, plasma frequency, and a density dependency of the interband transitions. The models proposed in this study permit an outstanding quantitative prediction of the ultrashort-pulse laser ablation’s final state and transient observables. Consequently, it is anticipated that in the future, these models will provide a quantitative understanding of the dynamics and behavior of laser ablation.


Introduction
After the invention of ultrafast chirped-pulse-amplification lasers in the mid-1980s [1], ablation with ultrashort-pulse (USP) lasers has been the subject of extensive investigation.Understanding time-resolved optical pump-probe measurements is of utmost relevance in this regard, since these approaches constitute the key experimental tools for resolving the ultrafast transient material response [2,3].In addition to tackling industry-relevant concerns connected to laser micromachining, opportunities are provided for basic science to explore high-energy, nonequilibrium physical phenomena on picoseconds time scales [3].
In metals, free (valence) electrons absorb USP laserpulse radiation within the optical penetration depth [4].Excited states undergo relaxation via electron-electron collisions within a few fs to hundreds of fs [5], resulting in thermalization of the electronic subsystem and following a Fermi distribution.Consequently, the moderate heat capacity of valence electrons in metals gives rise to extremely high electron temperatures in the range of approximately 10•10 3 K, so that for pulse durations exceeding the electron-electron thermalization time, optical properties are influenced by varying electron collision frequencies.This is commonly described semi-classically for metals in terms of a temperature-dependent collision frequency within the pure Drude model [6], the more advanced Drude-critical point (DCP) [7], or the Drude-Lorentz (DL) model [8].Likewise, several attempts have been made to simulate these influences using density functional theory [9,10].
Once an electron temperature has been determined, energy dissipation within the electronic subsystem occurs via thermal diffusion into the bulk [11] and direct heating of the lattice by electron-phonon collisions.This phenomena are described by the two-temperature model (TTM) [12,13], which accounts for material-dependent parameters such as the thermal conductivity of electrons k e , the heat capacity of the electrons and the lattice c V,EL , c V,IO , and the electron-phonon coupling factor g [14].It is common to neglect the thermal conductivity of phonons in metals because of its indirect proportionality to the lattice temperature T IO [15] and its generally low contribution for good conducting pure metals.
Energy transfer from hot electrons to the cold lattice causes a rise in lattice temperature near the surface, in a depth which lies usually in the range of the optical penetration depth and the electronic diffusion length.Depending upon the pulse duration and absorbed energy density, a number of dynamic phenomena occur, such as isochoric heating if stress confinement is fulfilled, accompanied by photomechanical ablation (spallation) after mechanical relaxation of compressive stresses [16].At higher fluences or pulse durations, however, the ablation process becomes predominantly photothermal, with prominent liquid-vapor formation in a phase explosion, which occurs at lattice temperatures exceeding 90 % of the critical temperature [17].Rapid phase transitions, such as superheated solid liquefaction by heterogeneous melting, are also typically observable on timescales of a few picoseconds, even regardless of the applied fluence [18].These physical patterns have been identified by numerous simulation as well as experimental techniques and are well established [19,20].
Macroscopic representations of lattice dynamics are often derived from a hydrodynamic (HD) approach [21,22], albeit sometimes in the context of solid mechanics (SM) [23].TTM-HD expresses the TTM as a function of internal energy and extends it by the energy densities associated with a pressure increase.In addition, the system of equations incorporates two further continuity equations, particularly for the conservation of mass and momentum [22].Traditionally, in TTM-HD, multiphase equations of state (EOS) are utilized to calculate the appropriate thermodynamic state variables from the given temperatures and pressures [24].The key drawback of this macroscopic technique is that phase transitions must be compensated energetically and spallation must be artificially triggered.For the latter, simulation cells are segregated according to predetermined criteria [25][26][27].Depending upon the underlying material models, TTM-HD simulations closely reflect reality [22], and owing to the modest computational effort, TTM-HD is an excellent tool for drawing conclusions on the dynamics of the whole ablation process.
A plethora of optical pump-probe methods have been developed to study the USP ablation dynamics, including time-resolved reflectometry [28,29] or microscopy [30], surface ellipsometry [31,32], dual-angle reflectometry [33], as well as interferometry [34].The experimental observations lead to a greater comprehension of laser-induced phenomena, such as the nonequilibrium dynamics of electrons in bulk metals [5], semiconductors, or nanoparticles [35], structural alterations (melt transition) [36], and density variations owing to surface expansion [37].Nonetheless, the diversity of physical processes and fast nonequilibrium dynamics pose a tremendous challenge for modeling optical properties within these states.
Among this spectrum of techniques, pump-probe ellipsometry (PPE) stands out for its ability to directly measure the transient complex refractive index n & k (without assuming material-specific models) within the optical penetration depth of the sample [38,39].Thus, two quantities reflection and absorption are accessible at every instant, which, unlike pump-probe reflectometry, enable broad and quantitative conclusions regarding the free system variables to be derived (electron and lattice temperature as well as density).
Although an ultrashort laser pulse may induce temperature and density gradients near the materials surface, resulting in a non-homogeneous complex refractive index where Fresnel's equations are not strictly valid, the resulting optical response can be understood as an average within the optical penetration depth [40][41][42].
A previous work combining PPE experiments and simulations of copper revealed the underlying physical principles as well as the most significant influences on the transient optical response [43].Three unique results were attainable: (1) pulse absorption causes thermal excitation of d-band electrons, resulting in a steep increase in the transient refractive index; (2) material rarefaction directly affects the plasma frequency and is reflected in the transient extinction coefficient k, while the refractive index n directly correlates with the (temperature-dependent) collision frequency; and (3) the thermal electron-phonon equilibrium time was deduced from the transient minimum of the extinction coefficient.Prior studies involving time-resolved reflectometry (PPR) and optical response simulations provide comparable findings; however, their results remain in some aspects more qualitative with respect to the overall dynamics [28,44,45].Clearly, the synergy of experiments and physical modeling is essential for an experimentally validated and more complete description of USP laser ablation, especially for fluences approaching the ablation threshold.
Understanding the electronic properties is just as crucial as describing the motion of matter.Thermal conductivity and electron-phonon coupling have a significant influence on the energy deposition within the lattice, which initiates the aforementioned dynamics.Consequently, broad range models of the optical, transport, and thermodynamic properties of matter are required to adequately capture the properties of the irradiated target material as well as its time evolution.The majority of models account for certain mechanisms of physical phenomena, often missing additional effects that occur during USP laser-matter interaction.For instance, the electron-electron collision frequency ee is often employed to be proportional to the electron temperature via ee ∝ T 2 e [46].In contrast is the electron-ion collision frequency directly proportional to the lattice temperature T i according to ei ∝ T i [47].Both dependencies are in fact connected to the electrical conductivity rather than thermal conductivity and are classically limited to the low-temperature regime.Therefore, saturation effects related to screening effects at high temperatures are not considered.Also, these models lack a density dependency, especially important when the material suffers substantial dilution during USP laser ablation.
Therefore in this work, we propose a closed-form theory of electron-phonon coupling, electronic thermal conductivity, and optical properties of aluminum (Al), as prevailing in USP laser ablation, which incorporate the dependency on the phase, electron, as well as lattice temperatures andmost important-the density.We compare results from a one-dimensional TTM-HD simulation including our derived material models to PPE experiments conducted at fluences close to the ablation threshold, which serve as a benchmark.In this way, the electron temperature, lattice temperature, and density can be resolved from the measured optical response during ultrafast laser ablation, which yields quantitative understanding on the dynamic behavior of the studied system, and extends the up-to-date qualitative picture.

Sample preparation and experimental setup
High-purity (99.999%) bulk Al samples with an area of 1 cm × 2 cm and a thickness of 500 µm were used as ablation target, like in our previous works [48].The initial rough surface was sanded and subsequently polished with a polycrystalline diamond suspension with decreasing grain sizes from 6 µm to 1 µm.A colloidal SiO 2 solution (20 nm particle size) was applied as finish, yielding a measured root-mean-square surface roughness of R a < 14 nm.In this way, phenomena such as plasmonic excitation due to rough surfaces with subsequent ripple formation can be avoided even in the singlepulse regime [49], while also enabling precise measurements of crater depths and volumes.Laser ablation experiments were performed using a femtosecond-laser source (FemtoRegen, Spectra Physics) with a central wavelength of 1056 nm (FWHM = 5 nm) and a pulse duration of τ p = 525 fs.The ultrashort-laser pulses were emitted with a repetition rate of 500 Hz at a constant average power of 30 mW.The pulseto-pulse energy fluctuation was below 1 % and the beam quality was M 2 ≤ 1.3 (measured with BeamSquared, Ophir).The experimental ablation threshold and self-reflection of the pertinent data used for comparison with our simulations were taken from previous research [48].
The transient optical response was gauged using the same pump-probe-ellipsometry setup and laser system as in our previous study [38].Initially, the laser pulse was split into a pump and a probe pulse, which was frequency doubled (528 nm and 425 fs).The pump pulse was focused at perpendicular incidence on the sample with a plano-convex lens of f = 100 mm.The laser peak fluence F 0 on the sample was calculated via F 0 = (2•E p ) ⁄ (πw 0 2 ) and set by attenuating the pulse energy E p with a combination of a half wave plate and polarization beam splitter.The beam-waist radius (1/e 2 intensity decay) was measured to w 0 = 15.0(3)µm (Micros-potMonitor, PRIMES).The probe pulse was guided through an optical delay line introducing a temporal delay between pump and probe pulse on the substrate.The complex refractive index ni•k was measured using a rotating analyzer 665 Page 4 of 14 ellipsometry setup.The complex refractive index was calculated in the spatial center of the pump pulse, where the local fluence equals the peak fluence F 0 .A maximum delay time of ∆t = 20 ps was obtained for the PPE measurements with a temporal resolution close to the probe pulse duration.A detailed description of the experimental setup, the measurement procedure, and the image and data post-processing can be found in Ref. [38].To validate the results obtained from PPE, pump-probe reflectometry experiments were conducted using our setup, as detailed in Ref. [30,50], employing the same laser parameters and fluences at multiple of ablation thresholds.Equally to PPE, the relative reflection change ΔR∕R 0 from PPR was assessed at the center of the Gaussian beam profile.
The steady-state optical indices of Al as a function of the photon wavelength at ambient conditions were acquired with a Sentech SE850 ellipsometer at three angles of incident of 60°, 65° and 70°.The pseudo-dielectric function was calculated from the ellipsometric angles tan Ψ and cos Δ, providing the complex refractive index [51].A deviation from the initial steady-state values at the probe wavelength to the PPE measurements at negative delay times was less than 10% and might be related to the sensitive calibration of the PPE setup [39].

Modeling of thermal transport and optical properties
Compared to the used underlying theory, the material models for electronic transport and the dielectric function presented in this work contain numerous modifications.All electronic properties were calculated within the framework of semi-classical methods.First, electronic interactions were included by means of a mean-field theory, using the Thomas-Fermi screening potential [52], in both electronic thermal conductivity and electron-phonon coupling calculations.Also, the electron-phonon coupling was calculated in dependence of the metal density.For the calculation of thermal conductivity, the liquid phase was subjected to the same rigorous ansatz as the solid phase.The collision frequencies derived from the thermal conductivity are incorporated into the dielectric function calculation.Consequently, a self-contained coherent formalism for electron transport and optics is guaranteed.The contributions to the temperatureand density-dependent collision frequencies were weighted based on the results of spectrally resolved ellipsometric measurements at room temperature.In addition, a density dependence for the interband contribution in the underlying Drude-critical point model was introduced.In the following, each step is discussed separately.In the supplementary materials, a comprehensive derivation of the models used in this work is provided.

Electron-phonon coupling
The electron-phonon coupling factor g 0 in this work is calculated based on an approach proposed by Petrov and Inogamov et al. [53,54] and is depicted in Fig. 1a (solid black line).Hereby, the electronic density of states (DOS) and the probability for electrons interacting with longitudinal phonons were taken into account in the form of a Fourier transform of the screened Coulomb potential (see Suppl.Mat.S1).The obtained result of g 0 as a function of the electron temperature T e is compared to various references.At first glance (solid black line), there is an almost constant value for g 0 , with a slight increase toward higher temperatures.The constant behavior at T e < 5000 K is related to the type of the used Fermi screening potential, which gives different results in contrast to the Lindhard shielding used in Ref. [54].For higher temperatures, the curve is in a good agreement with various reference data.However, over the entire temperature range, a drastic difference is seen with Brown et al. [55] and Gorbunov et al. [56].This is attributed to the deviating calculation of the electron-phonon mass enhancement parameter λ and the second moment of the phonon spectrum <ω> 2 within the Eliashberg formalism [57], which are strongly sensitive on the used density functional theory methods [58].Finally, we implemented the density dependence of the electron-phonon coupling in the form of g(T e , x) = x g 0 , where the density ratio x is defined as x = ρ/ρ 0 with ρ 0 as the density at standard conditions [59].

Electronic thermal conductivity
The thermal conductivity of electrons due to the electron-electron scattering k ee and electrons with ions k ei is given by 1/k tot = 1/k ee + 1/k ei according to the Matthiessen rule [60].The electron-electron contribution is primarily not affected by the liquid or solid phase, while scattering of electrons by ions must be considered separately for both states [61].First, we calculated the electronic contribution k ee based on the solution of the Boltzmann transport equation in the relaxation time approximation using the efficient algorithm from Ref. [62].Hereby, the exchange interaction between electrons was neglected, while Thomas-Fermi shielding was taken into account.The electron-ion contribution k ei , on the other hand, comprises the mean-free-path of electron-ion collisions in the solid and liquid states.The variation of the mean-free-path under compression and rarefaction was derived from the density dependence of the cold pressure, corresponding to previous works [61,63].A guided derivation is given in the supplementary materials (Section S2). Figure 1b shows the total electronic thermal conductivity for the solid and liquid phases for selected ion temperatures T i as a function of the electron temperature T e .For the solid phase at T i = 293 K, a strong maximum is observed at the electron temperature of 7500 K, followed by a decrease and eventually passing into a slightly increasing plateau beginning at 20•10 3 K.For higher ion temperatures, the pronounced maximum smears out, accompanied by a general decrease of k tot (S Fig. 1b red dashes and blue dashdots).In the liquid phase, the thermal conductivity at equal electron temperatures is roughly by a factor of two lower than in the solid phase.A smooth continuous increase with electron temperature is found.Here, the strong influence of ion temperatures, especially for T e < 20•10 3 K, is not given, whereas the slight rise with increasing T e > 20•10 3 K is present, as it is the case for the solid phase.

Dielectric function and complex refraction index
Subsequently, we modeled the dielectric function of Al with the DCP model based on steady-state ellipsometric measurements at room temperature and the calculated refractive index n and the extinction coefficient k (s.Fig. 1c and supplement S3 for details on the DCP modeling).The parameters were obtained using a differential evolution algorithm with a nonlinear least-squares routine [68].The used optical model and resulting fit parameter are given in the supplementary materials (S Table 6).The spectral interval of our ellipsometry measurement spans from 350 nm (0.6 eV) to 2250 nm (3.1 eV) (Fig. 1c, open black squares).The fit to the measured spectrum is displayed in Fig. 1c with a blue line and is closely matching the measurement.It can be seen that the interband transition peak at 1.48 eV is well represented within the DCP model.It should be noted that, due to the asymmetric behavior of the given peak, two critical points around the central peak wavelength were assumed [69] (s. also supplementary materials, S3).
To model the optical properties at the higher intensity during laser ablation within the DCP model, we included a density dependence to the plasma frequency ω p and the interband contribution ϵ b (ω) of the dielectric function [70].In addition, the dependency on electron and lattice temperature, T e and T i , respectively, is taken into account for the effective collision frequency ν eff [6].For ω p and ϵ b (ω), the density dependence was considered via ω p 2 = x•ω p,0 2 and ϵ b (ω,x) = x•ϵ b,0 (ω).This model implies the linear polarizability of the atomic charge density, in analogy to the simple Drude-Lorentz oscillator.The effective collision frequency ν eff comprises contributions from electron-electron and electron-phonon collisions, ν ee and ν ei , respectively, which were calculated by means of the Drude relation from the heat conduction via the electrical resistivity, which is dependent on electron and lattice temperatures as well as the density [60].The electron-phonon collisions are thus given separately for the solid as well as the liquid phase, in analogy to the thermal conductivity.Since not every photon-electron collision leads to absorption, the electron-electron and electron-phonon collision frequencies in the DCP model are multiplied by an effective weighting factor η [62,71,72].We validated the modeling of the dielectric function using different experimental techniques: spectral resolved steadystate ellipsometry and time-integrated self-reflectance measurements at pump-pulse wavelength as well as PPE measurements at probe pulse wavelength.

Collision frequencies
Figure 1d shows the effective collision frequencies ν eff as a function of electron temperature for different lattice temperatures.According to the Fermi liquid theory, a dominating quadratic dependence is observed at low T e [73], whereas, at higher T e , a gradual increase due to shielding effects may be observed [54].In contrast, the contribution of the electronion collision frequency v ei in Fig. 1d (arrow) only yields a nearly constant upward shift.For the solid-to-liquid transition an abrupt increase of v ei is given, which reflects the major decrease of thermal conductivity during the solid-liquid transition (see Fig. 1b).The final dielectric function used in the TTM-HD simulation was calculated for the solid and liquid phases, respectively.The solid-liquid mixture region within the EOS was modeled using the Landau expression [74] for an equal solid-liquid volume fraction.

1D TTM-HD simulation model
Simulations were conducted using a one-dimensional twotemperature hydrodynamics model (TTM-HD) as provided by Povarnitsyn et al. [22].Further, a multiphase equation of state (EOS) was implemented to realistically describe phase states and spallation due to cell separation when the spallation condition is reached [24].This model allows for the simulation of the most crucial phenomena, from spallation to phase explosion via melting and vaporization, with accurate pulse absorption.
Multiple simulations were carried out for the 525 fs pulse duration and corresponding wavelengths for the pump and probe pulse.The maximum simulation time extends up to 20 ps following the pump impact, utilizing an internal time increment of 1 fs and a logarithmic spatial grid with a 1 nm resolution near the surface.
The optics in terms of absorption and reflection of the laser radiation were determined in the TTM-HD with methods of classical electrodynamics using Helmholtz equations [75].The dielectric constant is space-dependent (the target surface is at z = 0), so that the Helmholtz equations must be solved numerically by means of the transfer matrix method (TTM) [44].
We have implemented the transport and optical material parameters provided in the supplementary materials into the hydrodynamic code.To enable comparison of the simulation with the PPE measurement, a continuous probe pulse with a wavelength of 528 nm was applied, making it possible to evaluate the spatially and temporally resolved complex refractive index.These values were then averaged within the optical penetration depth (d opt = λ/4πk) to ensure an appropriate comparison with the experimental results of the PPE.In this manner, potential issues due to refractive index gradients within the optical penetration depth and, therefore, a possible non-Fresnel-like reflection are circumvented.A detailed verification of the methodology can be found in Supplementary Materials, Section S4.

Results
The simulations conducted in this study were performed with a pulse duration of τ p = 525 fs and a wavelength of λ = 1056 nm on a semi-infinite Al target with a thickness of 1 µm to facilitate comparisons with previously published experimental outcomes.Detailed analysis of the ablation dynamics is performed for two irradiated fluences, below (0.75 F thr ) as well as above the ablation threshold (1.5 F thr ).

Ablation energetics including time-integrated observables
The ablation threshold F thr is defined in the TTM-HD simulation as the minimum amount of energy density at which delamination of a molten layer can be observed, which is commonly referred to as spallation.In the TTM-HD simulation, the spallation criterion is fulfilled, when a negative pressure of -2 GPa is exceeded [76] in a metastable liquid phase.Once this criterion is satisfied, the hydrodynamic cells are separated in the code at the appropriate position, taking into account energy conservation due to formation of two new surfaces [22].In our simulations, we find an ablation threshold for Al at an irradiated fluence of 0.62 J/cm 2 , which is in perfect agreement with the previously reported experimental value of 0.63(3) J/cm 2 [48,77].
Page 7 of 14 665 In a next step, we used temporally integrated self-reflectance measurements of the pump pulse from Ref. [48] to validate the calculated laser energy absorption in the simulation.In this regard, experimental data have been previously reported for 0.75 F thr as well as 1.5 F thr with values of 94.1(1.7)% and 91.5(6) %, respectively, using our chosen laser parameter and sample [48].From the simulation, values of 93 % and 91.1 % are determined, which are in excellent agreement with the measurement.This confirms that the modeled dielectric function ϵ adequately predicts the laser energy deposition of the pump pulse in the Al sample.

Time-dependent observation of ultrafast changes in the complex refractive index
To understand the change in the dielectric function of the laser-irradiated Al, the measured transient complex refractive index (from ∆t = -5 ps to 20 ps) is meticulously compared to the TTM-HD simulations.Two fluences below and above the ablation threshold are examined in the following, namely 0.75 F thr and 1.5 F thr , respectively.The optical response of the probe pulse in the simulation is calculated for a continuous wave laser radiation with a constant intensity of 10 5 W/cm 2 (F abs ≈ 150 µJ/cm 2 ), which is several orders of magnitude below the pump-pulse intensity of 10 11 W/cm 2 (F abs ≈ 62 mJ/cm 2 ).Therefore, negligible thermal excitation due to the probe pulse is expected.Furthermore, this implies that the optical response in the simulation is instantaneous.Contrary, experimentally measured signals represent an averaged value over the probe pulse duration for each specific delay time ∆t.This results in a temporal broadening of the experimentally determined signals, which can be described as the convolution of temporal envelope of the probe pulse with the transient dielectric function.This leads to intrinsically expected slight deviations from the simulated data.For a better comparability, the change of the complex refractive index of the probe pulse with respect to its value at ambient condition in dependency of the delay time is highlighted in Fig. 2b and c) as ∆n = n-n 0 and ∆k = k-k 0 , where n 0 = 0.69 and k 0 = 5.16.Experimental PPE values are given as empty squares and filled circles (green) for fluences of 0.75 F thr (left column) and 1.5 F thr (right column), respectively.The pump pulse is indicated as brown filled area at the delay time zero.The almost constant experimental values of Δn ≈ 0 for negative delay times increase upon pulse irradiation.This rise corresponds approximately to an error function with the width of the pulse duration, which is indicating a linear interaction.At the lower fluence of 0.75 F thr (Fig. 2b,) a maximum of the ∆n curve is reached at about ∆t ≈ τ p , followed by a plateau.Here, ∆n increases by 0.13, corresponding to a relative change ∆n/n 0 of 20 %.For the fluence above the threshold of 1.5 F thr (Fig. 2c), a characteristic maximum with ∆n = 0.26 (∆n/n 0 = 35%) is formed at ∆t ≈ 1 ps, which is followed by a slight decay to a plateau with ∆n = 0.23.After reaching the plateau or maximum, respectively, both experiments exhibit a gentle falling slope, with a decrease of less than 5 % to 10 % for delay times of 20 ps.
For the imaginary part of the complex refractive index, an opposite behavior is observed.During the pump pulse, the curve of ∆k begins to decrease until a plateau is reached, whereby the decay time of approximately 2 ps is larger than the rise time of ∆n for both fluences.Exceeding delay times of about ∆t ≈ 5 ps, ∆k saturates at ∆k ≈ -0.5 (0.75 F thr ) and at ∆k ≈ -1 (1.5 F thr ).This corresponds to a relative fall-off of about 10 % and 28 %, respectively.
Figure 2a presents PPE and PPR measurements as ΔR/ R 0 transients in a direct comparison.At first glance, measurements for both 0.75 F thr and 1.5 F thr fluences exhibit an excellent correspondence, which confirms that both PPE and PPR inspect the optical properties within the optical penetration depth.In line with the PPE curves, the ΔR/ R 0 drops after pulse impact by about 5 % for 0.75 F thr and around 15 % for 1.5 F thr .Both curves reach saturation at about 5 ps and 10 ps for 0.75 F thr and 1.5 F thr , respectively.The PPR trend aligns well with measurements in Refs.[44,45] and similarly exhibits the same qualitative indications of ablation dynamics on a comparable time scale.These curves demonstrate a pronounced interaction between temperature effects on both electrons and the lattice, along with onset of early mechanical expansion.The quantification of this phenomenon will be performed using forthcoming simulation outcomes.
The calculated transient curves from the TTM-HD simulations are displayed in Fig. 2b for 0.75 F thr and Fig. 2c for 1.5 F thr as black solid lines.At first glance, there is an excellent agreement between simulation and experiments.Solely for the fluence of 1.5 F thr and delay times exceeding 2 ps, ∆k oscillates slightly around the experimental points.However, the simulated n & k transient curves around the zero delay time indicate faster changes in comparison to the experimental values.This is, as stated before, attributed to the finite temporal resolution of our pump-probe experiment.
An interesting feature in the measurement as well simulation becomes visible for 1.5 F thr with the real part of the refractive index forming a characteristic maximum in Δn at ∆t ≈ 1 ps (Fig. 2c).This maximum is well reproduced within the TTM-HD simulation, also when compared to the absence of a maximum at 0.75 F thr (Fig. 2b).
Changes in the complex refractive index within the DCP model are caused by the complex interplay of the plasma frequency ω P , effective collision frequency ν eff , and further the density dependence of interband transitions [43].However, this contributions may be separated and investigated in a more detail fashion by separately considering the real and complex part of the dielectric function, ϵ r and ϵ i , which are calculated from the measured complex refraction index as ϵ r = n 2 − k 2 and ϵ i = 2nk, respectively.

Discussion
As shown in the Results section, by comparing simulation and experiments in terms of energetic observables and transient complex refractive index, a realistic description of the ablation process is feasible within the presented framework.Henceforth, the optical response of laser-irradiated Al due to transient changes of the thermodynamic quantities T e , T i , and ρ is discussed for the fluences of 0.75 F thr and 1.5 F thr .A sophisticated understanding of the mutual influences of these quantities on the contributions to the dielectric function is of crucial interest for a quantitative understanding of the dynamics of ultrafast laser ablation.To this end, in analogy to the complex refractive index, the transient change of the real and imaginary part of the dielectric function is investigated, depicted as ∆ϵ r and ∆ϵ i for both fluences in Fig. 3a, b.
Along with the experimental data (open black squares and filled green circles), different solutions of the TTM-HD optical response are exhibited: i) neglecting the overall density dependence (ω P = ω P,0 , ϵ b = ϵ b,0 and ν eff = ν eff (T e , T i ), red dash-dotted line), ii) modeling of a density-dependent collision and plasma frequency (ω P = ω P (x), ϵ b = ϵ b,0 and ν eff = ν eff (x, T e , T i ), blue dashed line), and iii) density dependence of all parameters in the dielectric function ϵ (ω P = ω P (x), ϵ b = x ϵ b,0 and ν eff = ν eff (x, T e , T i ), black solid line).It should be noted that all three cases include the dependencies on the electron as well as ion temperatures in the collision frequency.
Figure 3 further depicts simulation outcomes for the electron and lattice temperatures calculated with TTM-HD within the optical penetration depth (s.Fig. 3c and d), as well as the density and surface displacement (s.Fig. 3e  and f).The experimental values for density and surface displacement were derived from PPE measurements using the referenced approximation [77].Figure 3c and d portrays temperature curves for typical TTM solutions with a slightly fluence-dependent electron-phonon coupling time ranging from 4 ps to 6 ps at 0.75 F thr and 1.5 F thr , respectively.The maximum electron temperature, which occurs approximately after half of the pulse duration τ p /2, is dependent on the fluence, with corresponding values of 9•10 3 K and 18 •10 3 K (s.Fig. 3c and d).In Fig. 3e and f, the densities within the optical penetration depth decrease at a slower rate than the increasing lattice temperature.This is due to the fulfilled condition of stress confinement, and therefore partial isochoric heating of the lattice with subsequent pressure relaxation.The density decrease after thermal equilibrium at about 5 ps is likewise proportional to lattice temperature and fluence, with relative decreases of approximately 20% and 40% for 0.75 F thr and 1.5 F thr , respectively.
For a better understanding of the ablation dynamics and the associated outcomes, three successive time regimes as known from literature (1-3 in Fig. 3c and d) are examined separately below.
Time regime 1 (∆t < τ p /2 )-heating of electrons and lattice: Considering the material response in the early stage of pulse irradiation at delay times below a half of the pulse duration (∆t < τ p /2), optical properties are dominated by the fast heating of the electronic subsystem.Hereby, electron temperatures up to 9•10 3 K (0.75 F thr ) and 18•10 3 K (1.5F thr ) are reached within the pump pulse maximum (see Fig. 3c).Comparing these values to the findings of the thinfilm experiments performed in Ref. [45] using a similar PPR setup in conjunction with simulations, it was found that the maximum T e for fluences of 0.75 F thr and 1.5 F thr is about 6.5•10 3 K and 16•10 3 K, respectively.This results are in robust agreement with our observations.With such high electron temperatures, the total collision frequency ν eff increases (see Fig. 1d), being dominated by electron-electron collisions ν ee .Contrary, the plasma frequency ω P as well as the bound part of the dielectric function ϵ b remain nearly unchanged, since the density does not change considerably in this time regime.Thus, the initial growth of ϵ r , driven by the temperature-dependent increase of ν eff , may be approximated within the DCP model by ∆ϵ r ≈ ω P,0 2 /ν eff,0 2 ω P,0 2 /ν eff 2 (T e , T i ) (with ω L 2 < ν eff 2 , s. supplementary materials S3).Note that the fraction ω P,0 2 /ν eff 2 is larger than unity at normal densities and electron temperatures up to about 10•10 3 K, and therefore, ∆ϵ r begins to rise with increasing ν eff .As a consequence, the maximum electron temperatures as well as the magnitudes of ∆ϵ r are directly proportional to the irradiated fluence and double for 1.5 F thr when compared to 0.75 F thr .
On the other hand, the imaginary part ϵ i of the dielectric function for the given conditions is indirectly proportional to ν eff with ∆ϵ i ≈ ω P,0 2 /ν eff,0ω P,0 2 /ν eff (T e , T i ).Thus, the curve of ∆ϵ i within this time regime is likewise increasing.As it is the case for ∆ϵ r , the change of ∆ϵ i is again nearly proportional to the electron temperature and applied fluence, respectively.Double electron temperature results in doubled change of real part ∆ϵ i in the pulse maximum at τ p /2.The absence of a density change for times Δt < τ p /2 is reflected in the coincidence of all three simulated PPE curves for the cases i) to iii), where different levels of density dependence in the DCP model were taken into account (s.Fig. 3a and b, red dot-dashed, blue dashed, and black solid lines).
Time regime 2 (τ p /2 < ∆t < 5 ps)-onset of the early mechanical motion: At τ p /2, in the simulation, the lattice temperature reaches about 1200 K, which is 1.3 times higher than the equilibrium melting temperature of T melt = 933 K. Subsequently, the overheated Al surface undergoes a rapid phase transition from the solid to liquid phase.The phasetransition timescale is noticeably shorter than the lattice thermalization time τ eq , which is in the range of 3 ps to 6.5 ps (see Fig. 3d and c).Guo et al. [79] reported that in such a case, Al may undergo a nonthermal, electronically induced, phase transition.However, this conclusion was later refuted by Kandyla et al. [36], who found that the melting transition in ultrafast laser-irradiated Al is purely thermal in nature.Ultrafast phase transitions with time scales with 1 ps to 2 ps were found by these authors.Since in the TTM-HD model, phase transitions are governed inter alia by the equation of state, rapid isochoric heating of the lattice will yield in a superheated solid phase with possible temperatures above 1.3 T melt .Subsequent ultrafast thermal melting by homogeneous nucleation is found in our case for delay times of less than 1.5 ps [18].
For USP processes, fulfilling the stress confinement condition, a quasi-isochoric heating of the lattice is induced, yielding a maximized pressure build-up within the heated region.The characteristic time for the mechanical relaxation of compressive stresses τ mech is typically in the order of 3 ps to 5 ps for liquid metals [48,80].However, an onset of mechanical motion may be observed at lower delay times as pump-probe experiments show [37,45].This transition region, where the dielectric function is affected by a complex interplay of rising lattice temperature, decreasing electron temperature, decreasing density, and solid-to-liquid-phase transitions, is marked as second time regime in Fig. 3c and d.Starting with the given temperatures at τ p /2, the increase of the lattice temperature T i due to the sustained electron-phonon coupling is continuous until the thermal equilibrium τ eq at 3 ps-4 ps for the fluence of 0.75 F thr and 6.5 ps for 1.5 F thr (s. Figure 3c,d).We find that for 0.75 F thr , τ eq is comparable with the value of 2 ps from Ref. [78].For 1.5 F thr , τ eq is in the range of 6.5 ps, which is nearly 100% longer than the sub-threshold value.This difference may be explained by an analytic approximation of τ eq = C e (T e )/g(T e ).Since C e (T e ) is almost linear in T e within the range of simulated electron temperatures, and g(T e ) remains mostly constant, τ eq should approximately increase proportional to the maximum T e , and thus laser fluences, respectively.This consideration is in excellent agreement with the simulated equilibration times: For 1.5 F thr , a T e maximum of 18•10 3 K corresponds to 6.5 ps, and for 0.75 F thr , the maximum T e of 9•10 3 K corresponds to τ eq = 3 ps -4 ps.
The decrease of the plasma frequency ω p due to the resulting density decrease is faster and stronger than the rise of ν eff , where at this point, ν eff follows mainly the lattice temperature.The influence of T i on ν eff is four-to-five times greater than that of T e (s.Fig 1d).At this point, the different cases of successively decreasing density influence in the optical response (given at the beginning of the discussion) exhibit a diverging progression.
First, the real part Δϵ r appears virtually unaffected by the density dependence of the bound contributions (cases ii) and iii) in Fig. 3a and b but depends primarily on the density dependence of the collision frequency and plasma frequency Page 11 of 14 665 (cases i) and iii)).Moreover, since ν eff 2 <<ω L 2 regains validity (decreasing electron temperatures and eventually electron-phonon equilibration), Δϵ r can be approximated by ω p 2 (x)/ω L 2 .Therefore, for delay times of a few picoseconds, the response of Δϵ r is exclusively governed by the rarefaction of the molten material and thus the decreasing plasma frequency ω p .
In our previous work, we have shown for Al (if n 2 ≪ k 2 is satisfied) that the changes in the real part of the dielectric function Δϵ r can be related to the density variation via the relation Δϵ r ≈ -Δk 2 ∝ ∆ρ [77].Based on the density, the transient surface displacement d surf is derived in terms of a surface velocity u disp with u disp ≈ c liq Δρ⁄ρ and d surf = ∫dt u surf [81], where c liq is the speed of sound in the liquid phase.In this way, we can indirectly determine the transient density through the measurement of the transient Δϵ r with PPE.
Figure 3e and f shows the so calculated density ρ averaged over the optical penetration depth of the probe pulse and the surface displacement d surf for irradiated fluences of 0.75 F thr (left column) and 1.5 F thr (right column).Solutions of the TTM-HD simulation are plotted as a black solid line and colored symbols refer to mentioned approximations from the experimental PPE data.We find an outstanding agreement between the transient density approximated from the experiments and the simulations.From this, we may conclude that starting at approximately ∆t = τ p /2, the material gradually expands due to previous isochoric heating of the lattice.As the time evolution continues, the material density reaches a local minimum with a decrease to about 80% (0.75 F thr ) and 60% (1.5 F thr ) at delay times of 5 ps (s.Fig. 3 e, f, upper row).This decrease in density results in a continuous surface expansion with a velocity of about 600 m/s (0.75 F thr ) and 1500 m/s (1.5 F thr ) due to the release of the laser-induced pressure build-up near the surface (s.Fig. 3 e, f, lower row, red highlighted).Here, we can also clearly observe that the n & k values measured with PPE represent the average value within the optical penetration depth, as the average density corresponds to that of the simulation.
For the imaginary part of the dielectric function change Δϵ i , however, a contrary effect is observed: First, taking into account the temperature dependence of Δϵ i on its own (case i), red dot-dashed line), there is a steady increase in Δϵ i , which enters saturation at about 2 ps.Taking density effects into account when calculating the collision and plasma frequency (case ii), blue dashed line), the behavior of Δϵ i alters slightly, although the same trend remains.For both of these cases i) and ii), no pronounced peak shortly after delay time zero is found, which is in significant contrast to the measurements.By adding a density dependence to the bound contribution via ϵ b = x ϵ b,0 within the optical model (black line, complete modeling as described for case iii), the experimental data are clearly reproduced.Here, the peak at approximately τ p /2 is very well captured for both fluences.The density influence at 1.5 F thr is more pronounced as expected from the previous analysis of the time-resolved density data.With this, it becomes apparent that the peak in Δϵ i at τ p is a result of the onset of a density decrease in the bound part of the dielectric function.It is also shown that the DCP model reproduces very well the strong opposite influence of density on the bound part in the real and imaginary part of the dielectric function.
The linear impact of the relative density on the bound contribution of the DCP model is particularly clearly illustrated in the Δϵ i measurements.Here we can observe that as the density is reduced to 80% (0.75 F thr ) or 60% (1.5 F thr ) of the original value, the Δϵ i value decreases by about 1 or 2, respectively.
Time regime 3 (∆t > 5 ps)-continuous surface bulging: After a delay time of 5 ps, no significant changes in the dielectric function may be observed.According to the discussion of time regime 2, the density dependence in the plasma frequency plays the central role in the real part of the dielectric function.Contrary, the bound part of the critical points affects preliminary the imaginary part of Δϵ.This observation is consistent with the finding in Ref. [43].An analysis of the temporal evolution of the surface density, given in Fig. 3e and f, shows that a local minimum density of 80% (0.75 F thr ) and 60% (1.5 F thr ) is reached at approximately 5 ps -6 ps, respectively.At this point, the electronic and ionic subsystems are in thermal equilibrium with maximum values of 1250 K at 3 ps -5 ps (0.75 F thr ) and 3000 K at about 6 ps (1.5 F thr ).A robust agreement to reported values in Refs.[45,82] is found, with 1500 K for F < F thr and 3000 K to 4000 K for above threshold fluences within the spallation regime.
From this times on, electronic heat conduction distributes the absorbed energy deeper into the material and the surface slowly begins to cool.Both processes, the propagation of the tensile stress wave and the surface cooling, cause the density of the material to gradually increase near the surface.At the end of our time scale, the surface displacement is up to 6 nm (0.75 F thr ) and 16 nm (1.5 F thr ), respectively.Our simulations show that the total displacement can add up to a few tens of nm for delay times of a few tens of ps, when the ablation process is initiated as removal of molten material by photomechanical spallation.Therefore, there is a strong interplay between an expanding surface and the material removal process.This results, inter alia, in efficient ablation when the pulse duration is shorter than the mechanical relaxation time, and a reduction in ablation efficiency when the laser pulse is applied on an already expanding surface [48,80,83].

Summary and conclusion
In this work, an already existing TTM-HD simulation has been extended by a new and consistent temperature and-more important-density-dependent model for the electronic thermal transport and optical response of bulk aluminum.This model has been extensively evaluated on experimental final state as well as time-resolved observables and has been found to be sufficient for a realistic prediction of the USP laser ablation process.
Having deciphered the individual contributions to the dielectric function, we conclude that the real part of the complex dielectric function ϵ r directly represents the density of the material within the optical penetration depth.On the other hand, we find that the imaginary part ϵ i exhibits a complex interplay of transient changing collision frequency, plasma frequency, as well the density-dependent interband transitions.The linear density dependence of the bound contribution is crucial for the interpretation of the time-resolved PPE experiments.
In addition, we show that the introduction of an effective weighting factor η in the summation of collision frequencies leads to an appropriate representation of the optical response over the entire visible spectrum.Without the factor η, the total collision frequency ν eff , derived from thermal conductivity and the electrical resistivity, would be overestimated due to the partial contribution of photon-electron scattering processes during absorption.Other research groups have previously suggested this factor [71,72], which is again supported by the results presented here.
Furthermore, it was found that conventional modeling of scattering rates based on Fermi liquid theory (T e and T i ) is insufficient and leads to an underestimation of thermal conductivity, especially at elevated temperatures.This would result in lower ablation thresholds in USP laser simulations when compared to experimental reference data.
Finally, our PPE experiments and simulations corroborate the previously described interpretation of the ablation dynamics within the stress confinement regime in three subsequent time regimes: namely ultrafast heating of electrons within the pulse duration, lattice heating during electron-phonon coupling within 2 ps -3 ps delay time, and at last, mechanical relaxation after 3 ps -5 ps, leading to the onset of mechanical motion and continuous surface expansion.The aluminum sample studied here displays somewhat attenuated stress confinement compared to stainless steel and other alloys, which indicate a stronger electron-phonon coupling and electron localization, reaching thermal equilibrium within 1 ps to 2 ps [48,80,84].However, the somewhat weaker electron-phonon coupling in Al here allows us to temporally separate all three components of stress confinement-electronic heating, lattice heating, and mechanical relaxation-by time-resolved ellipsometry as well as in the simulations.
Our findings and the proposed material models provide a great leap toward an experimentally validated model of ultrashort-pulse laser ablation, which in future can tackle important open questions such as pulse duration-dependent ablation efficiency.Furthermore, such a validated model will allow resolving temperature and density within the optical penetration depth of metals with ultrafast temporal resolution.

Fig. 1
Fig. 1 Summary of calculated parameters in this work.a Electron temperature dependence of the electron-phonon coupling parameter g 0[14,55,56,64,65] together with values at room temperature[66,67].b Total thermal conductivity k tot of electrons for solid and liquid phases at different ion temperatures as a function of the electron tem-

Fig. 2 a
Fig. 2 a Transient ΔR/R 0 for fluences of 0.75 F thr and 1.5 F thr obtained with PPE (black empty squares) and PPR (green filled circles) setups.b, c Transient real part of the complex refractive index ∆n and the extinction coefficient ∆k together with the correspond-

Fig. 3
Fig. 3 The experimentally measured real and imaginary parts of dielectric function for fluences of 0.75 F thr (black empty squares) in (a) and 1.5 F thr (green filled circles) in (b) and assuming various dependencies are illustrated as colored lines together with (c, d) electron and ion temperatures calculated by the TTM-HD simulation, and (e, f) density and transient surface displacement.Simulation and exper-