Ultrashort pulsed laser induced complex surface structures generated by tailoring the melt hydrodynamics

Fotis Fraggelakis; George D. Tsibidis; Emmanuel Stratakis

doi:10.29026/oea.2022.210052

Abstract

We present a novel approach for tailoring the laser induced surface topography upon femtosecond (fs) pulsed laser irradiation. The method employs spatially controlled double fs laser pulses to actively regulate the hydrodynamic microfluidic motion of the melted layer that gives rise to the structures formation. The pulse train used, in particular, consists of a previously unexplored spatiotemporal intensity combination including one pulse with Gaussian and another with periodically modulated intensity distribution created by Direct Laser Interference Patterning (DLIP). The interpulse delay is appropriately chosen to reveal the contribution of the microfluidic melt flow, while it is found that the sequence of the Gaussian and DLIP pulses remarkably influences the surface profile attained. Results also demonstrate that both the spatial intensity of the double pulse and the effective number of pulses per irradiation spot can further be modulated to control the formation of complex surface morphologies. The underlying physical processes behind the complex patterns’ generation were interpreted in terms of a multiscale model combining electron excitation with melt hydrodynamics. We believe that this work can constitute a significant step forward towards producing laser induced surface structures on demand by tailoring the melt microfluidic phenomena.

Introduction

During the past decades, materials’ micro and nano fabrication using ultrashort laser pulses have emerged as a key technology, which has contributed to major advances in science, technology and industry^1-7. Compared to expensive and time consuming techniques, such as photo- or electron beam-lithography^{8, 9} that do not provide the necessary specificity and scalability for industrial applications, laser-based processing constitutes an efficient and precise technology to fabricate novel functional materials and devices^10-12.

Laser surface processing is a rapid and scalable chemical-free technique for surface functionalization that can be applied on almost any kind of material. Laser Induced Periodic Surface Structures (LIPSS) is the most common type of morphology that can be produced via irradiation of solids’ surfaces by Gaussian beams^{2, 13, 14}. Depending on the laser wavelength, λ_L, used to fabricate LIPSS, their period Λ varies from (i) deep/shallow subwavelength (Λ<<λ_L/2), termed as High Spatial Frequency LIPSS (HSFL)^{13, 15} to (ii) subwavelength (λ_L/2 < Λ< ~λ_L), coined as Low Spatial Frequency LIPSS (LSFL)^{16, 17}, to (iii) suprawavelength (Λ>λ_L), called grooves^18-21 or spikes^{20, 22}, as well as other complex structures^{23, 24}. The formation mechanisms of both the LSFL and HSFL structures is attributed to a spatially modulated deposition of the laser energy as a result of the interference of the incident beam and electromagnetic waves scattered on a rough surface in the far- and near-field respectively^{2, 17, 25, 26}.

Apart from LIPSS formation, an alternative laser processing method, namely the Direct Laser Interference Patterning (DLIP), has recently been introduced to fabricate periodic structures on solid surfaces. DLIP is based on the use of two or more coherent laser beams that interfere and the resulting laser intensity profile irradiates the material to fabricate periodic structures on large areas^27-30. The imprinted pattern is determined by the angle, the number and the coherence of the interfering beams^{31, 32}. The capability to tailor the features of a surface topography by controlling the angle of incidence of the coherent beams offers the capability to form a broad range of surface structures for potential applications^27-30. Nonetheless, despite the ongoing research and investigation of the various surface patterns, a quantitative analysis of the underlying physical mechanisms behind the DLIP structures’ formation is still missing. In some recent works, the role of the induced thermal effects has been investigated^{33, 34}. Recently, we have introduced a multiscale theoretical model taking into account the impact of the microfluidic factor as well to interpret the DLIP pattern formation upon double pulse irradiation (DPI)³⁵. Theoretical simulations of the periodic patterns formed by two- and four-beam DLIP double pulses, confirmed by experiments, demonstrate the significant role of melt hydrodynamics in the surface structuring process³⁵.

Despite the significant progress to date, the desired level of control of laser induced structures’ characteristics including symmetry, size and hierarchical features formation is far from being accomplished. LSFL and HSFL with multiple axes of symmetry have been only reported as a result of polarisation modification or DPI^{22, 24}. Matching the extraordinary variety of structures and functionalities found in nature requires an in-depth understanding of the formation mechanism along with developing the suitable irradiation techniques.

Given the impact of the intensity profile and polarisation in the produced structures, it would be crucial to explore whether appropriate control of the spatiotemporal energy deposition profile would lead to more complex surface morphologies with enhanced application-based functionalities. To the best of our knowledge, no previous study has been conducted on the investigation of the type of topographies that can be produced if double pulses of different spatial intensity are used. In principle, double pulses of Gaussian^{36, 37} or DLIP³⁵ shapes have been previously reported, however the combined effect of Gaussian and DLIP profiles should be important to explore. As pointed out above, the Gaussian beam leads to the formation of LSFL structures, while DLIP yields morphologies that depend on the angle between the coherent beams and the polarisation. Therefore, a thorough investigation of the combined effect of such pulses and inter-pulse delay on surface structures’ formation would be important to study. Of further importance is the investigation of the resultant microfluidic movement derived from the phase transition caused under laser irradiation. The application of delayed pulses is crucial to reveal the impact of melt flow, considering that sufficient time is required for hydrodynamical effects to dominate, prior to the action of the second pulse that arrives at the materials’ surface and influences the melt flow. For this purpose, the pulse separation here is kept on the order of ~500 ps, which has been proven to be suitable to highlight the role of the melt microfluidic movement³⁸.

Experimental and simulation protocols

Laser processing

