## 2 Theoretical Model
### Fluid Transport. Ripples and Grooves
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 \(NP\geq 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.