Modelling ultrafast laser ablation

Free electrons are present in the conduction band of a metal. Electrons below the Fermi edge can absorb photons and end at energies above the Fermi edge. The average energy increase due to intraband absorption can be described in terms of Drude's approach [17]. In this picture, electrons move in the electric field (here the field of the laser light, alternating with frequency ${{\omega}_{\text{L}}}$) and collisions, occurring with a frequency ${{\nu}_{\text{D}}}$, hamper the movement of the electrons. This leads to a dielectric function of

$$\begin{eqnarray}&&\varepsilon \left({{\omega}_{\text{L}}}\right)={{\varepsilon}_{r}}-{{\left(\frac{{{\omega}_{\text{P}}}}{{{\omega}_{\text{L}}}}\right)}^{2}}\frac{1}{1-\text{i}\,{{\nu}_{\text{D}}}/{{\omega}_{\text{L}}}},\end{eqnarray} \tag{ 4 }$$

where ${{\varepsilon}_{r}}$ is the dielectric constant of the unperturbed material and ${{\omega}_{\text{P}}}=\sqrt{{{n}_{\text{e}}}{{e}^{2}}/{{m}_{\text{e}}}{{\varepsilon}_{0}}}$ the plasma frequency, determined by the density of the free electrons n_e and their effective mass ${{m}_{\text{e}}}$. Note that the choice of the collision frequency ${{\nu}_{\text{D}}}$ varies in literature. Matthiesen's rule can be applied to include several types of collisions in this so-called Drude frequency. This frequency ${{\nu}_{\text{D}}}$ is roughly connected with the transport time τ_tr, appearing in equations (21) and (28); however, when interested in absorption, collisions among free electrons of the same effective mass should not be included since the collective momentum of this electronic system is not changed in such collision. Often a constant ${{\nu}_{\text{D}}}$ is assumed [18–21] which lies in the range of ${{10}^{15}}~\text{Hz}$, corresponding to a collision time in the femtosecond range. However, ${{\nu}_{\text{D}}}$ is likely to depend on material parameters which may change during irradiation. For instance, the electron–phonon collision time ${{\nu}_{\text{ep}}}$ depends on the lattice temperature T_i. For further increasing temperature, the limit of hot plasma may be reached. Here, the collision frequency is given by Spitzer's formula [22], where the collision frequency decreases with increasing electron temperature. A transition from the initially cold metal to a hot plasma was discussed in [23]. Moreover, since the Drude model also applies roughly for free electrons in semiconductors and dielectrics, the dependence of ${{\nu}_{\text{D}}}$ on the electron density n_e may become crucial [24]. Generally, it should be noted that v_D is rather a phenomenological than a microscopical parameter [26].

The square root of the dielectric function (4) equals the complex refraction index $\hat{n}$,

$$\begin{eqnarray}&&\hat{n}=n+\text{i}~k=\sqrt{\varepsilon \left({{\omega}_{\text{L}}}\right)},\end{eqnarray} \tag{ 5 }$$

where n is responsible for the refraction of light and k determines the attenuation of an electromagnetic wave in material. The reflectivity R and the free carrier absorption coefficient a_fca follow as

$$\begin{eqnarray}&&R=\frac{|\hat{n}-1{{|}^{2}}}{|\hat{n}+1{{|}^{2}}}\quad \text{and}\quad {{a}_{\text{fca}}}=2k{{\omega}_{\text{L}}}/c,\end{eqnarray} \tag{ 6 }$$

with c denoting the velocity of light in vacuum. In weakly excited metals, linear absorption can be assumed, where the absorbed energy per unit volume and time is proportional to the photon flux, i.e. the intensity I. This results in an exponential decay of laser energy in the material, which is described by the Lambert–Beer extinction law,

$$\begin{eqnarray}&&\begin{array}{*{20}{l}}I(z,t)={{I}_{\text{las}}}(t)(1-R)\,\exp \left(-z\,{{a}_{\text{abs}}}\right),\\ \quad \quad \ \;={{I}_{S}}(t)\,\exp \left(-z\,{{a}_{\text{abs}}}\right)\end{array}\end{eqnarray} \tag{ 7 }$$

where I_las is the laser intensity irradiating the material, I_S = I_las (1 – R) is the resulting intensity at the surface after reflection, and ${{a}_{\text{abs}}}={{a}_{\text{fca}}}$ the absorption coefficient. The source term in equation (2) is then given by

$$\begin{eqnarray}&&Q(z,t)={{a}_{\text{abs}}}\,{{I}_{S}}(t)\,\exp \left(-z\,{{a}_{\text{abs}}}\right)\end{eqnarray} \tag{ 8 }$$

Note that interband absorption mechanisms may also play a role in metals. In polyvalent materials like aluminum, nearly-parallel bands lead to an increased absorption at certain laser frequencies [25–28]. An even more complicated picture appears for noble metals, where the d band lies energetically in the s band of free electrons. Here, d electrons may be excited into the s band above the Fermi level with sufficiently high photon energies. Such transitions cannot be described in the frame of the Drude model and more sophisticated approaches have to be applied [29].

In practice, the strength of absorption, i.e. the parameter ${{a}_{\text{abs}}}$ in equations (7) and (8) is often taken from tabulated data [30], in particular when the behavior of a real, non-idealised material shall be modelled. Note, however, that even in metals this parameter may change considerably under strong electronic excitation [28, 29, 31].

Experimental works have shown an effective penetration of laser energy considerably larger than the optical light penetration: time-resolved front pump–rear probe experiments revealed that a thin layer of gold is heated homogeneously up to a thickness of approximately 100 nm [14]. This effect was attributed to ballistically moving electrons, heating a layer with a thickness of their free path ${{\lambda}_{\text{ball}}}\approx 100~\text{nm}$. The authors proposed to consider this effect in the TTM, applying in the source term (8) an effective penetration depth

$$\begin{eqnarray}&&{{a}_{\text{abs}}}={{\left(1/{{a}_{\text{fca}}}+{{\lambda}_{\text{ball}}}\right)}^{-1}}.\end{eqnarray} \tag{ 9 }$$

In [32] this expression was incorporated into the TTM for copper, silver and tungsten. The numerical results on ablation depths were compared with experimental observations. The authors found a ballistic depth of about 15 nm for copper and 50 nm for silver (and none for tungsten). The inclusion of an effective absorption coefficient according to equation (9) has been found to be important only for fluences close to ablation threshold. This finding is consistent with an estimation connecting the ballistic range with the time of free flight of the electrons (entering transport, see equation (28)), leading to a decreasing ${{\lambda}_{\text{ball}}}$ for increasing temperatures [33]. Generally, the concept of ballistic electrons is not directly compatible with a temperature-based description. Therefore, expression (9) seems to be practical—however, it is questionable in the kinetic view.

The description of laser light absorption in semiconductors has to take into account changes in the density of free electrons in the conduction band of the material. Since conduction band and valence band are separated by a bandgap which is typically in the range of the photon energy of visible light, the generation of an electron–hole plasma upon irradiation plays a crucial role [18]. The absorption of laser energy by free electrons in the conduction band can be described with the Drude formalism as sketched in the previous section. Also absorption by holes can be described in this framework [18]. Note that the plasma frequency ${{\omega}_{\text{P}}}$ in equation (4) depends on the density n_e of free carriers (i.e. electrons and holes).

The attenuation of the laser intensity I in the material is strongly connected with the excitation of the free electron density n_e (or, equivalently, the density of the electron–hole plasma ${{n}_{\text{e}-\text{h}}}$). In many cases, single photon absorption across the band gap has to be considered, as well as two photon absorption. For instance in silicon irradiated with visible light, two photon absorption may be a direct transition while single photon absorption is indirect and thus accompanied by a momentum change. Depending on photon energy, three photon absorption may also be of considerable importance, but is neglected in the formula presented in the following. The attenuation of laser intensity by single- and two photon absorption as well as free carrier absorption is given by [18, 24, 34, 35]

$$\begin{eqnarray}&&\frac{\partial I(z,t)}{\partial z}=-({{a}_{0}}+{{a}_{\text{fca}}}(z,t))~I(z,t)-{{b}_{0}}~I{{(z,t)}^{2}},\end{eqnarray} \tag{ 10 }$$

where a₀ denotes the single photon absorption coefficient and b₀ the two-photon absorption coefficient. Free-carrier absorption can be identified with Drude absorption, compare equations (4) to (6). However, note that the dielectric function $\varepsilon ({{\omega}_{\text{L}}})$ may consist of more terms than equation (4) for excited semiconductors [18, 24]. In case two-photon absorption can be neglected (i.e. b₀ is small in equation (10)), the Lambert–Beer extinction law (7) applies with ${{a}_{\text{abs}}}={{a}_{0}}+{{a}_{\text{fca}}}$. Note that ${{\alpha}_{\text{fca}}}$ may strongly change during irradiation [24].

The increase of free electron density due to these direct laser-induced excitation processes is given by [18]

$$\begin{eqnarray}&&\frac{\partial {{n}_{\text{e}}}(z,t)}{\partial t}=\frac{{{a}_{0}}}{\hbar {{\omega}_{\text{L}}}}~I(z,t)+\frac{{{b}_{0}}}{2~\hbar {{\omega}_{\text{L}}}}~I{{(z,t)}^{2}}.\end{eqnarray} \tag{ 11 }$$

For the complete description of the transient free electron density in laser-excited semiconductors, additional excitation and recombination processes have to be included. In [34], electron impact ionisation and Auger recombination have entered the balance, thus

$$\begin{eqnarray}&&\frac{\partial {{n}_{\text{e}}}(z,t)}{\partial t}=\frac{{{a}_{0}}}{\hbar {{\omega}_{\text{L}}}}\,I(z,t)+\frac{{{b}_{0}}}{2\,\hbar {{\omega}_{\text{L}}}}\,I{{(z,t)}^{2}}+{{\delta}_{\text{imp}}}{{n}_{\text{e}}}(z,t)-{{\gamma}_{A}}{{n}_{\text{e}}}{{(z,t)}^{3}},\end{eqnarray} \tag{ 12 }$$

where ${{\delta}_{\text{imp}}}$ denotes the impact ionisation coefficient and ${{\gamma}_{A}}$ the coefficient of three-body Auger recombination. When particle transport comes into play, further processes influence the transient local density of the electron–hole plasma [24, 34, 35]. The resulting modified description of the energy density similar to the TTM, equations (2) and (3), will be discussed in section 4.4. Note that the description of the optical parameters according to equation (6) leads to a very good comparability of calculated melting thresholds with experimental values for a large range of pulse durations [24], once a density- and temperature dependent Drude collision frequency is taken into account.

A modified description of free electron density evolution in semiconductors tracing also nonequilibrium electron distributions has been presented in [36]. It is based on the multiple rate equation for dielectrics [37, 38], which will be described in the following section. In the extended multiple rate equation [36], Auger recombination processes and electron–phonon coupling have been implemented additionally. The results have shown good agreement with experimental findings [18].

Strong excitation of semiconductors may considerably change their optical properties. In [39–41], it was reported that GaAs turns into a metallic state after irradiation with a femtosecond laser pulse. The underlying ultrafast phase transitions are discussed in section 5.2 of this review, with a focus on the loss of crystalline order. Note that the accompanying change of the optical properties may influence the energy input for sufficiently long laser irradiation.

In wide-band gap dielectrics, absorption of visible light is possible only when high intensities are applied, allowing for nonlinear optical processes. The studies of laser-irradiation of dielectrics in the context of material ablation focus mainly on the time-resolved description of the free electron density. The common argumentation includes that for increasing electron density the real part of the dielectric function (4) may vanish and further become negative, leading to an imaginary refractive index (5) accompanied by a strong metallic-like absorption, compare section 2.1. The critical density n_crit of this transition can be calculated under the assumption of ${{\omega}_{\text{L}}}/{{\nu}_{\text{D}}}\gg 1$ to

$$\begin{eqnarray}&&{{n}_{\text{e},\text{crit}}}={{\varepsilon}_{0}}~{{\varepsilon}_{r}}\frac{{{m}_{\text{e}}}\omega _{\text{L}}^{2}}{{{e}^{2}}}.\end{eqnarray} \tag{ 13 }$$

Figure 2 shows the dependence of optical parameters on the free electron density n_e [42]. A laser with a wavelength of $\lambda =500~\text{nm}$ (corresponding to a frequency of ${{\omega}_{\text{L}}}=3.77\times {{10}^{15}}~\text{Hz}$) and intensity ${{I}_{0}}=9.7\times {{10}^{13}}~\text{W}~\text{c}{{\text{m}}^{-2}}$ has been assumed to irradiate on a quartz-like material with ${{\varepsilon}_{r}}=2.4$. The Drude collision frequency was taken as ${{\nu}_{\text{D}}}={{(1~\text{fs})}^{-1}}$ and the effective mass m_e was assumed to equal the free electron mass. Note that in this case the assumption ${{\omega}_{\text{L}}}/{{\nu}_{\text{D}}}\gg 1$ is not fulfilled and the considerable increase of absorption, given by $k=\text{Im}(\hat{n})$, is a smooth function rather than resembling a threshold behaviour.

Zoom In Zoom Out Reset image size

Note that besides the free-electron response following from its dielectric function (4) other interesting optical phenomena occur in laser-irradiated dielectrics. For instance, the light attenuation (as described for semiconductors with equation (10)) includes self-focussing and defocussing, which are nonlinear optical processes becoming important for high intensities and high free electron densities [3, 43]. Propagation of the laser in certain parameter ranges may therefore result in long-distance filamentation, observed in various materials [44–47].