In the experiments, we utilize a DPI irradiation scheme employing two temporally delayed fs pulses with different intensity distribution. Namely, a pulse with Gaussian distribution (Fig. 1(d), indicated as G) is combined with a pulse with a 1D or 2D intensity distribution produced by DLIP (Fig. 1(d), indicated as V, H or D). For the experiments, fs pulses, emitted by a Pharos (Light Conversion) laser source with duration of τ_p 170 fs at a wavelength λ_L = 1026 nm have been used. As shown in Fig. 1 a delay line is built to generate double pulses. In particular, the primary pulse is initially divided into two pulses with perpendicular polarisation orientation by a polarizing beam splitter. The energy distribution between the pulses is regulated by a half waveplate (HWP) placed before the polarizing beam splitter (PB). The two generated pulses pass through the arms A and B corresponding to optical paths that differ by 150 ± 1 mm acquiring an interpulse delay of Δτ = 500 ± 3 ps. The two pulses are guided to a computer controlled programmable Spatial Light Modulator (SLM) and via a suitable optical system, comprising two lenses, and interfere onto the sample. The two focusing lenses f₁ = 400 mm and f₂ = 30 mm are placed on the appropriate distances recombining the beams on the sample giving a spot diameter of 190 μm radius (2ω). The working principle of the SLM screen is described in Fig. 1(b). Although the s-polarized pulse is reflected from the SLM and maintains its Gaussian shape (Fig. 1, G), the p-polarized one interacts with the SLM screen and is divided into two or four beams acquiring on the sample a profile comprising vertical (Fig. 1, V), or horizontal (Fig. 1, H) lines, or a periodic array of dots (Fig. 1, D). The periods of DLIP patterns used, Λ_DLIP, were comparable to the laser wavelength. Owing to a rotating HWP placed before the DLIP part, the initial order of the s- and p- pulse can be inverted enabling one to control whether the Gaussian or the DLIP pulse reaches the sample first. The total energy, E_tot, for the pair of pulses, the total number of irradiation shots, NP, as well as the order of the pulses was varied, as indicated in the experimental section. The energy distribution was 2/3 of E_tot for the Gaussian and 1/3 E_tot for the DLIP pulse respectively, tuned to maintain a similar effective fluence in the irradiated areas. NP was varied in the range of 10 and 500 pulses. Images of the processed surfaces were acquired via Scanning Electron Microscopy (SEM) microscopy, while a Fast Fourier Transformation (FFT) was performed to calculate the periodicities of the induced structures, through the use of the open source software Gwyddion. For all the experiments, mirror polished commercially available 316-stainless steel (provided by RS) substrates have been used.

Figure 1.(a) Experimental setup. Abbreviations: half-wave plate (HWP), linear polarizing cube (LPC), beam splitter (BS), spatial light modulator (SLM), focusing lenses (f₁,f₂). (b) SLM function. (c) Spot profile distribution at the sample: Gaussian (G) and DLIP (V, H, D) profiles. (d) The red arrow indicates the polarisation vector.

Theoretical model

DLIP Intensity profiles

DLIP is based on the superposition of the electric field vectors of coherent beams. The electric field of a light wave n is equal to ${\overset{⇀}{E}}_{n} = A_{n} {\overset{⇀}{p}}_{n} e^{i (ω t - {\overset{⇀}{k}}_{n} \cdot \overset{⇀}{r})}$ where $A_{n}$ is the amplitude, ${\overset{⇀}{p}}_{n}$ is the polarisation vector, ${\overset{⇀}{k}}_{n}$ (i.e. $| {\overset{⇀}{k}}_{n} | \equiv k = \frac{2 π}{λ_{L}}$ ) stands for the wavevector, $\overset{⇀}{r}$ is the position vector, ω is the angular frequency, and t is the time. It is noted that the resulting interference pattern for N beams is calculated through the expression ${\overset{⇀}{E}}_{total} = \sum_{n = 1}^{N} {\overset{⇀}{E}}_{n}$ . Each laser beam irradiates the material at an incident angle with the vertical axis equal to $θ_{n}$ and azimuthal angle $φ_{n}$ . It is assumed that all laser beams are linearly polarised and the polarisation direction is the same for all DLIP beams (i.e. perpendicular to that of the Gaussian). The total spatial intensity distribution is, then, provided by the formula $I_{0} = c ε_{0} {| {\overset{⇀}{E}}_{total} |}^{2}$ /2, where c and $ε_{0}$ stand for the speed of light and dielectric vacuum permittivity, respectively. For interference with two ( $φ_{1} = φ_{2} = 0$ , and $θ_{1} = θ_{2} = θ$ ) and four beams ( $φ_{1} = φ_{2} = 0$ , $φ_{3} = φ_{4}$ = π/2, $θ_{1} = θ_{2} = θ_{3} = θ_{4} = θ$ ), the following total intensities are $I_{0}^{(2)}$ and $I_{0}^{(4)}$ on the surface of the material, respectively (by calculating the Poynting vector and by taking the average value $〈 I_{0} 〉$ over the laser period $λ_{L} / c$ ).

1 $\begin{array}{l} I_{0}^{(2)} = I_{1} [1 + cos (2 k x sin θ)], \\ I_{0}^{(4)} = I_{1} {[1 + \frac{1}{2} cos (2 k x sin θ) + \frac{1}{2} cos (2 k y sin θ) \\ + 2 cos (k x sin θ) cos (k y sin θ)]}, \end{array}$

