Laser-induced periodic surface structures: Arbitrary angles of incidence and polarization states

We demonstrate that the finite-difference time-domain (FDTD) method can be used to study the formation of laser-induced periodic surface structures (LIPSSs) under oblique incidence and arbitrary polarization states. The inhomogeneous energy absorption below a rough surface at oblique incident angle is studied by the FDTD method and compared to the analytical efficacy theory developed by Sipe et al. [Phys. Rev. B 27, 1141 (1983)]. The two approaches show excellent agreement. The surface topographies of both low-spatial-frequency LIPSSs (LSFLs) and high-spatial-frequency LIPSSs (HSFLs) at oblique incidence and various polarization states (linear, circular, radial, and azimuthal) are simulated by taking into account topography-driven interpulse feedback effects. For sand mixed sand p-polarized beams at large incident angles, the simulated LSFLs show orientations that are neither perpendicular nor parallel to the laser polarization. Their physical origin is explained by a nonzero angle between the in-plane wave vector of the incident light and the scattered light. By using circularly-polarized beams, the simulated surface topographies are further shown to be in excellent agreement with the triangular LSFLs and the fine-dot nanoscaled structures reported in the literature. In addition, the simulation of HSFLs suggests that their periods have almost no dependence on the incident angle nor the polarization angle, while their surface topographies resemble that of blazed gratings at large incident angles for p-polarized and mixed sand p-polarized beams. The origin is explained by the appearance of asymmetry in the scattered near-field light. The extension to oblique incidence and arbitrary polarization states demonstrates that the FDTD approach on LIPSSs is a highly versatile and powerful technique which has the potential to be one of the foundations for a complete modeling and understanding of LIPSSs under various irradiation conditions.