The increase of free electron density in the initially empty conduction band of a dielectric is triggered first by photoionisation from lower bands. During multiphoton ionisation, many photons are absorbed in one step. For very high intensities this process transfers to tunnel ionisation [48]. The transition is characterised by the so-called Keldysh parameter ${{\gamma}_{\text{K}}}={{\omega}_{\text{L}}}~\sqrt{{{m}_{\text{red}}}{{\epsilon}_{\text{gap}}}}/e{{E}_{\text{L}}}$, where m_red is the reduced mass of electrons and holes, ${{\epsilon}_{\text{gap}}}$ the band gap and ${{E}_{\text{L}}}$ the electric laser field amplitude. The Keldysh parameter is essentially determined by the ratio of the band gap energy ${{\epsilon}_{\text{gap}}}$ to the ponderomotive energy of the electron oscillating in the electric laser field, ${{U}_{\text{P}}}={{e}^{2}}E_{\text{L}}^{2}/\left(4{{m}_{\text{red}}}\omega _{\text{L}}^{2}\right)$. For a large band gap and small electrical field amplitude (i.e. small ponderomotive energy) many photons have to be absorbed. This multiphoton regime is characterised by ${{\gamma}_{\text{K}}}\gg 1$. When the ponderomotive energy exceeds by far the band gap, ${{\gamma}_{\text{K}}}\ll 1$, the latter becomes unimportant and electrons can appear in the conduction band. This process may be interpreted as a tunnelling process across potentials bent by the electric laser field. The ionisation rate for the process of photoionisation in strong laser fields was calculated by Keldysh in [48]. For various material and laser parameters, it is plotted for instance in [4, 43, 49]. Approximative expressions describe either the multiphoton- or the tunnelling regime [48]. In [49], a simplified complete description as a weighted transition between both regimes was proposed. Modified derivations have been proposed in [50, 51]. Experimentally determined values of the multiphoton ionisation cross sections are also available [52, 53].

When electrons are present in the conduction band of the dielectric, they may enable further electrons to overcome the band gap [54, 55]. This is the case for free electrons with a high energy above a certain critical energy ${{\epsilon}_{\text{crit}}}$, which is on the order of the band gap energy [49]. The corresponding impact ionisation rate may be assumed to directly depend on the free electron density, compare also equation (12). Considering photoionisation with a rate ${{\dot{{n}}}_{\text{pi}}}$ and impact ionisation with a probability ${{\delta}_{\text{imp}}}$, the increase of free electron density in dielectrics is given by the addition of both contributions [55, 56],

$$\begin{eqnarray}&&\frac{\partial {{n}_{\text{e}}}}{\partial t}={{\dot{{n}}}_{\text{pi}}}\left({{E}_{\text{L}}}\right)+{{\delta}_{\text{imp}}}\left({{E}_{\text{L}}}\right)~{{n}_{\text{e}}}.\end{eqnarray} \tag{ 14 }$$

However, experimental studies applying equation (14) have led to contradictory results [19, 52, 53, 56, 57]. Moreover, theoretical investigations state fundamental doubts as to whether this standard rate equation is applicable in general in the subpicosecond time regime [49, 58, 59]. One basic assumption of equation (14) is that impact ionisation depends directly on the total density of the free electrons. However, since impact ionisation needs a certain critical energy of the ionizing electron, this process also depends on the energy of a particular electron in the conduction band. While photoionisation generates electrons with low kinetic energy in the conduction band, impact ionisation requires electrons of high kinetic energy. This additional energy is absorbed from the laser light by intraband absorption. If this absorption process takes time comparable to the laser pulse duration, it is obvious that equation (14) is oversimplified. Such discrepancy will be reflected in a non-stationary shape of the free electron energy distribution. Examples of such distributions on ultrashort timescales, obtained with kinetic approaches [49, 60], are shown in the following section 2.4. Note that in [61] a description similar to equation (14) was proposed, where the intensity-dependence of the avalanche parameter was discussed in more detail. The idea of the authors is that impact ionization may be assisted by further photon absorption. For sufficiently high intensities such a process seems to be likely; however, its inclusion into the parameter ${{\delta}_{\text{imp}}}$ still does not capture the time dependence of the electrons' energy increase. Another qualitative way to include such time dependence into the description of impact ionization has been proposed in [3]. It is based on an empirical delay time and applies the electron density at earlier times to calculate the current rate of density increase by impact ionization.

In a microscopic picture, the energy distribution of the electrons has to be traced. The possibility to keep track of this distribution without applying a complex full kinetic approach has been introduced with the multiple rate equation (MRE) in [37, 38]. Here, the rate of intraband absorption W_1pt and thus the time of energy increase by one photon energy has been included in a description maintaining the conceptual simplicity of the rate equation (14). The MRE introduces virtual intermediate levels of energy in the conduction band and counts the density of electrons at these levels. Photoionisation provides electrons at the lowest energy, increasing the density n₀, while the density of electrons with energies above critical energy for impact ionisation ${{\epsilon}_{\text{crit}}}$ is denoted as n_k. The number of virtual levels k is given by the integer above ${{\epsilon}_{\text{crit}}}/\hbar {{\omega}_{\text{L}}}+1$. The multiple rate equation then reads

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{\overset{\cdot}{{n}}\,}_{0}} & = & {{\overset{\cdot}{{n}}\,}_{\text{pi}}}\left({{E}_{\text{L}}}\right)+2{{\tilde{{\delta}}\,}_{\text{imp}}}\,{{n}_{k}}-{{W}_{1\text{pt}}}\left({{E}_{\text{L}}}\right){{n}_{0}} \\ {{\overset{\cdot}{{n}}\,}_{1}} & = & {{W}_{1\text{pt}}}\left({{E}_{\text{L}}}\right){{n}_{0}}-{{W}_{1\text{pt}}}\left({{E}_{\text{L}}}\right){{n}_{1}} \\ {} & \vdots & \\ {{\overset{\cdot}{{n}}\,}_{k-1}} & = & {{W}_{1\text{pt}}}\left({{E}_{\text{L}}}\right){{n}_{k-2}}-{{W}_{1\text{pt}}}\left({{E}_{\text{L}}}\right){{n}_{k-1}} \\ {{\overset{\cdot}{{n}}\,}_{k}} & = & {{W}_{1\text{pt}}}\left({{E}_{\text{L}}}\right){{n}_{k-1}}-{{\tilde{{\delta}}\,}_{\text{imp}}}\,{{n}_{k}}. \\\end{array}\end{eqnarray} \tag{ 15 }$$

Here, ${{\tilde{\delta}}_{\text{imp}}}$ denotes the pure probability of impact ionisation in case an electron with energy $\epsilon >{{\epsilon}_{\text{crit}}}$ is present. The absorption probability W_1pt depends on the laser electric field amplitude E_L. The inclusion of a possible dependence of ${{W}_{1\text{pt}}}$ on the electrons' energy was discussed in [38]. The free electron absorption probability W_1pt can be roughly described with the Drude model, compare equations (4)–(6). Note again that the choice of the collision frequency ${{\nu}_{\text{D}}}$ is crucial for quantitative comparisons [24], and that different expressions and dependencies can be found in literature [18–21, 23, 24, 28].

Summing up the MRE, we obtain the rate of increase of the total free electron density n_e, which is given by

$$\begin{eqnarray}&&\frac{\partial {{n}_{\text{e}}}}{\partial t}={{\dot{{n}}}_{\text{pi}}}\left({{E}_{\text{L}}}\right)~+~{{\tilde{\delta}}_{\text{imp}}}~{{n}_{k}}.\end{eqnarray} \tag{ 16 }$$

On timescales larger than the timescale of the intraband absorption process, the shape of the electron distribution does not change significantly and the fraction of high-energy electrons ${{n}_{k}}/{{n}_{\text{e}}}$ is temporally constant. Then, equation (16) reduces to the standard rate equation (14) with ${{\delta}_{\text{imp}}}\left({{E}_{\text{L}}}\right)={{\tilde{\delta}}_{\text{imp}}}~{{n}_{k}}/{{n}_{\text{e}}}$. For the non-stationary case, the fraction of high-energy electrons changes with time and the difference in the last term of equation (16) as compared to equation (14) is substantial. The transition from the non-stationary regime on ultrashort timescales to the asymptotic avalanche regime for longer timescales was calculated in [38]. In figure 3, it is marked for constant laser intensity I of pulse duration ${{\tau}_{\text{L}}}$ and photon energy $\hbar {{\omega}_{\text{L}}}=2.48~\text{eV}$ as the line ${{t}_{\text{MRE}}}$ in the $I-{{\tau}_{\text{L}}}$ plane. The intensity refers to the intensity inside the material. Material parameters taken from SiO₂ have been applied [38]. The line t_MRE separates the parameter region of ultrashort irradiation where multiphoton ionisation dominates from that for longer pulses, where the collisional avalanche dominates the free electron generation.

Zoom In Zoom Out Reset image size

Several improvements and applications of the MRE (15) have been published since its introduction. It has successfully been applied to interpret experimental observations on laser damage and ablation with shaped laser pulses [62, 63]. An intensive discussion of the model was published by Christensen and Balling in [64]. In this work, energy conservation in the process of impact ionisation has been improved, the changing optical parameters during irradiation have been included and a depth-dependent analysis of the dielectric breakdown and ablation threshold has been performed. Also the limitations of the model in direct comparison with experiments have been discussed. For instance, the inclusion of defects would open the way to the description of more realistic material. For several applications also multishot experiments would be of interest. Other refinements mentioned in [64] as being necessary have meanwhile been performed, in particular the inclusion of fast recombination processes [42] as well as of nonlinear light propagation [43, 65] and relaxation processes [36, 43]. Note that when considering longer timescales up to the picosecond regime, recombination processes are of increasing importance. Particularly in SiO₂, self-trapped excitons (STE) lead to a decrease of the free electron density on a timescale of only 150 fs [19, 66, 67]. Such recombination processes can be included in the rate equation with an additional term $-{{n}_{\text{e}}}/{{\tau}_{\text{rec}}}$, where ${{n}_{\text{e}}}$ is the density of free electrons (or of the virtual level in the MRE) and ${{\tau}_{\text{rec}}}$ is the characteristic recombination time [42, 44, 68, 69].

It is important to note that apart from the large amount of publications considering the free electron density as the crucial parameter for laser damage of dielectrics, finally the absorbed energy is responsible for the damage of the material. Both thresholds are connected [70, 71], however, situations may occur where n_e exceeds n_crit at the end of the pulse duration and no damage will be observed [60]. In [72] two criteria have been applied to calculated ablation depths, i.e. a criterion of electron density as well as of energy density. Both yield indistinguishable results for a broad range of incident fluences. In [72] it was shown that a model based on the MRE (15) is capable to reproduce all experimentally obtained quantities like change of phase shift, absorption, reflectivity and final depth of the ablation crater with one single set of model parameters.

As indicated in the previous section, the energy distribution of laser-excited electrons may determine the appropriate description of the absorption. It may also determine further energy dissipation processes like electron–phonon coupling and transport. The description in terms of a temperature, like in equation (2), might not be justified on ultrashort timescales, where strong nonequilibrium distributions can occur. The knowledge of the particular distribution may also be relevant for the interpretation of experimental data [73], especially on ultrashort timescales, where equilibrium distributions should not be assumed from the first. In this section, we show some examples of nonequilibrium distributions, review their influence on measurable quantities and the timescales of their thermalisation.

Kinetic approaches like the Boltzmann equation or asymptotic trajectory Monte Carlo simulations trace the electronic distribution function after ultrafast excitation in detail [49, 60, 73–77]. Here, the collision probability at each electron's energy enters the approach and no assumption on the energy distribution has to be made. Kinetic approaches are therefore applicable for the description of absorption and relaxation processes far from thermodynamic equilibrium, when no temperature can be defined. Figure 4 shows the energy distribution of electrons in the conduction band of a SiO₂-like material during irradiation with a laser pulse of constant intensity. Details of the calculation are given in [49]. It can be seen that after a few tens of femtoseconds a smooth shape of the distribution is reached. This coincides with the calculations in [60] finding thermalisation times for electron distributions in dielectrics in the range of several tens of femtoseconds. Note that the thermalisation time reflects the smoothing of the distribution function towards a Fermi distribution (or, for sufficiently low densities, towards a Maxwell–Boltzmann distribution), whereas the transition time in figure 3 refers to the influence of two different ionisation processes, increasing the amount of electrons in this distribution. Ionisation and thermalisation are independent processes; for the given intensity in figure 4 the transition time is larger than the thermalisation time.

Zoom In Zoom Out Reset image size

In [78], the excitation and relaxation of different metals have been studied with help of Boltzmann collision integrals. Details of the density of states have been implemented in the calculation through an effective one-band model; they are reflected in the excited electron distributions [78]. The entropy of the excited distribution has been monitored after excitation. It increases to an asymptotic maximum value; the timescale of this increase can be identified with the thermalisation time ${{\tau}_{\text{therm}}}$. Figure 5 shows the results for aluminum, gold and nickel. For nickel two curves are shown, which refer to two calculations with different screening parameters. The screening of the Coulomb potential by free electrons is important for the collision probability and thus the interaction strength of electrons with each other. The Coulomb potential decreases exponentially with a characteristic length $1/\kappa$. The so-called screening parameter κ can be directly calculated from the nonequilibrium energy distribution [79, 80]. For a thermalised distribution, i.e. in metals a Fermi distribution, the screening is known as Thomas–Fermi-screening, here denoted by the parameter ${{\kappa}_{\text{therm}}}$. Figure 5 shows that the nonequilibrium screening remarkably increases the resulting thermalisation time. The influence of the screening length on ${{\tau}_{\text{therm}}}$ is roughly a fourth power law, as the inset of the figure shows [78].

Zoom In Zoom Out Reset image size

For thermalised distributions, the electrons can be described in terms of temperature-based equations like (2). For strongly disturbed electron distributions such temperature ${{T}_{\text{e}}}$ is not defined from the first. Only when the electrons thermalise to a new equilibrium distribution, the description in terms of an electronic temperature is justified.