where I₁ is the intensity of each of the constituent laser beams of the DLIP in which a Gaussian envelope is included (~ $e^{- 2 (\frac{x^{2} + y^{2}}{R_{0}^{2}})}$ ) with an e^–2 Gaussian spot radius equal to R₀; similarly, I₁ includes the Gaussian distribution in time, $e^{- 4 ln 2 (\frac{{(t - 3 τ_{p})}^{2}}{τ_{p}^{2}})}$ . It is noted that a Cartesian coordinate system is used (x,y,z). It is evident that the choice of $θ_{n}$ can be used to define the periodicities of the interference pattern. More specifically, for two-beam DLIP, a sinusoidal (of periodicity equal to $Λ_{DLIP} = \frac{λ_{L}}{2 sin θ}$ along x-axis) or for a four-beam DLIP, a dot-type intensity distribution (of periodicity equal to $Λ_{DLIP} = \frac{λ_{L}}{sin θ}$ along x-axis and y-axis and a secondary $Λ_{DLIP} = \frac{1}{\sqrt{2}} \cdot \frac{λ_{L}}{sin θ}$ along the direction defined by $\hat{x}$ + $\hat{y}$ ( $\hat{x}, \hat{y}$ are unit vectors) are derived as illustrated in Fig. 1.

Modelling of laser-matter interaction

To simulate the surface modification processes and reveal the impact of excitation levels and hydrodynamical effects following irradiation with fs pulses, a multiscale description of the underlying physical mechanisms is required to account for the response of the material. Given that the intensity profile changes spatially with the combination of temporarily separated Gaussian and DLIP pulses, a detailed investigation of the excitation, relaxation processes and propagation of the produced hydrothermal waves can be performed through a theoretical framework that comprises modules that model the processes in various temporal regimes. A detailed description of the theoretical model that simulates the physical mechanisms that characterize irradiation of stainless steel with DLIP pulses is presented in ref.³⁵. In that report, the thermophysical properties and optical parameters of the material are also provided.

In summary, a three-dimensional Two Temperature Model (TTM) represents the standard theoretical framework to investigate laser-matter interaction and the energy transfer between the electron and lattice subsystems

2 $\begin{array}{l} C_{e} \frac{\partial T_{e}}{\partial t} = \overset{\leftarrow}{\nabla} (k_{e} \overset{\leftarrow}{\nabla} T_{e}) - g (T_{e} - T_{L}) + W, \\ C_{L} \frac{\partial T_{L}}{\partial t} = \overset{\leftarrow}{\nabla} (k_{L} \overset{\leftarrow}{\nabla} T_{L}) + g (T_{e} - T_{L}), \end{array}$

where $C_{e}$ and $C_{L}$ are the heat capacities of the electron and lattice subsystems, respectively while $T_{e}$ and $T_{L}$ are the temperatures of the two systems. By contrast, $k_{e}$ ( $k_{L} \sim 0.01 k_{e}$ ) stand for the electron (lattice) conductivity, $g$ is the electron-phonon coupling parameter while W corresponds to the absorbed laser power density that is related to the laser intensity. It is noted that the wavelength of all laser beams is equal to λ_L= 1026 nm. As the laser intensity comprises two (temporarily separated with a delay Δτ) pulses, W includes the variation of the total intensity of the laser beam (i.e. αΙ_total, where α is the absorption coefficient of the material at λ_L= 1026 nm). Due to the presence of two pulses, one Gaussian and one DLIP (with either two or four beams), the total laser intensity on the material surface is provided from the following expression (see derivation in the Supplementary information)

3 $\begin{array}{l} I_{total} (t, x, y, s u r f a c e) = \\ P_{1} e^{- 4 ln 2 {(\frac{t - 3 τ_{p} - G_{1} Δ τ}{τ_{p}})}^{2}} + P_{2} e^{- 4 ln 2 {(\frac{t - 3 τ_{p} - G_{2} Δ τ}{τ_{p}})}^{2}}, \end{array}$

where τ_p is the pulse duration of each laser beam and P₁ and P₂ are the two temporarily separated³⁹ pulses with a Δτ delay. With respect to the laser beams: (i) the Gaussian beam irradiates the material first if G₁= 1, G₂= 0, otherwise, G₁= 0, G₂= 1 if the DLIP irradiation comes first; (ii) P₁ (or P₂) is the (spatial) profile of the Gaussian (or the DLIP) beam. More specifically, for the DLIP beam, P₂= $I_{0}^{(2)}, I_{0}^{(4)}$ which depends on the number of beams that comprise the total beam. In the experiments performed in this work, the total laser energy E_tot that is distributed in the system is split into two parts, 2/3 E_tot are provided from the Gaussian beam while the rest is attributed to the DLIP beam.

It is noted that the values of the energy per pulse used in the experiments (i.e. 40 μJ, 60 μJ, 80 μJ) are found, by the simulations, resulting in mass removal (i.e. ablation)³⁵. To account for the experimental results and considering that mass removal is observed to be minimal in the experiments, all simulations have been carried out using lower energy values of the order of ~25 μJ.

It is emphasised that special attention is required to determine the impact of the (absorbed) laser energy given that the pulse separation in the double pulse experiment (Δτ =500 ps) indicates that the second pulse irradiates a material, which is in its molten phase. Results for the maximum lattice temperature $T_{L}^{max}$ evolution with time at various depths are illustrated in Fig. 2 following irradiation of a flat surface with E_tot= 25 μJ for G₁= 1, G₂= 0 (Fig. 2(a)) and G₁= 0, G₂= 1 (Fig. 2(b)); results show that this energy value leads to a phase change in the second case. Interestingly, when the Gaussian pulse irradiates the material first, the impact of the second pulse leads also to minimal mass removal that is simulated with the removal of material lattice points with temperatures higher than 0.9T_cr [T_cr stands for the thermodynamic critical temperature which for stainless steel is ~8500 Κ]^{35, 40, 41} (Fig. 2(c)). Although the volume of the ablated material is minimal at E_tot= 25 μJ, a substantially larger ablation is predicted at higher energies (results are not shown).

Figure 2. $T_{L}^{max}$ evolution at different depths for (a) G₁= 1, G₂= 0 and (b) G₁= 0, P₂=1, after irradiating a flat profile (NP= 1) with E= 25 μJ, Δτ = 500 ps. (c) $T_{L}^{max}$ on the sample surface as a function of time in the cases shown in (a) and (b).

Interestingly, both the values and evolution of $T_{L}^{max}$ are dependent on which pulse irradiates the material first, since the energies of the two pulses are different (Fig. 2(c)). Furthermore, the material response is related to what phase the material is in at the time of irradiation; more specifically, it is shown that if the second pulse irradiates a molten volume, a significant reflectivity variation³⁵ leads to a distinct thermal response of the material causing different: (i) absorption levels produced by the sequence of the double pulses and (ii) temperature gradients. The evolution of $T_{L}^{max}$ , illustrated in Fig. 2(c) is calculated for a flat surface. Although such $T_{L}^{max}$ evolutions illustrate the thermal response of the system for NP=1, similar conclusions are also deduced for non-flat profiles.

On the other hand, to interpret the induced surface modification, the development and evolution of hydrothermal waves and the dynamics of the produced fluid movement as a result for the phase transformation is mathematically described by the Navier-Stokes equations (NSE)³⁹

4 $\begin{array}{l} ρ_{0} (\frac{\partial \overset{\leftarrow}{u}}{\partial t} + \overset{\leftarrow}{u} \cdot \overset{\leftarrow}{\nabla} \overset{\leftarrow}{u}) = \\ \overset{\leftarrow}{\nabla} \cdot (- P 1 + μ (\overset{\leftarrow}{\nabla} \overset{\leftarrow}{u}) + μ {(\overset{\leftarrow}{\nabla} \overset{\leftarrow}{u})}^{T}), \end{array}$

where $μ$ and $ρ_{0}$ stand for the viscosity and density, respectively, of the molten (uncompressed) material, while P and $\overset{\leftarrow}{u}$ are the pressure and velocity of the fluid and superscript T denotes the transpose of the vector $\overset{\leftarrow}{\nabla} \overset{\leftarrow}{u}$ , and 1 is a 3×3 identity matrix. It is noted that velocity gradient tensor is a 3×3 matrix ${(\overset{\leftarrow}{\nabla} \overset{\leftarrow}{u})}_{i j} = \frac{\partial u_{i}}{\partial x_{j}}$ (see ref.³⁹ and Supplementary information) where (u,v,w) are the components of $\overset{\leftarrow}{u}$ in Cartesian coordinates (x, y, z). The solution of NSE is conducted through the employment of appropriate thermocapillary boundary conditions. More specifically, to describe the Marangoni flow, a shear stress balance on the free surface is described by the equation $\frac{\partial u_{τ}}{\partial n} = - \frac{1}{μ} \frac{\partial σ}{\partial T_{L}} \frac{\partial T_{L}}{\partial τ}$ at the liquid free surface where $u_{T}$ corresponds to the component of the velocity on the surface tangent direction τ while n indicates the component normal to the free surface. This expression turns into two equations for NP=1 (i.e. flat surface): $\frac{\partial u}{\partial z} = - \frac{1}{μ} \frac{\partial σ}{\partial T_{L}} \frac{\partial T_{L}}{\partial x}$ and $\frac{\partial v}{\partial z} = - \frac{1}{μ} \frac{\partial σ}{\partial T_{L}} \frac{\partial T_{L}}{\partial y}$ . In Eq. (4), σ stands for the surface tension of the material (see ref³⁵ for parameter values and for a more detailed description of the fluid dynamics module). As discussed in previous reports^{17, 35, 42}, inhomogeneous energy deposition leads to the development of Marangoni effects (i.e. surface tension-driven molten material flow) and displacement of material from regions of high to low temperatures. In addition, it is assumed that non-slipping conditions ( $\overset{\leftarrow}{u} = 0$ ) are applied on the interface between solid-and liquid interface. Finally, the molten material is considered to be an incompressible fluid ( $\overset{\leftarrow}{\nabla} \cdot \overset{\leftarrow}{u} = 0$ ). To simulate the fluid movement and dynamics of the molten material, the affected region has been divided in two subregions, one that contains a material in solid phase and another in liquid phase. The hydrodynamic equations are solved in both regions. To include the ‘hydrodynamic’ effect of the solid domain, material in the solid phase is modelled as an extremely viscous liquid (μ_solid=10⁵μ), which will result in velocity fields that are infinitesimally small. An apparent viscosity is then defined with a smooth switch function of Gaussian form to emulate the step of viscosity at the melting temperature¹⁷. A moving boundary for the two phases (solid-liquid interface) is based on position of the isothermal T_L=T_melting, where T_melting stands for the melting point of the material.

As the theoretical model aims to describe both the thermal response and dynamics of the molten material that leads eventually to morphological changes, topography variations should be considered due to consecutive irradiation as a result of increasing the energy dose. Therefore, apart from the aforementioned boundary conditions used to solve NSE (i.e. shear stress balance on a free but nonflat surface), appropriate conditions need to be introduced for the electron and lattice temperatures. More specifically, negligible heat loss from the free surface should be considered while it is also assumed that at large distances from the affected region both the electron and lattice temperatures are ~300 K. Special attention is required to evaluate the lattice temperature resulting from multiple irradiations. While the initial conditions for the velocity and lattice/electron temperatures before the delayed pulse irradiates the material are determined from the combined solution of NSE and TTM at the moment the second pulse heats the material, it is assumed that every consecutive double pulse in the train always irradiates a material at room temperature.

Results and discussion

Single pulse irradiation

The morphology obtained following irradiation with trains of single pulses are illustrated in Fig. 3, assuming various intensity profile types: Gaussian distribution (G), vertical DLIP (V), horizontal DLIP (H) and a DLIP profile produced with four beams (D) (Fig. 3(a)). Different pulse energies were used for the Gaussian and the DLIP distribution in order to maintain similar effective fluence value on the irradiated areas; the energy per pulse is 80 μJ for G, and 57 μJ for V, H and D. The total number of pulses is NP = 50 in all cases. For the Gaussian distribution subwavelength LSFL are formed throughout the irradiated area with orientation perpendicular to the laser beam polarisation (Fig. 3(b)) as reported in previous works^{2, 17, 43}. For the vertical DLIP some small corrugation has been observed in-between the DLIP pattern (Fig. 3, V). The type of the additional corrugation is dependent on the DLIP periodicity. More specifically, small protrusions inside the crater are produced if the DLIP periodicity is comparable to λ_L, otherwise, LIPSS are formed for $Λ_{DLIP} > \sim 2$ λ_L³⁵. The formation of these structures have been reported in a previous work³⁵. On the other hand, for the horizontal DLIP, LIPSS are formed unperturbed among the DLIP grooves (Fig. 3, H). It is possible that a surface plasmon wave-assisted mechanism is capable to explain the fabrication mechanism of those structures, however, further investigation is required³⁰. Finally, for a four beam-based DLIP, formation of ripples inside holes produced in regions of high intensity is observed³⁵ (Fig. 3, D).

Figure 3.(a) Spot intensity distribution. (b) SEM images of processed surface. (c) Fourier transform diagrams, produced by trains of single pulses. The total number of pulses is NP = 50, the energy per pulse is E_tot = 80 μJ for G, and E_tot = 57 μJ for the DLIP distribution V, H and D. The red arrow in A indicates the polarisation direction.

The FFT diagrams corresponding to the aforementioned topographies are illustrated in Fig. 3(c). In these FFT images, the spatial frequency values of LSFL structures correspond to the green dotted circles, while those of DLIP structures to the blue ones. More specifically, for G, the calculated LIPSS periodicity is equal to ~818 ± 7 nm. On the other hand, for V, the competition between DLIP and the impact of G leads to a more complex profile that is indicated by the presence of an additional protrusion corresponding to a larger periodicity. As a result, the period of the DLIP, for V, is Λ_DLIP = 1660 ± 4 nm marked with blue circles in the corresponding FFT diagram. In the case of H, the LSFL formation is not disturbed by the DLIP pattern due to the configuration of the polarisation with the DLIP. Results indicate that the induced LSFL have a period of 907 ± 15 nm for Λ_DLIP = 1666 ± 5 nm. The areas corresponding to LSFL frequencies are marked in the FFT diagram, nonetheless the zone does not feature a clear peak, compared to the other cases, due to the fact that ripples formation occurs only inside the DLIP groove. In the case of D, the period obtained from the FFT diagram in the horizontal and vertical axes are Λ_DLIP = 3230 ± 100 nm and Λ_DLIP = 3114 ± 60 nm, respectively, which is equivalent to the distance between the craters. Finally, the period of the LIPSS formed inside the craters is measured to be equal to 818 ± 20 nm, (i.e. similar to that of the G beam).

Combining a Gaussian with a DLIP pulse in DPI

In this section, the result of irradiation of stainless steel with two subsequent pulses, each having a Gaussian or DLIP profile is investigated both experimentally and theoretically. The experimental results together with theoretical simulations are shown in Fig. 4, where the order of the pulses is varied as indicated. The interpulse delay was always Δτ = 500 ps, while the total energy per pair of pulses was E_tot = 80 μJ and NP = 50. The label of each figure indicates the type of the distribution and the order of the pulses; for example, “G $↓$ + D $\leftrightarrow$ ” in Fig. 4(e) corresponds to a morphology that was obtained after irradiation of the sample with 50 double pulses where the Gaussian pulse precedes while the second 2D DLIP pulse irradiates the material after 500 ps. To account for the induced topography, a parametric study has been performed to describe both the interpulse and the intrapulse surface modification and how a corrugated profile influences the energy absorption and thermal response of the irradiated material^{17, 35}. Theoretical predictions for the induced structures are derived in all cases and compared with the experimental findings for the laser conditions used. At the same time, appropriate experimental protocols were applied, aimed to emphasise the role of hydrodynamic melt flow due to the cumulative action of the two pulses. All the results obtained for double pulses at the different sequences tested, i.e. G+ V, V+ G, G+ H, H+ G, G+ D, D+ G are presented in Fig. 4 and discussed below. It should be noted that the size of the arrows in the third and fourth columns of Fig. 4 present an estimate of the magnitude of the speed of the material flow at each spatial point; these results are derived from the solution of the Navier-Stokes equations providing the dynamics of the molten material at those time points and constitute a guide to eye. A more detailed and accurate representation is presented in Section Underlying mechanism, in which the spatial distribution of the magnitude and direction of the calculated fluid velocity vectors are illustrated for a representative G+ V case. Simulation results for the rest of the irradiation schemes applied are provided in the supplementary material. Surface patterns illustrated in the fifth column of Fig.4 depict contour plots that indicate the depth variation for NP = 50 (see also Supplementary information). The maximum depth attained for the pattern is ~230 nm.

Figure 4.SEM images of the double pulse irradiation schemes used (1^st column) together with the respective FFT analyses (2^nd column) and simulation results (columns 3−5).The pulse profile and the order are indicated in the legends of the 1^st column, while the superscript arrows indicate the pulse polarization in each case. The total energy of the double pulse pair is 80 μJ/pulse and NP = 50 pulses. Simulation results are shown in each case, just before (3^rd column, Δ $τ$ = 490 ps) and right after (4^th column, Δ $τ$ = 520 ps) the arrival of the second pulse. The final morphology predicted by the simulation experiment is shown in the 5^th column (contour plots to illustrate the height). The arrows in the 3^rd and 4^th column provide a guide to eye representation of the fluid transport direction. The color bars at the bottom of the figure provide the temperature (3^rd and 4^th columns) and height (5^th column) values range.

G+ V and V+ G

Experimental results illustrate (Fig. 4) that the G+ V combination leads to the formation of LSFL structures that are affected from the impact of the DLIP pulse on the irradiated area. In particular, LIPSS originating from G, with Λ_LIPSS = 946 ± 17 nm, and DLIP structures from the V, with Λ_DLIP = 1700 ± 55 nm, pulse smoothly overlap. Theoretical predictions show that when the V pulse (i.e. $Λ_{DLIP} = 1700$ nm, for θ=18⁰) irradiates the material after Δτ, LIPSS that have been produced from the Gaussian pulse G are not destroyed from V in the region where the DLIP intensity is low; by contrast, shallow ripples are formed along the region where the intensity of V is higher. To explain the produced topography, a theoretical investigation is conducted to evaluate the role of electrodynamics, thermal effects and fluid dynamics that yield those particular morphological features. As the increase of NP affects the morphology, a feedback mechanism is incorporated into the model to derive the absorbed energy and response of the material. To evaluate the influence of the second beam (Fig 4), the lattice temperature on the transverse plane (i.e. x-y plane) is illustrated at time t= 490 ps and at t = 520 ps (the third and fourth column images in Fig 4(a), respectively), for NP= 50. Simulations indicate that there are regions in which sufficiently high values of laser energy is absorbed along the wells of the rippled zone that leads to phase transition and therefore, V irradiates predominantly a molten material. The Gaussian beam dictates the development of the horizontally orientated ripples through the excitation of surface plasmons and the frequency of LIPSS is determined by the effect of the interference of SP with the incident beam⁴⁴. For NP= 50, the LIPSS periodicity is calculated to be equal to Λ_LIPSS~915 nm that is comparable to the experimental value $Λ_{LIPSS}^{exp}$ = 946 ± 17 nm. What is also observed is that relaxation processes and temperature gradients of the induced molten material lead to melting movement away from the valleys of the rippled zone and always along the wells of the rippled zone (i.e. at t= 490 ps). On the other hand, the V pulse increases the temperature at locations on the topography where the energy deposition from the DLIP is higher (i.e. at t= 520 ps). As a result of the enhanced temperature gradient along the x-axis, melt movement occurs, indicated by respective vectors in the temperature profiles presented in Fig. 4(b). The combined effect leads to the formation of depressed zones at the valleys of the ripples (depicted by ‘darker regions’ of high temperatures in Fig. 4) and to a periodic topography that agrees with the experimental results (see the left column in Fig. 4)). A more precise illustration of the hydrothermal melt movement at each point in a particular region that illustrates the effect of the hydrodynamical process is presented in the Supplementary information. In a V+ G scheme, the DLIP ( $Λ_{DLIP} = 1700$ nm) first irradiates the material surface; as a result, in locations where the DLIP intensity is higher, the enhanced energy absorption influences the response of the material when the Gaussian beam arrives with a delay of 500 ps. In that case, the main pattern observed is the DLIP with period of Λ_DLIP = 1701 ± 9 nm. Moreover, on top of the DLIP grooves, dots that originate from a spatially modulated periodic energy absorption, induced from the Gaussian beam, lead to ‘dot-like’ periodic structures of period equal to Λ_LIPPS = 885 ± 32 nm. More specifically, G first excites charge carriers and subsequently, through relaxation, the temperature rises in the regions among the DLIP intensity maxima, i.e. the DLIP valleys, as surface depression will take place there due to hydrothermal melt movement. At the same time, the induced surface plasmon excitation in the region outside such valleys and coupling with the incident radiation leads again to inhomogeneous energy deposition.

