Introduction
Laser irradiation of solid surface by ns and fs pulses at medium and high power densities generates a variety of structures like surface undulations, ripples, periodic waves, wrinkles, surface roughening, etc. The study of their evolution elucidates the mechanism of lasermatter interaction depending on the laser and material parameters and the interaction conditions (Bonse et al., Reference Bonse, Baudach, Kruger, Kautek and Lenzner2002; Tran et al., Reference Tran, Lam, Zheng, Murukeshan, Chai and Hardt2005; Wang and Guo, Reference Wang and Guo2005; Liu et al., Reference Liu, Jiang, Yang, Guan and Dai2011; Reif et al., Reference Reif, Varlamova, Varlamov and Beestehorn2011; Graf et al., Reference Graf, Kunz, Engel, Derrien and Muller2018). It was shown that ripples of the wavelength comparable to the laser wavelength on semiconductors and dielectrics induced by circularly polarized light originate from the excitation of surface plasmons/polaritons (Varlamova et al., Reference Varlamova, Costache, Ratzke and Reif2007). The subwavelength ripples on silicon formed by the IR laser pulse result from the synergy of the electron excitation and the capillary wave formation – basically the hydrodynamic mechanism (Tsibidis et al., Reference Tsibidis, Barberoglou, Loukakos, Stratakis and Fotakis2012). nteraction between the oscillating plasma (actually the ion acoustic waves in plasma) and molten surface was identified as the origin of the ripple formation on metals (Huang et al., Reference Huang, Liu and Zhu2011). However, the superwavelength ripples on silicon formed by the UV laser pulses originate from the solitary waves induced by the hydrodynamic instability (Lugomer et al., Reference Lugomer, Maksimovic, Geretovzsky and Sorenyi2013).
Regarding the hydrodynamic mechanism, the formation of surface structures may be associated or induced by the RayleighTaylor, RichtmayerMeshkov, or other type of instabilities (Lugomer Reference Lugomer2016, Reference Lugomer2017). Here we show that laser-generated RayleighTaylor Instability (RTI) induces formation of the new paradigm of surface structures.
The RTI evolves from the acceleration of a lowdensity fluid ρΗ into highdensity fluid ρL (or vice versa), as a consequence of baroclinic generation of vorticity at the perturbed interface. The evolution of RTI depends on the Atwood number, A = (ρH ρL)/(ρH + ρL) and on the initial conditions. It follows the exponential law in time with the growth rate that depends on the perturbation amplitude and the spatial period of the perturbation mode.
Initial conditions
The initial perturbation serves as the seed of the interface instability. Shock interaction with a narrowband combination of high-frequency modes triggers growth of a turbulent mixing layer in the linear regime. Later on the evolution enters the nonlinear regime when the growth rate slows down. This occurs purely via mode coupling of the wavenumbers when modes stay in some defined phase relations (Zhou et al., Reference Zhou, Cabot and Thornber2016). Abarzhi and collaborators (Abarzhi et al., Reference Abarzhi, Cadjan and Fedorov2007; Abarzhi and Rosner, Reference Abarzhi and Rosner2008; Pandian et al., Reference Pandian, Stellingwerf and Abarzhi2017) have shown that in such multimode perturbation growth of coherent periodic structures occurs as a result of interference of perturbation modes for the specific phase and amplitude relations. However, when a broad range of perturbation wavelengths is used, it happens that short wavelength modes are grouped around some short common wavelength, while the long wavelength ones are grouped around larger common wavelength and cause formation of structures in the band of certain width (Lugomer, Reference Lugomer2017). The density interface transforms into the structures which can be periodic or aperiodic, small but also the large ones (Lugomer, Reference Lugomer2017) or into complex wavevortex structures (Zabusky et al., Reference Zabusky, Lugomer and Zhang2005; Lugomer, Reference Lugomer2007).
Mixing width
The growth of structures that follows after initial perturbation is triggered by formation of a turbulent mixing layer and the mixing process which relates to the mass, momentum and energy. In RTI, the ratio of the spike/bubble growth is sensitive to the particular definitions of the mixing mass and mixing width (e.g. concentration threshold, product width, or integrated width and their evolution in time (Zhou and Cabot, Reference Zhou and Cabot2019). The mixing width can be also defined as visual width of the mixing zone – based on the identification of structures that are separated from diffuse areas of strong fluctuations and «noise», those above the identification threshold.
Mixing mass and Atwood number
An important issue in the mixing process is the A number because it affects the mixing mass. Namely, in multimode perturbation, different initial modes cause different growth rates of spikes and bubbles as well as different mixing mass. The mixed mass M, defined as 
$M = \smallint 4{\rm \rho} Y_1Y_2dV$ (Y 1 and Y 2 are the mass fractions and ρ is the mixed density), is lower fo higher Atwood number (Zhou and Cabot, Reference Zhou, Clark, Clark, Glendinning, Skiner, Hantington, Huricane, Dimits and Remington2019). The study of dependence of the mixing mass on the Atwood number between A = 0.2 and 0.8 by Zhou et al. (Reference Zhou, Cabot and Thornber2016) has shown that M nonlinearly increases with time, being higher for the lower A. It is possible, the mixed mass of two fluids of densities ρ1 and ρ2 may be taken as appropriate measure in multimode perturbation experiments.
Examples of the physical systems showing the RTI extend from the micro- to meso- and to the astrophysical scale (Dimonte, Reference Dimonte1999). The collection of research and review papers relating to turbulent mixing, non-equilibrium processes from atomistic to astrophysical scales by Abarzhi et al. (Reference Abarzhi and Sreenivasan2010, Reference Abarzhi, Gauthier and Sreenivasan2013a, Reference Abarzhi, Gauthier and Sreenivasan2013b) and by Zhou (Reference Zhou2017a, Reference Zhou2017b) give a detail view into the models, simulations, experiments and interpretation.
We consider the RTI and formation of wrinkle structures in planar geometry by multipulse laser irradiation of metal targets. Generated in the pulse overlapping area due to the decrease of the melting threshold they appear after N ≥ 5 pulses (silicon) (Lugomer et al., Reference Lugomer, Maksimovic, Farkas, Geretovzsky, Toth, Zolnai and Barsonyi2011), after N = 20 (fused silica and borosilicate glass) (Ben-Yakar and Byer, Reference Ben-Yakar and Byer2004), and even after N ≥ 2 (indium). Since the melting threshold of indium is decreased after the first pulse, the RTI evolution is generated after only two pulses. For infinitesimally thin layer, the RT interface has hypocycloidal shape. Multimodal perturbation of the target surface causes pulsations of the shockmomentum and affects the evolution of wavy structures in the C-zone of the spot. Solidified after the pulse termination they become wrinkles – henceforth called the wavewrinkles.
This paper describes the evolution of wavewrinkle structures associated with the RTI in the plane of the target surface at the liquid/vaporplume interface for medium laser energy density. The target surface irregularity (in the overlapping area of two pulses) causes radial and angular variation of the fluid layer thickness and velocity so that the wavewrinkles are different in various spatial domains. Their heterogeneous morphology ranges from wrinkles formed by the solitary waves to the solitary waves without exponential tail called compactons and to the hierarchical series of smaller wrinkles called wrinklons. We formulate the conceptual frame assuming that the nonlinear (2 + 1) KadomtsevPetviashvili (KP) equation – and the equations formulated in the KP sense – describe all local structures. These nonlinear evolution equations (NLEE) are different with respect to the terms like dissipation, dispersion, and their relation (whether nonlinearity and dispersion are balanced or not).
Motivation for this study is the elucidation of the wavewrinkle generation in the C-zone of the RTI which evolves into the new paradigm of the nanoscale wavewrinkle structures in the plane of target surface. We also show that the wavewrinkle structures in the environment of RTI bubbles are different from that of RTI spikes. The wavewrinkles in the bubble environment are discussed in this paper (paper I), while those in the spike environment will be discussed in paper II.
Outlines of the experiment
The experiments were performed in the open configuration in which the target is directly irradiated. The sample was situated in the gas chamber and irradiated in the presence of air at the pressure of P 0 = 1 atm. Irradiation (unpolarized) was performed by two pulses of Gaussian-like power profile by Q-switched ruby laser (τ = 30 ns; λ = 628–693 nm; spot size r ~850–900 µm, E~150 mJ; Es ~6.5–7.0 J/cm2). Indium plates of 1 cm × 1 cm × 0.1 cm (the melting point T M = 429 K and boiling point T B = 2345 K) were used as targets. chematic of the experimental setup can be found elsewhere (Lugomer, Reference Lugomer2016).
The phase explosion causes the formation of vapor/plasma plume spheroid which starts fast expansion and acceleration of molten metal. The horizontal shock wave generates the RTI in the plane of the target surface with a wavylike shape of spikes and bubbles. The lateral vapor/plasmaplume expansion above the molten layer in the ambient gas is shown in Lugomer (Reference Lugomer2016), and details are described in its Supplementary. The twolayers of fluids which move horizontally at the same velocity are formed: a molten metal layer of density ρ1 and vaporplasma plume layer of density ρ2. It is usually assumed that the density interface between ρL and ρH layers is just a zero-thickness plane, , the mathematical interface. Although this is true for the «low temperature system», for the high temperature system (such as generated by lasermatter interaction), the situation is more complex. We assume that the boiling surface layer of thickness h 1 is separated from the vapor/plasma layer by the density interlayer of thickness 
${h}^{\prime}_1$ (Fig. 1a). Thus, the thickness of fluid layer in which the waves are excited is not h 1 but h; h = (h 1 + 
${h}^{\prime}_1$). The displacement of the fluid layer from z = 0 is η* (Fig. 1b).
Fig. 1. Twolayers of fluids which move horizontally with the same velocity. (a) olten metal layer of density ρ1 and the vaporplasma plume layer of density ρ2 formed during laser interaction. The boiling surface layer of thickness h 1 is separated from the vapor/plasma layer by the density interlayer of thickness 
${h}^{\prime}_1$. The thickness of fluid layer in which the waves are excited is not h 1 but h; h = (h 1 + 
${h}^{\prime}_1$). (b) The parameter h is the thickness of the molten layer and η* is the displacement of the layer from the level z = 0.
The velocity of accelerated fluid in the C-zone is (sub)sonic U ≲ 1500 m/s and the kinematic viscosity of liquid In is ν ≤ 5 × 10−6 (m2/s). The thickness of the fluid layer which is the highest near the RTIfront (h ≲ 10 µm), decreases with distance. The corresponding Re number, Re = hU/ν decreases from Re ~3 × 103 to the values of <103 in distant zones. The Re number is below the critical value (Recrit ~104) for the wave rollup, and the KelvinHelmholtz (KH) rolls (vortex filaments) do not form. However due to inhomogeneous RTI it happens that in the vicinity of some RTI bubbles/spikes the KH may be formed.
The RTI front and the wavewrinkle structure become solidified by the end of the laser pulse in the fast cooling process so that a posteriori study by the scanning electron microscope (SEM) JEOL reveals their final morphology.
Conceptual frame
Generation of the wavewrinkles occurs due to the perturbation of the interface of twolayers of fluids by the gravitational perturbation in front of the RTI spikes and bubbles which have nonlinear evolution. For the irregular target surface (uneven bottom after the first pulse) and variable fluid layer thickness, the C-zone is inhomogeneous. Uneven bottom in nonlinear regime and the small amplitude assumption play significant role in the evolution of the velocity field. Variable bottom topography requires one to take into account not only the amplitude a and wavelength of the waves, λ, but also the order of amplitude of the variations of the bottom topography, the wavelength of the bottom variations, the reference depth and the nonlinearity parameter (Israwi, Reference Israwi2010). escription of such system requires introduction of the new class of differential equations – with necessity to confirm relevance of the steps in the solution procedure– what becomes mathematically challenging problem.
To simplify the problem – we assume as the first approach – that the variable surface topography does not change much over small areas (local domains) of about few tens of micrometers in size. We look for the evolution of the surface instability from the perturbation amplitude, d, of the fluid layer of thickness, h, (h = h 1 + 
${h}^{\prime}_1$) (Fig. 1a). For a weak nonlinearity (d < h), the wave equation can be derived from general equation for the molten layer bounded by a solid base from below and by the vapor/plasma layer from above. The wave equation follows from the continuity equation and boundary conditions at the interface for viscous liquid (Yang, Reference Yang2012; Lugomer et al., Reference Lugomer, Maksimovic, Geretovzsky and Sorenyi2013). For the small thickness of the molten layer – the “shallow fluid layer” which gives rise to the nonlinear waves is the most adequate approach.
Consider a layer of twodimensional fluid over an even bottom (Fig. 1a). The surface elevation is η*(x, t*) measured form the undisturbed level z* = 0. The domain of the fluid D(η*) is defined as D(η*) = {(x*,z*) : −h < z* < η*(x*, t*)} (Fig. 1b). Assuming that the particle velocity can be expressed as the gradient of velocity potential ϕ(x*, z*), u(x*, z*) = Δϕ(x*, z*) with Δϕ = 0 in D(η*) and with the boundary condition at the bottom ∂ϕ/∂z* = 0 at z* = −h, where t*, x*, y*, and z* are the dimensional time the space variables. As the prototype, for the shallow fluid layer thickness h, of the surface elevation η* above the zero-level, the equation for the waves propagating in the x-direction can be written (Infeld et al., Reference Infeld, Senatorski and Skorupski1994, Reference Infeld, Senatorski and Skorupski1995; Berger and Milewski, Reference Berger and Milewski2000; Oikawa and Tsuji, Reference Oikawa and Tsuji2006):
$$\equalign{&{\rm \partial} /{\rm \partial} x{^\ast} [{\rm \partial} {\rm \eta} {^\ast}/{\rm \partial} t{^\ast} + c_0(1 + 3{\rm \eta} {^\ast}/2h){\rm \partial} {\rm \eta} {^\ast}/{\rm \partial} x{^\ast} } \cr&+{ c_0h^2/6 (1-3S/{\rm \rho} gh^2){\rm \partial} ^3{\rm \eta} {^\ast}/{\rm \partial} x{^\ast}^3] + c_0/2\cdot {\rm \partial} ^2{\rm \eta} {^\ast}/{\rm \partial} y{^\ast}^2 = 0.$$Using the transformation
$$\eqalign{& u = \pm 3{\rm \eta} {^\ast}/2h,\quad x = \pm {\rm \beta} /h\lpar {x{^\ast}-c_0t{^\ast}} \rpar, \cr & t = {\rm \beta} c_0t{^\ast}/6h,\quad {\rm \beta} = [ \pm (1\ndash 3S/{\rm \rho} gh^2]^{-1/2},} $$where ρ fluid density, S surface tension, g acceleration of gravity, while t, u, and x are dimensionless time and space variables, one gets the equation (Oikawa and Tsuji, Reference Oikawa and Tsuji2006):
$$(u_t + u_{xxx} + {\rm 6}uu_x{\rm )}_x + {\rm 3}{\rm \sigma} ^2u_{{\rm yy}} = 0.$$ The last term, u yy, is added to the above 1D equation and gives the (2 + 1) KadomtsevPetviashvili (KP) equation for the waves on a shallow fluid layer. Assuming – for the domain D(x,y) – periodic boundary conditions u(x, y, t) = u(x + L x, y + L y, t), the solution of the above equation can be expressed as a summation of Fourier components in the form 
$F{(k_x, k_y) {\rm exp} ({ik_x}x + {k_y}y-{\rm \varpi} t)}$. The plane–wave solutions for the phase variable kx +my – ωt satisfy the dispersion relation (Kao and Kodama, Reference Kao and Kodama2012):
$${\rm \varpi}(k_x,k_y) = - (3{k_x}^2 / {k_y} + {k_x}^3 ) $$where 
$k = \sqrt{{kx}^2 + {ky}^2 }$ is the wave number.
Depending on the signature of dispersion, σ2, the solutions of this equation can be (1) stable nonlinear waves for the negative dispersion (σ2 = −1), described by the KP-II equation, and (2) unstable nonlinear waves for the positive dispersion (σ2 = +1), described by the KP-I equation. The KP-II equation describes the dynamics when gravity effects dominate over surface tension and gives stable wave solutions (Berger and Milewski, Reference Berger and Milewski2000). The KP-I equation describes dynamics when surface tension dominates over gravity effects causing the instability and decay of waves. Which of the above equations is relevant for description of the waves in particular case is determined by the criterion (Berger and Milewski, Reference Berger and Milewski2000; Oikawa and Tsuji, Reference Oikawa and Tsuji2006):
$$S/(\Delta rgh^2) ⋚ 1/3 $$where S is the surface tension at the interface of the liquid layer and the vapor layer, h is the thickness of the fluid layer, g is the acceleration of gravity. The fluid density ρ is actually the fluid density at the interface and written as Δρ is the difference in density of the two phases (the liquid molten layer and the vapor layer).
In the first case, the expression (5) is <1/3, and the gravity effects dominate over surface tension so that stable cnoidal and line-solitary wave solutions exist. In the last case, the expression (5) is >1/3, and the surface tension dominates over gravity effects so that cnoidal and line solitons become unstable. In the intermediate case (=1/3), the onset of complex higher order hydrodynamic instability is possible (Ablowitz and Clarkson, Reference Ablowitz and Clarkson1992; Berger and Milewski, Reference Berger and Milewski2000; Chakravarty and Kodama, Reference Chakravarty and Kodama2009; Song and Liu, Reference Song and Liu2012).
Assuming that the surface tension S at the interfacial layer (between the vapor/plasma plume and hot liquid metal) is very low, S ~10−6 N/m, the fluid layer thickness h ≤ 10 µm g ~10 m/s2, and Δρ (of the interfacial layer) ~6.5 × 103 kg/m3, one finds S/(Δρgh 2) = 0.153 < 1/3, meaning that the KP equation gives stable solution.
The irregular target surface and variable fluid layer thickness affect the fluid layer velocity, nonlinearity and dispersion in various domains of inhomogeneous fluid layer. The SEM analysis reveals that domains are well separated (for a distance l ≥ λ – mostly having no common boundaries) and selforganized in radial direction (along the extension of the bubblefront). Since the fluid layer thickness decreases from ~10 µm (near the RTI front) to h ≤ 1 µm (at the periphery), the conditions for the wavewrinkle evolution change in the radial direction. Assuming the shallow fluid layer, weak nonlinearity, and weak dispersion – the (2 + 1)-dimensional KadomtsevPetviashvili equation and the equations formulated in the KP sense – are appropriate for the simulation of structures in spatially separated domains. In this respect we relay on the solutions of NLEEs derived by various mathematical groups. For a shallow fluid layer a 2 + 1 description based on KP equation is the most appropriate. Owing to the fact that the flow conditions vary from one domain to another – which are well separated and show different wave structures – the basic KP equation modified by taking the new terms into account (characteristic for particular domain) may describe the solitary wavewrinkles of various characteristics.
Before continuing, we should shortly consider the question how these waves become solitary waves and how solitary waves become standing waves that are frozen at the end of laser pulse.
Traveling waves become solitary waves
The wavewrinkles appear as traveling waves and become solitary waves. Although it is known that transition of traveling into solitary waves occurs when the traveling waves are dominated by group velocity dispersion (rather than diffusion) and by nonlinear frequency shift (rather than nonlinear saturation) (Dudley et al., Reference Dudley, Peacock and Millot2001; Kuriakose and Porsezian, Reference Kuriakose and Porsezian2010; Liu and Dodin, Reference Liu and Dodin2015), this transition is not well understood in dynamical systems.
Solitary waves become standing solitary waves
The intriguing phenomenon in dynamical systems is the transition of solitons into stationary solitary waves. Generally, when two solitary waves of the same frequency and amplitude traveling through a medium with the same speed but in opposite direction superimpose on each otherthey give rise to a standing solitary wave (Wai et al., Reference Wai, Chen and Lee1989; Aziz, Reference Aziz2011; Tofeldt and Ryden, Reference Tofeldt and Ryden2017 Supplementary A. Such case occurs when the incoming solitary wave is reflected from the wall, but in the fluid system the reflection from inhomogeneity has the same effect on the incoming solitary wave. We assume that inhomogeneities of a fluid layer coincide with the boundaries of local domains and have important role in the transition of solitary waves into the stationary wavewrinkle soliton patterns.
These experimentally obtained, solidified, stationary solitary wavewrinkle patterns in particular domains are juxtapositioned with the solitary wavewrinkles generated by mathematical simulation.
Results and discussion
The surface morphology formed by two overlapping pulses on indium target is shown in Figure 2a. The SEM micrograph of Gaussian-like spot reveals scalloped circumference at the leftside and the heataffected zone (HAZ) with small chaotic structures in the C-zone (Fig. 2b). The rightside (overlapped area with strong melting) shows RTIfront and the wavewrinkles in the C-zone (Fig. 2c). For infinitesimally thin layer the RT interface has hypocycloidal shape. For the thicker layer it is 3D but one can assume it a quasi-2D because the RTIfront still resembles the hypocycloidal curve which in the case of multimodal perturbation appears as two or more intersecting hypocycloidal-like RTI fronts (Fig. 2c). The evolution of the RTI bubble-front (red) is associated with formation of different wavewrinkle structures in the columns A and B due to inhomogeneity of the C-zone of the spot (Fig. 3).
Fig. 2. SEM micrograph of the ruby laser spot on Indium target irradiated in the air at the atmospheric pressure. (a) Two laser shots have being fired forming an overlapping zone. (b) The leftside of the second shot showsirregular scalloped edge and the heataffected zone (HAZ) with small irregular structures of indium oxide particles. (c) The rightside of the second shot forms the overlapping zone with the first shot. The RayleighTaylor instability (RTI) with smooth wavy-like circumference of prominences (spikes) and walleyes (bubbles) is formed by horizontal fluid acceleration. The wavewrinkle structures formed in the circumferential zone (C-zone) follow the shape of the RTI spikes and bubbles in the plane of target surface.
Fig. 3. Schematic presentation of the wavewrinkles that evolve in the C-zone in front of the RTI. The enlarged segment shows the wavewrinkle morphology which changes with distance from the RTI front: the C-zone is divided into ear field (NF) edium Field (MF) and the far field (FF). Red arrows: The outflow of a hot low-density fluid and wavewrinkle structures formed in front of the bubbles. Blue arrows: The inflow of a cold high-density fluid in front of the spikes. The C-zone is inhomogeneous due to irregular target surface, variation of the fluid velocity and radial and angular variation of the fluid layer thickness between few microns and few hundreds of nanometers. The column (A) of the wavewrinkles is formed on a «thick» (micron thick) molten metal layer, while the column (B) is formed on a «thin» nanofluid layer. Notice that wavewrinkle structures are different in these two columns. The amplitude and wavelength of the wavewrinkle structures decrease with distance from the RTIfront from NF to MF and to FF.
Column A: The wavewrinkles are formed on relatively thick fluid layer. The structures in the NF, MF and FF zones are similar to the structures observed in various laser experiments on solid targets (Gorodetsky et al., Reference Gorodetsky, Kanicki, Kazyaka and Meicher1985; Zabusky et al., Reference Zabusky, Lugomer and Zhang2005; Lugomer, Reference Lugomer2007; Trtica et al., Reference Trtica, Gakovic, Radak, Batani, Desai and Bussoli2007) and shall not be discussed in this paper.
Column B: The structures formed on a thin fluid layer are different from those reported in the literature and reveal the new wavewrinkle paradigm. This new paradigm evolves from the inhomogeneous environment that afects the formation of nonlinear wavewrinkles. Not only their wavelength, amplitude and the profile but also their nature change with distance from the RTIfront; the C-zone may be divided into the ear field (NF) (close to the front of RTI), edium field (MF) and the far field (FF) waverinkle morphology (enlarged in Fig. 4).
Fig. 4. SEM micrographs of the wavewrinkle structures in the «thin fluid layer» of few hundreds of nanometer. (a) The NF wavewrinkle morphology shows rectangular horseshoe structures in front of the RTI bubble are bended following the curvature of the emerging bubble. (b) The MF morphology reveals few sets of wavewrinkle structures.
Near-field zone: Figure 4a shows sharphorseshue parabolic-like rectangular waves formed in front of the RTI bubble and the sharpoblique rectangular waves.
Medium-field zone: Figure 4b shows sharp rectangular linewaves, smooth rectangular aperiodic and the smooth rectangular periodic linewavewrinkles.
Far-field zone: Figure 4c shows periodic wavewrinkles. They mostly appear as the small set(s) of waves which correspond to the early phase of evolution.
By the end of the pulse an ultrafast cooling wave is formed which moves from the periphery of the C-zone toward the RTI boundary, freezing the wave wrinkles in various phases of their evolution. The waves become solidified after pulse termination in the time between ~100 ns and few μs stay as permanently frozen wrinkles with characteristics of the former waves.
Wave wrinkles on nanofluid layer
Near-field zone
Sharp horseshue and sharp oblique rectangular waves. The SEM micrograph of the NF zone (Fig. 5) reveals sharp rectangular «horseshue» waves. The asymmetry of the RTI bubble and the irregular target surface cause the asymmetry of the «horseshue» waves that make robust compact pattern. The fact that such wave pattern emerges in the parabolic-like form leads to conjecture that analogous dispersive systems support localized patterns (Liu and Dodin, Reference Liu and Dodin2015). They move without change of shape with different velocity: the first «horseshue» wave is followed by the next one of the higher amplitude which moves faster what indicates solitary waves. Their rear-side collision leads to jammed configuration. At some distance from the jamming zone the sharp rectangular waves become more stright and oblique with respect on the symmetry axis of the RTI bubble. The angle of the oblique waves increases with distance until they become transversal to the axis of the bubble. They interact forming the «X»- and «Y»-configuration characteristic for the nonresonant and resonant soliton interaction, respectively (Fig. 5).
Fig. 5. SEM micrograph of the NF zone showing the formation of the wavewrinkles in front of the asymmetric RTI bubble. The «horseshoe» solitary wavewrinkles of rectangular profile generated in front of the RTI bubble show jamming and the rear-side collision. The two sets of the oblique rectangular wavewrinkles collide and merge. The yellow frame (down left) comprises the set of rectangular unidirectional waves outside the jamming and collision zone. The rectangular singlelinesoliton, as well as the interaction of the solitary waves at the leftside that form «Y»type configuration called «resonant» and other one that form the «X»type configuration called «nonresonant», can be seen.
Medium-field zone
The MF zone (Fig. 4b) shows wavewrinkles divided into the left-side and the right-side systems. By the laser pulse termination the cooling wave causes solidification of a molten layer. Becoming thin elastic sheet it starts to shrink causing the lateral tension and breakup of wavewrinkles. Figure 6 shows the left and the rightside of different wavewrinkles localized in the rectangular domains – (A), (B) and (C). They form large-scale coherent structures (waves) inside the cells (local domains) – without the small scale structures. The flow is well ordered in the domains close to the RTI front and disordered in distant domains. The rectangular domains are selforganized into a “domain network” pattern. ormation of domains with the long range order should be associated with spontaneous symmetry breaking that is usually accompanied by topological defects such as domain boundaries (Park et al., Reference Park, Chai, Chai, Lee, Jia, Park, Lee and Park2014). Tentative explanation may be that fluid inhomogeneities form the rectangular cells that tend to organize into 2D “domain network” with p2mm symmetry [Fig. 7 (inset)]. However, the symmetry breaking transforms it into an irregular network which still resembles the short range p2mm organization.
Fig. 6. SEM micrograph of the MF zone showing various kinds of the wavewrinkle structures in local domains (yellow frames) A, B, and C. Domains are well separated (for a distance l ≥ λ – having no common boundaries) and selforganized in radial direction – along the extension of the bubblefront and comprise different kinds of waves: (A) Sharpperiodic rectangular waves – the compacton-like solitons; (B) periodic rectangular waves with rounded top surface; (C) eriodic rectangular waves with rounded top surface.ormation of the selforganized cascade of wrinklons (pper left).
Fig. 7. Selforganization of rectangular domains (Voroni cells) of inhomogeneous fluid in a domain network of p2mm symmetry. The formation of domains with long range order should be associated with spontaneous symmetry breaking that is usually accompanied by topological defects such as domain boundaries. Wavewrinkle direction in all domains (x-direction) is transversal to the direction of their propagation that coincides with direction of the fluid acceleration and direction of the RTI bubble-front propagation (red arrow; y-direction).
Constant wrinkle profiles are caused by homogeneously well-organized waves inside the domain. Wavewrinkle direction in all domains (x-axis in Fig. 7) is transversal to the direction of their propagation that coincides with direction of the fluid acceleration and direction of the RTI bubble-front propagation (y-axis in Fig. 7). Notice that wavewrinkles do not meet at the domain boundaries because every domain is well separated from the other one what justifies separate treatment of wavewrinkles in the domains A, B and C.
Sharp rectangular wave wrinkles: compacton-like structures
Domain A: The frame A in Figure 6 – enlarged in Figure 8a – shows robust unidirectional sharp rectangular line-solitary waves of equal width with sharp vertical front that suddenly decreases to zero. Their profile (nset to Fig. 8a) indicates the compacton-like solitons which are defined as solitons with the absence of infinite wings (Rosenau and Hyman Reference Rosenau and Hyman1993; Rosenau, Reference Rosenau2000, Reference Rosenau2005; Wazwaz, Reference Wazwaz2005; Li and Song, Reference Li and Song2014). If the field u(x) becomes equal to zero at some point x = x 0, the dispersion mechanism shuts off and the soliton remains stable, preserving its shape from spreading out. Really, Figure 8 shows the equal width solitons without exponentially decaying tail; instead the field vanishes identically outside some finite interval. For such solitons, the function is zero outside of a compact set.
Fig. 8. Domain A: Compactonlike soliton wavewrinkles. (a) SEM micrograph of the sharp rectangular wavewrinkles (in the yellow frame (a) in). Inset: Profile of the compacton-like wavewrinkles along the section line. (b) Sharp rectangular traveling waves (compacton-like) obtained as the explicit solution of the generalized KadomtsevPetviashvili Modified Evolution Equation (KP-MEW). These nonlinear waves can be juxtapositoned with the rectangular waves in Figure 6a. (Courtesy of M. Khalique; Reprinted with permission of Hindawi Publishing Corp., from Adem and Khalique (Reference Adem and Khalique2015); Copyright Hindai Publishing Corp, 2015).
From the mathematical aspect (asuming conditions given in the Conceptual frame), the equal-width compacton solitary waves may be generated by the equation formulated in the KP sense (Wazwaz, Reference Wazwaz2005; Li and Song, Reference Li and Song2014; Adem and Khalique, Reference Adem and Khalique2015)
$$(u_t + {\rm \alpha} (u^n)_x + {\rm \beta} u_{xxt})_x + {\rm \gamma} u_{yy} = 0$$where α, β, γ and are real-valued constants. Taking n = 2 the above equation includes 
$u_x^2$ (quadratic) nonlinear convection term and the term u xxx that relates to the dispersion effects. Eq (6) is known as generalized (2 + 1) KadomtsevPetviasvili Modified EqualWidth Equation (KP-MEW).
For detail mathematical procedure see Adem and Khalique (Reference Adem and Khalique2014, Reference Adem and Khalique2015). The mathematical procedure is schematically presented in Supplementary B and C.
The solution u(x,y,t) which describes sharp rectangular solitary waves is obtained in the form of the trigonometric function (Adem and Khalique, Reference Adem and Khalique2015)
$$\eqalignb{u(x,y,t)& = A_0 + A_1\lsqb {-{\rm \lambda} /2 + {\rm \delta}_2(-C_1\sin ({\rm \delta}_2z } \cr&\quad + C_2\cos ({\rm \delta}_2z))/(C_1\cos ({\rm \delta}_2z) + C_2\sin ({\rm \delta}_2z)) \rsqb \cr&\quad + A_2\lsqb {-{\rm \lambda} /2 + {\rm \delta}_2(-C_1\sin ({\rm \delta}_2z)} \cr&\quad+ C_2\cos ({\rm \delta}_2z))/(C_1\cos ({\rm \delta}_2z) + C_2\sin ({\rm \delta}_2z) \rsqb } $$where z = t − x − ((a + c)/b)y with (a = b = c = 1); 
${\rm \delta} _2 = (1/2)\sqrt {4{\rm \mu} -{\rm \lambda} ^2}$; λ and μ are constants and C 1 and C 2 are arbitrary constants. The solution (7) shows a set of compacton-like waves (Fig. 8b). The generalized KP-MEW equation is solved under boundary condition 
${\lim_{{\rm \xi \to} \mp \infty}} \,u(\rm \xi ) = a$ (amplitude). A regular compacton solution corresponds to the case a = 0 (Zhong et al., Reference Zhong, Tang, Li and Zhao2014). In that case the dispersion operator becomes degenerate at the front, and its degeneracy generates sharp fronts supporting the formation of robust compact localized patterns.
The waves are developed between two boundaries of the omain A (fluid inhomogeneities) that caused their reflection under an angle ~180° thus making posible interference with the original solitary wave. The size of the domain, l ≥ 800 nm, establishes the interference conditions for the formation of only three compacton-like waves of the wavelength λ ~200–250 nm. The compacton-like waves obtained as the exact solution of the KP-MEW equation in Figure 8b can be succesfully juxtapositioned with the standing solitary periodic waves of the omain A in Figure 8a.
Rectangular waves with the rounded top surface
Aperiodic rectangular waves with kinematic dispersion
Domain B: The aperiodic wavewrinkles in the frame B (Fig. 6) – enlarged in Figure 9a – are the smooth rectangular waves with rounded top surface and the amplitude a ~220 nm ee the inset to Fig. 9a. They may be assumed the long waves with the wavelength which gradually decreases from λ ~350 to ~250 nm. The profile of the waves indicates that characteristics of flow are somewhat different in the omain A and that besides dispersion associated with the term u xxx, the other dispersion mechanisms are present. Generally, the parameters like velocity and dispersion coefficient vary from one domain to another, namelymorphologic (target surface inhomogeneity) and hydrodynamic dispersion play a role in the evolution of flow structures. However, for the Domain B we assume that the target surface morphology affects the flow of the fluid parcels of the molten layer that travel the same distance by different pathways, by different velocities, in different times that give rise to the kinematic dispersion (Saco and Kumar, Reference Saco and Kumar2002). Therefore, the total dispersion arises from the contribution of the linear hydrodynamic (D LH), morphologic (D M) and nonlinear kinematic (D NK) dispersion: = D NK + D LH + D M (Saco and Kumar, Reference Saco and Kumar2002). (For details see Supplementary D) The contribution to total dispersion <D> from a kinematic dispersion overcomes the contribution of other two components. One can define the nondimensional kinematic dispersion D NK+ = D NK/D M ≡ δ (Saco and Kumar, Reference Saco and Kumar2002; Snell et al., Reference Snell, Sivapalan and Bates2004), which has two components, δ1 and δ2 for the x- and y-directions in the flow basin (Domain B), respectively.
Fig. 9. Domain B: Aperiodic rectangular wavewrinkles with rounded top surface. (a) SEM micrograph of aperiodic rectangular wavewrinkles (in the yellow frame (b) in Fig. 6. Inset: Profile of the wavewrinkles with kinematic dispersion along the section line. (b) Aperiodic traveling unidirectional waves with rectangular profile and rounded top surface obtained as the exact solution of the generalized KadomtsevPetviasviliBenjaminBonaMahony (KP-BBM) equation with presence of kinematic dispersion (δ1 = 0.1, δ2 = 0.2), expressed in the Weierstrass ℘-elliptic functions for v = 1.5, with α = β = γ = 1, B 1 = B 2 = 0.5 at t = 0. (Courtesy of A. Das; Reprinted with permission of Elsevier from Ganguly and Das (Reference Ganguly and Das2015). Copyright, Elsevier, 2015).
The longwaves with rounded top surface in the Domain B can be generated by the nonlinear equation formulated in the KP sense which incorporates kinematic dispersion coefficients δ1 and δ2 (Ganguly and Das, Reference Ganguly and Das2015)
$$(u_t + u_x-{\rm \alpha} (u^2)_x-{\rm \beta} u_{xxt})_x + {\rm \gamma} u_{yy} = {\rm \delta} _1u_{xxx} + {\rm \delta} _2u_{yyy}$$where terms u t + u x describe the propagation of wave and (u 2)x is nonlinear term. Equation (8) is known as the generalized (2 + 1) KadomtsevPetviasviliBenjaminBonaMahony (KP-BBM) equation.
Assuming periodic boundary conditions, the solution of Eq. (8) can be expressed as a summation of Fourier components in the form 
$F(k_x,k_y) = \exp (ik_xx + ik_yy-{\rm \varpi} t)$ and searched in the form of the traveling wave
$$u(x,y,t) = u({\rm \varsigma} );\quad {\rm \varsigma} = B_1x + B_2y-{\rm \nu} t$$The parameters B 1 and B 2 are related to the inverse width of the soliton in x- and y-directions, respectively, and v represents the traveling wave velocity.
After few transformations and integrations one finds the solution of the KP-BBM equation with dispersion δ1, δ2 > 0. The exact traveling wave solution is derived in terms of the elliptic Weierstrass We-function (Ganguly and Das, Reference Ganguly and Das2015)
$$\eqalign{u(x,y,t)& = \displaystyle{1 \over {2{\rm \alpha} B_1^2}} \left[ {6\displaystyle{{{(B_1^3 {\rm \delta}_1 + B_2^3 {\rm \delta}_2)}^2} \over {25B_1^3 {\rm \nu}}} \pm \sqrt {\displaystyle{{{(6{(B_1^3 {\rm \delta}_1 + B_2^3 {\rm \delta}_2)}^2)}^2} \over {{(25B_1^3 {\rm \nu} )}^2}}-4{\rm \alpha} B_1^2 \Lambda ({\rm \delta}_1,{\rm \delta}_2,{\rm \nu} )}} \right] \cr &-\exp \left( {\displaystyle{{2(B_1^3 {\rm \delta}_1-B_2^3 {\rm \delta}_2)} \over {5B_1^3 {\rm \nu}}} \varsigma} \right)\cdot {\rm We}\left[ {C_4 + \displaystyle{{5B_1\sqrt {{\rm \alpha \nu}}} \over {\sqrt {6B_1} (B_1^3 {\rm \delta}_1 + B_2^3 {\rm \delta}_2)}}} \right]\cdot \exp \left( {\displaystyle{{B_1^3 {\rm \delta}_1-B_2^3 {\rm \delta}_2} \over {5B_1^3 {\rm \nu}}} {\rm \varsigma}, 0g_3} \right)} $$where g 3 and c 4 are constants of integration; g 3 is constant of integration that depends on the boundary conditions, and c 4 = c 2 + c 3 (c 2 is also some constant of integration). The solution obtained for δ1 = 0.1, δ2 = 0.2, v = 1.5 with α = β = γ = 1, B 1 = B 2 = 0.5, in t = 0, shown in Figure 9b are the waves with variable wavelength which have the rectangular-like profile and rounded top surface.
The reflection of solitary waves from both «wall»-boundaries (inhomogeneities) in the fluid layer under an angle ~180° makes possible interference with the original wave. The interference of the solitary waves in this domain of the size l ~1500 nm establishes the conditions for the formation of four stationary waves. These waves (Fig. 9a) – with the wavelength which increases in radial direction of the spot from λ ~350 to ~500 nm – are developed between two boundaries (fluid inhomogeneities) of the Domain B and can be successfully juxtapositioned with the simulated waves in Figure 9b.
Periodic rectangular waves without kinematic dispersion
Domain C: The periodic rectangular wavewrinkles with rounded top surface in the frame C (Fig. 6) – enlarged in Figure 10a – have the profile shown on the inset to Figure 10a. Such morphology of periodic waves indicates even different flow conditions from the Domain B. They may be attributed to the kinematic dispersion which becomes small or zero in the omain C. Therefore, the generation of such waves can be described by Eq. (8) taking δ1 ~δ2 ≃ 0
$$(u_t + u_x-{\rm \alpha} (u^2)_x-{\rm \beta} (u_{xxt})_x + {\rm \gamma} u_{yy} = 0$$With assumption of periodic boundary conditions the solution can be expressed as a summation of Fourier components 
$F(k_x,k_y) = \exp (ik_xx + ik_yy-{\rm \varpi} t)$ (Ganguly and Das, Reference Ganguly and Das2015) and also searched in the form of traveling wave, 
$u(x,y,t) = u(x + y-{\rm \nu} t) = u({\rm \varsigma} )$. One takes 
${\rm \varsigma} = B_1x + B_2y-{\rm \nu} t$ and ν = B 1 + (γB 2)2/B 1 with parameters B 1 and B 2 which are related to the inverse widths of the soliton in the x- and y-directions, respectively. To find the class of solutions that travel with the velocity v, requires some additional assumptions to be made. Then, the KP-BBM equation is transformed into ODE the solution of which can be expressed by the Jacobi elliptic functions. The exact traveling wave solution of Eq. (12) may be written in the explicit form For detail mathematical procedure see Ganguly and Das (Reference Ganguly and Das2015):
$$u(x,y,t) = B\left[ {1-\sqrt 3 \displaystyle{{1\mp cn\left[ {\displaystyle{{\sqrt {2B{\rm \alpha}}} \over {\root 4 \of 3 \sqrt {B_1^2 + {\rm \gamma} B_2^2}}} \left( {B_1x + B_2y-\displaystyle{{B_1^2 + {\rm \gamma} B_2^2} \over {B_1}}t} \right),k} \right]} \over {1 \pm cn\left[ {\displaystyle{{\sqrt {2B{\rm \alpha}}} \over {\root 4 \of 3 \sqrt {B_1^2 + {\rm \gamma} B_2^2}}} \left( {B_1x + B_2y-\displaystyle{{B_1^2 + {\rm \gamma} B_2^2} \over {B_1}}t} \right),k} \right]}}} \right]$$where B is integration constant. The solution for v = 0.5, with α = β = γ = 1, B 1 = B 2 = 0.5 at t = 0 is shown in Figure 10b that can be successfully juxtapositioned with the experimental unidirectional standing solitary periodic waves of the omain C in Figure 10a. As in the previous cases, the perpendicular reflection of solitary waves from the «wall»-domain boundaries (inhomogeneities of the fluid layer) makes possible interference with the original solitary wave. For the waves with the wavelength of λ ~300–400 nm, this domain of the size l ~1000–1200 nm establishes the interference conditions that form fourstanding solitary waves.
Fig. 10. Domain C: Periodic rectangular wavewrinkles with rounded top surface. (a) SEM micrograph of periodic rectangular wavewrinkles (in the yellow frame (c) in Fig. 6. Inset: rofile of the wave wrinkles without kinematic dispersion along the section line. (b) Periodic traveling waves with rectangular profile and rounded top surface obtained as the exact solution of the generalized KadomtsevPetviasviliBenjaminBonaMahony (KP-BBM) equation without kinematic dispersion (δ1 = δ2 = 0), expressed in the elliptic Jacobi functions for v = 0.5, with α = β = γ = 1, B 1 = B 2 = 0.5 at t = 0. (Courtesy of Prof. A. Das; Reprinted with permission of Elsevier from Ganguly and Das (Reference Ganguly and Das2015). Copyright, Elsevier, 2015).
Formation of wrinklons
By the end of pulse the ultrafast cooling which starts at the periphery of the C-zone moves toward the RTI boundary. Ultrafast solidification transforms the liquid layer into an elastic sheet which starts to shrink at the liquid/solid phase transition. Lateral component of surface tension (tension stress ɛ) on the wavewrinkles causes their breakup and formation of the left and rightside of the solitary wavewrinkles [Fig. 6 (upper left)]. The ends of broken wrinkles on the elastic sheet became fixed by solidification at the point (line) that represents constraint boundary for the lateral tension field. A strong tension on the elastic sheet causes focusing of energy and formation of the cascade of wrinkles at the smaller wavelengths called wrinklons (Fig. 11). This is similar to wrinklons on hanged curtains which are fixed at the top point (line) and exposed to tension by gravitational field. The formation of selfsmilar wrinklon structures has been observed on thin elastic sheets of various materials including graphene (Vandeparre et al., Reference Vandeparre, Pinerua, Brau, Roman, Bico, Gaz, Bao, Lau, Reis and Damman2011; Meng et al., Reference Meng, Geng, Yu, Dou, Nie and He2013).
Fig. 11. Formation of wrinklons as selfsimilar hierarchical selforganized structures in the cooling phase after laser pulse termination at the constraint boundary of the solidified sheet. SEM micrograph of wrinklons from the yellow frame in Figure 6 (upper left) shows few scales of the wrinklon cascade. Inset: Schematic representation of the wrinklon formation with merging of two wrinkles of the wavelength λ into the wrinklon with 2λ.
A single wrinklon corresponds to the localized transition zone needed for transformation of two wrinkles of wavelength λ into a larger one of width 2λ [Fig. 11 (inset)]. This transition requires distorsion of the layer sheet which relaxes over distance L, meaning that wrinklon can be characterized by the length L. A series of wrinklons of different L but selfsimilar are organized in the same way being mostly parallel establish the wrinklon hierarchy. However, their parallelism is disturbed due to the variation of the local fluid velocity, dispersion, variation of the layer thickness etc (Fig. 11). Generally, the wrinklons follow the scaling laws different for thin and thick sheets analogous to the «light» or the «heavy» curtains.
The analysis of Figure 11 shows linear dependence of the wrinklon amplitude A on the wavelength λ, A = Qλ (Fig. 12a), similar to 
$A = {\rm \lambda} \sqrt \Delta$ (where Δ is the relative change of the elastic sheet) of Vandeparre et al. (Reference Vandeparre, Pinerua, Brau, Roman, Bico, Gaz, Bao, Lau, Reis and Damman2011), Deng and Berry (Reference Deng and Berry2015) and Meng et al. (Reference Meng, Geng, Yu, Dou, Nie and He2013). Dependence of the wrinklon's relaxation length L on the wavelength which follows the power law, L = Kλ3/2 (K ≃ 0.4 nm2/3) (Fig. 12b), is also in agreement with the behavior of wrinklons on the elastic sheets. Dependence of the length, L, on the amplitude follows the power law L = CA 1/2 (C ~ 81 nm1/2), characteristic for the “heavy curtains” (Fig. 12c). Such interdependence of the wrinklon parameters A, L, and λ is a mimic of the heavy curtain behavior under tension generating the hierarchy of wrinklons.
Fig. 12. Diagrams showing relation between characteristic wrinklon parameters. (a) The wrinklon amplitude versuswavelength for the cascade of wrinklons show linear dependence: A = Qλ. (b) The length L of the wrinklon energy equilibration versus wrinklon wavelength: L (nm) = Cλ3/2 (C ~0.4 nm2/3). (c) The length L of the wrinklon energy equilibration versus wrinklon amplitude: L (nm) = KA 1/2 (K ~81 nm1/2).
Relative change of the length, Δ, of the elastic sheet along the y-axis occurring in the cooling process (which is equivalent to our case) can be estimated by Meng et al. (Reference Meng, Geng, Yu, Dou, Nie and He2013)
$$\Delta y = \int_{T_0}^{T_1} {{\rm \alpha} _G(T)dT} $$where αG is the thermal expansion coefficient of material, T 0 and T 1 are initial and final temperature during the cooling process, respectively. Here we use the argument of Meng et al. (Reference Meng, Geng, Yu, Dou, Nie and He2013) that the temperaturedependent expansion coefficient of material generates an effective force per unit length of the sheet
$$F = Eh\int_{T_0}^{T_1} {{\rm \alpha} _G(T)dT} $$where E is the Young's module and h is the thickness. Assuming T 1 as the temperature at which the indium fluid layer solidifies and becomes an elastic sheet (T 1 ~430 K), and T 0 the room temperature (T 0 = 298 K) at which the cooling ends, one finds the temperature interval in which the wrinklons are formed, T 1 –T 0 = 132. During ultrafast cooling (between ~100 ns and ~1 µs) this temperature interval is passed in few hundreds of nanoseconds.
The expression for total energy after minimization gives the length L and its dependence on the wavelength λ, and the thickness, h, of the elastic sheet (Meng et al., Reference Meng, Geng, Yu, Dou, Nie and He2013)
$$L = \lsqb {3\Delta y(1-{\rm \nu}^2)} \rsqb ^{1/2}{\rm \lambda} ^2/h{\rm \pi} $$where ν is the Poisson ratio (ranging for great number of materials between 0 and 0.5). Estimating the average layer thickness h ~300 nm, the shrink distance in the ydirection Δy ~250 nm, and taking the wrinklon wavelength from the micrograph (Fig. 11), we use the relation (15) to estimate L for various wrinklon wavelengths. For the wrinklons with the wavelength λ ~160 nm, h ~300 nm and ν ~0.3, one finds L ~809 nm, while for λ ~220 nm, one finds L ~1326 nm. These values are in very good agreement with the values of L for the corresponding wrinklon wavelength in the diagram (Fig. 12b), what confirms that the wrinklon generation occurred in the cooling process after pulse termination.
Conclusion
We have shown that multipulse laser-induced RTI in the plane of target surface is associated with formation of the new paradigm of the wavewrinkle structures in the circumferential zone of Gaussian-like spot. Multipulse laser irradiation causes inhomogeneity of the surface fluid layer, fluid velocity, dispersion and the target surface variation so that formation of the nonlinear wavewrinkles become different in various local domains. The superposition and growth of the density interlayer perturbation with the background fluid motion cause the formation of the large-scale coherent structures (waves) inside the cells (domains) – without occurrence of the smallscale structures. The flow is well ordered (coherent) in the domains close to the RTI front and is disordered (incoherent) in distant domains. Domains are organized into a “domain network” with p2mm symmetry and oriented in radial direction of Gaussian-like spot, in the direction of the fluid motion. Wrinkles which appear as the nonlinear traveling waves become transformed into solitary waves due to balance of the group velocity dispersion and nonlinear frequency shift. The inhomogeneity of the fluid layer – which establishes the domain boundaries – causes the reflection of solitary waves and interference inside the domains. Solitary waves become transformed into stationary soliton wavewrinkle patterns – different in various domains. The heterogeneous solitary wavewrinkle morphology varies ranging from the solitary waves, compacton-like solitons, to the aperiodic rectangular waves (withrounded top surface) and to the periodic ones. The structure of the velocity fields and the flow symmetry determine the family of solutions with distinct local properties characteristic for the local domains. Consequently, all types of waves may be successfully juxtapositioned with the exact solutions of the nonlinear differential equations formulated in the KadomtsevPetviashvili sense taking into account the conditions in local domains of inhomogeneous fluid layer.
The wavewrinkle structures experience ultrafast cooling after laser pulse termination. The cooling wave that starts at the periphery travels toward the center causing sudden solidification and transformation of a thin molten layer into an elastic sheet which starts to shrink generating lateral tension on the wrinkles. This causes their breakup and formation of the cascade of smaller wrinkles called wrinklons. Wrinklons, called the elementary excitations of condensed matter in the form of elastic sheets and membranes that appear under tension stress, manifest some unique characteristics regarding the relation between their amplitude, equilibration length and the wavelength.
Regarding the experimental conditions and eventual control of the traveling and solitary wavewrinkles formation in lasersolid interactions, it may be said that they can be formed at the medium- and higher-energy densities, Es (below the plasma detonation threshold), but not on the very high ones. Namely, at very high Es the plasma detonation causes strong horizontal fluid acceleration and jetting along the target surface instead of the formation of traveling and solitary waves.
In the future work we will try to shed more light on the A dependence of the RTI growth.
Supplementary material
The supplementary material for this article can be found for this article can be found at https://doi.org/10.1017/S0263034620000105.
Acknowledgments
grateful to Dr. G. Baranovic, “Rudjer Boskovic” Institute, Zagreb, for critical reading of the manuscript and to Dr. Ye Zhou, Lawrence Livermore National Laboratory, Livermore, USA, for valuable suggestions. The author is specially grateful to Prof. S.I. Abarzhi, Carnegie Mellon University, USA, and The University of Western Australia, Perth, Australia, for very important comments and suggestions.
Financial support
This work has been supported by the European Union through the European Regional Development Fund the Competitiveness and Cohesion Operational Programme (KK.01.1.1.01.0001).