However, kinetic descriptions as mentioned in this section are feasible only when spatial transport can be neglected. This is the case for homogeneously heated thin films or if only timescales below the picosecond regime are relevant. Since thermalisation is fast for high-energy excitation (below 100 fs for the considered cases, see figure 5), it is reasonable to assume that the TTM, equations (2) and (3), is applicable to describe ultrafast laser ablation. In case the nonequilibrium of the electron system leads to modifications of the energy absorption or dissipation processes, they may often be included in material parameters. Examples are, for instance, the modified absorption coefficient, equation (9), to capture ballistic transport, or a time-dependent electron–phonon coupling strength, which will be discussed in more detail in section 3.3.

If a thin film is irradiated with an ultrashort laser pulse, the energy is distributed by ballistic electrons on a femtosecond timescale and the film is homogeneously heated [14]. In this case thermal conduction can be neglected and the equation of heat exchange, equation (1), can be applied. Here, we assume that thermalisation is a fast process and the electronic system can be described by a temperature T_e. It follows from equation (1) that the characteristic time of electron gas cooling due to energy exchange with the lattice is ${{\tau}_{\text{e}}}={{c}_{\text{e}}}/\alpha$, and the characteristic time of lattice heating is ${{\tau}_{i}}={{c}_{\text{i}}}/\alpha$. The value of ${{\tau}_{\text{e}}}$ depends on electron temperature and ranges typically from 0.01 to 1 ps. Both characteristic times apply for the case that the temperature of the other subsystem is fixed. In case of electron cooling in laser-irradiated metals, it is ${{c}_{\text{e}}}\ll {{c}_{\text{i}}}$ and the assumption of a constant lattice temperature is reasonable as long as ${{T}_{\text{e}}}\gg {{T}_{\text{i}}}$. If no further energy input is assumed and for ${{T}_{\text{e}}}\gg {{T}_{\text{i}}}$ and constant parameters c_e, c_i and α, equation (1) leads to a characteristic relaxation time of the temperature equilibration of [82]

$$\begin{eqnarray}&&{{\tau}_{\text{relax}}}=\frac{{{c}_{\text{e}}}\,{{c}_{\text{i}}}}{\left({{c}_{\text{e}}}+{{c}_{\text{i}}}\right)\,\alpha}.\end{eqnarray} \tag{ 17 }$$

Since usually ${{c}_{\text{e}}}\ll {{c}_{\text{i}}}$ this time reduces to ${{\tau}_{\text{relax}}}\simeq {{c}_{\text{e}}}/\alpha$.

Assuming that the lattice temperature is initially equal to room temperature ${{T}_{\text{i},\text{init}}}={{T}_{0}}$ while the initial electron temperature T_e,init is much larger, the solution of equation (1) reads

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{T}_{\text{e}}}(t) & = & {{\Theta}_{\text{e}}}\exp \left(-\frac{\alpha \left({{c}_{\text{i}}}+{{c}_{\text{e}}}\right)}{{{c}_{\text{e}}}{{c}_{\text{i}}}}\cdot t\right)+\mathcal{T} \end{array}\end{eqnarray} \tag{ 18 }$$

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{T}_{\text{i}}}(t) & = & {{\Theta}_{\text{i}}}\exp \left(-\frac{\alpha \left({{c}_{\text{i}}}+{{c}_{\text{e}}}\right)}{{{c}_{\text{e}}}{{c}_{\text{i}}}}\cdot t\right)+\mathcal{T}, \end{array}\end{eqnarray} \tag{ 19 }$$

where $\mathcal{T}={{c}_{\text{i}}}/\left({{c}_{\text{i}}}+{{c}_{\text{e}}}\right)\cdot {{T}_{\text{i},\text{init}}}+{{c}_{\text{e}}}/\left({{c}_{\text{i}}}+{{c}_{\text{e}}}\right)\cdot {{T}_{\text{e},\text{init}}}\,\gtrsim \,{{T}_{0}}\,$ is the joint asymptotic temperature of both systems. Here, ${{ \Theta }_{\text{e}}}=$ ${{T}_{\text{e},\text{init}}}-\mathcal{T}\approx {{T}_{\text{e},\text{init}}}\gg |{{ \Theta }_{\text{i}}}|$ and ${{ \Theta }_{\text{i}}}=\left(1-\left({{c}_{\text{i}}}+{{c}_{\text{e}}}\right)/{{c}_{\text{i}}}\right)\,{{ \Theta }_{\text{e}}}\lesssim 0$ [83]. The exponential function has thus a much higher influence on the electron cooling, while the change of lattice temperature remains small. For times  ∼${{c}_{\text{e}}}/\alpha$, the neglection of the change in lattice temperature is thus a reasonable approximation. With the heat capacity of a Fermi distributed electron gas given as ${{c}_{\text{e}}}=\gamma {{T}_{\text{e}}}$, we end at

$$\begin{eqnarray}&&{{c}_{\text{e}}}\frac{\partial {{T}_{\text{e}}}}{\partial t}=\gamma \,{{T}_{\text{e}}}\frac{\partial {{T}_{\text{e}}}}{\partial t}=-\alpha {{T}_{\text{e}}}+\alpha {{T}_{\text{i}}}\simeq -\alpha {{T}_{\text{e}}},\end{eqnarray} \tag{ 20 }$$

with the simple solution ${{T}_{\text{e}}}(t)\simeq -\frac{\alpha}{\gamma}~t.$ We thus expect a linear temperature decrease in case the TTM is valid and γ and α are constants. This relation was utilised for the interpretation of ultrafast reflectivity measurements on laser-irradiated gold films of variable thickness in [14].

The first attempts to determine the relaxation between electrons and the crystalline lattice has been published by Kaganov, Lifshitz and Tanatarov in the context of the excitation of metals under ion impact [9]. Its relevance for ultrafast laser excitation was later pointed out by Anisimov et al in [7, 8]. Here, the energy exchange was described by equation (1). The creation of phonons by high-energy electrons was explained to correspond to the classical picture of Cherenkov radiation, emitted by a supersonic electron moving in the crystal. The integration over the probability of this process yields the electron–phonon coupling parameter α entering equation (1), for temperatures ${{T}_{\text{e}}}\gg {{T}_{\text{i}}}$ and both temperatures much greater than the Debye temperature,

$$\begin{eqnarray}&&\alpha =\frac{{{\pi}^{2}}}{6}\frac{{{m}_{\text{e}}}{{n}_{\text{e}}}c_{\text{s}}^{2}}{{{\tau}_{\text{tr}}}{{T}_{\text{e}}}}.\end{eqnarray} \tag{ 21 }$$

This expression contains the sound velocity in the material ${{c}_{\text{s}}}$ and a characteristic time ${{\tau}_{\text{tr}}}$, which is the time of free flight of electrons (occurring in transport equations, e.g. ${{\tau}_{\text{tr}}}=1/\nu$ in equation (28)) under the conditions of equal temperatures. For sufficiently low excitation strengths, in metals ${{\tau}_{\text{tr}}}\propto T_{\text{e}}^{-1}$ and thus the electron–phonon coupling parameter can be assumed as constant.

Allen [10] applies a more general expression based on Boltzmann collision integrals. His approach is valid for metals with arbitrary density of states, thus not bound to the assumption of free electrons. He connects the electron–phonon coupling parameter to the Eliashberg spectral function for electron–phonon coupling, commonly applied in the theory of superconductivity [10, 15]. The derivation of a coupling parameter near room temperature yields

$$\begin{eqnarray}&&\alpha =\frac{3~\hbar ~\lambda \langle {{\omega}^{2}}\rangle}{\pi {{k}_{\text{B}}}{{T}_{\text{e}}}}.\end{eqnarray} \tag{ 22 }$$

The moment of the spectral function, $\lambda \langle {{\omega}^{2}}\rangle$, can be obtained experimentally [15, 84]; the parameter λ is known in the theory of superconductivity as the electron–phonon coupling constant.

A large number of experimental studies have been performed to determine the electron–phonon coupling strength experimentally. However, significant deviations resulting from different measurements have been reported [85, 86]. This might be due to neglection of nonequilibrium effects as described in the next section [78]. A more simple explanation is that these measurements have been performed for different excitation strengths, yielding different values of the temperature dependent coupling parameter.

In [15], Lin et al derived a further expression, particularly considering possible temperature dependencies of the coupling parameter. The authors apply an expression based on the theory of Allen [10] and include the particular energy distribution of the electrons $f\left(\epsilon \right)$, depending on their kinetic energy . It is valid for temperatures in the range of electronVolts and reads

$$\begin{eqnarray}&&\alpha =\frac{\pi \hbar {{k}_{\text{B}}}\lambda \langle {{\omega}^{2}}\rangle}{D\left({{\epsilon}_{\text{F}}}\right)}{\int}_{-\infty}^{\infty}D{{\left(\epsilon \right)}^{2}}\left(-\frac{\partial f}{\partial \epsilon}\right)d\epsilon,\end{eqnarray} \tag{ 23 }$$

where $D\left(\epsilon \right)$ denotes the density of states and ${{\epsilon}_{\text{F}}}$ the Fermi energy. Lin et al [15] applied experimentally determined $\lambda \langle {{\omega}^{2}}\rangle$ together with density of states $D\left(\epsilon \right)$ obtained with a simulation package based on density functional theory (VASP, [87–89]). By that, the authors calculated electron–phonon coupling parameters for a wide range of electron temperatures in various metals. The resulting data base [90] is often applied when parameters for the two temperature model, equations (2) and (3), are needed.

Recently, Petrov et al [84] have shown that considerable corrections of the electron–phonon coupling parameter can be calculated when the assumption of a single electron band is dropped and, for instance, different properties of d- and s-electrons are considered. Waldecker et al [91] have studied the energy transfer to different phonon modes in aluminum and found considerable differences of the corresponding coupling strengths. Consequently, electron–phonon coupling constants obtained by fits to experimental results depend on the chosen model.

Above, expressions for the electron–phonon coupling parameter in metals have been reviewed. However, when semiconductors and dielectrics are excited, the energy is also initially absorbed by the electrons. The number of free electrons in such materials strongly changes during irradiation and depends for ultrashort pulse irradiation on a complicated interplay of changes of optical parameters and nonlinear excitation processes, see sections 2.2 and 2.3. The density of free electrons n_e is in turn a crucial parameter for the electron–phonon coupling strength, since it determines how many electrons are able to interact with the phonon subsystem at all.

Experimental measurements of the electron–phonon coupling parameter are difficult in semiconductors and dielectrics, since n_e and therefore many connected parameters are rapidly changing during and after irradiation. An indirect possibility is to measure the equilibration time between electronic and phononic temperature, ${{\tau}_{\text{relax}}}$. As equation (17) shows, this time is directly connected with α. For ${{c}_{\text{e}}}\ll {{c}_{\text{i}}}$, we obtain $\alpha ={{c}_{\text{e}}}/{{\tau}_{\text{relax}}}$ or, when carriers, i.e. electrons and holes, all share a joint heat capacity c_e−h [34], $\alpha ={{c}_{e-h}}/{{\tau}_{\text{relax}}}$. Such an expression enters the modified two temperature description ('nTTM' [24]), which will be discussed in section 4.4. Values for ${{\tau}_{\text{relax}}}$ are given between 240 fs and 500 fs for various semiconductors [34, 35, 92–95]. The electron–hole heat capacity can be calculated with the assumption of a Maxwell–Boltzmann distribution to $3{{n}_{\text{e}}}{{k}_{\text{B}}}$, thus it increases linearly with n_e. Degeneration leads to an asymptotic behavior. Assuming a constant relaxation time and applying the above approximation for the electron–phonon coupling parameter α, the dependence of α on the free electron density n_e resembles a linear increasing function with a transition to a constant asymptotic value.

However, the approximation $\alpha ={{c}_{\text{e}}}/{{\tau}_{\text{relax}}}$ actually holds only when the heat capacity is assumed as constant. Nevertheless, the resulting behavior of a coupling parameter linearly increasing with density and asymptotically resembling a constant value was confirmed in more detailed investigations [24]. In [96], calculations on the basis of Boltzmann collision integrals have been performed for SiO₂. Deviations from a constant asymptote leading to a slightly decreasing coupling parameter are attributed to screening effects becoming important for higher electron densities.

The difference in the temperatures of electrons and lattice is often referred to as 'laser-induced nonequilibrium'. However, as described in the beginning of this chapter, the electrons themselves can be in a nonequilibrium state, prohibiting description in terms of temperature. Such non-thermalised electrons may show a delayed cooling, thus a smaller electron–phonon coupling than a thermalised electron gas of the same energy density [74]. Depending on the density of states, an increased coupling strength may also be expected, in particular for noble metals where excited d-electrons contribute more strongly to the coupling than would be expected for thermalised electronic systems [78].

The coupling strength under electronic nonequilibrium may also depend on laser parameters like photon energy [97]. However, thermalisation is fast for high excitation strengths [78], thus the coupling parameter quickly tends towards its equilibrium value. A possibility to include electronic nonequilibrium effects in the frame of the TTM (2) and (3) is a time-dependent coupling parameter, as proposed in [97]:

$$\begin{eqnarray}&&\begin{array}{*{20}{l}}\alpha (t)=\alpha \left({{T}_{\text{e}}}(t),{{T}_{\text{i}}}(t)\right)\\ \qquad +\left[{{\alpha}_{\text{noneq}}}-\alpha \left({{T}_{\text{e}}}(t),{{T}_{\text{i}}}(t)\right)\right]\cdot \exp \left(-t/{{\tau}_{\text{therm}}}\right),\end{array}\end{eqnarray} \tag{ 24 }$$

where ${{\alpha}_{\text{noneq}}}$ denotes the coupling parameter in nonequilibrium and ${{\tau}_{\text{therm}}}$ the corresponding thermalisation time.

Also for strongly heated materials such as fluids or warm dense matter, strong deviations from values of the electron–ion coupling parameter predicted by the according theories have been reported [86, 98–101]. The transition from electron–phonon coupling of solids to electron–ion collisions as known in plasma theory has been proposed in [23], however, a sophisticated theory connecting these regimes is still outstanding. In case phase transitions are expected during the electron–phonon relaxation time, modified expressions for the energy transfer might become important.