As a result of the temperature gradients induced within the molten material the melt is pushed away from the DLIP valleys along the x-axis (shown for t= 490 ps in Fig. 4(b)). On the other hand, the G pulse increases the temperature outside the DLIP valleys driving the produced hydrothermal waves along the y-axis (shown for t= 520 ps in Fig. 4(b)). The combined effect of the above phenomena gives rise to the formation of surface patterns that comply well with the experimental observations (Fig. 4(b)). The theoretical prediction of the period of LIPSS is estimated to be ~890 nm that is similar to the experimentally measured value of ~885 ± 32 nm.

To emphasise the special role of hydrothermal movement and provide a more detailed picture of the impact of the induced hydrothermal waves, the melt fluid velocity vectors at every lattice point in a region ~3×3 μm³ are illustrated for G+V at 490 ps and 520 ps respectively (Fig. 4). These vectors present an estimate of the magnitude of the speed of the material flow at each spatial point, providing the dynamics of the molten material at those time points.

G+ H and H+ G

In the case of a double pulse experiment in which the contributing pulses to the DLIP yield an interference pattern parallel to the x-axis, the impact of the DLIP pulse leads to a different surface topography. More specifically, as in the G+ V case, it is expected that the first pulse (G) excites surface plasmons and leads to subwavelength LIPSS oriented perpendicularly to the laser polarisation that will lead to a rippled profile. Nonetheless, analysis of the experimental results indicate that only one ripple frequency is observed, that of the DLIP, with a period of Λ_DLIP = 1630 ± 14 nm. In Fig. 4(c), the respective temperature profile at t= 490 ps is illustrated, i.e. before the arrival of the H pulse onto the hot material. Upon the incidence of the H pulse (i.e. t= 520 ps in Fig. 4(c)), the temperature in the valleys of the periodic regions formed by G will further rise above that attained with the initial DLIP pulse. As shown by the respective simulations shown in the same figure, hydrodynamical movement drives the melt flow parallel to the y-axis before solidification, giving rise to a topography formation that closely resembles the experimental one. It should be noted that $Λ_{DLIP}$ is a multiple of the predicted LIPSS period (Λ_LIPSS~850 nm) and therefore there is no overlap between the patterns produced from the combined effect of G and H.

