The Role of Crystalline Orientation in the Formation of Surface Patterns on Solids Irradiated with Femtosecond Laser Double Pulses

Featured Application

Laser processing appears to be a scalable method to fabricate highly ordered and complex biomimetic surface structures which can be used in an abundance of industrial applications. On the other hand, a multiscale theoretical approach represents a powerful tool that can provide significant information for designing materials. The approach is aimed to offer a cost-effective method towards minimising the use of trial and error approaches to reduce a large number of experiments towards fabricating application-based patterns.

Abstract

A theoretical investigation of the underlying ultrafast processes upon irradiation of rutile TiO₂ of (001) and (100) surface orientation with femtosecond (fs) double pulsed lasers was performed in ablation conditions, for which, apart from mass removal, phase transformation and surface modification of the heated solid were induced. A parametric study was followed to correlate the transient carrier density and the produced lattice temperature with the laser fluence, pulse separation and the induced damage. The simulations showed that both temporal separation and crystal orientation influence the surface pattern, while both the carrier density and temperature drop gradually to a minimum value at temporal separation equal to twice the pulse separation that remain constant at long delays. Carrier dynamics, interference of the laser beam with the excited surface waves, thermal response and fluid transport at various pulse delays explained the formation of either subwavelength or suprawavelength structures. The significant role of the crystalline anisotropy is illustrated through the presentation of representative experimental results correlated with the theoretical predictions.

1. Introduction

The impact of the employment of ultra-short pulsed laser sources in material processing has received considerable attention due to its important applications, in particular in industry and medicine [1,2,3,4,5,6,7,8,9,10]. Due to its unique advantages, the laser-assisted processing method (LAPM) emerged as a significant tool for precise fabrication at nano- and micro-scales.

Laser-induced periodic surface structures (LIPSS) on solids is one type of surface patterns that has been explored extensively, becausethe features of the produced structures provide impressive properties that can be used in many applications including microfluidics [1,11], tribology [12,13,14], tissue engineering [11,15] and advanced optics [16,17]. There exists a wealth of reports related to the elucidation of the underlying physical processes for the formation of these structures in a variety of laser conditions (i.e., fluence, pulse duration, laser wavelength, pulse separation, polarisation state, energy dose) and in different materials [17,18,19,20,21]. Some of the most representative types of LIPSS are the coined low-spatial-frequency (LSFL) ripples [22,23], high-spatial-frequency (HSFL) ripples [18,24], grooves [20,25] and spikes [25]. A series of multiphysical phenomena were analysed to explain the physical origins of the LIPSS formation [26,27,28,29,30]. More specifically, various theoretical models have been proposed to interpret the production of periodic structures: interference of the incident wave with an induced scattered wave [31,32,33], or with a surface plasmon wave (SP) [21,22,34,35,36,37], or due to self-organisation mechanisms [38].

Although, the formation of LIPSS at long time scales is closely related to a phase transformation, one process that undoubtedly has to be investigated is the response of the material at short time scales. More specifically, as stated above, the impact of electromagnetic effects and plasma contribution [39], surface plasmon excitation [22,34], wave interference effects [40] are some of the processes that need to be explored when the material is in a non-equilibrium state. Therefore, the ultrafast dynamics of the excited carriers is an important ingredient that greatly influences the LIPSS periodicities. To provide a comprehensive elucidation of the ultrafast dynamics of the charge carriers, it is important to elaborate on the influence of the laser parameters on the excited plasma characteristics and then, on the thermal response of a material. The existing multiscale models and experimental studies have successfully provided a detailed analysis of the physical mechanisms that lead to a plethora of structural effects [22,41,42,43,44,45]. In principle, the laser beam parameters determine the first steps of the processes because energy absorption and relaxation processes influence significantly the material heating [43].

Nevertheless, although for many semiconducting materials, the physical mechanism that describes ultrafast dynamics is well-established and the theoretical model works effectively in various conditions, there is still a missing picture for some types of semiconductors such as TiO₂ and more specifically for materials with different crystal orientations. In particular, in a previous report, a parametric study of the laser conditions and the induced structures revealed that there is a difference in the response of the material if crystals with orientation (001) and (100) are irradiated [46,47,48]. In that study, it was emphasised that a key difference between the two crystal orientations is the effective electron mass (apart from the fact that interactions between atoms, thermomechanical response and melting are strongly affected by the interatomic distance which varies with the crystal orientation [49]); the discrepancy of the optical mass (measured or calculated in precious reports [47,48,50,51,52,53,54]) affects optical parameters of the material, amount of the absorbed energy and, eventually, thermal response of the system that is also projected on the type of LIPSS that are formed [46]. Nevertheless, despite these early results and conclusions, a more elaborated approach is required to evaluate more explicitly the impact of crystal orientation in the energy absorption; it is noted that TiO₂ is a uniaxial birefringent material with optical axis parallel to direction (001). Therefore, in the case of an irradiation under normal incidence of the face (001), the electromagnetic wave of laser propagates according to the normal refractive index $\left(n_{\perp }, k_{\perp }\right)$. On the other hand, in the case of irradiation under normal incidence of the face (100) or (110), the optical axis lies in the plane of the surface. Thus, there are waves propagating according to refractive indexes $\left(n_{\perp }, k_{\perp }\right)$ and $\left(n_{\parallel }, k_{\parallel }\right)$ respectively. Furthermore, although, the difference in the optical response of the system induced by the distinct crystal orientation for irradiation with femtosecond pulses has been investigated [46], no previous study has been conducted to explore the impact of the shape of the laser pulse and more specifically of temporarily separated pulses. In previous reports, the effect of the irradiation of various types of materials with temporarily separated pulses has been investigated which revealed a significant variation in the ultrafast dynamics, thermal response of the material and finally, changes in the LIPSS frequencies [29,36,55]. Therefore, and given the role of the crystal orientation, a detailed description of the ultrafast dynamics resulting from irradiating semiconductors with different crystal orientation with double pulses of variable temporal separation would provide significant information. More specifically, elucidation of this issue is of paramount importance, firstly, to elucidate the underlying multiscale physical mechanisms of laser-matter interaction and secondly to associate the resulting thermal effects with the induced surface patterns.