Note that strong electronic excitation may lead to coherent phonon excitation or nonthermally induced phase transitions, which will be considered in sections 5.2 and 5.3. These processes occur on the timescale of phonon vibrations, i.e. about 100 fs, and are thus considerably faster than the gradual lattice heating by incoherent electron–phonon interaction discussed in the present chapter, which for most materials last for several picoseconds.

A considerable influence of the description of heat transport on surface temperature, thus on melting and ablation thresholds, was studied in [107, 108] and will be reviewed here. We assume a laser pulse duration τ_L with ${{c}_{\text{e}}}/\alpha \ll {{\tau}_{\text{L}}}\ll {{c}_{\text{i}}}/\alpha$ and neglect the heating of the lattice when focussing on the dynamics of the electrons, thus ${{T}_{\text{i}}}={{T}_{0}}$. We also neglect the thermal inertia of the electrons described by ${{c}_{\text{e}}}~\partial {{T}_{\text{e}}}/\partial t$. This is justified for low and moderate intensities approximately up to melting threshold [107, 108]. In case the spatial scale of thermal conduction is much larger than the laser penetration depth, we can consider the heating as a surface effect and take the energy source with the intensity I_las irradiating on the material surface as a boundary condition

$$\begin{eqnarray}&&{{I}_{S}}(t)=(1-R)\,{{I}_{\text{las}}}(t)=-{{\kappa}_{\text{e}}}\,\frac{\partial {{T}_{\text{e}}}}{\partial z}{{|}_{z=0}}\end{eqnarray} \tag{ 25 }$$

of an equation, following from the electron heat conduction equation (2) as

$$\begin{eqnarray}&&0=\frac{\partial}{\partial z}~{{\kappa}_{\text{e}}}~\frac{\partial {{T}_{\text{e}}}}{\partial z}~-~\alpha ~\left({{T}_{\text{e}}}-{{T}_{0}}\right).\end{eqnarray} \tag{ 26 }$$

Substituting ${{\left(\frac{\partial {{T}_{\text{e}}}}{\partial z}\right)}^{2}}=2u\left({{T}_{\text{e}}}\right)$ into the above equation, we obtain an implicit solution

$$\begin{eqnarray}&&{{I}_{S}}{{(t)}^{2}}=\,2\,{\int}_{{{T}_{0}}}^{{{T}_{\text{S}}}(t)}{{\kappa}_{\text{e}}}(T)\,\alpha (T)\,\left(T-{{T}_{0}}\right)\,\text{d}T,\end{eqnarray} \tag{ 27 }$$

determining the electron surface temperature ${{T}_{S}}(t)={{T}_{\text{e}}}(0,t)$, as a function of laser intensity I_S(t). Equation (27) is valid for any dependence of the parameters κ_e and α on temperature. Considering heat transport we find a number of different approaches for the parameter of thermal conduction of free electrons in a metal.

The electron thermal conductivity in solids can be written as

$$\begin{eqnarray}&&{{\kappa}_{\text{e}}}=\frac{1}{3}\cdot \frac{{{c}_{\text{e}}}~{{v}^{2}}}{\nu},\end{eqnarray} \tag{ 28 }$$

where v² is the electron mean square velocity and $\nu ={{\nu}_{\text{ee}}}+{{\nu}_{\text{ep}}}$. Here, v_ee is the collision frequency of electron–electron collisions, whereas v_ep denotes the electron–phonon collision frequency. The so-called transport time ${{\tau}_{\text{tr}}}=1/\nu$ is the time of free flight of electrons also entering the electron phonon coupling parameter in equation (21). For Fermi-distributed electrons we can assume ${{c}_{\text{e}}}=\gamma {{T}_{\text{e}}}$, v = v_F, ${{\nu}_{\text{ee}}}\propto T_{\text{e}}^{2}$ and ${{\nu}_{\text{ep}}}\propto {{T}_{\text{i}}}$, thus

$$\begin{eqnarray}&&{{\kappa}_{\text{e}}}\left({{T}_{\text{e}}},{{T}_{\text{i}}}\right)=\frac{1}{3}\,v_{\text{F}}^{2}\,\gamma \,\frac{{{T}_{\text{e}}}}{AT_{\text{e}}^{2}+B{{T}_{\text{i}}}},\end{eqnarray} \tag{ 29 }$$

where A and B are constants. Note that the dependence of ${{\nu}_{\text{ee}}}$ on T_e and its consequences for the thermal conductivity of metals have been intensively investigated in [84, 109].

At ${{T}_{\text{e}}},{{T}_{\text{i}}}\ll 1~\text{eV}/{{k}_{\text{B}}}$ the collision frequency is determined by electron–phonon collisions only, ${{\nu}_{\text{ee}}}\ll {{\nu}_{\text{ep}}}$, and equation (29) reduces to

$$\begin{eqnarray}&&{{\kappa}_{\text{e}}}=\frac{1}{3}~v_{\text{F}}^{2}~\gamma ~\frac{{{T}_{\text{e}}}}{B{{T}_{\text{i}}}}=:{{\kappa}_{0}}~\frac{{{T}_{\text{e}}}}{{{T}_{\text{i}}}}.\end{eqnarray} \tag{ 30 }$$

For temperatures considerably larger than the Fermi temperature, equation (29) looses its validity. More general expressions for α and ${{\kappa}_{\text{e}}}$ that are valid in a wide range of electron and lattice temperatures were derived in [110, 111]. Calculations similar to those performed in [110, 111] result in the following approximate formula for the electron thermal conductivity [107]:

$$\begin{eqnarray}&&{{\kappa}_{\text{e}}}=K\cdot \frac{{{\left({{\vartheta}^{2}}+0.16\right)}^{5/4}}~\left({{\vartheta}^{2}}+0.44\right)~\vartheta}{{{\left({{\vartheta}^{2}}+0.092\right)}^{1/2}}~\left({{\vartheta}^{2}}+b~{{\vartheta}_{i}}\right)},\end{eqnarray} \tag{ 31 }$$

where $\vartheta ={{k}_{\text{B}}}{{T}_{\text{e}}}/{{\epsilon}_{\text{F}}}$ and ${{\vartheta}_{i}}={{k}_{\text{B}}}{{T}_{\text{i}}}/{{\epsilon}_{\text{F}}}$. The parameters K and b are constants depending on the material. For low temperatures with $\vartheta \ll 1$, equation (31) reduces to equation (29). The constants K and b can then be determined with the identities $v_{\text{F}}^{2}~\gamma /3=K\cdot 0.147~A~{{\epsilon}_{\text{F}}}/{{k}_{\text{B}}}$ and $b=B~{{k}_{\text{B}}}/A~{{\epsilon}_{\text{F}}}$. At high temperatures, when $\vartheta \gg 1$ and the electron gas becomes non-degenerate, equation (31) results in the well known dependence ${{\kappa}_{\text{e}}}\sim T_{\text{e}}^{5/2}$ which is characteristic for low-density plasmas [22, 23].

Figure 7 shows the dependence on electron temperature of these different expressions of the electron thermal conductivity for the case of gold. The lattice temperature was assumed as ${{T}_{\text{i}}}={{T}_{0}}=300~\text{K}$. Other parameters are the Fermi velocity ${{v}_{\text{F}}}=1.39\cdot {{10}^{6}}~\text{m}~{{\text{s}}^{-1}}$, the heat capacity proportional to $\gamma =67.6\,\text{J}\,{{\text{m}}^{-3}}\,\text{K}{{\,}^{2}}$, and a constant value of ${{\kappa}_{0}}=315\,\text{W}\,{{\text{m}}^{-1}}{{\text{K}}^{-1}}$ (all from [112]). The constants A and B determining the collision frequencies entering equation (29) are $A=1.2\cdot {{10}^{7}}~{{\text{s}}^{-1}}{{\text{K}}^{-2}}$ and $B=1.23\cdot {{10}^{11}}~{{\text{s}}^{-1}}{{\text{K}}^{-1}}$, taken from [113]. The constants in equation (31) then read K = 353 W m⁻¹K⁻¹ and b = 0.16.

Zoom In Zoom Out Reset image size

One easily sees the different ranges of validity for the different approximations: equation (30) strongly overestimates the heat conduction for temperatures larger than ${{T}_{\text{e}}}\gtrsim 4\cdot {{T}_{0}}=1200~\text{K}$. At temperatures around Fermi temperature ${{T}_{\text{e}}}\simeq {{T}_{\text{Fermi}}}\simeq 64~000\text{K}\simeq 200~{{T}_{0}}$, equation (29) also ceases to be correct and leads to a strong underestimation of ${{\kappa}_{\text{e}}}$.

The temporal evolution of the surface temperature strongly depends on the choice of the parameter for thermal heat conduction, ${{\kappa}_{\text{e}}}$. One can see from equation (26) that in its solutions the electron temperature follows the shape of the laser pulse instant by instant. Assuming a constant coefficient of electron thermal conduction, the electron temperature at the surface ${{T}_{\text{e}}}\left(z=0,t\right)={{T}_{S}}(t)$ is proportional to the absorbed intensity ${{T}_{S}}(t)\propto {{I}_{S}}(t)$; assuming expression (30) with ${{\kappa}_{\text{e}}}\propto {{T}_{\text{e}}}$ as valid for low excitations, we obtain ${{T}_{S}}(t)\propto {{I}_{\text{abs}}}{{(t)}^{2/3}}$ [108]. Generally, the solution (27) of equation (26) can be applied for arbitrary expressions of ${{\kappa}_{\text{e}}}$ and α.

Figure 8 shows the time history of the electron surface temperature ${{T}_{S}}(t)$ for gold, irradiated with a Gaussian laser pulse of 1 ps duration (full width at half maximum, FWHM) with a maximum surface intensity of ${{I}_{S,\text{max}}}=4\cdot {{10}^{10}}\,\text{W}\,\text{c}{{\text{m}}^{-2}}$, corresponding to a fluence of ${{F}_{S}}=35\,\text{mJ}\,\text{c}{{\text{m}}^{-2}}$. The initial temperature was chosen equal to the constant lattice temperature of ${{T}_{0}}=300~\text{K}$. The additionally shown increase of lattice temperature has been calculated through equation (1), with a constant electron–lattice coupling parameter of $\alpha =3.5\cdot {{10}^{16}}\,\text{W}\,\text{K}{{\text{m}}^{-3}}$, after the surface temperature of electrons ${{T}_{S}}(t)$ has been determined for ${{T}_{\text{i}}}={{T}_{0}}$. The figure shows the heating of the electrons, following the Gaussian shape of the laser pulse. The solid line was calculated with the general expression for the thermal heat conduction, equation (31), valid in the shown temperature regime. The dashed–dotted curve, calculated with equation (29), leads to an overestimation of the surface temperature since the thermal conduction approaches zero for high temperatures, compare figure 7. Equation (30), instead, underestimates the surface temperature of electrons, since the thermal conduction is far too large even for moderate electron temperatures. For laser intensities under consideration, the time dependence of the surface temperature ${{T}_{S}}(t)$ calculated with a constant ${{\kappa}_{\text{e}}}$ is close to that obtained with the full expression, equation (31). To understand this result, we note that, though the thermal conduction coefficient determined by equation (31) varies several-fold in the temperature range from 300 K to 20 000 K, it has in this range two extrema, maximum and minimum. As the result, ${{T}_{S}}(t)$, defined by equation (27), appeared to be only slightly dependent on the particular form of ${{\kappa}_{\text{e}}}\left({{T}_{\text{e}}}\right)$. Also the subsequently calculated maximum surface temperature of the lattice, reaching slightly above melting temperature when the full expression (31) for the electronic heat conductivity is applied, is reproduced best with a constant heat conductivity. Since the approximation of a constant ${{\kappa}_{\text{e}}}$ provides a good accuracy and simplifies the analysis, it can be employed to obtain analytical formulas to estimate melting and evaporation thresholds.

Zoom In Zoom Out Reset image size

Note that the simplifications made above, i.e. neglection of electronic heat capacity and cooling down of electrons to initial temperature, are not expected to have significant influence on the eventual lattice temperature. The reason is that in the full solution of equations (2) and (3) (see next section) the equilibrium temperature of electrons and lattice is mainly determined by the heat capacity of the lattice.

The TTM given by the equations (2) and (3) can be readily solved numerically to analyse the heat penetration into the bulk of the material. Figures 9 and 10 show the results of such a run with parameters close to the ablation threshold. Parameters for gold have been applied, as given in the previous section, with a heat conductivity according to equation (31). Figure 9 presents the transient surface temperatures of electrons and lattice, which values are comparable with the analytical approximate solution in figure 8, applying the same heat conductivity given in equation (31). The temperature profiles of the lattice at different instants of time are shown in figure 10. We see that after a few picoseconds the lattice temperature strongly exceeds the melting temperature, also in the bulk of the material up to several tens of nanometers.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

A contour plot of such a calculation for the case of copper is shown in figure 11 for an absorbed fluence just above ablation threshold [32]. The upper panel shows the electron temperature, whereas the lower panel shows the lattice temperature. The black line in the plot of the electron temperature represents the laser intensity as a function of time. The black line in the plot of the lattice temperature is drawn at an 'ablation temperature' defined and calculated in [32]. Here, the sum of both energy densities, of lattice as well as of electrons, was taken for equilibrated temperatures and compared to the latent heat of evaporation. The resulting evaporation temperature of copper was 9100 K. Note that the latent heat of phase transitions was taken as a criterion here but was not explicitly included in the calculations underlying figure 11.

Zoom In Zoom Out Reset image size

A numerical result of the TTM, equations (2) and (3), including an effect of phase transition on resulting temperatures by considering the latent heat of melting is shown in figure 12. Here, after reaching the equilibrium melting temperature (dashed line), the further energy increase was attributed to the melting process until the latent heat of melting was provided. Only after that, further energy transfer from the electrons further increases the temperature of the lattice. Note that the combination with molecular dynamic (MD) simulations leads to similar temperature values [114]. The abovementioned consideration of latent heat of melting does not explicitly account for melting kinetics, which are inherently included in MD simulations, see sections 7 and 8.