A different profile is induced for the H+ G pulse sequence (Fig. 4(d)). In that case, the DLIP pulse leads to the formation of long stripes (tops and valleys) periodically situated at distances defined by $Λ_{DLIP}$ , however, the polarisation beam direction of the DLIP pulse (along the x-axis) yields excitation of surface plasmon waves of ~880 nm periodicity and hydrothermal waves that propagate along the x-axis. Similar structures have been observed and interpreted in a previous report³⁵. On the other hand, irradiation with a G pulse with a 500 ps delay, enable excitation of surface plasmon along the y-axis, with hydrothermal waves propagating along this axis, as well and yielding LIPSS parallel to the x-axis. As noted in the previous paragraph, the fact that the $Λ_{DLIP}$ is a multiple of the predicted period of the LIPSS, produced from G, leads to surface patterns similar to experimental results (Fig. 4(d)). By contrast, a different relation between the periodicity due to the $Λ_{DLIP}$ and G would lead to a distinctly different profile with a substantial overlapping of the respective periodicities. In Fig. 4(d), the temperature profile and melt movement are illustrated at t= 490 ps and t= 520 ps, respectively. According to the simulation results, vertical ripples produced by H are much more pronounced in H+ G case because the melt flow is already developed upon arrival of G, in contrast to the G+ H pulse sequence case.