To this end, we present a detailed theoretical approach that describes the ultrafast dynamics in rutile TiO₂ [22,26,30,56,57,58], to account firstly, for excitation and electron-phonon relaxation upon irradiation of TiO₂ in two crystal orientations with periodic ultrashort laser double-pulses separated by a temporal delay from zero to tens of picoseconds (Section 2). The details of an experimental protocol that has been developed is described in Section 3, while a thorough analysis of results and validation of the theoretical model are presented in Section 4 by discussing the ultrafast dynamics, thermal response and induced surface modification predictions. Concluding remarks follow in Section 5.

2. Theoretical Model

2.1. Ultrafast Dynamics

It is known that during laser irradiation of a material, various physical processes occur on a femtosecond timescale. In this work, excitation of a semiconducting material, rutile TiO₂ with (001) and (100) surface orientation was considered through a laser beam of wavelength λ_L = 248 nm that corresponds to photon energies equal to ~5 eV which is higher than the band gap of the material (~3.0 eV, at 300 K); thus, it is assumed that predominantly one- and (secondly and with a substantially smaller probability) two-photon absorption mechanisms contribute to excitation of the charge carriers in the valence and conduction bands, while higher order photon processes are less likely to occur. Furthermore, linear free carrier photon absorption on the intraband transitions are possibly increase the carrier energy (but not the carrier density), while Auger recombination and impact ionization processes lead to a decrease and increase, respectively, of the number of carriers in the conduction band.

The relaxation time approximation to Boltzmann’s transport equation [22,26,30,56,57,58] is employed to determine the spatial ($\overrightarrow{r} = \left(x,y,z\right)$, in Cartesian coordinates) and temporal dependence (t) of the carrier density number, carrier energy and lattice energy; more specifically, the following set of coupled energy and particle balance equations are used to describe the temporal variation of the carrier density number N_c, carrier temperature T_c and lattice temperature T_L:

$$\begin{array}{l}{C}_c\frac{\partial T_c}{\partial t}=-\frac{C_c}{{\tau }_{c-l}}\left(T_c-T_L\right)+S(\overrightarrow{r},t)\\ C_L\frac{\partial T_L}{\partial t}=\overrightarrow{\nabla }\cdot \left(K_L\overrightarrow{\nabla }{T}_L\right)+\frac{C_c}{{\tau }_{c-l}}\left(T_c-T_L\right)\\ \frac{\partial N_c}{\partial t}=\frac{{\alpha }_{\mathrm{S}\mathrm{P}\mathrm{A}}}{\hslash {\omega }_L}I(\overrightarrow{r},t)+\frac{{\beta }_{TPA}}{2\hslash {\omega }_L}{I}^2(\overrightarrow{r},t)-\gamma N_c{}^3+\theta N_c\end{array}$$

(1)