Zoom In Zoom Out Reset image size

Heat transport, as well as electron–phonon coupling, may depend strongly on the transient temperatures of electrons and lattice. In [115] a systematic study on the solution of the TTM was performed for varying pulse durations. It was shown that the maximum lattice temperatures are reached for pulse durations close to the electron–phonon relaxation time (17).

With the help of hybrid TTM-MD simulations (see section 8), the quality and efficiency of ablation has been studied additionally [116]. Here, the functional dependencies of heat conduction and electron–phonon coupling parameter on electron temperature have been particularly important. It was predicted that their interplay determines the character of the processing and therefore influences the quality of the final morphological structure.

In the case of laser absorption in dielectrics or semiconductors, the number density of free carriers, excited to the conduction band of the solid (compare sections 2.2 and 2.3), is spatially not constant. This is due to a large but finite penetration depth of the laser in the solid. The excitation of the carriers itself determines the weakening of the laser intensity in space.

As a result, density gradients of free carriers occur in the solid. In [34], an extension of the TTM (2) and (3) was proposed, which takes into account not only Fourier's law describing the heat conduction due to temperature gradients but also Ohm's law (electrical current due to electrical fields, i.e. gradient in charge density) as well as the cross-over processes of the Seebeck's and Peltier's effect. Electrons and holes have been assumed to propagate together, a process known as ambipolar diffusion which leads to equal densities of both kinds of carriers all over the solid. This density-dependent two temperature model (nTTM, [24]) allows to study excitation and transport of electrons and holes as well as carrier-lattice coupling for semiconductors on equal footing as the classical TTM for metals. The numerical implementation is, however, more complex. After the proposal by van Driel in 1987 [34], the model was recalled by Chen et al in 2005 [35]. In recent years, further implementations have been realised [21, 24, 117–119], some of them in a hybrid description with molecular dynamic methods as will be mentioned in section 8.

Another thorough discussion of transport in all kinds of materials (metals, semiconductors and dielectrics) has been given by Bulgakova et al in [120]. The two-temperature transport description for metals was completed by the rate of free electron generation, which was implemented in the description for semiconductors and dielectrics. In contrast to [34], the assumption of ambipolar diffusion was abandoned here. The resulting electrical field was calculated self-consistently, balancing the currents due to both the field and density gradients. Consequences of that work will be discussed in section 5.3.

Usually melting starts at the surface, where the energy barrier for heterogeneous nucleation of a liquid layer at the solid–vapor interface vanishes. Then, a melting front proceeds from the surface into the material with a velocity ultimately limited by the speed of sound. Experimentally, melt front velocities of many hundreds of m s⁻¹ were observed [121, 122]. The lower limit of the melting time is the time an acoustic wave travels through the heated layer. Typical estimated times for melting by this mechanism are in the range of  ≈$100~\text{nm}/1000~\text{m}{{\text{s}}^{-1}}=100~\text{ps}$.

In [123] the possibility of laser-induced melting of crystals due to homogeneous nucleation was considered. Figure 13 shows that a sufficiently superheated bulk crystal (about 1.5 times the melting temperature) would melt completely in less than one picosecond assuming homogeneous nucleation of the molten phase. Here, the statistical rate of the appearance of critical nuclei was exploited to calculate the time after which the complete volume is filled with critical nuclei.

Zoom In Zoom Out Reset image size

The resulting melting time strongly depends on the ratio $T/{{T}_{\text{m}}}$ of overheating and on the applied material parameters and can therefore not be calculated exactly. However, the reverse statement is reliable: melting times less than a picosecond are easily reached. The time for homogeneous melting is thus limited by the time for lattice heating. In the case of high power laser irradiation with sufficient intensity to initiate ablation, melting temperature as well as strong overheating can be reached on a timescale of a few picoseconds, see section 4. Thus, for sufficient overheating the total time of lattice melting is governed by the time of electron–lattice heating, typically in the range of a few picoseconds. Such melting times of metals have been confirmed experimentally [124, 125].

An atomic picture of nucleation kinetics can been obtained with molecular dynamic (MD) approaches, which will be described in sections 7 and 8. In [114, 126] a combined simulation of two-temperature heat conduction (2) and MD approach giving access to the transient atomic structure has been performed. Figure 14 shows snapshots of such simulation: The atoms are colored according to their local order parameter, red atoms have crystalline surrounding while blue regions are molten. Irradiation of a 50 nm Nickel film was assumed with a laser of 200 fs duration and a fluence considerably above melting threshold but below ablation threshold. Figure 14 represents locations close to the irradiated surface where the temperature is clearly above equilibrium melting temperature. Note that the MD simulation implicitly accounts for nonequilibrium melting temperatures. In particular the pressure dependence of melting dynamics has been studied in [126]. The snapshots start 14 ps after irradiation with a 200 fs laser pulse, each panel is taken one picosecond after the preceding one. Small nuclei of molten phase can be observed in the volume after 15 and 16 ps, which grow fast and appear as a large molten region within the next picosecond.

Zoom In Zoom Out Reset image size

In contrast, figure 15 shows a cutout at the backside of the film where no large overheating was observed. These snapshots are separated each by 20 ps. A melt front proceeds gradually from both sides (at the left-hand side the fast melting process discussed above has reached this depths, at the right-hand side a melt front has built at the free surface). The melt front velocity is much slower than the effective velocity of the ultrafast process discussed above.

Zoom In Zoom Out Reset image size

Thus, both types of thermal melting have been observed in one simulation. Since heterogeneous nucleation of the molten phase appears without any energy barrier at the solid-vapor interface, it occurs on both surfaces. For strong overheating above melting temperature, the gradually proceeding melt front is superimposed by molten nuclei which build homogeneously due to fluctuations. Once a phase transition has started, the melt front further proceeds into regions of lower overheating with temperatures around melting temperature. Recently, both thermal melting processes have also been observed in a simulation of a laser-irradiated semiconductor film [119].

When a material is excited by high-fluence femtosecond laser pulses, an extreme nonequilibrium state is created, which can live for a couple of picoseconds, see section 1. This state is characterised by the presence of very hot electrons that, after absorbing a large amount of photons from the laser pulse and thermalizing through electron–electron collisions, acquire an extremely high temperature. On the other hand, the lattice still remains cold, since the laser excitation and electron thermalisation occur on a time scale that is generally much shorter than lattice heating through incoherent electron–phonon coupling, which has been described in section 3. Although ions remain at room temperature, they experience strong forces due to the dramatic changes suffered by the potential energy surface (PES), i.e. the effective potential acting on the ions, which depends on the electronic state, after the laser pulse excited a considerable fraction of electron–hole pairs [127, 128]. Thus, ions 'feel' a sudden change of the potential, which leads to the appearance of forces, and start moving. This motion is in most cases cooperative (i.e. the ions do not move independently from each other) and can sometimes even be coherent [129, 130]. Such mechanisms have been reviewed in [131–133].

If the forces arising from the above mentioned changes in the potential energy surface are strong enough, then the material disorders very quickly, even before incoherent electron–phonon heating becomes effective. This ultrafast disordering can be interpreted as melting, which occurs at room temperature. Such phenomenon is called nonthermal melting, and has been studied both experimentally and theoretically since the early eighties [134, 135]. In case sufficient energy is provided by the exciting laser pulse, the material stays in the new phase also after recombination of the electron–hole pairs and/or electron cooling.

The original theory formulated by Stampfli and Bennemann [127, 128] states that the functional form of the potential energy surface (PES) strongly depends on the occupation of the electronic levels. As an example, Stampfli and Bennemann computed the PES along a particular transversal and a particular longitudinal phonon mode at the L-point of the Brillouin zone of silicon. When the electrons are in the ground state, the PES looks like a two dimensional harmonic potential, in which both phonon modes oscillate independently of each other. Now, if 15$\%$ of the valence electrons become excited, the PES acquires a form which can no longer be described as a two dimensional harmonic potential. Moreover, in one of the phonon directions it even becomes repulsive (negative curvature), which means that large atomic displacements have to be expected along this direction [128]. Stampfli and Bennemann used for their theory the tight-binding model, which relies on a simple single-electron Hamiltonian. Later, different density functional theory approaches were used to account for the changes in the PES upon electron–hole-pair excitation, which were based on the ideas first developed by Stampfli and Bennemann [136–143]. The basic procedure to obtain the so called 'laser excited ab initio PES' is first to assume that electrons have a particular (usually very high) temperature as a consequence of the excitation and subsequent thermalisation. Thus, electronic occupations are considered to follow a Fermi distribution at the resulting electronic temperature T_e. Then, the Kohn–Sham equations are solved self-consistently under this constraint for the electronic occupations. In addition, the electronic entropy at T_e is determined. Then, by applying a generalised Born–Oppenheimer approximation, the PES is obtained as $V\left(\left\{{{\mathbf{r}}_{i}}\right\},\,\left\{{{n}_{\mathbf{k}}}\left({{T}_{\text{e}}}\right)\right\}\right)={{E}_{e}}-{{T}_{\text{e}}}{{S}_{\text{e}}}$, where $\left\{{{\mathbf{r}}_{i}}\right\}$ are the atomic coordinates, $\left\{{{n}_{\mathbf{k}}}\left({{T}_{\text{e}}}\right)\right\}$ the electronic occupations at T_e, E_e is the total energy of electrons (which depends on $\left\{{{n}_{\mathbf{k}}}\left({{T}_{\text{e}}}\right)\right\}$), and S_e is the electronic entropy, also depending on the electron occupations. ab initio MD simulations were applied to study laser induced structural changes, in semiconductors like Si [136, 137, 139, 141, 142], Te [144], InSb [140, 145], As [146], semimetals like Bi [138, 147], and metals like Al [139], Au [139], Ag [143, 148] and Cu [143].

In semiconductors, laser pulses create electrons in the conduction band (mostly in antibonding states) and holes in the valence band (mostly in bonding states). As a consequence, the overall bond strength of the material is reduced. This phenomenon is known as laser induced bond softening and has been confirmed by a wide range of experimental studies [149, 150]. While the first experiments were based on optical methods [134, 135], later time resolved x-ray diffraction or x-ray spectroscopy measurements allowed a direct analysis of the lattice order [130, 151–156]. All these results show clear evidence of ultrafast structural changes in different elements, mostly semiconductors and semimetals. Experiments have also revealed rather short melting times of metals after femtosecond laser excitation [124, 125]. However, they could be explained with (homogeneous) thermal melting [123] or other ultrafast destruction mechanisms [157].

The main feature of non thermal melting is the fact that it occurs when atoms are still either at room temperature or at least well below the melting temperature. In the framework of the theory of Stampfli and Bennemann, nonthermal melting should occur when the presence of excited carriers leads to lattice instabilities, represented by imaginary phonon frequencies. All static ab initio calculations agree in the occurrence of lattice instabilities for high enough electronic temperatures. Note that static calculations, based on the calculation of phonon spectra in the laser excited state, cannot predict which kind of lattice instabilities will lead to nonthermal melting. In semiconductors, the first phonons to become unstable are the transversal acoustic ones.

In order to confirm that lattice instabilities lead to nonthermal melting, simulations had to be performed. The first ab initio simulations to describe nonthermal melting of semiconductors were performed by Parrinello and co-workers [136, 137]. Those simulations on Si at constant volume and using a basis set of plane waves yielded a rapid disordering of the lattice structure and convergence of the pair correlation function to that of the liquid state.

It is important to mention that it is experimentally very difficult to characterise nonthermal melting and to distinguish it from ultrafast thermal melting. Therefore, experimental efforts in recent years have been aimed at obtaining information about the trajectories of the atoms immediately after laser excitation. Ultrafast diffraction techniques, based either on x-rays or on electrons, are the most suitable to follow the average atomic trajectories. However, different experiments on semiconductors have led to different and contradictory results. For instance, laser pump x-ray probe diffraction experiments on InSb were interpreted in terms of an inertial model [158]. According to it, the laser excitation leads to a completely flat PES, so that atoms move after excitation with constant thermal velocity corresponding to room temperature. On the other hand, ultrafast electron crystallography experiments performed on Si were compatible with a diffusive motion of the atoms after laser excitation, which is governed by stochastic interatomic collisions. In both experimental studies the time resolution was not small enough to capture possible acceleration of atoms due to lattice instabilities. Recent ab initio simulations performed on larger cells and using the code CHIVES [141, 142], based on Gaussian wave functions, clarified this controversy. According to those simulations, atoms first undergo an acceleration and are then decelerated for a time between half and one picosecond. This subdiffusive motion could be characterised by being fractionally diffusive [142]. Then, after enough interatomic collisions have taken place, the motion becomes purely diffusive, satisfying Einstein's diffusion equation. The onset of diffusion corresponds to the presence of a liquid state [142].

Density functional theory calculations have shown that, in contrast to semiconductors with covalent bonding, ultrafast laser excitation induces bond hardening in metals like gold, copper and silver [139, 143, 148]. However, bond softening or no effects on phonon frequencies have also been obtained for various laser-excited metals [139, 159]. In [160] it was mentioned that the question on bond hardening or softening may depend also on boundary conditions.

Apart from nonthermal melting discussed in the previous section, another ultrafast phase transition is under discussion, namely Coulomb explosion. This process is known for large molecules or clusters, when highly excited electrons leave the structure. The resulting repulsive forces between the remaining ionised atoms lead directly to destruction of the cluster.

It has been intensively debated whether this process may also be responsible for ablation from large solid samples. Molecular dynamic approaches were used to study the effect of such repulsive Coulomb forces on the crystal lattice [161]. Experimental time-of-flight measurements have found monoenergetic ion beams emitted from irradiated surfaces, which were interpreted as stemming from Coulomb explosion [162]. It is however questionable whether in solids strong electric fields at the surface can be compensated by conduction band electrons from the bulk of the material. The large reservoir of mobile electrons is the main difference from the case of small solid-like clusters.