G+ D and D+ G

Interesting morphologies are obtained when a combined G plus a D pulse with $Λ_{DLIP}$ = 3100 nm is used (Fig. 4(e)). Theoretical results indicate that the action of the G beam leads to an expected surface plasmon excitation and hydrothermal waves, while the D pulse increases further the temperature at points where the DLIP intensity is high. Thus, a periodic profile is produced ( $Λ_{LIPSS} = 824 \pm 46 nm$ ) with the ripples oriented parallel to the x-axis that are interrupted with periodically situated subsided regions due to the impact of the D pulse (Fig. 4(e)). The latter regions exhibit the ripple periodicity attained by G. The respective simulated temperature profile as well as the propagation of the induced hydrothermal waves are presented in Fig. 4(d) for t= 490 ps and t= 520 ps respectively, for NP= 50.

By contrast, when D precedes the irradiation to G, a quite different profile, comprising high-aspect ratio asperities, is obtained. The horizontal spacing among those structures is 3163 ± 55 nm, while the vertical one is 3072 ± 67 nm, coinciding with the distance between the craters in the D case of Fig. 3. Among the peaks, traces of LIPSS with a period in the range of 620–720 nm are observed both in horizontal as well as in vertical direction. More specifically, the intensity profile of D (Fig. 2) leads to a temperature rise in the locations where the energy deposition is high. Furthermore, repetitive irradiation with D (at NP > 2) will enable excitation of surface plasmons and the subsequent formation of ripples oriented parallel to y-axis. Following irradiation with a G pulse, i.e 500 ps after the impact of the first pulse, the material in the region which is not affected by the D pulse will be excited and surface plasmon waves will be developed, before leading to a rippled periodic pattern similar to the one observed in SEM images (Fig. 4(f)). The respective temperature profile and hydrothermal melt movement at t= 490 ps and t= 520 ps, respectively, are illustrated in Fig. 4(f).

To summarise, the experimental observations and simulations indicate the following: comparing G+ V and V+ G, LIPSS originating from G are more pronounced when G comes first. Comparing G+ H and H+ G, LIPSS resulting from the action of H are observed only in the case that H irradiates the material first. Similarly, for G+ D and D+ G, LIPSS are more pronounced when the G precedes the action of D. LIPSS are always formed perpendicular to the polarisation of the first pulse. In addition, the above discussion shows that promotion of SP excitation-based mechanisms for the formation of LIPSS are influenced a lot on the following: (i) G always leads to excitation of SP waves. (ii) the DLIP polarisation and periodicity (as also shown in ref.³⁵) account on whether SP excitation is possible or fluid transport will determine surface pattern formation. (iii) The order of irradiation with Gaussian and DLIP pulses can impact on the conditions for periodic structure formation through SP excitation or fluid movement along prepatterned topographies.

Generation of complex patterns

In the previous sections, the effect of different pulse intensity profiles and pulse sequences was investigated and evaluated. It is evident that the spatio-temporal intensity distribution significantly influences the hydrothermal melt waves and the subsequent superposition of fluid transport derived from each pulse separately. It was shown that the combined effect of different intensity profiles predominantly influences the behaviour of the produced molten material through variation of its vorticity; this effect has proven to influence the final topography attained, giving rise to patterns of higher geometrical complexity. Two parameters that are expected to tune the surface morphology are the laser energy and the number of pulses; the increase of both parameters leads to higher temperature gradients and therefore to large vorticity changes. A characteristic example of the evolution of surface topography upon increasing NP at specific pulse energy is presented in Fig. 5 (rows), for the D+ G pulse sequence.

Figure 5.Stainless steel surface morphology following variation in pulse energy and NP, as indicated in the respective legends, for the D+ G pulse sequence case. The corresponding result of theoretical simulations for 20 μJ pulse energy is shown as an inset defined by the dashed rectangle, illustrating (i) the temperature profile for NP = 50 and t = 520 ps and (ii) the final surface morphology for NP = 50.