where C_c (C_L) is the carrier (lattice) specific heat, K_L corresponds to the lattice heat conductivity, $\hslash {\omega }_L$ stands for the photon energy and ${\alpha }_{SPA}$ and ${\beta }_{TPA}$ correspond to the single and two-photon absorption coefficients, respectively. C_c is provided by the expression $C_c = 3N_c{k}_B$, where $k_B$ stands for the Boltzmann constant. The lattice specific heat and heat conductivity for TiO₂ are C_L = 690 J Kgr⁻¹ K⁻¹ [59] and K_L = 13 Wm⁻¹ K⁻¹ [60], respectively. At λ_L = 248 nm, ${\alpha }_{SPA}\sim {10}^6 \times \exp \left(T_L/300\right)$ cm⁻¹ [61] while ${\beta }_{TPA}\sim$37 cm/GW [46,62]. In Equation (1), τ_c-l is the carrier-phonon energy scattering time which is derived from $1/{\tau }_{c-l}\simeq {10}^4{N}_c^{1/3}$ s⁻¹ [50]. Finally, γ is the coefficient for Auger recombination (γ~10⁻³² cm⁶/s [63]) and θ is the impact ionization rate coefficient. In this work, an intensity dependent expression is used for the impact ionization parameter θ [33,64] (unlike an approximating formula used for other semiconductors and at higher laser wavelengths, i.e., $~e^{-E_g/k_B{T}_c}$ s⁻¹, for Silicon at 800 nm [22,57].

$$\begin{array}{l}\mathit{\theta }=\frac{e^2{\tau }_{c}I}{{c\epsilon }_0{nm}_r{E}_g\left({\omega }_L{}^2{\left({\tau }_c\right)}^2+1\right)\left(2{-m}_r{/m}_e\right)}\\ \end{array}$$

(2)

In Equation (2), c is the speed of light, e is the electron charge, ε₀ stands for the vacuum permittivity, m_r (i.e., m_r = 0.5m_e) is the reduced electron mass and n is the refractive index of the material, while I is the peak intensity of the laser beam, and τ_c is the electron collision time. Based on previous reports [50], scattering rates for electron-hole and electron-electron collisions scale as $N_c$ and $N_c^{1/3}$, respectively, yielding approximate value τ_c~1 fs for TiO₂ regardless of the orientation. This is the value that was used in this work. $E_g$ stands for the energy band gap of rutile TiO₂ (~3 eV) [46]. The ‘source term’ that appears in the first equation in Equation (1) [22,56,57] is given by

$$S(\overrightarrow{r},t)=({\alpha }_{SPA}+{\alpha }_{FCA})I(\overrightarrow{r},t)+{\beta }_{TPA}{I}^2(\overrightarrow{r},t)-\frac{\partial N_c}{\partial t}\left(E_g+3k_B{T}_c\right)$$

(3)

In Equation (1), it is evident that carrier transport and carrier heat diffusion have been ignored. Although more rigorous expressions have been presented in previous works in which the contribution of those effects have been incorporated, reports showed that estimation of damage thresholds are not influenced significantly if those terms are ignored [56,57]. Similarly, for the sake of simplicity, it is assumed that the energy band gap does not vary significantly within the pulse.

The laser intensity $I\left(\overrightarrow{r},t\right)$ in Equations (1)–(3) is obtained by considering the propagation loss due to one-photon, two-photon and free carrier absorption, respectively [22,26,30,36,57,65],

$$\frac{\partial I(\overrightarrow{r},t)}{\partial z}=-({\alpha }_{SPA}+{\alpha }_{FCA})I(\overrightarrow{r},t)-{\beta }_{TPA}{I}^2(\overrightarrow{r},t)$$

(4)

for a laser beam that is Gaussian, both temporally and spatially, the transmitted laser intensity on a flat surface for double pulses is provided by the following expression

$$I(\mathrm{x},\mathrm{y},\mathrm{z}=0,t)=\frac{\sqrt{\mathrm{l}\mathrm{n}2}{E}_p(1-R(\mathrm{z}=0,\mathrm{t}))}{\sqrt{\pi }{\tau }_p}{e}^{-\left(\frac{2({\mathrm{x}}^2+{\mathrm{y}}^2)}{R_0^2}\right)}\left[e^{-4\ln 2{\left(\frac{t-3{\tau }_p}{{\tau }_p}\right)}^2}+e^{-4\ln 2{\left(\frac{t-3{\tau }_p-t_{delay}}{{\tau }_p}\right)}^2}\right]$$

(5)

where 2E_p is the total fluence of the laser beam and τ_p is the pulse duration (i.e., full width at half maximum), R₀ is the irradiation spot-radius (distance from the centre at which the intensity drops to 1/e² of the maximum intensity, and R is the reflectivity while irradiation under normal incidence was assumed. The value of 3τ_p in the exponent in Equation (5) has been chosen to ensure that the intensity switches off at t = 0. The objective of the current work is to illustrate the impact of temporarily separated pulses, therefore, Equation (5) provides the total intensity of a double pulse of pulse delay equal to t_delay. In Equation (5), Cartesian coordinates have been used.

The calculation of the free carrier absorption coefficient and the reflectivity are derived from the dielectric parameter of the material (assuming also corrections due to band and state filling [66]), ${\epsilon }^{\prime }$

$${\epsilon }^{\prime }=1+\left({\epsilon }_{un}-1\right)\left(1-\frac{N_c}{N_v}\right)-\frac{e^2{N}_c}{{\epsilon }_0{\omega }_L{}^2}\frac{1}{m_{opt}^{*}\left(1+i\frac{1}{{\omega }_L{\tau }_c }\right)}$$

(6)

where N_v stands for the valence band carrier density (~5 × 10²³ cm⁻³ for TiO₂) and ε_un is the dielectric parameter of the unexcited material at λ_L = 248 nm. In previous reports, the values of ε_un [67,68,69,70,71,72] and the optical mass $m_{opt}^{*}$ [47,48,50,51,52,53,54] were measured and estimated for both (001) and (100), respectively. As noted above, the dielectric parameter is expected to affect the energy absorption and therefore precise parameter values for $m_{opt}^{*}$ and ε_un are required. Based on the aforementioned works, the values ${\epsilon }_{un}^{\left(100\right)} = 1.3+7.61i$, ${\epsilon }_{un}^{\left(001\right)} = 3.09+6.95i$, $m_{opt}^{*\left(100\right)} = 6m_e$, $m_{opt}^{*\left(001\right)} = 2m_e$ were used in the current study. For a flat surface, the reflectivity and free carrier absorption coefficients are computed through the following expressions

$$\begin{array}{l}{\alpha }_{FCA}(\mathrm{x},\mathrm{y},\mathrm{z},t)=\frac{2{\omega }_L{k}_{ext}}{c}\\ R\left(\mathrm{x},\mathrm{y},\mathrm{z}=0,t\right)=\frac{{(1-n)}^2+k_{ext}^2}{{(1+n)}^2+k_{ext}^2}\end{array}$$

(7)

where k_ext is the extinction coefficient of the material that is calculated from ${\epsilon }^{\prime }$.

2.2. LSFL Formation

Among the physical processes that are aimed to account for the formation of LSFL structures, the most dominant scenarios are based on: (a) the interference of excited Surface Plasmon Polariton waves (SPP) with the incident beam [22,30,34,73], (b) the interference of the incident laser beam with induced surface scattered waves [40] or (c) the role of near field effects [39]. A key ingredient, though, for the excitation of SPP or surface scattered waves is the presence of an initial corrugation (for example, surface defects [22]). In regard to the first scenario that will be adopted in the current work, the interference of the incident and SPP leads to the formation of periodic ripples (LSFL) that are oriented perpendicularly to the electric field of the laser beam [18,21,22,31,32,33,74,75,76,77,78,79,80]; this is generated by a combination of multiscale processes such as the spatially periodic energy distribution (as a result of the interference of the incident beam and a surface wave), an induced characteristic spatial modulation of the lattice system through electron-phonon scattering and effects related to fluid dynamics (for conditions that allow the material to undergo a phase transition) leading eventually to a rippled profile when the resolidification process ends [22,78,79,81].

To describe the excitation of SPP, dispersion conditions are derived through the elaboration of the boundary conditions (continuity of the electric and magnetic fields at the interface between a metallic and dielectric material) (ε_d = 1) for a flat surface (number of pulses, NP = 1). Therefore, a requirement for a semiconductor to obey the above relation and conditions [22,82] is that $Re\left({\epsilon }^{\prime }\right)<-1$ and the computed SPP wavelength Λ_SPP is given by

$${\Lambda }_{\mathrm{S}\mathrm{P}\mathrm{P}}=\frac{{\lambda }_L}{{}_{Re}\sqrt{\left(\frac{{\epsilon }^{\prime }{\epsilon }_d}{{\epsilon }^{\prime }+{\epsilon }_d}\right)}}$$

(8)

The condition $Re\left({\epsilon }^{\prime }\right)<-1$ and Equations (6) and (8) can be used to derive the range of values of the excited carrier densities that lead to SPP excitation. According to the above condition, only carrier densities larger than ~13 × 10²² cm⁻³ (for (100)) are capable to lead to the excitation of SPP (Figure 1a) for nearly flat profiles. Although, Equation (8) can be used to calculate the SPP for nearly flat surfaces, there exists a discrepancy between the experimentally measured ripple periodicities and predictions due to the above expression. This is due to the fact that in order to observe periodic profiles experimentally, a sufficient number of pulses are required. On the other hand, repetitive irradiation is accompanied with a continuously varying depth of the periodic profile which results from a shift of the SPP wavelength to smaller values [21,78,83,84] that leads to SPP frequencies that deviate from the value resulting from Equation (8).

On the other hand, for nearly flat profiles, it is evident that due to the different dielectric parameters and optical masses of TiO₂ for (001 and (100), a larger range of Λ_SPP values for (100) is predicted (~205 nm) compared to the computed value (~232 nm) for (100). The minimum values of Λ_SPP occur for ~27 × 10²² cm⁻³ (001) and ~44 × 10²² cm⁻³ (100), respectively.

To describe the formation of LSFL with scenario (b), the electron density values and surface corrugation are used to determine the orientation and the periodicity of LIPSS. To correlate N_c with the formation of laser-induced surface structures, the inhomogeneous energy deposition into the irradiated material is computed through the computation of the product $\eta \left(\stackrel{\rightharpoonup }{K}, \overrightarrow{k_i}\right) \times \left|b\left(\stackrel{\rightharpoonup }{K}\right)\right|$ [40]. In this expression, $\eta$ stands for the efficacy with which a surface roughness at the wave vector $\stackrel{\rightharpoonup }{K}$ (i.e., $\left|\stackrel{\rightharpoonup }{K}\right|$ = ${\lambda }_L/\mathsf{\Lambda }$, where $\mathsf{\Lambda }$ is the predicted structural periodicity along the direction in which the periodic structures are formed) induces inhomogeneous radiation absorption, $\overrightarrow{k_i}$ stands for the component of the wave vector of the incident laser beam on the material’s surface plane and b is a measure of the amplitude of the surface roughness at $\stackrel{\rightharpoonup }{K}$. In principle, surface roughness is represented by spherically shaped islands and standard values from Sipe’s theory for the ‘shape’ (s = 0.4) and the ‘filling’ (f = 0.1) factors [25,40,85] are considered. In the current work, linearly p-polarised laser beams are used and based on Sipe’s theory, LIPSS are formed where $\eta$-maps show sharp features. The position of maximum value of $\eta$ along K_x (Figure 1b,c) for 27 × 10²² cm⁻³ (001) and ~44 × 10²² cm⁻³ (100) yields $\mathsf{\Lambda } =$233 nm and $\mathsf{\Lambda } =$206 nm, respectively, which are close to the one derived from Equation (8). According to Figure 1b,c, the produced structures are oriented perpendicularly to the polarisation vector of the incident beam.

2.3. Fluid Transport. Ripples and Grooves

To model the surface modification of TiO₂ upon irradiation of material with repetitive pulses, the laser fluence is assumed to be sufficiently high that results in a phase transformation. The melting temperature of the material ($T_{melting}$ = 2110 K [86]) is considered as the threshold at which phase transformation occurs. Following a slow resolidification process, a non-flat relief is induced on the material. Moreover, the fluence values used for the simulations and experiments in the current work are sufficiently high to cause ablation. To simulate mass ablation, it is assumed that it occurs if the material is heated above a critical temperature (~$T_L$>$T_L^{boiling}$ where $T_L^{boiling}$ ~3245 K for TiO₂). It is noted that the choice of the critical temperature T_critical above which mass removal occurs has been investigated in previous studies. According to the discussion in reference [87], a typical value of T_critical is 1–2 times higher than the boiling temperature. In other studies, it was noted that a solid material which is heated with ultrashort pulsed lasers at sufficiently high fluences undergoes a phase transition to a superheated liquid whose temperature reaches values ~$0.90 \times T_{cr}$ ($T_{cr}$ stands for the thermodynamic critical point [88]) and evaporation due to dynamics of the Knudsen layer (adjacent to the liquid-vapor interface [22,88]). To argue the choice of the ablation temperature threshold, it would be useful to know an estimate of $T_{cr}$, however, to the best of our knowledge, a relevant value for $T_{cr}$ for TiO₂ does not exist in the literature (except from a very high value for Ti (~15,500 K [89])). Based on the above discussion, for the sake of simplicity, the boiling temperature was selected to represent the threshold for mass removal.

To describe the dynamics of an uncompressible fluid a Navier-Stokes equation (NSE) is used [83]

$${\rho }_0\left(\frac{\partial \overrightarrow{u}}{\partial t}+\overrightarrow{u}·\overrightarrow{\nabla }\overrightarrow{u}\right) = \overrightarrow{\nabla }·\left(-P+\mu \left(\overrightarrow{\nabla }\overrightarrow{u}\right)+\mu {\left(\overrightarrow{\nabla }\overrightarrow{u}\right)}^T\right)$$

(9)

The boundary conditions imposed at the liquid free surface of the material are the following: $\frac{\partial u}{\partial Z} = -\sigma /\mu \frac{\partial T_L}{\partial X}$ and $\frac{\partial v}{\partial Z} = -\sigma /\mu \frac{\partial T_L}{\partial Y}$, where (u,v,w) are the components of $\overrightarrow{u}$ in Cartesian coordinates. ${\rho }_0$ (${\rho }_0 = 3.8375-2.8 \times {10}^{-4} T_L$ gr/cm³) [90] and $\mu$ stand for the density and viscosity of molten TiO₂, while P and $\overrightarrow{u}$ stand for the pressure and velocity of the fluid, respectively. On the other hand, σ corresponds to the surface tension of the material (see [9], for a more detailed description of the fluid dynamics module; similarly, temperature dependent values for μ, σ are used [91]). It is noted that in this work, Cartesian coordinates are used and (x,y,z) converts to (X,Y,Z) assuming a profile change with increasing NP. In Equation (9), superscript T denotes the transpose of the vector $\overrightarrow{\nabla }\overrightarrow{u}$.

A usual approach followed to simulate the multiscale dynamics and solve Equations (1)–(9) is the employment of a finite difference method on a staggered grid which represents an effective approach towards suppressing numerical oscillations. More specifically, T_c, T_L, N_c, and P are computed at the centre of each element; on the other hand, the time derivatives of the displacements and first-order spatial derivative terms are evaluated at locations midway between consecutive grid points. A common technique to solve the NSE is the projection method and the velocity and pressure fields are calculated on a staggered grid using fully implicit formulations [92,93]. On the other hand, the horizontal and vertical velocities are defined in the centres of the horizontal and vertical cells faces, respectively. It is noted that a multiple pulse irradiation scheme is required to derive the formation of a periodic relief [9]. More specifically, LIPSS are formed in the following steps:

• The first pulse irradiates a flat surface which leads to the formation of a crater and small protrusions (humps) at the edges on the surface of the heated zone due to mass displacement [22,33]. Moreover, due to the high fluence value, some ablation also occurs. The first pulse irradiates a flat surface with no corrugations, therefore the formation of periodic structures is not expected to happen. It is noted that due to the axial symmetry of a Gaussian beam, for NP = 1, Equations (1)–(9) can be solved in 2D.
• The second pulse, then, irradiates the attained pattern and therefore the spatial symmetry breaks; as a result, 2D modelling can no longer be used. The coupling of the electric field of the incident beam with the induced surface-scattered wave produces a nonuniform, periodic distribution of the absorbed energy. The periodic variation of the absorbed energy, in turn, leads to a periodic excited electron density distribution [9]. It is noted, however, that the computation of the amount of the absorbed energy at each position requires the evaluation of the energy deposition on a curved surface (i.e., Equation (7) for reflectivity is valid for flat profiles). Therefore, appropriate computational schemes are used to compute the absorbed energy on each point of the curved surface [9]. The calculated spatially modulated electron energy distribution is transferred to the lattice system (through the second equation of Equation (1)) and subsequently, upon phase transition fluid transport and resolidification processes, LIPSS are formed.
• The above methodology is used to describe the formation of LIPSS for N$P\ge$2 (including a correction to the surface plasmon wavelength shift to smaller values with increasing depth of the profile following an increase in dose [13,32,39]); however, there is a resonance at which further excitation of surface plasmons stops being the driving force behind the induced surface profile and suprawavelength structures are produced. In a previous report, it was shown that if the surface profile becomes sufficiently deep (at large NP) normal thermocapillary waves which lead to regular LIPSS are not produced [8,23]. By contrast, another solution of NSE dominates, namely, hydrothermal waves that propagate between the wells of the ripples in a perpendicular direction to the laser beam polarisation [8]. Another important feature of these solutions is that, only, waves of a certain periodicity (i.e., larger than the laser wavelength) lead to stable structures upon solidification, which are orientated perpendicularly to the beam polarisation and they are termed grooves.

3. Experimental Protocol

The experiments were performed at the Ultraviolet Laser Facility operating at IESL-FORTH (Heraklion, Greece). The experimental installation, samples and irradiation conditions have been described in another of our papers [46]. In brief, (001) and (100) oriented polished surfaces of rutile TiO₂ single crystals (Crystal-GmbH), maintained at room temperature, were irradiated by a hybrid distributed feedback dye laser/KrF excimer system delivering linear polarised 450 fs pulses at 248 nm at repetition rate varied of 10 Hz. The double-pulse irradiation was performed with a temporal delay between pulses varying between 0.5 and 450 ps using a Michelson-type interferometer. After appropriate spatial filtering and focusing, the laser beam on the sample was of quasi-Gaussian of size 19 μm measured at 1/e² of the maximum intensity. The surface morphology after the irradiation was examined by Jeol Scanning Microscope JSM 7000F, which employs a Schottky type field-emission (T-FE) gun for the electron source operating at pressure 10⁻⁸ Pa, with a function for high-resolution image observation. The spatial resolution achieved was 1.2 nm at 30 kV (max accelerating voltage).

4. Discussion

To solve the Equations (1)–(9), a numerical scheme based on the use of a finite difference method is followed; the discretization of time and space has been chosen, appropriately, to satisfy the Neumann stability criterion. Von Neumann boundary conditions are considered and heat losses at the front and back surfaces of the material are negligible. The initial conditions are T_c(t = 0) = T_L(t = 0) = 300 K, and N_c = 10¹² cm⁻³ at t = 0. On the other hand, for NSE (Equation (9)), it is assumed that the velocities of the fluid at t = 0 are zero. The values of the laser beam used in the simulation are: the (peak) fluence of $E_p\left(\equiv \sqrt{\pi }{\tau }_p{I}_0/\left(2\sqrt{ln2}\right)\right)$, where I₀ stands for the peak value of the intensity while R₀ in Equation (5) is taken to be equal to 20 μm and τ_p = 450 fs. Double pulses of equal constituent pulses are used with temporal separation in the range [0 ps, 20 ps]. The ultrafast dynamics and the thermal response of the heated material are investigated to take into consideration the role of fluence, pulse separation, and crystal orientation.

4.1. Single Pulse Excitation (t_delay = 0)

Firstly, to highlight the role of crystal orientation, single pulses are used. More specifically, double pulses of zero delay (t_delay = 0) with no interference effect are considered and the ultrafast dynamics are investigated. Two sets of the equal laser pulses with E_p = 50 mJ/cm² and E_p = 125 mJ/cm², respectively, lead to a total pulse fluence on the solid sample of 100 mJ/cm² and 250 mJ/cm², respectively. The evolution of the carrier density, carrier temperature and the lattice temperatures on the surface of the material (at the centre of the Gaussian pulse where the energy absorption is maximal) are illustrated in Figure 2a–c, respectively, for the aforementioned fluences and (100) and (001) crystal orientation. Results illustrate that the maximum carrier densities are higher for (100) than for (001) which highlights the important role of the value of the optical mass. A similar monotonicity is followed by the maximum carrier temperature and simulations show that the maximum value is acquired shortly after the time the carrier density reaches the peak value which agrees with results in previous reports for other semiconductors [26,56,57,94,95]. We note there is an initial rise of the carrier temperature (Figure 2b) during the first moments of irradiation due to the significantly large heat capacity of the carrier system with respect to that of the lattice.

The main processes of energy increase from the electron system are through a (linear) one-photon and free-carrier energy absorption. After a short period, T_c rises rapidly with an increase in the pulse energy. The initial increase can be explained through the examination of the contribution of the competing mechanisms indicated by the various components in the “source term” (see an enlarged region in Figure A1 in Appendix A). Similarly, this behaviour is followed by a small decrease. At large fluences, larger amounts of carriers are produced and (given the fact that one-photon absorption and Auger recombination are the two predominant factors that alter carrier density) Auger recombination becomes significant because it varies as $N_c^3$. Hence, the enhanced Auger recombination [94] converts carrier ionization energy into kinetic energy that results in an increase in the carrier temperature. It is noted that T_c attains very high values and this could potentially influence the energy band-gap response. Certainly, a more rigorous approach (e.g., based on first principles [96]) should also be considered to estimate a possible E_g variation for large T_c. Nevertheless, without loss of generality, it is assumed that substantial qualitative and quantitative changes (that yields different morphological changes) do not occur if E_g is considered constant. Results indicate that differences of the maximum values of T_c and N_c for (001) and (100) at E_p = 125 mJ/cm² are not very distinct compared to those at lower fluences due to the fact that this energy value rapidly yields maximum value of N_c close to N_v (thus, the upper value of available carriers for excitation is nearly reached). Furthermore, it is noted that although longer wavelengths and other semiconductors both impact ionization and two-photon absorption play an important role in the carrier dynamics and lattice response, for TiO₂ and Λ_L = 248 nm (that yield energies well above E_g) the latter processes appear to yield insignificant changes in N_c and T_c (see Appendix B). Nevertheless, for the sake of completeness and aiming to present a general theoretical framework that describes processes that also take place in other semiconductors, even at higher photon energies, a full model is considered assuming a nonzero ${\beta }_{TPA}$ and θ (Equation (3)).

Simulation results that represent the maximum surface lattice temperature also show that T_L evolution rises faster for (001) and therefore, the crystal orientation plays an important role (Figure 2c). To avoid confusion, it is noted that the maximum value of T_L (~3245 K) occurs at different depths of the irradiated material because ablation is more intense at higher fluences and therefore a larger volume is ablated. More specifically, the size of the ablated volume should be larger at higher fluences as the absorbed energy is distributed at larger depths. Nevertheless, the approach followed in the current analysis takes into account the fact that the temperatures attained are high enough to induce evaporation of the material and therefore, to avoid an overestimation of the material surface temperature in the calculations, a correction to the simulation scheme is required. Hence, becausethe boiling temperature is assumed to be the condition for evaporation, in the simulations, all points which exceed that temperature are removed. Thus, the region which is simulated at each timepoint varies dynamically and upon removal of all points with $T_L>T_L^{boiling}$, the depth of the induced morphology increases. It is noted that the constant temperature ($T_L = T_L^{boiling}$) in Figure 1c is due to the fact that energy relaxation between the lattice and carrier subsystems is not been reached for some picoseconds; this indicates that until the moment of equilibration between $T_L$ and $T_c$, the carrier subsystem will continue to give energy to the lattice system which might rise above the $T_L^{boiling}$ even if the surface points dynamically change. Hence, the surface continues to reach the boiling temperature and it evaporates leading to ablation while this elevating high temperature remains constant until relaxation processes start to become stronger. By contrast, a cutoff in the temporal decay of the carrier density or the reflectivity is not observed in the simulations. This results from the fact that the carrier density drops slowly after the end of the pulse due to an Auger recombination. Other processes that influence the variation of carrier density (see third equation in Equation (1)) such as the single photon absorption process that is dependent on the lattice temperature has ended before the lattice reaches $T_L^{boiling}$. Similarly, for the same E_p, the enhanced T_c for (001) yields a larger T_L which, in turn, implies a larger ablation compared to the mass removal for (100).

Special attention is also required to the behaviour of the reflectivity evolution. In previous studies, it was reported that at longer wavelengths and other semiconductors (i.e., silicon), the reflectivity initially drops below the value at 300 K before it starts rising towards reaching a peak as the carrier density increases [26,66,94]. Then, as the Auger recombination becomes important and the pulse switches off, the number of carrier drops which forces the reflectivity to fall before relaxing to the initial value. A similar behaviour is also present for both orientations of rutile TiO₂ (Figure 2d)_. Notably, although the shape of the reflectivity curves is similar for both crystal orientations, the different value of the initial reflectivity for (100) and (001) causes variation in the energy absorption. On the other hand, for the same crystal orientation, irradiation with higher fluences increases the metallisation process of the material. This leads to a rapidly enhanced carrier density followed by a larger peak value of the reflectivity of the excited material. The latter is attributed the generation of an increased number of carriers at higher fluences (Figure 2a).

4.2. Double Pulse Excitation (t_delay ≠ 0)

To investigate the influence of irradiation with double pulses on the ultrafast dynamics Equations (1)–(7) are solved for pulse delays in the range [0,−20 ps]. In Figure 3, Figure 4 and Figure 5, the surface carrier density, reflectivity and carrier temperature evolution, respectively, are illustrated for three representative values of t_delay = 0.5 ps, 0.9 ps, and 1.5 ps for E_p = 50 mJ/cm² and 125 mJ/cm² and both crystal orientations. Results are shown for a temporal delay between pulses smaller than 2τ_p (0.9 ps) where the two pulses overlap as well as for pulses with a larger separation (1.5 ps). According to simulations (Figure 6b,c), the maximum surface carrier density and temperature decrease gradually for short delays (t_delay < 2τ_p), reaching a minimum value that remains unchanged even at longer t_delay. To explain this behaviour, it is important to note that the first of the constituent pulses leads to ablation; hence, the enhanced surface lattice temperature (~$T_L^{boiling}$) of the material yields similar ${\alpha }_{SPA}$ (i.e., $\sim {10}^6 \times \exp \left(T_L^{boiling}/300\right)$ cm⁻¹) before the second pulse irradiates the material for t_delay > ~2τ_p. Figure 6a, illustrates that $T_L~T_L^{boiling}$ at, approximately, t_delay = 1 ps which justifies the appearance of the minimum $T_L$ at t_delay > 2τ_p. This behaviour is different from the one exhibited in other laser conditions (i.e., nonablation conditions at Λ_L = 800 nm) for which the second pulse leads to higher carrier densities and temperatures [36,97]. Due to the substantially larger contribution of one-photon based excitation at low wavelengths and the negligent contribution of impact ionization and two-photon ionization, the fact that ${\alpha }_{SPA}$ remains constant at t_delay > ~2τ_p leads to unchanged values for T_c and N_c.

The above investigation is important because it also reveals both the ultrafast dynamics and temperature evolution inside the volume of the material in ablation conditions. A detailed picture of the induced thermal effects below the surface of the material will help to model the morphological changes and surface patterning procedure. However, results show that T_c, T_L and N_c can provide a detailed description of the laser energy distribution and thermal effects on the affected layer of the material, exploration of dynamics through analysis of the thermal changes on the surface (in ablation conditions at E_p = 50 mJ/cm² and 125 mJ/cm²) is hindered by the fact that the material remains at the isothermal $T_L = T_L^{boiling}$ for a long time and therefore, important details of the thermal effects below the surface of the material are not easily manifested. It is evident that one way to illustrate how pulse separation influences the energy absorption and thermal effects at larger depths in ablation conditions is by presenting the affected zone characteristics (i.e., temperature distribution, volume of affected zone, position of surface after ablation). By contrast, to show the impact of the laser parameters on the induced thermal effects below the ablated region (and at depths characterized by lower energy distribution), an alternative investigation is performed; more specifically, simulations are conducted assuming excitation of the material at lower fluences (i.e., that do not cause ablation; this is equivalent to exploration of dynamical effects in places where material is excited but is not removed). To emphasise on the significant role of the one-photon excitation and at conditions that do not induce ablation (i.e., $T_L<T_L^{boiling}$) before the second pulse irradiates the material, simulations have been performed for irradiation with E_p = 10 mJ/cm² for t_delay = 0.7 ps, 1.5 ps, 4 ps and 10 ps (Figure 7). Interestingly, the carrier density is higher for (001) than for (100) (Figure 7a–d) unlike the behaviour for ablation conditions. This can be attributed to the almost insignificant variation of reflectivity at low fluences and the fact that energy absorption (due to lower reflectivity) is higher for (001) than for (100) at low carrier densities (Figure 4). Results for the carrier density and lattice temperature will help to understand the thermal response of the material at larger depths when high temperatures (but not equal to the boiling point are generally reached). Another interesting issue that is revealed by comparing the carrier density evolution in the two regimes (i.e., ablative and nonablative regimes) is the drop of N_c as a function of time. In contrast to the pronounced smoother gradual decrease in the carrier density value at lower fluences (part of the curve in Figure 7b–d due to the first pulse) for which the maximum attained value for N_c is low, there is a substantially steeper decrease when a substantially peak value of N_c has been reached (Figure 2a, Figure 3 and Figure 7a–d (after exposure to the second pulse)). It is noted that the smooth decrease in N_c for small carrier densities has also been reported in previous works [26,56,57,94,95].

On the other hand, given that the single photon absorption coefficient is essentially unchanged within the first of the double pulse sequence (because T_L has not started to increase), any change in the carrier density is predominantly caused from the contribution of I and secondly to the Auger recombination. At low fluences (~10 mJ/cm²), no sharp increase in the carrier density is computed because it occurs at higher fluences (Figure 2a and Figure 3). By contrast, when the second constituent pulse irradiates the material, the energy of the lattice system is sufficiently high (Figure 7e–h) to lead to high levels of excitation as ${\alpha }_{SPA}$ increases rapidly. Thus, a large number of excited carriers is produced.

Results show that for subablation conditions, the maximum T_L increases as the temporal delay becomes larger. This can be attributed to the enhanced one-photon absorption (${\alpha }_{SPA}$ attains large values during the first pulse, but it is about two orders higher when the second pulse irradiates the material due to its T_L dependence (Figure 7e–h)) and free carrier absorption while the influence of impact ionisation appears to be minimal. According to Figure 7, the second pulse irradiates an already hot material (at increasing T_L) that further increases the single photon absorption that leads to a higher carrier density which affects eventually T_L through electron-phonon scattering. This is also manifested at a higher t_delay where ablation starts to occur (Figure 7h). It is noted that the t_delay value (Figure 7h) for which the temperature reaches ablation conditions, there is large jump of the carrier density to the increase in the energy absorption and production of a large number of energetic carriers. It is noted that at these fluences, reflectivity variation during the pulse is insignificant and it does not induce changes in energy absorption. To summarise, Figure 7 illustrates that in subablation conditions, an increase in the maximum values of T_L occurs and due to heat diffusion and heat transfer it is more enhanced at higher depths inside the material volume as t_delay increases. These results indicate that in regions where the energy deposition is not high enough to induce material removal (at higher depths inside the material) a larger volume of the material will be affected with increasing t_delay. This conclusion is of paramount importance to interpret the distinct surface patterning at longer delays (i.e., ripples and grooves).

4.3. LIPSS Formation

The above investigation provided a detailed description of the ultrafast dynamics after the heating of rutile TiO₂ of different crystal orientation with ultrashort double laser pulses.

As noted in the introductory section, a question of paramount importance is how the pulse separation and the consideration of irradiation of rutile TiO₂ (001) and (100) influences the surface morphological features (i.e., types of structures and periodicities). The images in Figure 8 show the TiO₂ (001) surface after irradiation with two pulses of (peak) fluence (for each pulse) E_exp = 125 mJ/cm² and different temporal delays between laser pulses. According to results in a previous report [46], this fluence corresponds to the LSFL domain of morphological surface maps (MSM). In agreement, ripples were formed at a small delay of t_delay ≤ 4 ps (Figure 8a,b). However, at a longer delay of ~10 ps (Figure 8c), suprawavelength grooves were formed; they coarsened at even longer delay of ~30 ps (Figure 8d) and disappeared at ~60 ps, being replaced by spike-like columns (not shown).

To relate experimental data with theoretical results, it is important to estimate the dependence of the produced carrier densities with the induced surface structures as an increase in the irradiation dose (i.e., increase in NP) modifies the energy absorption and characteristics of the morphological features. For the investigation of the LIPSS formation, between the two proposed mechanisms stated in the previous section based on Sipe’s theory and the excitation of SPP, the former suffers from fundamental inconsistencies: more specifically, it is already known that one of the limitations of the efficacy factor-based theory is the neglect of the so-called ‘feedback mechanism’ which is very important to calculate the evolution of the periodicity of the induced periodic structures [34,40].

On the other hand, the SPP-based approach is used to compute the periodicities considering the produced carrier density values following the surface morphology that is induced for every pulse. In contrast to more precise electrodynamics simulations based on Finite Difference Finite Domain Schemes (FDTD) used in previous reports to correlate the induced LIPSS periodicities with increasing NP [21,24,39,84,98], an alternative and approximating methodology was employed, in this work, to relate the SP wavelength with the produced maximum depth of the corrugated profile [78,79] (i.e., which is directly linked with NP). The methodology is based on the spatial distribution of the electric field on a corrugated surface of particular periodicity and height and how continuity of the electromagnetic fields influences the features of the associated SPP. Results for SPP wavelength as a function of NP are illustrated in Figure 9 for E_p = 125 mJ/cm² for (100) and (001) TiO₂, respectively (t_delay = 4 ps). As mentioned in the previous section, the decrease in SPP wavelength with the increase in the irradiation dose NP is related to the expected shift of SPP wavelength to smaller values as the profile becomes deeper [21,78]. Furthermore, according to Figure 7a–h, the larger carrier densities (and thermal response) produced for the orientation (001) lead to a more rapidly attained deeper profile with respect to that for (100). Hence, the curve that shows the decrease in the ripple periodicity as a function of NP should be steeper for (001) while SPP wavelengths also are expected to expand to a larger range of values (Figure 9). Furthermore, simulation results predict a substantially smaller SPP wavelength for (001) orientation which should be also projected on the periodicity of the induced ripples.

Apart from the comparison of the LIPSS periodicities following irradiation with the same fluence, pulse separation (t_delay = 4 ps), and NP, it is also very important to evaluate potential morphological variations for different pulse separation. In Figure 10, simulations and experimental results are shown after a double pulse of E_p = 125 mJ/cm² irradiated (001) TiO₂ for two delays, t_delay = 4 ps and t_delay = 10 ps. An interesting point regarding the two values of the pulse separation is that they are situated before and past an optimum pulse separation value (t_delay~6 ps) for which the depth is critical to lead to a transition of structures from ripples (i.e., subwavelength structures oriented perpendicularly to the laser polarization direction) to grooves (i.e., suprawavelength structures that are aligned with laser polarization direction); therefore it is important to evaluate the impact of small and large pulse separation. Due to the fact that for a smaller pulse separation (t_delay = 4 ps), damage related temperatures will be developed into smaller regions, a shallower profile will be produced compared to that for t_delay = 10 ps. Hence, the condition for promoting the production of grooves will be satisfied more easily at longer delays. Simulation results show rippled structures (Figure 10a) for NP = 80 for the above conditions at t_delay = 4 ps of periodicity ~220 nm and grooves (Figure 10b) for NP = 80 at t_delay = 10 ps of periodicity ~1.8 μm. Experimental results on (001) TiO₂ rutile surface that validate the simulations are illustrated in Figure 8b,c, which show the production of ripples and grooves, respectively, for NP = 100. Interestingly, irradiation of (100) rutile TiO₂ surface does not lead to the formation of grooves [46]. Instead of grooves, a featureless flat area replaces the ripples when the delay between laser pluses increases (see Appendix C). The formation of this unusual structure is not yet explained by simulations. This emphasises that crystalline orientation of the irradiated surfaces is a key parameter in the formation processes of suprawavelength periodic structures under femtosecond laser irradiation.

The multiscale model used in this work described the physical processes that account for the formation of LSFL on TiO₂ while a transition from ripples to grooves was also derived. It is evident that more accurate conclusions will be drawn if more appropriately developed experimental (for example time-resolved) protocols are also introduced to evaluate ultrafast dynamics and relevant parameter valuessuch as electron-electron and electron-hole scattering frequencies, electron-phonon relaxation times, and reflectivity changes. It is noted that although the theoretical predictions appear to agree with the evolution of the surface structures, simulation results manifested a significant ablation that takes place for the fluence values used in the experiment which to some extent deviates from experimental observations. This can also be attributed to the choice of various parameters such as the electron effective mass that affects the dielectric constant and eventually the optical properties of the material and energy absorption; this is a critical issue to both the induced carrier density and thermal response of the material. Values for $m_{opt}^{*}$ used in the current work and at low wavelengths might need to be revised appropriately through a more rigorous theoretical approach and suitable experimental protocols. Similarly, in extreme laser conditions that induce severe ablation, it is very important to estimate more accurately the role of evaporation and employ more advanced theoretical schemes such as FDTD algorithms to evaluate electrodynamical effects that are related to the underlying mechanisms of surface modification [24]. Nevertheless, the predominant aim of the study focused on the exploration of the underlying physical mechanisms that correlate laser parameters with the material’s crystal orientation and the induced ultrafast dynamics and produced surface patterns. The qualitative and quantitative results presented in the current investigation were aimed to provide a first insight into the patterning techniques of rutile TiO₂ with femtosecond pulses in ablation conditions and highlight the significant role of the material characteristics of the material related to the crystal orientation in the response of the solid. The detailed analysis emphasized the role of the different direction dependence of the electron effective masses for the two crystal orientations as the distinct two cases lead to different excitation processes, energy absorption and anisotropic laser surface patterning.

5. Conclusions

A detailed theoretical framework was presented that describes both the ultrafast dynamics and surface modification physical mechanisms after the heating of rutile TiO₂ of different crystal orientations with ultrashort double laser pulses in ablation conditions. Results manifest that the crystal orientations and interpulse delays play an important role in the onset of surface pattern formation because they influence both the carrier dynamics and thermal response of the irradiated structure. Simulations showed that although both the maximum carrier temperature and density drop with increasing pulse separation, these parameters reach a minimum value at about t_delay > 2τ_p after which no variation occurs. Interestingly, irradiation of TiO₂ with (100) orientation lead to higher T_c and N_c values than for (001) which are well explained by the underlying physical processes and the excitation conditions. To illustrate the response of the material deeper in the surface of the irradiated solid, simulations indicated an increase in the maximum values of T_L due to heat diffusion and heat transfer that is more enhanced at higher depths inside the material volume as t_delay increases. These results are very important because they are also associated with the induced surface patterning. Furthermore, surface modification is greatly affected by the pulse separation and both simulation results and experimental observations explain the occurrence of grooves on (001) compared to the formation of subwavelength structures at smaller delays. For values of the pulse separation past an threshold value (~ 6 ps), the depth of the surface pattern allow the formation of grooves instead of ripples, which indicates that the modulation of the frequencies of the induced LIPSS on the surfaces of solids can be tailored by a control of t_delay and the number of pulses. Furthermore, simulation results indicate that the crystal orientation influences the periodicity of the LSFL structures, leading to patterns of lower periodicity for (001) TiO₂.

The approach is aimed to provide a systematic laser-based processing strategy of materials and tailor the morphology of an irradiated surface according to the demand of exciting applications ranging from biomedical engineering to photovoltaics and nanoelectronics.