Bulgakova et al [120] have studied the possibility of electrostatic disintegration of surface layers in detail, by modelling photoelectron emission and electronic transport due to the resulting fields for metals, semiconductors and dielectrics. Typical representations of these material classes have been chosen with gold, silicon and sapphire. It has been assumed that the surface of a bulk material was irradiated with a 100 fs (FWHM) laser pulse. The fluence has been chosen slightly above the respective threshold for ion emission, namely 1.2 J cm⁻² for gold, 0.8 J cm⁻² for silicon and 4 J cm⁻² for sapphire. The resulting electric field in a thin surface layer is plotted over time in figure 16. The calculated field strength for gold and silicon are enlarged by factors of 100 and 50 respectively. However, for sapphire rather large transient electric fields are observed. In [120] a critical field strength necessary to break atomic bonds in sapphire was estimated with help of the latent heat of sublimation for a single atom. The resulting critical electric field is shown in figure 16 as a dashed line. For a short time the calculated electric field at the surface of the dielectric material exceeds this value. Bulgakova et al conclude that Coulomb explosion is possible in dielectrics. Note that a similar amount of emitted electrons have also been obtained in the calculations for metals and semiconductors. However, the resulting space charge was neutralised fast enough, due to the high amount of high-mobility charge carriers in these materials [120]. In [163] experimental results were interpreted with the process of Coulomb explosion also for a semiconductor, namely silicon, which have later been under discussion in [164, 165].

Zoom In Zoom Out Reset image size

The theoretical results of [120] have been compared with experiments on sapphire [166] concluding that the sharp peak of electrical field at the surface (see figure 16) may be responsible for the observed fast ion emission from the surface marking the onset of ablation. However, the main part of ablation occurs later and is accompanied by a much larger amount of removed particles and the formation of a crater on the target surface [166–168]. The latter ablation regime, also occurring in metals and semiconductors, is the main mechanism determining ultrafast laser ablation. The dynamics of matter expansion will be discussed in the following sections of this review.

The initial rapid laser excitation can be viewed as isochoric heating. Equilibration of electron–lattice temperature is generally faster than expansion of the material, which, in turn, is usually fast in comparison with the time scale of heat conduction. Therefore, we consider here adiabatic expansion of matter with a density equal to normal solid density ${{\rho}_{0}}$. We further assume a temperature on the order of the critical temperature of the material and study the expansion of a uniformly heated, semi-infinite layer.

The flow of expanding matter is described by the equations of gas dynamics

$$\begin{eqnarray}\begin{array}{*{35}{l}} \frac{\partial \rho}{\partial t}+\frac{\partial}{\partial x}~\left(\rho u\right) & = & 0, \end{array}\end{eqnarray} \tag{ 32 }$$

$$\begin{eqnarray}\begin{array}{*{35}{l}} \frac{\partial u}{\partial t}+u~\frac{\partial}{\partial x}u & = & -\frac{1}{\rho}~\frac{\partial p}{\partial x}, \end{array}\end{eqnarray} \tag{ 33 }$$

where ρ is the density, and p is the pressure of the material which expands in the x direction with the mass velocity u. The set of equations (32) and (33) is completed in the case of adiabatic expansion by the equation of the isentrope $S\left(\,p,\rho \right)=\text{const}.$, resulting in a direct dependence of pressure on density, $p={{p}_{S}}\left(\rho \right)$. The flow of matter then depends only on the self-similar coordinate x/t [173, 179] and we obtain an implicit equation for the density profile as

$$\begin{eqnarray}&&\frac{x}{t}\left(\rho \right)=-{{c}_{\text{s}}}\left(\rho \right)+{\int}_{\rho}^{{{\rho}_{0}}}{{c}_{\text{s}}}\left({{\rho}^{\prime}}\right)\frac{\text{d}{{\rho}^{\prime}}}{{{\rho}^{\prime}}}.\end{eqnarray} \tag{ 34 }$$

Here, ${{c}_{\text{s}}}\left(\rho \right)=\sqrt{{{\left(\partial p/\partial \rho \right)}_{S}}}$ is the local sound velocity, which is the key parameter in the following discussion.

Figure 17 shows the phase diagram of aluminum with four isentropes of different entropy and the isochore of solid state density ${{\rho}_{0}}$. They were taken from the tabulated equation-of-state (EOS) data of aluminum which were constructed by the methods developed in [180, 181]. After the initial rapid isochoric heating (vertical dotted line) by laser excitation, the material is in a state marked by the letter A. Varying excitation strengths lead to different temperatures and pressures within the material, therefore the entropy S at that stage and during the following adiabatic expansion depends on the excitation strength. Expansion occurs along an isentrope, four examples are marked in the figure 17. The solution of equation (34) for these four isentropes is plotted as a density profile depending on the self-similar coordinate x/t in figure 18.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The square root of the slope of the isentrope is the sound velocity which enters equation (34). Outside the two-phase region the sound velocity of the liquid phase has a comparably large, approximately constant value (note the logarithmic scale of the pressure axis in figure 17). For a constant sound velocity, the self-similar solution (34) yields an exponential expansion profile, this can be compared with the head and the initial part of the rarefaction waves in figure 18. However, the sound velocity decreases sharply when the isentrope enters the two-phase-region [179]. For the isentrope with the smallest entropy, this point at the binodal is marked with the letter B in figure 17. The strong drop of the sound velocity by several orders of magnitude causes a jump in $x/t~\left(\rho \right)$ in equation (34), since the term ${{c}_{\text{s}}}\left(\rho \right)$ changes rapidly at the moment the density ρ drops below the density at the binodal, ${{\rho}_{\text{B}}}$. Therefore, a plateau of constant density $\rho (x/t)={{\rho}_{\text{B}}}$ appears in the self-similar solution, compare equation (34) and figure 18. More importantly, due to the low value of sound velocity for densities smaller than ${{\rho}_{\text{B}}}$, the self-similar coordinate x/t in equation (34) stays nearly constant with further decreasing ρ. For the $\rho (x/t)$ profile in figure 18, this leads to a steep, step-like density gradient. Note that for higher excitations, a layer of very low density (corresponding to the gas phase) can be observed in front of this step-like density gradient. For the lower entropies considered here, corresponding to lower initial temperatures, the density decreases on a scale of $\lesssim$10 ms⁻¹. Such a profile is, also after several nanoseconds, steep enough to act as an optically sharp interface.

Here, we have discussed the self-similar profile for the special case of expanding aluminum. In [173, 174] a general qualitative case was discussed, considering a simplified isentrope with two different slopes. The assumed analytical expression for the pressure dependence on such an isentrope was $p\propto {{\rho}^{\gamma}}$ with $\gamma ={{\gamma}_{\text{liq}}}\gg 1$ for the case of expansion in the liquid phase and $\gamma ={{\gamma}_{2\text{ph}}}\ll 1$ for the two phase region, reflecting its small sound velocity. Both branches are matched at a certain point (${{\rho}_{\text{B}}},{{p}_{\text{B}}}$) in the phase diagram, namely on the binodal, compare point B in figure 17. This consideration directly shows that the point (${{\rho}_{\text{B}}},{{p}_{\text{B}}}$) corresponds to a finite segment in the self-similar representation, building the plateau segment. Its width is given by the difference of the sound velocities ${{c}_{\text{liq}}}\left({{\rho}_{\text{B}}}\right)-{{c}_{2\text{ph}}}\left({{\rho}_{\text{B}}}\right)\approx {{c}_{\text{liq}}}\left({{\rho}_{\text{B}}}\right)$ taken at both respective sides of the binodal point B, while its density equals ${{\rho}_{\text{B}}}$ [1, 173, 174].

Note that another mechanism to explain an optically sharp interface was discussed by Bulgakova in [182]. Here, heating and expansion far above the critical point was assumed and a certain curvature of the isentrope led to a rarefaction shock appearing as a sharp front in the density profile. However, as will be discussed in the following section, an optically sharp interface seems to exist in particular close to ablation threshold, thus for lower excitations than assumed in [182].

Up to now, the expansion of a semi-infinite, uniformly heated material into vacuum has been analyzed. Such condition of a uniform temperature can be fulfilled in a thin film of heated material with thickness $d\sim 100$ nm [14], compare discussion in section 4. However, when the rarefaction wave hits the unperturbed substrate, i.e. at a time of ${{t}_{r}}=d/{{c}_{\text{s}}}\left({{\rho}_{0}}\right)$, it is reflected and the self-similar solution loses its validity. Reflection of the rarefaction wave will also occur in bulk material at the non-ablating, shock-compressed internal part of the target.

Figure 19 sketches qualitatively different stages of the evolution of the density profile. Figure 19(a) shows the unperturbed self-similar rarefaction wave with a head travelling with sound velocity ${{c}_{0}}={{c}_{\text{s}}}\left({{\rho}_{0}}\right)$ towards the substrate, as shown also in the density profiles in figure 18. After reflection at the substrate, the head of the rarefaction wave travels with the local sound velocity $c\left(\rho \right)$ in the same direction as the expanding fluid, see figure 19(b). However, when reaching the plateau, the local sound velocity is $c\left(\rho \right)={{c}_{\text{liq}}}\left({{\rho}_{\text{B}}}\right)$ while at the other end of the plateau it is given by the much smaller sound velocity on the two-phase side of the corresponding point on the binodal ${{c}_{2\text{ph}}}\left({{\rho}_{\text{B}}}\right)$. Therefore, the rarefaction wave will not overtake the plateau and the latter remains nearly unperturbed when moving away from the substrate. It follows directly from mass conservation that the density in the region between the substrate and the plateau must decrease. The resulting density profile at that stage, expected to survive for comparably long time, is sketched in figure 19(c). Numerical gas dynamic calculations [173] as well as molecular dynamic simulations [178, 183] confirm this picture. The latter simulations will be discussed in more detail in section 7.

Zoom In Zoom Out Reset image size

The low density region in figure 19(c) consists of two-phase matter, while the high density region is in condensed phase. With increasing fluence the thickness of the two-phase matter increases, while the thickness of the condensed shell decreases. Since the intensity distribution over the focal spot is not uniform but in good approximation follows a Gaussian distribution, the expanding plume takes the shape of a cupola [175–178]. A laser beam of visible light can pass through the shell building this cupola. It can also pass through the low-density region between shell and substrate. This picture therefore represents an explanation for the observation of Newton fringes during ultrafast microscopy of an ablating surface [6].

For high laser fluences, the expansion isentrope enters the two-phase region near the critical point. Here, the drop in sound velocity is less pronounced and the resulting density profile is less steep, compare figures 17 and 18. Therefore, the outer condensed shell of the expansion plume will smear out as time proceeds; first in the center of the Gaussian laser profile, where the largest intensity has excited the target. Figure 20 shows a schematic picture of the fluence profile of the irradiating laser beam (a) and the resulting profile of the expansion plume (b). The point 'ev' in figure 20 marks the rim of the open aperture of the cupola, where the fluence exceeds a critical value F_ev for the appearance of a shell of condensed matter. Note however, that a comparably sharp density gradient still exists along the curve '1'. The region in front of the sharp density gradient is filled with vapor 'v'; its thickness increases with increasing fluence [175].

Zoom In Zoom Out Reset image size

In [184] the limits for the observation of Newton fringes were estimated. The lower limit was identified with an isentrope reaching the binodal at zero or normal pressure, marking roughly the ablation threshold. The upper limit was defined through the sharpness of the density profile of the self-similar solution. It turned out that the ratio of the starting temperatures (at solid density, point A in figure 17) approximately equals the ratio of critical temperature to the boiling temperature at normal conditions, ${{T}_{\text{crit}}}/{{T}_{\text{boil},\text{n}}}$, which is the relative extension of the two-phase region.

In case of linear absorption two pulses of different fluences heat the target to different temperatures of approximately the same ratio. In order to verify whether our estimations on the thresholds for the observation of the Newton fringes are reasonable, we therefore compare the ratio of the experimentally observed fluence limits of appearance of Newton fringes with the temperature ratio ${{T}_{\text{crit}}}/{{T}_{\text{boil},\text{n}}}$ extracted from EOS data. Table 1 shows these values for various materials. In case the disappearing of Newton fringes was not observed in the experiment, the largest ratio experimentally reached is given.

It can be seen from table 1 that the fluence range for the observation of Newton fringes corresponds to the relative size of the two-phase region. The denoted ratios do not exactly coincide; however, the comparison is satisfactory and supports our explanations for the observation of Newton fringes in near-threshold ablation experiments of different materials. We conclude that the upper limit for the observation of Newton fringes is given by the sharpness of the expansion profile, while the ablation threshold is connected with an expansion isentrope crossing the binodal at approximately zero pressure. The latter assumption has also been verified through molecular dynamic simulations, see section 7.

Note that this picture of ablation close to threshold is universal not only for different kinds of material, metals as well as semiconductors [184, 185], but also independent of the wavelength of excitation, up to the regime of extreme ultraviolet (XUV) [186].

There exist a number of hydrodynamic approaches studying distinct features of ultrafast laser ablation, like progression of a pressure wave into the target [178], spinodal decomposition [187], influence of a critical tension on spallation [188], formation of nanoparticles in the ablation plume [189], and many more.

There are also several hydrodynamic simulations aiming to follow the complete ablation process. They incorporate descriptions of the initial excitation and energy dissipation to hydrodynamic approaches in order to include all relevant mechanisms [190–192].