Results demonstrate the pronounced influence of the first pulse in the sequence in the feedback mechanism, while each subsequent dose of the series enhances the depth of the induced topography. Similar results (not shown here) are derived for the rest of the pulse sequences. On top of that, the role of hydrothermal phenomena and the impact of the topography in the produced temperature gradients that lead to a dynamic evolution of the surface patterns at increasing dose, are very critical. Simulation results indicate that deeper profiles promote larger energy deposition locally that influences the melt movement.

In Fig. 5, indicative results on the combined effect of pulse energy and energy dose (ΝP) on the induced topography, upon irradiation with the D+ G sequence, are also shown. Notably, significant variations in the surface morphology attained, is observed, and an ensemble of complex structures can be realized. In particular, for E_tot = 80 μJ and for low NP random LIPSS structures are formed, while the fingerprint of the DLIP pattern becomes evident at elevated NP. Similarly, for E_tot = 60 μJ and NP = 10, the DLIP undulation is not visible and the surface comprised LIPSS patterns oriented in different directions. For low NP, neighbouring ripples create square patterns, while at the higher NP= 20 a supra-wavelength morphology resulting from DLIP takes place. Further irradiation (e.g. NP= 50) reveals the impact of the DLIP, as the DLIP pattern is clearly visible on the surface surrounded by dotted structures (see also Fig. 5). Finally, at even lower energies (E_tot = 40 μJ and NP = 50), LIPSS-like complex structures, resulting from the combined effect of the two pulses, become present. These results indicate that hydrothermal waves possibly coupled with mass removal effects should be the predominant factors that explain the induced topography. Comparing the three images of high magnification (Fig. 5, zoomed region for NP = 50) the impact of fluence in the generation of high-complexity structures, is unveiled. When the fluence is low and is close to the ripple formation threshold (E_tot = 40 μJ) ripples are only formed in perpendicular directions due to DPI. When the fluence is higher (E_tot = 60 μJ) the temperature gradient applied by the DLIP pattern is strong enough to displace material and form valleys. Blue arrows indicate the trajectory of the flow and the orange dotted circle indicates the area that the material is accumulated. In even higher pulse energy (E_tot = 80 μJ) the temperature gradient is stronger, and the volume of the displaced material is higher (blue arrows in Fig. 5, zoomed region for NP = 50 and at 80 μJ). In that case, the displaced material is accumulated away from the areas irradiated with higher intensity and forms the observed protrusions marked by the orange dotted circle. Simulations have been performed for E_tot = 20 μJ to illustrate the topographies attained at low energies and the role of hydrothermal phenomena. Results show that the ultrafast dynamics and the induced hydrothermal effects lead to a different surface profile than the one produced for 25 μJ (shown in Fig. 4). More specifically, lower energies produce smaller temperature gradients therefore the fingerprint of either D or G pulses will be more pronounced.

The ability to simultaneously control the two key structure formation mechanisms, namely the inhomogeneous absorption and the microfluidic melt motion, offers potentially endless possibilities in controlling laser induced morphology in the submicron scale. Apart from the complex structures presented above, in the following we show that suitable combination of DPI parameters enables the formation of hierarchical complex structures as well, which in some cases are similar to those observed in natural systems⁴⁵).

For example, Fig. 6 presents the generation of a multiple-length-scale morphology upon irradiation with E_tot = 70 μJ and NP = 100 for the D+ G irradiation scheme. The features of the three prominent periodicities, visualized via a reciprocal FFT algorithm are presented in the same figure (Reconstructed Surface Profile). In particular, the first with smaller period (~500 nm) corresponds to LIPSS which are oriented perpendicularly in adjacent areas, similarly to E_tot = 40 μJ and NP = 50 (see Fig. 5). The highest one corresponds to the effect of D pulse of DLIP, whilst the middle one corresponds to the 1D DLIP period, when the V and H schemes are combined. The advantage of such complex surface morphology attained, compared to the previously reported DLIP morphologies, is that it combines multiscale structures (DLIP and LIPSS) over the whole processed area in a single step.

Figure 6.Complex structures generated with E_tot = 70 μJ, NP = 100 with the D+G pulse sequence irradiation scheme. The SEM image is analysed with FFT and reconstructed surfaces comprising structures at different length scales are illustrated. The areas corresponding to the different periodicities are highlighted in FFT diagrams.

Underlying mechanism

In the previous sections, the characteristics of complex structures’ formation upon irradiation of stainless steel with double pulses of Gaussian and DLIP profiles, were investigated, both experimentally and theoretically. The impact of parameters such as the pulse profile, order of pulses, interpulse delay, pulse energy and energy dose was presented and discussed. A key element of the formation processes is the dynamic variation of the hydrothermal melt movement due to the action of the two irradiation pulses, which depends on the interpulse delay. In particular, at t<500 ps a melt flow is developed by hydrothermal waves created due to induced temperature gradients. Upon the impact of the second pulse the temperature profile changes drastically, as indicated by the pronounced decrease of the reflectivity in areas that are in molten phase³⁵; indeed such pre-heated regions enable high energy absorption levels locally when the second pulse irradiates the material. As a consequence of the absorbed energy following the second pulse, together with the pre-existing surface temperature due to the first pulse, a complex temperature profile is generated, that is unfeasible to be created by single pulse irradiation. To further shed light on the effect of such temperature profile on hydrothermal melt movement, we present in Fig. 7 simulation results, for the G+ V case, for very low NP values, i.e. during the first steps of structures’ formation. When NP = 2, the molten material at t = 490 ps (Fig. 7(a)) develops a complex flow trajectory that differs substantially from that observed at t= 520 ps (Fig. 7(b)), as a response to the applied temperature profile. Considering that resolidification occurs in the order of nanoseconds, for the range of energies tested⁴⁶, a significant volume of the material has been already displaced upon the impact of the second pulse. Therefore, the overall trajectory followed by the molten material upon a single pair of double pulse irradiation can be considered as the product of the two separate trajectories imposed by each individual pulse. This demonstrates that the shape of the induced structures can be controlled by imposing different pulse profiles with suitable interpulse delays.

Figure 7.Simulation results showing the surface temperature profile obtained for the G+V case at different time moments and NP, as indicated. The arrows in the figures in the zoomed regions illustrate the spatial distribution of the magnitude and direction of the calculated fluid velocity vectors.

A key element is the value of the pulse delay (~500 ps). The choice of such a large value for the pulse separation between the two constituent pulses of the train of pulses was to illustrate the impact of the molten material when the second pulse arrives³⁵. Although electron-phonon relaxation processes are completed within some tens of picoseconds and molten material is produced during that timescale, the lattice temperature of the material is quite high (>6000 K as seenFig. 3(c)). Furthermore, high temperature gradients (that occur till up to some hundreds of ps) is set to influence greatly the response of the material leading to complex dynamics. Thus, any effort to highlight the impact of the hydrothermal profile that is generated before the second pulse irradiates the material will be hindered. By contrast, if we allow longer pulse separations, the temperature of the fluid relaxes to substantially lower (but still above the melting point) values and therefore it is easier to reveal the variation in the fluid dynamics that is caused as a result of exposure to the second pulse.