I. INTRODUCTION
Laser-induced periodic surface structures (LIPSSs, also termed ripples) were first observed by Birnbaum [1] in 1965 on semiconductors irradiated by a ruby laser. Since then, LIPSSs attracted growing interest from both the fundamental and application points of view and are under intense research. LIPSSs usually emerge as gratinglike arrangements of matter composed of (quasi) periodic lines on the irradiated surface and they are found to form on almost any material with a huge span of pulse durations ranging from continuous wave to femtosecond [2]. This observation implies that regular patterning of surfaces using LIPSSs is a universal manifestation of an intriguing matter reorganization potential. Their control, increased uniformity, and high resolution have led to advanced applications in functional surface engineering [3][4][5] from wettability control to optical color marking and surface blackening, microfluidics, and optical data storage [6]. Despite these successful applications, a complete understanding of their physical origin under various material properties and irradiation conditions is still unestablished and remains to be a subject of intense investigations and debates. Theoretical framework in the electromagnetic theory [7][8][9][10][11][12][13] and mechanisms based on self-organization and surface instability [14][15][16][17] were equally developed to explain the origin of LIPSSs. There is plenty of experimental evidence that the period and orientation of LIPSSs show a clear correlation to laser irradiation parameters such as polarization, angle of incidence, and fluence [8,[18][19][20][21][22], favoring the electromagnetic nature of LIPSSs, at least in their earliest formation stage. The most widely found type of ripple, the low-spatial-frequency LIPSS (LSFL) with a period comparable to the illuminating wavelength [2], is well-explained in certain aspects by a pioneering electromagnetic theory developed in the 1980s by Sipe et al. [7], also known as the efficacy theory. This theory, with a wide acceptance today, attributes LSFLs to inhomogeneous laser energy absorption due to interference between surface scattered light and the incident light. The efficacy theory, a pioneering theory on LIPSSs [18], is successful in explaining the period and orientation of LSFLs. With the advancement of ultrafast laser techniques, another type of LIPSS, the high-spatial-frequency LIPSS (HSFL) [2] with a period smaller than half of the illuminating wavelength oriented either orthogonal [8,[23][24][25][26][27][28][29] or parallel [9,29] to the polarization, was found to form after irradiations only with subpicosecond pulses. The HSFLs show real potential for fabricating highly uniform regular nanopatterns, renewing interest in the topic. The experimental evidence that HSFLs only form upon irradiation by ultrashort laser pulses and that their orientations are strictly correlated to the direction of laser polarization indicates a preference on ultrafast energy deposition mechanisms as the origin of HSFLs. However, it should be noted that the efficacy theory fails to describe HSFLs [30][31][32], mainly due to a lack of feedback mechanisms [10,11].
To overcome theses limitations, a numerical approach [10][11][12][13][33][34][35][36][37][38] adopting the finite-difference time-domain (FDTD) method was developed within the last decade to simulate the formation of LIPSSs. The main advantages of the FDTD approach over the pioneering efficacy theory are the flexibility to model light-matter interactions for any geometries directly from Maxwell's equations, a natural and straightforward incorporation of both topography-driven interpulse feedback effect [11,12,35], and intra-pulse feedback of permittivity dynamics [13,33,34]. The FDTD approach has shown great potential and success in unravelling the origin of many types of LIPSSs including HSFLs. Skolski et al. predicated a variety of different types of LIPSSs by FDTD feedback simulations [10,11]. Rudenko et al. elucidated the respective role of quasicylindrical and surface plasmon waves in the formation of LSFLs by electromagnetic calculations adopting the FDTD method [35]. FDTD simulations with the incorporation of permittivity dynamics [13,33] were also used to understand the formation of volumetric nanograting in glass. It has been shown that the complexity of LIPSSs is intrinsically initiated by the coherence of the laser field via superposition between topography-driven scattered field and the incident field [12]. Very recently, multiphysical simulations combining Maxwell's equations and Navier-Stokes equations were proposed to unravel the origin of LIPSSs beyond the laser absorption stage, in which the first stage is to model the absorption stage using the FDTD method [39,40]. Despite these fruitful outcomes, to the authors' knowledge, all the published FDTD approaches on LIPSSs are restricted to normal incidence and linearly polarized light [10][11][12][13][33][34][35][36][37][38][39] with only very few exceptions on nonlinearly polarized beams at normal incidence [41]. Yet, it is widely known that the incident angle is an important parameter which determines the period [42][43][44][45][46][47][48] and even the orientation of LSFLs [18]. The polarization state of the laser beam is another important parameter in controlling the orientation, local curvature, and shape of the resulting LIPSSs. It has been shown that a variety of surface patterns such as triangular LSFLs [49,50], nanoscaled fine-dot structures [27], chains of spherical nanoparticles [51,52], and speckled structures [53,54] can be obtained by circularly and elliptically polarized beams. LIPSSs with radial and azimuthal polarization states were also obtained [55][56][57][58], which were shown to reveal the local polarization states of tightly focused beams [59]. In addition, the correlation between the HSFL period and incident angle is a very interesting problem.
It is the purpose of this paper to extend the FDTD approach on LIPSSs to arbitrary angles of incidence and polarization states. In an effort to validate our approach, we first compare the results of the FDTD approach with the efficacy theory. Second, we go beyond the efficacy theory by taking into account the topography-driven interpulse effect. We focus on the periodicity and orientation of the simulated LIPSSs and their relation with the surface scattered field. Anomaly LIPSSs with orientations that are neither perpendicular nor parallel to the laser polarization are found under particular conditions. Their physical origins are addressed. Third, the surface topography of HSFLs is studied as a function of incident and polarization angles. The dependency of its period and orientation as well as its general topography are discussed in detail. Finally, the surface topographies of LSFLs and HSFLs under circularly polarized beams are simulated and are in good agreement with the triangular LSFLs and the fine-dot HSFLs reported in the literature. LIPSSs simulated by beams of radial and azimuthal polarizations are also addressed.

II. EFFICACY THEORY AND THE FDTD APPROACH
The efficacy theory [7] is based on analytical solutions of Maxwell's equations in the presence of a randomly distributed rough surface and yields an expression for the inhomogeneous energy deposition into the material. The main assumptions of the theory are (1) the surface roughness is randomly distributed and (2) the roughness thickness must be much smaller than both the illuminating wavelength and the period of LIPSSs being studied. Assumption (1) restricts the theory from adding any topography-driven feedback effect as the second pulse sees a surface which has already been modified and thus cannot be seen as random. Assumption (2) prevents the theory from properly describing HSFLs, as pointed out by Skolski et al. [10], that their amplitude in Fourier space continues to increase with large values of k to infinity, which is apparently unphysical. Assumption (2) also limits the efficacy theory to distinguish the role of the surface roughness type on the energy depositions, which is easily accessible by the FDTD approach [60]. Other limitations of the efficacy theory include the omitting of intrapulse feedback and limited ability to model beams with complex polarizations. Nevertheless, its comparison with numerous experimental observations demonstrates that the efficacy theory is very successful in predicting the period and orientation of LSFLs.
Instead of using the analytical efficacy theory, we tackle the inhomogeneous energy deposition below a material's rough surface by numerically solving Maxwell's equations with the FDTD method. The FDTD method was first introduced by Yee in 1966 [61] and it is based on numerically solving the two coupled curl equations, where E is the electric field, H is the magnetic field, 0 is the vacuum permittivity, μ 0 is the vacuum permeability, and J is the current density. The current density can be written as the vector sum of bound current and free current: The bound current is calculated as where χ is the linear electric susceptibility. We use the Drude model to simulate materials with metallic optical properties. The equation of motion for the free electrons yields the following differential equation: where ν e is the collision frequency and n e is the free electron density. Upon irradiation of ultrafast laser pulses, the free electron number and collision frequency can be dynamically varying during the pulse, resulting in intrapulse feedback to the energy deposition. This effect is the most pronounced in semiconductors and dielectrics as the free electron density increases several orders of magnitude due to interband transitions. For metals, since there is already a large amount of free electrons before irradiation, the intrapulse feedback effect is expected to be small. It has been shown that including the permittivity dynamics does not strongly affect the predicted LIPSS period in metals [62]. We note that the objective of the paper is not to provide an accurate model of LIPSSs for a specific material but to extend the FDTD approach to oblique incidence and arbitrary polarization states, thus the material permittivity is kept as a constant in the simulations. The exclusion of any intra-pulse feedback also means that a direct comparison with the analytical efficacy theory is possible. Surface roughness in the simulations is introduced by following the method first proposed by Skolski et al. [10]: The roughness layer is represented by a randomly distributed binary function b(y, z). When b(y, z) = 1 the permittivity in the roughness layer is set to the material under study, while when b(y, z) = 0, the permittivity is set to the vacuum value. In the efficacy theory, there are two parameters which define the surface roughness, the filling factor f and the shape factor s. The filling factor is the averaged value of the function b(y, z), the shape factor is comparable to half of the aspect ratio of the roughness particle [7,10]. Figure 1 shows a schematic of the FDTD simulation scheme. The incident direction (k vector) of the laser beam is described by two angles, the angle θ with respect to the normal of the surface and the angle α which is the angle between plane of incidence and the y axis. These two angles completely determine all the possible incident directions. To compare the FDTD results with the efficacy theory, linear polarization is first assumed, other polarization states will be discussed later. The angle between the electric field vector and the normal of the incident plane ( n in Fig. 1) is defined as γ . The angle γ thus determines all the possible linear polarization states. Perfectly matched layers are implemented in the front and rear surfacesof the computation domain to absorb unphysical reflections due to the truncation of the domain. Periodic boundary conditions (PBCs) are implemented in the other four (lateral) boundaries.
The use of PBCs maximizes the interaction area and allows the implementation of oblique incidence. The inclusion of oblique incidence requires a special treatment at the periodic boundary. It can be shown that when writing the PBC at oblique incidence, field values at future time steps are needed, which brings a fundamental difficulty for a straightforward The light-blue ellipse illustrates the plane that is perpendicular to the wave vector of the incoming beam. The electric field vector is in the light-blue plane. For linear polarizations, the electric field vector forms the angle γ with respect to the surface normal n of the plane of incidence. For γ = 0 • , s polarization. For γ = 90 • , p polarization. The small cubes are the randomly generated roughness particles on the material's surface (not to scale, only for illustrative purposes). The periodic energy modulation is calculated in the plane that is immediately below the roughness (here it is shown only for illustrative purposes). The coordinate representation is used throughout the paper.
implementation. This difficulty is overcome by introducing the constant transverse wave number (CTW) method [63], which has been shown to be advantageous compared with other methods. In CTW, the PBC is written in the frequency domain and Fourier transformed to the time domain. Future field values are avoided by keeping the transverse wave number constant. The only drawback of this method is that complex field values are needed, thus twice the computer memory and computational workload are required. More details on CTW can be found in Ref. [63]. To implement laser beam incidence at oblique angles and arbitrary polarization states in the simulations, incident electric field components on the three cartesian axes need to be written as a function of the angles α, θ , and γ . For linear polarizations, it can be found by analytic geometry for the most general case-the projections of the incident field amplitude on the three cartesian axes.
In this paper, we only discuss the situation of α = 0, as other cases are equivalent by a coordinate transformation. For α = 0 and a linearly polarized beam with incident angle θ and polarization angle γ , the three cartesian components of the incident field amplitude are where E 0 is the electric field amplitude of the incident pulse. Equation (4) is used as a source to launch the incident wave in the simulations for linearly polarized beams.
We have developed a three-dimensional FDTD code that incorporates the CTW method to solve Eqs. (1)-(3). The lateral grid spacings ( y and z) are chosen between 10 nm and 40 nm, depending on the type of LIPSS being studied (HSFL or LSFL). The vertical grid spacing x used is 10 nm throughout the paper. The time step of the simulation is chosen to satisfy the Courant-Friedrichs-Lewy stability condition [64]. The time step is on the order of 0.02 fs, depending on the exact grid spacing used. The lateral dimension of the simulation domain varies between 7 μm × 7 μm and 30 μm × 30 μm to ensure that several dozens of ripple (either HSFL or LSFL) periods can be resolved. The thicknesses of the samples used in the simulations are 200-500 nm, depending on the number of pulses applied, to ensure a bulk optical response. The wavelength of the laser beam used in the simulations is 800 nm throughout the paper. The temporal envelope of the source is assumed to be a Gaussian pulse with full-width-half-maximum of 30 fs and a total simulation time spanning from −90 fs to 90 fs. The absorbed energy density at a certain depth below the rough surface E ab (y, z) is computed from the electric field amplitude, which itself is retrieved from the cycle-averaged electric field values [12,65]. Two-dimensional fast Fourier transforms F are applied to the absorbed energy densities below the rough surfaces to reveal the absorption features in the frequency domain and compared with the efficacy theory. In the efficacy theory calculations and the FDTD calculations, the same filling factor f = 0.1 as in Refs. [7,66] is used. The shape factor s used in the efficacy theory calculations is determined from the aspect ratio of the roughness particles used in the FDTD calculations, for an equal comparison between the two approaches.

III. COMPARISON WITH THE EFFICACY THEORY AT OBLIQUE INCIDENCE
The comparisons with efficacy theory are performed for silicon with different levels of excitations. We note that although the comparisons here are shown only for silicon, both methods support any values of the complex permittivity, and therefore are not limited to silicon. The complex permittivity of excited silicon is calculated via the Drude model, where un , n e , m e , m * , and ν e are, respectively, the unperturbed permittivity of silicon at 800 nm, the free-carrier density, the free electron mass, the optical effective mass of the electronhole pairs, and the Drude collision frequency. The values m * = 0.18 and ν = 1/1.1 fs −1 are used for femtosecondlaser-excited silicon [67,68]. shape and position are almost not affected by the angle of incidence. This is in agreement with reported experimental data [42], which implies that surface plasmon waves are not the origin of the type-d features. The spread of the typer features in the k y -direction seems to slightly increase at larger incident angles. In contrast, the type-s features show strong dependence on the incident angle. As the incident angle is increased, the center of the circle which encloses the type-s features move away from the origin (k y = 0, k z = 0) along the k y axis. Its center position in the k y axis can be well described by k ys = k 0 sin θ . This is to be expected, as at oblique incidence, the incoming photons already have momentum component k y = k 0 sin θ . The scattering due to surface roughness adds additional wave vectors with preferential directions depending on the nature of the scattered field (surface plasmon waves and/or the cylindrical scattered waves [35]). Through coherent interference with the incident wave, this scattered field gives rise to periodic energy absorption with the same preferential directions. Therefore, the wave vector components of the incoming photons determine the central position of the type-s features. As the incident angle is increased, the period of the type-s features along the z axis remains unchanged, yet both approaches show that their periods along the y axis are reduced. The overall effect is then to form ripples that are neither parallel nor perpendicular to the polarization direction. Their origins will be discussed further later in the multipulse results. Figures 3(a) and 3(b) show comparisons between the FDTD calculated |F (E ab )| and the results of the efficacy theory at various incident angles for p-polarized beams. Again, the FDTD results and the efficacy theory are in good agreement. At normal incidence, the p-polarized result is identical to the s-polarized result, except for a 90 • rotation. However, at oblique incidence, the p-polarized result is not a simple 90 • rotated version of the s-polarized result. The type-d and the type-r features remain largely unchanged (except for a 90 • rotation) as the incident angle increases although at large incident angles the type-d features become less pronounced. Similar to the s-polarized result, the circle center position of the type-s features follows the lateral wave vector component of the incident wave k y = k 0 sin θ . However, the position of the type-s features along its enclosure circle (the dotted circles) changes: The pattern is rotated by 90 • and the intensity of the two branches becomes asymmetric. Thus, the period of LSFLs for a p-polarized beam can be approximately written as If surface plasmon waves are the dominant origin of LSFLs, the above formula should be written as p = λ/[λ/λ sp ± sin θ ], where λ sp is the surface plasmon wavelength. This is in agreement with the ripple period obtained by the phase matching condition involving surface plasmon waves [44]. Both approaches show that the intensity of the solution with higher k y values is much higher than the lower branch, favoring a reduced ripple period at higher incident angles. This is due to the asymmetry in the scattered-light patterns (Fig. 7) and will be discussed later.
The agreement with the efficacy theory for s and p polarizations motivates us to study more complex cases: when the polarization angle γ is in between 0 • and 90 • . In these cases, the incident wave has both s and p components. Nevertheless, the efficacy theory describes s-and p-polarizations separately [7,30] and cannot be applied to these mixedpolarizations. We note that the overall solution here is not a simple linear combination of separate s and p solutions due to the coupling between the two modes in the presence of nanoscale surface roughness. Figures 4(a) and 4(b) show the FDTD calculated |F (E ab )| when the polarization angle is varied between 0 • and 90 • , at two oblique incident angles. From the figure, we see that the type-d and type-r features rotate through the same angle as the laser polarization directions around the (k z = 0, k y = 0) point. In Figs. 4(a) and 4(b), the intersected blue lines are plotted using the same angle as γ . We see that the position of the type-s features rotate through the same angle γ around the (k z = 0, k y = k 0 sin θ ) point. From Figs. 2-4 we can draw the following conclusions: (1) The shape and position of the type-d and type-r features (at all polarization angles) are not significantly affected by the angle of incidence. (2) The incident angle θ only affects the center of the circle in k space of the type-s features following (k z = 0, k y = k 0 sin θ ). (3) The position of all the features (type-d, type-r, and type-s) rotates through the same angle as the polarization angle γ varies, while the rotating axis is (k z = 0, k y = 0) for the type-d and the type-r features, and (k z = 0, k y = k 0 sin θ ) for the type-s features. We note that the insensitivity of type-d features as a function of incident angle has been demonstrated experimentally in Ref. [42].

A. Linear polarization at oblique incidence
In addition to the ability to model arbitrary polarization states and incidence angles, further advantages of the FDTD approach over the efficacy theory include the ability to model topography-driven interpulse feedback effects [11,39] and to study surface topography in the spatial domain. The interpulse feedback mechanism has been shown as the main process driving the final rippled topography [37,[69][70][71]. The interpulse feedback effect here is taken into account by the socalled holographic ablation model (HAM) first introduced by Skolski et al. [11]: After the simulation of the first pulse, the topography of the sample is updated according to the calculated energy absorption, implying a change of material permittivity to 1 (material removal) when the local energy density exceeds the ablation criterion. This newly updated topography is used for the simulation of the next pulse, and the cycle is repeated until the last pulse. The ablation criterion in HAM is determined beforehand as equal to the calculated energy density in the topmost layer of a flat sample under the same irradiation conditions [11]. In HAM, any intensityrelated effect such as the intrapulse feedback is not considered. However, since the ablation criterion used in HAM corresponds to the removal of the smallest layer numerically possible (the FDTD grid size x in the vertical direction), the incident fluence in HAM can be interpreted as the ablation threshold under the simulation conditions.
Material permittivities supporting a shorter decay length of the surface plasmon waves has been shown to favor the growth of more regular LSFLs [62]. For this reason, the Note that in the spatial-domain images, only a central part (4 μm × 4 μm) of the whole computational domain is shown for a better visibility of the ripples. In the spatial-domain images, a gray linear color scale is used, where the brighter area represents mountains and the darker area represents valleys. This color scale is used throughout the paper for the spatial-domain images. permittivity value of Cr ( = −1.98 + 21.9i) is used for the following multipulse simulations on LSFLs. Another consideration in choosing metals instead of Si for the multipulse simulations is to allow the omitting of intrapulse feedback mechanisms. Furthermore, we expect the results of Cr to be qualitatively similar to Si when the excitation level permits a dielectric to metallic transition. Figure 5 shows the simulated surface topographies after a consecutive illumination of N = 1, 3, 7 pulses at normal incidence. After the final pulse, three distinct features in k-space can be identified. The type-s features are LSFLs with a period around 700 nm oriented perpendicular to the laser polarization. The initial size of the surface roughness particle used in the simulation is 80 nm × 80 nm × 10 nm with a filling factor of 0.1. The FDTD feedback cycles show a steady growth of LSFLs, which is consistent with observations from Refs. [38,72] that LSFLs 245430-7 mainly result from the scattering by nanoparticles with sizes around 100 nm. Comparing the spatial-domain images, the positions of the ridges and valleys of the LSFLs generated by a subsequent pulse is always aligned with those of the previous pulses, indicating a phase-locking mechanism that stabilize and amplifies the periodic structures. This is in consistent with experimental results utilizing high-resolution pulse-to-pulse imaging techniques [37,[69][70][71]. The type-g features gradually grow after the illumination of three pulses. In the N = 7 image, the type-g features are formed with an orientation parallel to the laser polarization. These structures are often referred to as grooves, and they have periodicities that can be much larger than the laser wavelength. In Ref. [40], two possible origins of the grooves are proposed: an electromagnetic origin in the ablation regime and a hydrodynamic origin by melt flow. The grooves that we discuss in this paper are those with electromagnetic origins. Ripples with a period that equal to half of the LSFLs are also predicted in the simulations. In this paper, we name these features in the frequency domain as type 2s. All three of these features have been observed experimentally, such as in Ref. [30], and we find strikingly good agreement in the k-space patterns. We note that one should distinguish the type-2s features from the type-r features. Although they are both HSFLs, their physical origins are different: The type-2s features have a period that is strictly half of the LSFLs while the type-r feature can have various periods far below λ/2 [2,8,[23][24][25][26][27][28][29]. The type-r features are triggered by near-field scattering [12] of the incident pulse and are strongly influenced by the size of the initial surface roughness [24][25][26]36,38]. Our results show that the type-2s features start to grow only after the LSFLs are already formed and deepened, and their periods are strictly half of the LSFLs. A careful inspection of the surface topography at N = 7 suggests a start of ridge splitting into two equally smaller ridges, supporting the finding of Ref. [73] that the redistribution of the electric field due to the already formed LSFL valleys are responsible for the type-2s features. Figures 6(a) and 6(b) show the simulated surface topographies and their Fourier transform amplitudes after consecutive illumination of seven pulses by s-and p-polarized beams for a range of incident angles. For p-polarized beams, the orientation of the LSFLs is always orthogonal to the laser polarization at all incident angles. We find that their periods can be well described by p = λ/(1.12 + sin θ ). The surface  [74]. If surface plasmon waves are the only physical origin of LSFLs, it implies a ripple period of p = λ/(λ/λ sp + sin θ ) = λ/(1.43 + sin θ ) [44]. In Ref. [35], the role of scattered quasicylindrical waves are revealed to be equally important as the surface plasmon waves in defining the ripple period, with the former obeying a linear dispersion relation. If only scattered cylindrical waves are considered, it implies a ripple period of λ/[1 + sin θ ]. Therefore, both the quasicylindrical waves and surface plasmon waves are involved in defining the LSFL periods shown in Fig. 6(b). In contrast to the p-polarized beams, the results of s-polarized beams show a development of LSFLs that are neither orthogonal nor parallel to the polarization direction at large angles of incidence. As shown in Fig. 6(a), at oblique incident angles, the surface topographies consist of ripples that are orientated at two different directions which are neither perpendicular to the laser polarization. This can be understood by analyzing the direction of the surface scattering. Figure 7(a) shows the in-plane component of the scatteredfield ( E 2 y,s + E 2 z,s , where the subscript s represents scattered field) by a single square-shaped roughness particle placed in the center of the y-z plane. The size of the particle is 80 nm × 80 nm × 10 nm. Thanks to the superposition principle, the scattered field is obtained by the difference in the field values calculated with and without the roughness particle. The simulation shows a preference of scattering parallel to the laser polarization. For an s-polarized beam, the intensity of the scattered light is symmetric along the two scattering directions. The in-plane wave vector of the incident light is perpendicular to the scattering direction. According to the phase matching condition: k i = G ± k s , where k i is the in-plane wave vector of the incident light, G is the wave vector of the formed ripples, and ±k s are the wave vectors of the scattered light, we have G = k i ∓ k s . Therefore, the wave vector of the interference pattern defining the ripple orientation, forms a non-90 • angle with respect to the laser polarization. The symmetric scattering of an s-polarized beam means that the intensity of the energy modulation is equally strong in the two scattering directions, forming ripples with two orientations. It is worth noting that their features in the spatial-frequency domain are consistent with the efficacy theory results shown in Fig. 2(b). However, it is not straightforward to analyze the spatial-domain result in the efficacy theory. Anomalous ripples that are neither parallel nor perpendicular to the laser polarization have been observed experimentally for spolarized beams at oblique incidence [18]. For both s and p polarizations, Fig. 6 predicts that the orientation of the grooves (type g) is parallel to the laser polarizations at all incident angles. Their periods are reduced at larger incident angles for p-polarized beams, while for s-polarized beams, their periods are not strongly influenced. Figure 6 also shows the change of the type-2s features as the incident angle is increased. For p-polarized beams, the type-2s features always have the same orientation as the corresponding type-s features but with frequencies that are doubled with respect to the type-s features. For s-polarized beams, the results are more complicated: at larger incident angles (for example, θ = 70 • ), the type-2s features retain a branch which is the same as the type 2s of the normal incident case. It also contains an additional branch which doubles the spatial frequency of the type-s features in both y and z directions. This seems to suggest that its physical origin is the same as the case of normal incidence: the redistribution of the electric field due to the already formed LSFL valleys. Figures 8(a) and 8(b) show the simulated surface topographies at oblique incidence when the polarization angle is in-between 0 • and 90 • . At large incident angles [ Fig. 8(b)], the orientation of LSFLs forms a non-90 • angle with respect to the laser polarization. The corresponding scattered field (by a single roughness particle) and the wave vectors are illustrated in Figs. 7(b) and 7(d), respectively. The scattering is the strongest along the polarization direction. As discussed earlier, we write the phase matching condition as G = k i ∓ k s . The nonzero angle between k i and k s means that the orientation of the resulting LSFLs is neither perpendicular nor parallel to the laser polarization. At oblique incidence, for any linearly polarized beam with a nonzero p component, the distribution of the scattered-light intensity exhibits asymmetry, with a stronger scattering in the direction facing the incoming light. Owing to this asymmetry, the energy modulation is stronger in one of the two scattering directions and, through the amplification of interpulse feedback, results in LSFLs with a single orientation. In the frequency domain, the Fourier transforms of the surface topographies obey the same rule as the energy absorption features of the first pulse as shown in Fig. 4(b): The distribution of the type-s features rotates through the same angle γ around (k z = 0, k y = k 0 sin θ ) as the polarization angle varies. Figure 8 also predicts that the orientation of the type-g features in the frequency domain is always rotated by the same angle γ around (k y = 0, k z = 0) as the polarization angle varies, even at mixed s-and ppolarization states.
Another striking result of the FDTD approach is the prediction of the type-r features [11,12,38]. In the spatial domain, these are HSFLs with an orientation perpendicular to the laser polarization. Their formation has been revealed to be strongly correlated with the nature of the surface roughness [25,26,36,38], and the interference between the incident light and the scattered near-field light has been proposed to explain their origin [12]. Nevertheless, there are few studies concerning their orientation and periodicity at oblique incidence. Figures 9(a) and 9(b) show the simulated surface topographies of HSFLs for s-and p-polarized beams under various incident angles. The dominance of the type-r features is evidenced in their Fourier transforms. The HSFLs perpendicular to the laser polarization are usually found on dielectric and semiconductor materials under low levels of excitations. This is evidenced by the fact that HSFLs are usually found at the periphery of a Gaussian beam [8,23,25,26]. Therefore, the permittivity used in these simulations is taken from the value of Si with a low level of photoexcitation n e = 4 × 10 27 m −3 . It is worth mentioning that in our simulations with n e = 4 × 10 27 m −3 , the type-r features are easier to be found than those with a higher free-carrier density such as n e = 8 × 10 27 m −3 . The dimension of the initial roughness particle is kept small as well, as a nanoscaled surface topography has been found to favor the growth of HS-FLs [25,26,36,38]. In the Fourier transforms, we see that the brightest part of the type-r features for all the incident and polarization angles lies on the same circles, indicating the same HSFL periods for all cases (or at least a very weak dependence on the angles). This can be understood as a direct consequence of the origin of HSFLs as near-field scattering: The wave number of the scattered near-field light is so much higher than the wave number of the incident light (approximately seven times as in Fig. 10), such that after the vector sum, the in-plane wave number of the incident light only has minor contribution to the period of the energy depositions. The argument also explains the observed orientations for different polarization angles as shown in Figs. 10(a) and 10(b): the orientation of the type-r features is solely determined by the polarization angle γ at all incident angles. This is different from LSFLs, in which case the orientation depends on both γ and θ (Fig. 8). We should note that the above argument may not be strictly generalizable to all materials and experimental conditions: for HSFLs with a period that is not much smaller than the laser wavelength, both its orientation and period may show weak dependence on the incident angle.
Another interesting aspect is the appearance of asymmetry in the topographical profile of the resulting ripples by beams with pure p polarization or mixed s and p polarizations under   . The topographical profile of HSFLs under those conditions resembles that of a blazed grating: The slope of the grating lines tends to be steep on the side which is facing the incoming light and to be gentle on the side which is away from the incoming light. This can be understood as the asymmetry in the scattered-light intensity distribution showing in Fig. 7(b): The scattered-light intensity is stronger in the direction facing the incoming light, resulting in asymmetry in the energy depositions. In Ref. [46], the authors find similar blazed gratings on ripples with periods larger than the laser wavelength (grooves) at large oblique angles. A simple explanation by using geometrical optics is given in Ref. [46]. Our results suggest that blazed ripples may form even for subwavelength HSFLs produced by p-or mixed s-and p-polarized beams. We note that in such cases, geometrical optics does not hold and a full vectorial modeling of laser energy deposition is required.

B. Circular, radial, and azimuthal polarizations
The full numerical nature of the FDTD approach determines that the method is highly versatile. By defining the field components of the incident beam, it is in principle possible to implement any complex polarization state, including those with spatially and/or time-varying polarization directions. Among those beams, circularly polarized and cylindrical vector beams [75] are of particular interest. In Ref. [49], triangular LSFLs in hexagonal arrangements were fabricated by using circularly polarized fs beams, and the authors have demonstrated a superior superhydrophobic performance of the textured surface compared to the LIPSSs fabricated by linearly polarized beams. Other surface topographies such as chains of spherical nanoparticles [51,52], speckled and disordered structures [53,54], and nanoscaled fine-dot structures [27] by using circularly polarized beams are also reported. Nevertheless, whether those structures can be explained within the electromagnetic frame remains an unanswered question. Figure 11 shows the simulated LSFL and HSFL surface topographies after the illumination of 20 circularly polarized fs pulses at normal incidence. The LSFLs are obtained with the permittivity value of Cr and initial roughness particles of larger lateral dimensions (80 nm × 80 nm × 10 nm). The HSFLs are obtained with a permittivity value corresponding to excited Si with n e = 4 × 10 27 m −3 and roughness particles of smaller dimensions (40 nm × 40 nm × 10 nm). For LSFLs, the simulated surface topography is in strikingly good agreement with the scanning electron microscopy (SEM) images shown in Ref. [49] at a low number of applied circularly polarized pulses on stainless steel, which show a growing of triangular-LSFL in hexagonal arrangements. The simulated HSFLs topography shown in Fig. 11(c) is quite similar to the SEM image of nanoscaled fine-dot structures on diamondlike carbon films fabricated by circularly polarized fs pulses [27]. The agreement between the simulations and the experimental data strongly suggests an electromagnetic origin of those structures obtained with circularly-polarized beams. Figure 12 shows the simulated topographies of LSFLs and HSFLs for beams with radial and azimuthal polarizations. The LSFLs here belong to type-s features observed for the linearly polarized beams, while the HSFLs here belong to the type-r features. Both types of LIPSSs are shown to be with orientations that are perpendicular to the local electric field vector of the laser beams, thus revealing the local complex polarization states of the applied beams. We note that LIPSSs with such characteristics have been observed experimentally for both LSFLs [55][56][57][58] and HSFLs with subwavelength periods [59].
Circular polarization (including elliptical) is implemented in the FDTD approach by adding a phase delay between the electric-field components of the incident beam: where η determines the phase delay. When η = π/2, the polarization becomes circular. Radial and azimuthal polarizations (at normal incidence) are implemented by determining the position-dependent electric field components across a doughnutlike amplitude distribution [75]. For example, for radial polarization: E 0y = E 0 cos[arctan(z/y)], E 0z = E 0 sin[arctan(z/y)], E 0 = Ar exp(−r 2 /w 2 ), where r = z 2 + y 2 , A determines the amplitude of the incident beam, and w is the beam waist.

V. SUMMARY
In this paper, we have extended the FDTD modeling of LIPSSs to include oblique incidence and arbitrary polarization states. The inhomogeneous energy absorption below a material's rough surface at oblique incidence is calculated by the FDTD method and compared to the analytical efficacy theory. The two approaches show excellent agreement. The surface topographies of both LSFLs and HSFLs under oblique incidence and various polarization states are obtained in the same theoretical framework. Upon irradiation of p-polarized beams, the simulated LSFL period decreases as the incident angle increases and the orientation is always perpendicular to the laser polarization. While upon irradiation of s-polarized or mixed s-and p-polarized beams, LSFLs that are neither perpendicular nor parallel to the laser polarization are found at large incident angles. The origin is explained by the nonzero angle between the in-plane wave vector of the incident light and the scattered light in these cases. The simulations of HSFLs with linearly polarized beams show that their periods have almost no dependency either on the incident angle or the polarization angle, while their surface topographies resemble that of blazed gratings at large incident angles for p-polarized and mixed s-and p-polarized beams. The origin is attributed to the asymmetry in the scattered near-field light. Simulations with circularly polarized beams show the development of triangular LSFLs and nanoscaled fine-dot HSFLs which are both in very good agreement with the experimental results found in literature, strongly suggesting their electromagnetic origins. LSFLs and HSFLs with radial and azimuthal polarized beams are also obtained in the simulations, in which the orientation follows the local polarization direction of those beams. Our results suggest that the FDTD approach on LIPSSs is a highly versatile and powerful technique that goes well beyond the efficacy theory and, when combined with multiphysical models [39,40], may pave one of the foundations for a complete understanding and modeling of LIPSSs under various irradiation conditions.

ACKNOWLEDGMENTS
This research was supported by the European Research Council (ERC Starting Grant No. 637476) and the Netherlands Organization for Scientific Research (NWO). This work was conducted at the Advanced Research Center for Nanolithography, a public-private partnership between the University of Amsterdam, Vrije Universiteit Amsterdam, the Netherlands Organization for Scientific Research (NWO), and the semiconductor-equipment manufacturer ASML.