Among them, the work of Colombier et al [190] is particularly interesting, since the influence of certain parameters entering the simulation is systematically studied and the numerical results are directly compared with experimental measurements. In [190] the authors implement the two-temperature model, equations (2) and (3), into an earlier developed hydrocode [193]. They apply a constant coupling parameter α and a thermal conductivity ${{\kappa}_{\text{e}}}$ as given in equation (31). Moreover, they allow changing heat capacities due to changing densities of the solid. This feature is necessary because the authors assume that electron–lattice heating may proceed more slowly than expansion of matter. Thus, heating in this case is not assumed to be necessarily isochoric. In accordance with this assumption, an additional term for the electronic pressure was included in the hydrodynamic equations. The authors compare the calculated ablation depths for copper and aluminum with experimentally measured ablation depths and thresholds. The coincidence is good for fluences well above ablation threshold. However, the simulation overestimates the threshold fluence for copper, while it is underestimated in the case of aluminum. For the case of copper, a parametric variation was performed to analyse the influence of certain parameters entering the model. Namely, the electron–phonon coupling parameter was multiplied by a factor between 0.03 and 30 keeping all other model parameters at the original values. The same was done for the electron thermal conductivity. The results are repeated in figure 21, which shows the threshold fluence for ablation for both cases of parameter variations.

Zoom In Zoom Out Reset image size

An increased thermal conductivity leads to an increased threshold fluence, while an increased electron–phonon coupling decreases the predicted ablation threshold. The physical explanation is straightforward: a better energy confinement and a faster lattice heating leads to higher lattice temperatures at the surface and thus more efficient ablation.

This analysis gives an idea about the accuracy of calculated ablation thresholds. In each simulation, parameters of a distinct uncertainty enter. When the simulation is aimed to provide quantitative predictions of experimentally accessible observables, the influence of such parameters on the finally sought integrated quantity has to be analyzed carefully. Such variations are also a key to understanding the interplay of mutually influencing processes.

The model described in [190] was extended in [191], where the time evolution of electronic pressure and atomic pressure was included in the hydrodynamic model explicitly. It was shown that the electronic contribution to ablation is particularly strong in the first tens of nanometers close to the free surface. The same idea was followed by Chimier et al in [194]. Here, a scenario of instantaneous pressure relaxation was discussed, thus, expansion occurs not only before electron–lattice coupling is completed, but even during ultrafast laser irradiation. The resulting equation-of-state in dependence on electron temperature and lattice temperature was presented. While this EOS reveals a number of generally interesting features, conclusions based on the assumption of instantaneous pressure relaxation on subpicosecond timescales are questionable.

In [192], Povarnitsyn et al simulated ablation of aluminum. They assumed two temperatures in one fluid and connected the TTM with a semiempirical multiphase EOS with separate descriptions of electrons and heavy particles. Expansion and phase transitions far above ablation threshold were studied. For different layers in the expanding material, the thermodynamic trajectories in the phase diagram were followed and material's behavior was analysed. Two different mechanisms of liquid–gas phase transitions were studied, first, bubble growth in metastable liquids and second, fracture of material due to large strain rates. The latter has been shown to be the dominant ablation mechanism. Thus, processes at negative pressure have to be included to gain reasonable comparison with experimental results.

In [195] and [196] the model was further elaborated and applied for various materials. The method allows in particular analysis of the phases of the expanding matter. Figure 22 shows the time-versus-depth plot of expansion. Different colors refer to different phases. White areas are not occupied by any ablating matter. The ablation/movement of gas and liquid clusters to positions outside the inital surface and the formation of a crater at the depth of about 100 nm can be clearly seen. The oscillations on the crater surface between melting and metastable melting (see figure 22) are attributed to elastic (acoustic) oscillations of the thick liquid layer.

Zoom In Zoom Out Reset image size

The classical MD method is based on the solution of a number of Newton's equations, each of them corresponding to the motion of a particular atom. While the underlying equation is thus a rather simple one, the challenge lies in the choice of the interaction potential between the atoms and in the numerical complexity of simulating a huge amount of particles to reflect most realistic macroscopic properties and dynamics of matter.

The simplest example of a potential U used in MD simulations is the Lennard-Jones (LJ) potential,

$$\begin{eqnarray}&&U\left({{r}_{ij}}\right)=4\,\epsilon \left[{{\left(\frac{\sigma}{{{r}_{ij}}}\right)}^{12}}\,-\,{{\left(\frac{\sigma}{{{r}_{ij}}}\right)}^{6}}\right],\end{eqnarray} \tag{ 35 }$$

where ${{r}_{ij}}=|\left\{{{\mathbf{r}}_{i}}\right\}-{{\mathbf{r}}_{j}}|$ is the distance between the atoms at the respective positions $\left\{{{\mathbf{r}}_{i}}\right\}$ and ${{\mathbf{r}}_{j}}$. The energy and space parameters and σ can be determined through comparison with macroscopic properties of matter like the bulk modulus and the equilibrium volume. Note that this determination depends on the choice of these macroscopic parameters as well as on the fitting procedure [198]. The first term in equation (35) essentially corresponds to the repulsive force between atoms while the second term leads to an attractive force at longer distances. Such types of potentials, comprising interactions between pairs of atoms, are called pair potentials. Despite its simplicity, several important results on ablation dynamics have been obtained with large scale MD simulations assuming an LJ interaction potential. Some of these will be reviewed in the following sections.

An important advantage of MD based descriptions as compared for instance to hydrodynamic approaches (see section 6) is the a priori inclusion of nucleation kinetics and phase transitions. In contrast to hydrodynamic descriptions, where an instantaneous [178] or a delayed [192] transition to the two-phase mixture may be assumed, the MD simulations do not have to make any assumptions on nucleation kinetics. On the contrary, important information on the nucleation kinetics are gained from MD calculations. Molecular dynamic methods are therefore particularly suitable in studying fast non-equilibrium phase transitions. For instance when studying the expansion of overheated matter and the resulting density profile in sections 6.1 and 6.2, an instantaneous drop in sound velocity was the basis of the discussion. However, in reality this drop may proceed continuously, since the material proceeds through metastable states. It was in fact an important result of the first large scale molecular dynamic simulations in this context [183] that the general picture of a moving shell in front of a two phase region (see figure 19) also holds when metastable states are inherently included [175, 178, 183].

Quantities like the local density, temperature or pressure are averaged quantities. They can be extracted in MD simulations at each time step of the calculation in order to compare with hydrodynamic or other continuous approaches or to study nucleation kinetics in detail. The local density averaged over N particles, which are typically a small fraction of the total number of simulated particles, is straightforwardly defined through the volume V_N, these N particles occupy in the simulation region, thus

$$\begin{eqnarray}&&n=N/{{V}_{N}}.\end{eqnarray} \tag{ 36 }$$

The local temperature T_N is usually determined with the ideal gas approximation,

$$\begin{eqnarray}&&\frac{3}{2}n{{k}_{\text{B}}}{{T}_{N}}=\underset{j=1}{\overset{N}{\sum}}\frac{m}{2}{{(\mathbf{v}_{j}^{T})}^{2}},\end{eqnarray} \tag{ 37 }$$

where k_B is the Boltzmann constant, m the mass of the simulated particles and $\mathbf{v}_{j}^{T}=\text{d}{{\mathbf{r}}_{j}}\,/\text{d}t-{{\mathbf{v}}^{c}}$ the thermal velocity of the particle j, given by the change of its coordinate ${{\mathbf{r}}_{j}}$ and the center of mass velocity ${{\mathbf{v}}^{c}}$ of the corresponding computational cell of N atoms. The local pressure is often calculated as a sum of the ideal gas contribution and the virial term as

$$\begin{eqnarray}&&{{p}_{N}}=n{{k}_{\text{B}}}{{T}_{N}}+\frac{n}{6}\underset{j=1}{\overset{N}{\sum}}{{\mathbf{F}}_{j}}\cdot {{\mathbf{r}}_{j,N}},\end{eqnarray} \tag{ 38 }$$

where ${{\mathbf{F}}_{j}}$ is the force acting on particle j and ${{\mathbf{r}}_{j,N}}$ is the distance of particle j from the center of the considered volume [199].

The extraction of such quantities allows direct comparison of molecular dynamic results with hydrodynamic calculations, as performed for example in [178, 183, 197]. Some of these results will be reviewed in the following sections.

Large-scale simulations of the dynamics of 10–100 millions of particles have been performed assuming a Lennard-Jones potential and a semi-infinite target with periodic boundary conditions in two dimensions and free expansion in the third [176–178]. Here, an exponential heating profile has been assumed with its maximum at the free boundary of the simulation region.

Figure 23 shows the evolution of the ablating material, heated from the right-hand side and expanding towards it. From top to bottom different stages of the ablation process are shown, namely a melting front proceeding into the bulk matter, expansion of the surface, formation of a foam (two-phase region) consisting first of bubbles in liquid phase which later coalesce, leading to a droplet–vapor mixture.

Zoom In Zoom Out Reset image size

These dynamics reflect the results of hydrodynamic calculations as introduced in the preceding section. Qualitatively, figure 23 can be directly compared with figure 19. The upper picture in figure 23 shows the rarefaction of the material and a sharp density gradient at the expansion front as was found from the self-similar solution of expansion, compare figure 19(a). The gradient conserves as time proceeds and below the surface a phase transition to gaseous phase takes place, forming a foam structure behind a condensed shell moving further outwards. The low-density two-phase region expands and the density is decreasing further, as it was also concluded from the self-similar hydrodynamic calculations, compare figure 19(c). The shell is conserved during expansion. In the bottom picture of figure 23 it is marked with two vertical lines.

Note that in the molecular dynamic simulation the phase transition from solid to liquid (from ordered to disordered state) is also visible. A solid-density region, the boundary of which is also marked in the bottom picture of figure 23, will remain after ablation, forming the bottom of the ablation crater.

Such simulation results of ablation dynamics also compare quantitatively well with results of hydrodynamic approaches [178, 183]. A thorough discussion of the mentioned ablation dynamics comparing hydro- and molecular dynamic simulations qualitatively and quantitatively as well as analyzing them with respect to experimental observations with the help of ultrafast microscopy and interferometry was presented by Inogamov et al in [178].

A complete view of the ablation cupola is represented with expansion pictures for the same instants of time after initial excitation but for different excitation strengths. Figures 24 and 25 show the expansion for decreasing excitation strengths from left to right.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 24 represents an early stage of the formation of the liquid shell, which is visible for excitations marked with ${{T}_{0}}=3~\text{MDU}$ and ${{T}_{0}}=2~\text{MDU}$. The temperatures are given in Lennard-Jones units ('moleculardynamic units', MDU), i.e. they are normalised to the energy parameter of the Lennard-Jones potential , whereas the depth scales with σ, see equation (35). Here, $T=1.15~\text{MDU}$ corresponds to the ablation threshold. The pre-assumed heating profile is exponential with a characteristic length of $600~\sigma$, which corresponds for the applied simulation to a depth of d_T  =  204 nm. In figure 24 the undisturbed crystal phase can be identified as well as a molten region for lower excitations and a liquid–gas mixture for higher excitations. The threshold for the liquid shell in front of the two-phase mixture as discussed in section 6.2 lies in-between the two highest excitations; for ${{T}_{0}}=5~\text{MDU}$ no shell but a clear boundary of dense vapor is observed.

Figure 25 shows a later stage of expansion. At that time, droplets in the gas phase have formed, as can be clearly seen for the two highest excitations. On the other hand, in the panel for ${{T}_{0}}=1.25~\text{MDU}$, marking approximately the ablation threshold, bubbles have also appeared deep in the material. The physical processes close to ablation threshold are discussed in more detail in the following section. Note that similar viewgraphs of the fluence-dependent expansion profile have also been obtained for an interatomic potential of a metal [200, 201].

In figure 20 the ablation threshold is marked at a radius of r_cr. In contrast to this sketched profile, the material expands continuously in the real experiment. This is also visible in the MD simulation in figures 24 and 25 showing a slight excursion of the material surface also for ${{T}_{0}}=1~\text{MDU}$, which is slightly below ablation threshold.

The ablation threshold was investigated in detail with the help of MD simulations in [197]. Figure 26 shows the phase-space trajectory of expanding aluminum in the vicinity of ablation threshold. A homogeneously heated slab expands freely in one dimension with laterally periodic boundary conditions. The expansion in both directions leads to overlapping of the rarefaction waves corresponding to the image of a reflected rarefaction wave as discussed in section 6.2 [183]. While in [183] the complete picture of the shell structure could be confirmed for the first time, in [197] the thermodynamic pathway of the central two-phase region was at the focus of the study. To that end, density, pressure and temperature were evaluated for the central third of the simulation volume. The three parts of figure 26 represent different excitation strengths, increasing from top to bottom. In each figure, four isentropes are also plotted, obtained from an aluminum EOS [180]. Note that these isentropes also include metastable states, in contrast to those in figure 17, which are based on the same data.

Zoom In Zoom Out Reset image size

The interpretation of the thermodynamic pathways shown in figure 26 was possible through molecular dynamic snapshots and the relation to the corresponding thermodynamic phases [197]. The upper part of figure 26 shows expansion along an isentrope into a region of negative pressure, i.e. tensile stress, then pressure relaxation connected with contraction, and further acoustic oscillations around zero pressure. The density is reduced as compared to the initial situation, but the material is still in liquid phase. The central figure represents the expansion slightly below ablation threshold. The material expands along an isentrope, reaches metastable states at negative pressure and relaxes to nearly zero pressure. In that state it contains bubbles in a liquid phase. However, the initial energisation was not strong enough for a complete transition into gaseous phase and also, after the observed partial phase transition, the material is still under tensile stress—therefore it contracts again. The final state is thus in liquid phase and no ablation has occurred. The lowest part of figure 26 shows the thermodynamic pathway of aluminum slightly above ablation threshold: the material expands again to negative pressure, the gas phase nucleates while the material is still expanding. Zero pressure is reached, however, in this case it is slightly positive pressure. Therefore, after a short contraction, the material further expands into the gas phase—thus it is ablated.