The surface relief attained in the early stages of structures’ formation will eventually determine the topology of the structures obtained for elevated number of pulses, due to the optical and hydrodynamic feedback mechanisms. Indeed, upon increasing the dose to NP = 5, the temperature gradient of t= 490 ps and t = 520 ps differ mainly in amplitude, but not in topology. The stability of the topology of the structures is evident in the fact that the vector fields, observed in Fig. 7(b) and 7(c), are quite similar. Conclusively, the combined action of the two pulses has greater impact in the first steps of the structure formation, defining the pattern of the structures, whilst the extraordinary trajectory followed by the molten material as NP increases, enhances the structures’ aspect ratio. Following the formation of significantly deep valleys and high hills, i.e in NP ~ 50 (see Supplementary information), the surface morphology will prevail over the pulse profile in the determination of the temperature gradient. It is evident that the laser conditions and selection of NP (i.e. $N P \leq 50$ ) were appropriately chosen to correlate experimental observations with simulations results. Certainly, higher laser fluences or increasing NP are expected to produce deeper profiles and gradually alter the surface topography. This might be the objective of further investigation, however, in the current study we focused on how consecutive pulses of different intensity profiles influence the final topography. On the other hand, it is also important to correlate the value of Δτ with particular characteristics times such as those associated to electron-phonon relaxation and Marangoni flow. The pulse separation is much larger (about two orders of magnitude) than the electron-phonon relaxation time. By contrast, it is much smaller than the expression that has been used in some reports^{47, 48} as the characteristic time for Marangoni flow effect

5 $τ_{m} = \frac{η L^{2}}{| \frac{d σ}{d T_{L}} | T_{m} h},$

where η stands for the dynamic viscosity, T_m is the melting temperature, h corresponds to the thickness of the volume in liquid phase, $| \frac{d σ}{d T_{L}} |$ is the absolute value of the derivative of the surface tension with the lattice temperature and L is a radial dimension. Assuming the values for the parameters for stainless steel, and setting L (>10) μm,h=20 nm, $τ_{m}$ >50 μs. According to the use of $τ_{m}$ as a criterion that determines when thermocapillary effects become competent in the surface modification process, it appears that the Marangoni flow does not appear to play the predominant role in the surface structure formation at very small timescales. The fact, though, that pulse separations of Δτ=500 ps lead to significant surface modification and simulations demonstrate the impact of thermocapillary effects indicate that further exploration of Eq. (5) is required as it yields an underestimation of the fluid movement due to surface tension.

The analysis performed in this work demonstrated the impact of the polarisation-dependent process on the LSFL formation for either G or DLIP pulses. Experimental results showed that while periodic patterns were visible at small pulse energy levels or doses, the DLIP pattern imprint became more pronounced at higher pulse energies (Fig. 5). As shown in Fig. 5, experimental observations are significantly sensitive to the laser conditions and small changes in the parameters such as the energy values or dose are expected to lead to different topographies due to the involvement of more complex mechanisms (i.e. ablation processes, stronger temperature gradients, overlapping between DLIP and surface plasmon wave period, etc). Although the applied conditions used in this work were predominantly aimed to interpret the experimental observations in the absence of ablation effects (discussed previously in ref.³⁵), a generalized scenario can be addressed and investigated in more detail by incorporating appropriate modules in the theoretical model. For example, a more precise evaluation of the spatial distribution of the laser energy and the energy absorption from the material can be achieved through the employment of advanced computational electrodynamics based models²⁵. Furthermore, appropriate modules based on atomistic molecular dynamics could also be incorporated to provide an alternative description of microscopic processes related to kinetics and phase transition-related mechanisms⁴⁹. Another issue that requires further investigation is whether the large temperatures that are reached during the exposure of the material to a train of double pulses (>6000 K) are sufficient to cause ablation/evaporation. More specifically, the attained temperatures are substantially higher than the melting (1800 K)/boiling (3100 K) points. As emphasised in this work and previous reports^{35, 40, 41}, mass removal is associated with temperatures larger than 0.9T_cr; therefore more exploration is needed to evaluate if that threshold constitutes an overestimation of evaporation/mass removal. Nevertheless, while further development of the model towards providing a more complete description of the processes could be the objective of a future work, the present approach demonstrates in a consistent way the capability to control laser-matter interaction through tailoring the features of the spatial intensity profile. The multiscale model that was employed is aimed also to allow an optimisation of the pattern features assuming the influence of the underlying physical processes and more importantly the role of hydrodynamical fluid transport. The agreement of simulation and experimental data demonstrates that using double, spatially tailored, fs pulses it is possible to actively control the melt microfluidic flow and in principle to make possible the generation of laser induced structures by design.

Conclusions

In conclusion, the investigation performed in this work indicated that combining Gaussian beams with DLIP in double pulse trains enables the generation of unique sub-micron surface topographies with increased complexity. Results showed that the order of the Gaussian and DLIP pulses in train sequences influences the final surface profile attained. Furthermore, experimental results and multiscale modelling revealed that both the spatial intensity of the two pulses, as well as hydrodynamical effects have a significant impact on the pattern features. Our experiments and results presented here emphasise the capability to actively tailor the microfluidic melt motion that dominates the structure formation process, via controlling the applied temperature gradient’s temporal profile. More importantly, the unique irradiation schemes examined here lead to the generation of novel complex morphologies comprising features in multiple length scales. This demonstrates an unparalleled capacity towards tailoring laser-induced morphology and obtaining complex topographies for a variety of applications.

References

[6] Diels JC, Rudolph W. Ultrashort Laser Pulse Phenomena: Fundamentals, Techniques, and Applications on A Femtosecond Time Scale 2nd ed (Academic Press, Amsterdam, 2006).

[8] Cirelli RA, Watson GP, Nalamasu O. Optical lithography. In Buschow KHJ, Cahn RW, Flemings MC, Ilschner B, Kramer EJ et al. Encyclopedia of Materials: Science and Technology 6441–6448 (Elsevier, Oxford, 2001); https://doi.org/10.1016/B0-08-043152-6/01138-4.

[10] Sugioka KJ, Meunier M, Piqué A. Laser precision microfabrication (Springer, Berlin, 2010); http://www.loc.gov/catdir/enhancements/fy1316/2010931868-b.html.

[39] Pozrikidis C. Introduction to theoretical and computational fluid dynamics (Oxford University Press, New York, 1997); http://www.loc.gov/catdir/enhancements/fy0603/96026420-d.html.