Note the difference in the timescales of contraction in the two cases below ablation threshold. Each point on the trajectory was determined 0.2 picoseconds after the preceding one. Thus, close to ablation threshold the expansion recovers on a much longer time scale than further below ablation threshold. The timescale of expansion and contraction in the central picture of figure 26 by far exceeds that of acoustic oscillations in the liquid phase. Such long-scale surface excursions in the vicinity of ablation threshold were found experimentally [53, 202] and analyzed with help of hydrodynamic and molecular dynamic simulations [177, 197]. The reason can be attributed to the formation of gaseous bubbles, also below ablation threshold. Recently, a hybrid simulation including fast resolidification has shown that such bubbles can be conserved as a void structure below a re-crystalised silver surface [203].

The choice of a suitable potential is crucial in any MD simulation. Usually it is chosen with the aim of describing a real material with a most accurate interaction potential, reproducing known properties like crystal structure, melting temperature or sound velocity. It is likely that any chosen potential works best (meaning most realistically) only in a certain range of temperature, density or pressure. Finding a 'good potential' is a challenging task and beyond the scope of this review.

Instead, we briefly show in this section, which differences have been observed in the previously discussed works, appearing due to the particular choice of the MD potential. For instance, the central panels of figure 25 indicate that gaseous bubbles in expanding Lennard-Jones material appear in the solid phase. The reason lies in a much higher melting temperature of approximately 50% of the critical temperature, than obtained for metal potentials, e.g. aluminum where the melting temperature is only about 10% of the critical temperature [177, 197]. The threshold isentrope leading to ablation of the LJ material may therefore not enter the liquid state [197, 204].

The threshold energies for ablation, as well as those for melting and bubble formation, have been compared for different metals as well as LJ material in [197]. The main difference in the potentials is that the LJ potential is a pair potential whereas the metals were simulated with many body potentials reflecting the delocalised metallic binding. Figure 27 (top) shows the threshold energies for ablation (spallation), as well as those for melting and bubble formation (void) for LJ material and different metals [197]. The energies of initial excitation E₀ were normalised to the cohesive energy E_coh of the material, directly resulting from the applied respective potentials. It can be seen that breaking of the crystal structure is most expensive in terms of energy for the LJ material, whereas the energy range for melting and void formation are smaller than in all considered metals. Differences in melting thresholds were studied on atomistic scale in [177], where the potential barriers for relocation of a single atom in an fcc lattice were analyzed. The energetic costs for such displacement were much lower in aluminum than in LJ material and silicon.

Zoom In Zoom Out Reset image size

Figure 27 (bottom) shows also a comparison of a small simulation region (30 monolayers, as applied for the simulations presented in section 7.3) with a larger one (300 monolayers, used in the simulations presented in section 7.2). The increase in the thickness of the simulated film has a clear effect on the energy ranges for the observed phases. In particular for the thick film of LJ material the formation of voids appears directly from the solid phase. This behaviour was also observed in [204] and studied in dependence on the crystal's orientation.

It was pointed out in [175] that the transverse dimensions of the simulation region (perpendicular to the direction of expansion) are also important. A better agreement of MD simulations with hydrodynamic calculations has been achieved for larger lateral dimensions.

When heat conduction effects come into play, especially in metals, the electronic heat transport has to be taken into account. A possible way to model such conditions is the combination of the two temperature description with a molecular dynamic simulation. Such hybrid combination was introduced in [205] and further developed in [114, 126, 206]. In these simulations the heating of the electrons by the laser and the further energy dissipation by electron–phonon coupling as well as electron thermal transport is described in the frame of the two temperature model (TTM), while the movement of the lattice atoms is simulated with molecular dynamic (MD) methods. Here we briefly review the principal method as presented in [114].

To combine the two temperature model with a molecular dynamic simulation, the latter equation of the system (2) and (3) was substituted by an equation of motion for the lattice atoms and a coupling term accounting for the lattice heating by the electrons. The resulting equation system reads:

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{c}_{\text{e}}}\,\frac{\partial {{T}_{\text{e}}}}{\partial t} & = & \frac{\partial}{\partial z}\,\left({{\kappa}_{\text{e}}}\,\frac{\partial}{\partial z}{{T}_{\text{e}}}\right)\,-\alpha \,\left({{T}_{\text{e}}}-{{T}_{\text{i}}}\right)+Q(z,t), \end{array}\end{eqnarray} \tag{ 39 }$$

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{m}_{j}}\frac{{{\text{d}}^{2}}{{\mathbf{r}}_{j}}}{\text{d}{{t}^{2}}} & = & {{\mathbf{F}}_{j}}+\mathbf{F}_{j,\text{therm}}^{e-i}. \end{array}\end{eqnarray} \tag{ 40 }$$

Here, equation (40) has to be solved for each considered atom j with mass m_j at position ${{\mathbf{r}}_{j}}$. The force ${{\mathbf{F}}_{j}}$ acts on atom j and depends on the chosen interatomic potential. The force $\mathbf{F}_{j,\text{therm}}^{e-i}$ enters due to electron-lattice heating. The energy received by the lattice is transferred to motion of the lattice ions. This term is balanced with the loss of electron energy $\alpha \left({{T}_{\text{e}}}-{{T}_{\text{i}}}\right)$ to ensure energy conservation. $\mathbf{F}_{j,\text{therm}}^{e-i}$ is directly proportional to the momentum of particle j. In [114] it was shown that the thermal velocity $\mathbf{v}_{j}^{T}=\text{d}{{\mathbf{r}}_{j}}\,/\text{d}t-{{\mathbf{v}}^{c}}$ has to be considered here, thus $\mathbf{F}_{j,\text{therm}}^{e-i}=\xi {{m}_{j}}\mathbf{v}_{j}^{T}$, where ${{\mathbf{v}}^{c}}$ is the velocity of the center of mass of the computational cell to which atom j belongs, see section 7.1. The factor ξ which ensures energy conservation, is deduced in [114].

To solve equation system (39) and (40) numerically, the main challenge is to couple a finite difference scheme for the electronic heat conduction with a molecular dynamic (MD) simulation accounting for atomic movement of the surface atoms and a continuum description for the bulk of the material. Figure 28 shows schematically the parts of the model and the energy exchange between them. The MD method is used only in the very surface region of the target, where active processes of laser melting and ablation are taking place. The diffusion equation for the electron temperature is solved in a much wider region affected by the thermal conduction from the absorbing surface layer. In order to avoid reflection of the pressure waves propagating in the MD part from the irradiated surface, a dynamic boundary condition that mimics the interaction of the atoms in the boundary region with the outer 'infinite elastic medium' should be applied [206]. In the outer part the lattice temperature is calculated using equation (3), usually neglecting the ionic heat conduction ${{\kappa}_{\text{i}}}$.

Zoom In Zoom Out Reset image size

This scheme allows studying the influence of the electron dynamics on structural changes due to laser irradiation. For short laser pulses a prolonged lattice heating can thus be included and in contrast to the studies presented in the preceding section 7 the pre-assumption of a certain energisation of the atoms can be avoided.

In [114, 126, 206] the hybrid TTM–MD has been developed and kinetics of melting and disintegration of a metal film have been observed on an atomic scale.

In [206] a coherent emission of ablated material was found for laser irradiated copper. This picture is consistent with the observation of Newton fringes [6] and the results of large scale MD simulations applying the simplified LJ potential [178, 183] as reported in section 7. The application of a realistic metal potential in combination with a consideration of the electron dynamics of excitation and relaxation led to a good agreement of calculated ablation thresholds with experiment [206].

In [114] melting and disintegration of free standing nickel and gold films have been studied. Figure 29 shows contour plots of the resulting lattice temperature (a) and pressure (b) for irradiation of a 50 nm nickel film in the vicinity of but above the ablation threshold. A laser pulse of 200 fs duration was assumed to irradiate on the target from the bottom of the plots. On a timescale of a few picoseconds, the lattice temperature and the pressure within the material increase strongly. The onset and completion of melting are marked by black dotted and solid lines respectively. It starts from both free surfaces; however, the melting from the back of the film starts with delay, which results directly from the electronic heat transport assumed in this calculation. Melting kinetics on the atomic scale have been shown and discussed in section 5. Another effect of the implementation of the electron dynamics is the gradual lattice heating. Figure 29(a) shows a temperature increase at the irradiated surface of the film lasting about 20 ps. This timescale is connected with the electron–phonon coupling discussed in section 3. For the calculation of figure 29 a comparably large, constant coupling parameter was assumed. The strong expansion of the film observed on the same timescale leads finally to two voids in the center of the irradiated film. The rupture follows a strong tensile pressure, compare figure 29(b). After disintegration of the film, pressure waves, visible as blue and green zigzag lines, appear in both parts of the disintegrated film.

Zoom In Zoom Out Reset image size

Figure 30 shows a temperature plot for the similar case of the irradiation of a 50 nm film of gold, irradiated with a 200 fs laser pulse slightly above ablation threshold. We compare it to the case of nickel, figure 29. In gold, the lattice temperature increases for a duration up to 75 ps. The assumed coupling constant was about one order of magnitude lower than that for nickel. Another difference between the two materials appears in the temperature profile. The gold film shows a homogeneous temperature distribution, which is due to a large electronic heat conductivity. As a result, both parts of the disintegrated film are equally heated, while in nickel the disintegration of the film separates a hot front layer from a colder back layer. The effects described above can only be observed in the frame of a TTM–MD description when including a description of electron dynamics as heat transport and electron–phonon interaction in the MD simulation.

Zoom In Zoom Out Reset image size

Results for a bulk sample of nickel were presented in [207] and are shown in figure 31. The figure shows contour plots of lattice temperature (left plots) and pressure (right plots) for nickel irradiated with a 1 ps laser pulse of increasing fluence from top to bottom, i.e. from (a) to (d). The black line marks in each plot the border between liquid and crystal. For moderate excitations, inducing melting but no ablation, the figure (a) shows heating, heat transport, a pressure wave of positive and negative pressure proceeding into the material and expansion at the free surface. The phase border induces a slight discontinuity for most of these features. For slightly higher fluence, ablation of a single liquid layer is observed, as was also reported in the works reviewed in section 6. A pressure oscillation occurs in the ablated layer, visible in figure 31(b) right-hand side. It appears both in the ablated layer and in the remaining bulk. The surface of reflection of the bulk pressure wave is provided by the liquid–solid interface. Further above the ablation threshold (c) the ablated material consists of several layers or droplets, which are of decreasing size for even higher fluences (d).

Zoom In Zoom Out Reset image size

When semiconductors or dielectrics are to be described in the frame of a TTM–MD or similar scheme, the implementation of the electron dynamics becomes more complicated. In these materials, the density of conduction band electrons is not constant (compare section 2), and the transport of heat and also particles has to be considered (compare section 4).

Combinations of the density dependent TTM [34] with molecular dynamic simulations have recently been realised [21, 117, 119]. In [21, 117] the influence of certain parameters and approximations on the resulting dynamics have been discussed. The results of [119] show the interplay of pressure and temperature and their influence on the melting kinetics. Note that ablation has not been studied, since a changing potential energy surface (see section 5.2) has not yet been implemented in this large scale simulation (see sections 5.2 and 8.4).

Another way of implementing carrier dynamics into a MD simulation of semiconductors was chosen in [208]. Here, the authors applied a Monte Carlo approach to account for carrier heating on femtosecond timescales, together with a molecular dynamic calculation of the subsequent atomic displacement on pico- to nanosecond timescales. Picosecond heating followed by expansion up to the nanosecond regime was also investigated. The initial carrier dynamics were taken explicitly into account, including generation and transport of free carriers and phonon emission. The transition from the semiconducting solid silicon to a metallic-like liquid state has also been implemented. It drastically changes the reflectivity of the material (compare section 2). The authors attribute the onset of ablation to a direct conversion of translational mechanical energy into surface energy. The general view of ablating material as well as the pathway in the phase diagram look similar to those in the works reviewed above [178, 183, 197]. Figure 32 shows the near-threshold ablation after femtosecond excitation. Nucleation of a bubble can be observed—due to the small calculation region only one bubble is visible—as well as ablation of a liquid shell.

Zoom In Zoom Out Reset image size

Neither of the approaches presented in the previous section implementing electronic excitation in semiconductors to molecular dynamics consider changes of the potential due to electronic excitation. As we have discussed in section 5.2, such changes of the potential energy surface may, however, lead to strong changes of the atomic bonds and should thus be considered in a consistent molecular dynamic simulation, in particular in semiconductors exhibiting strong covalent bonds. The implementation of a transiently changing potential surface into such an approach is, however, challenging. Direct combinations of first principles calculations with MD simulations have been employed by different groups and are mentioned in Chapter 5.

Also, different molecular dynamics simulations including the dynamics of electrons within a tight-binding scheme have been performed, giving insight into the mechanism of laser induced lattice instabilities as a driving force for nonthermal melting. One of those approaches was based on the numerical solution of the time-dependent Schrödinger equation [209–211], whereas in the other approach an equation of motion for the density matrix was solved together with the molecular dynamics [212–216]. In spite of the model character of the tight-binding method, good agreement was obtained with experiments regarding pre-ablation and ablation thresholds in graphite [216] and regarding the dependence of the ablation and melting thresholds in Si as a function of the pulse duration [215]. Figures 33 and 34 show results of such hybrid approach for graphite combining a tight-binding description of electron dynamics with molecular dynamics simulations based on a potential energy surface calculated using the same tight binding Hamiltonian. Incoherent electron–phonon interactions and diffusion were also taken into account, but using a phenomenological approach [216]. In figure 33, a lower excitation has been assumed, which does not show melting of the graphite planes. However, a pre-ablation regime consisting in the successive ejection of intact and ordered planes has been obtained. This result was experimentally confirmed by Zewail and co-workers through ultrafast electron crystallography experiments [217]. The existence of this pre-ablation regime was exploited for producing ultra thin graphitic flakes by femtosecond laser irradiation of graphite [218].

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

For comparison to the pre-ablation regime, figure 34 shows the complete disorder of the graphite structure on a timescale of about 80 fs, solely due to electronic excitation. A liquid-like state and the ablation of single carbon atoms and chains is observed. Thus transient low density liquid carbon can be produced by laser excitation of graphite [216].