Solitary Nanostructures Produced by Ultrashort Laser Pulse

Laser-produced surface nanostructures show considerable promise for many applications while fundamental questions concerning the corresponding mechanisms of structuring are still debated. Here, we present a simple physical model describing those mechanisms happened in a thin metal film on dielectric substrate irradiated by a tightly focused ultrashort laser pulse. The main ingredients included into the model are (i) the film–substrate hydrodynamic interaction, melting and separation of the film from substrate with velocity increasing with increase of absorbed fluence; (ii) the capillary forces decelerating expansion of the expanding flying film; and (iii) rapid freezing into a solid state if the rate of solidification is comparable or larger than hydrodynamic velocities. The developed model and performed simulations explain appearance of microbump inside the focal spot on the film surface. The model follows experimental findings about gradual transformation of the bump from small parabolic to a conical shape and to the bump with a jet on its tip with increasing fluence. Disruption of the bump as a result of thinning down the liquid film to a few interatomic distances or due to mechanical break-off of solid film is described together with the jetting and formation of one or many droplets. Developed theory opens door for optimizing laser parameters for intended nanostructuring in applications.


Background
Fabrication of surface structures and/or nanodroplets is a popular direction in very promising laser applications. Such structures are used for surface-enhanced Raman spectroscopy (SERS) and tip-enhanced Raman spectroscopy (TERS) [1]. The enhancement of a Raman signal is based on the following. (i) The sub-wavelength structure/tip, like a frozen jet above the top of a bump, is produced. (ii) The structure is illuminated by weak light. (iii) There are spots of sub-wavelength volume (hot-spots) near the tip where orders of magnitude enhancement of a weak illuminating electromagnetic field takes place. (iv) The Raman effect is approximately proportional to the fourth power of the field. Thus, the Raman signal from a single molecule caught into the hot-spot can be enhanced *Correspondence: nailinogamov@gmail.com 1 Landau Institute for Theoretical Physics, Russian Academy of Sciences, Chernogolovka, Russia Full list of author information is available at the end of the article up to the detectable level. In the similar way, using nanostructures and hot-spots, the photoluminescence from a single molecule under illumination by weak light is measured [2]. In another limit of strong illuminating laser pulse, the surface nanostructures amplify (perhaps again as a result of the hot spot enhancement) electron and X-ray emission from the surface [3].
Two-four ray interference of an ultrashort pulse is applied to produce a two-dimensional lattice of solitary bumps or holes on a thin film [4,5]. Using UV laser pulse (second harmonics of femtosecond laser 785 nm), the lattice spacing between neighboring nanobumps on the thin gold film may be done less than a micrometer (spatial period of a lattice is 760 nm in [6]); exactly this type of bumps are considered in this paper. These plane two-dimensional "crystals" enhance absorption of light by dye [7]. This is important for the super sensor technique like the surface-enhanced infrared absorption (SEIRA) solving similar problems as SERS. The grating array is a fundamental component in nanophotonics. It can include metamaterials which show a negative refractive index and harmonic conversion [5,8].
The bumps and tips are used in nanophotonics/plasmonics as antennas for operations with light [1,9]. The bumps or solitary undersurface frozen bubbles can be utilized for production of two-dimensional photon crystals. Such bubbles were found and studied in [10][11][12]. Surface nanostructures are also employed for colorizing, tribology, and manipulations with wettability [13][14][15] as well as in bio-applications [16].
Thermomechanical separation of film from substrate described below is a real reason for the bump formation. Higher absorbed energy F c in Eq. (1) results in jetting from the top of bump. By increasing F c either a part of the jet, or whole jet, or whole jet together with a part of bump disintegrates into flying droplets. Thus, the production of droplets can gradually change from a single droplet to a few and then many droplets with increasing F c .
Droplets are characterized by distributions dn/dm, dn/ dE, dn/d throughout mass, kinetic energy, and an angle of escape. The simulations described below show how to adjust, e.g., the flying angle of droplet or to regulate the number of droplets. This is possible in the case of solitary nanobump with small size R L ∼ 1 μm in (1). For larger bump size, the conditions are less controllable. New observation given in the paper concerns the appearance of solid debris in an ejecta cloud; see Fig. 7 below. Their presence is caused by inertial stretching and mechanical break-off of solid or solidified part of shell forming a bump.
In this work, we study dynamics of a gold film of 40-100 nm thick on a silica substrate triggered by fast heating with ultrashort laser pulse. A sequence of involved physical processes leading to the formation of nanobump is discussed below. After several picoseconds of electron-ion relaxation, the pressure waves generated by such heating can result in spalling separation of film and its flying away from the substrate after the acoustic time t s = d f /c s determined by the film thickness d f and speed of sound in metal c s . Fluence of the heating laser beam in the plane of its cross-section can be represented by the spatial Gaussian distribution where in-plane radius r is measured from the beam axis. Therefore, the heating of film surface is nonuniform along the cylindrical radius r. As a result of such heating, the initial (here it means after separation from a substrate) velocity distribution v(r, t = 0) in material flying away from the substrate has a maximum v c (t = 0) = v(r = 0, t = 0) at the beam axis. Thus, the separated film has the dome-like shape which inflates with time, and volume of an empty cavity between the separated film and the substrate increases during inflation.
Typical initial flight velocities are in the range of v c (0) ∼ 30 − 200 m/s. The inflation stage can last from a few to several tens of nanoseconds if the diffraction-limited micrometer-sized laser focal spots (R L ∼ 1 μm) are considered. Capillary forces acting along the warped flying film decelerate inflation of a dome. Capillary deceleration of a bulging dome focuses mass flow along the dome shell in direction to its axis. This results in the formation of an axial jet and droplet on the top of the dome.
Here, our new simulation results and comparisons with experiments are presented. Those results in particular explain the formation of nanocrowns and appearance of debris in a form of frozen droplets lying on surface of irradiated spot. We demonstrate that the former is a consequence of break-off of the liquid part of the dome, and the latter is a result of the capillary return of droplet, respectively.

Previous theoretical suggestions
Main difficulty in understanding (solved recently in [22]), how the jet from the dome top appears, was around the focusing problem. It was obvious from the experiments with jetting that accumulation of mass of a heated spot takes place [17,23,24]. The material flows radially into the small central region forming a jet. This inflow causes strong thinning of a film outside to the central region. The thinning has been measured experimentally [24]. The scheme of experiments with a film/substrate target is shown in Fig. 1. A laser beam distributed according to (1) hits a target normally to its surface. The illumination produces a hot focal spot having a form of thin film disk marked out by the red rectangle in Fig. 1.
According to the previous suggestions, the hot disk expands laterally along a contact surface with velocities v rad as a result of heating. This expansion motion collides with and reflects from a cold surrounding film indicated by blue color in Fig. 1. The expansion and reflection velocities correspond to the red and blue arrows in Fig. 1. The reflected wave propagates back to the center. Then, the cumulation of the reflected wave in the center strongly amplifies the wave. Collapse of the wave in the center produces vertical velocities. Thus, the central jet appears according to this consideration.
Amplification indeed can have a place in a cylinder thick enough. If we replace in Fig. 1 the thin red disk by the infinite vertical red cylinder placed into an infinite rigid blue tube, then amplification will happen, because the expansion flow induced by an initial pressure maximum at the axis can propagate outward, reflect from the rigid tube, and converge back to the axis. But for the thin films of interest, the aspect ratio R L /d f is large, where R L ∼ 1000 nm is a beam radius, d f ∼ 40 − 100 nm is the thickness of a film. For such condition, the reflection is weak, and the pressure wave cannot propagate long distance through shallow liquid with a free boundary condition p = 0 on the upper film boundary.
In our simple model, the initial (i.e., after film separation) radial velocities v rad , see Fig. 1, are small relative to velocities normal to the contact v norm /v rad ∼ R L /d f . Origin of the velocity field v norm (r, t = 0) is connected with a vertical film/substrate hydrodynamic interaction through the underlying contact [22]; see Fig. 1. Velocity distribution v norm (r, t = 0) mainly follows the fluence distribution (1) if the film/substrate adhesion is weak. It is known that a gold film used in experiments [17,23,24] is weakly coupled to a glass substrate. We have v norm (r, 0) ∝ F(r) (1) everywhere inside the circle r < r sep with the exception of the vicinity of ring r = r sep . Velocity v norm (r sep , 0) = 0 at the ring for the small but finite adhesion, while the distribution F(r) (1) exponentially damps out but remains non-zero.
Our model presented in [22] considers the following processes and phenomena: (i) two-temperature electronion relaxation stage, (ii) melting, (iii) film/substrate hydrodynamic interaction, (iv) formation of vertical velocity field (shown by the black arrows in Fig. 1) as a result of this interaction and separation of film from substrate (contact rupture) due to tensile stress coming to the contact with an acoustic wave, (v) early stage of inflation of bump h(r, t) = v norm (r, 0) t following the initial vertical velocity distribution v norm (r, 0), (vi) deceleration of vertical velocity v norm (r, t) < v norm (r, 0) and focusing to the axis v rad (r, t) < 0 of the liquid material, which forms a curved shell of bump as a consequence of the capillary forces acting tangentially along the curved surface, (vii) formation of an outward jet and an inward counter-jet, and decay of the outward jet into droplets.
Here, we omit discussion of above items and instead present development of the model based on Monte Carlo modeling of electron heat conductivity in material simulated by a classical molecular dynamics code. This allows us to take into account the strong electronic thermal transport in the case of metal and thus to describe both the spreading of heat along a film and the very fast cooling/crystallization of hot melt in the bump shell. The cooling solidifies molten metal and arrests the further destruction of a bump.

Classification of experiments and corresponding simulations
There is a list of experimental parameters: F c , R L (they characterize a laser beam (1)), d f , Z f , Z s , σ , χ, T m (they characterize a target), here Z f , Z s are acoustic impedances Z = ρ c s of a film and a substrate, σ , χ are coefficients of surface tension and electron thermal diffusivity of metal χ = κ/c, κ and c are heat conduction and capacity, T m is a melting temperature of metal. Let us create important combinations of these values to narrow the parametric space. The processes of capillary deceleration and focusing from the one side and the cooling and freezing from the other side are of primary significance. Therefore we introduce the capillary and thermal numbers here v c (t = 0) is velocity in the tip of a future bump at the instant of separation of a film from a contact, see Fig. 1, ρ is density of metal. The non-dimensional numbers v 0σ , v 0χ (2) give strengths of the capillary and freezing effects relative to initial inflation velocity v c (0). The films used in experiments are thinner than thickness of a heat-affected zone d T , d T ≈ 140 nm for Au [25]. Thus, a film is uniformly heated along a normal direction at the stage of separation from a contact. This is the initial stage relative to the much later stage of deceleration of a bump. Initial increase of temperature is  (3) should be omitted. It is important that the pulse is ultrashort and that spread of absorbed heat from a skin layer is supersonic at the two-temperature stage [22,26]: duration t T of heating of a film in the normal direction is shorter than acoustic time scale t s = d f /c s . The latter defines the duration of mechanical film/substrate separation process. In these conditions, the pressure distribution created by a sudden temperature rise (3) is here ≈ 2 is Gruneisen parameter and T in [kK] is given by (3). This is an isochoric pressure increase.
Running acoustic waves are Expansion of heated (3) and suddenly pressurized (4) gold into the glass substrate produces a weak shock in glass and a rarefaction wave in gold. Pressure and velocity at the contact film/substrate are continuous, therefore where u cb and p cb are values at the contact boundary and p is given by expression (4). Pressure p cb (5) is positive, the contact moves in a substrate direction. This continues up to the instant 2t s when a wave sent from a contact at t ≈ 0 will reflect from a vacuum boundary of a film and will return back to the contact. Pressure We take values ρ f = 19.3 g/cm 3 , ρ s = 2.2 g/cm 3 , c s | film = 3.3 km/s, c s | substr = 3.9 km/s for gold (f ) and glass (s) when we calculate impedances and coefficients in these expressions.
At the instant t = 2t s , the expansion half-period 0 < t < 2t s finishes and the contraction half-period begins. In the linear acoustic approximation (valid if pressures are much less than a bulk module and motion is subsonic) the contact pressure p cb and velocity u cb given in (6) change their sign to the opposite sign during the contraction half-period. If contact adhesion is weak then a contact break-off at the beginning of a tensile stage, i.e., at the beginning of the second half-period, is done. In this case, a momentum of a film accumulated during the first halfperiod defines separation velocity (3); compare velocities v = 2p cb /Z f and u cb = p cb /Z s (5), v ≈ v cb /3.7 for the gold/glass film/substrate target. Comparing this approximate expression for v with the two-temperature hydrodynamics results, we see that the expression for v underestimates separation velocity ≈2 times. This follows from the table presented in [22], p. 24. The corrected expression for the separation velocity is Using velocity (7) as v c (0), we obtain values for the nondimensional parameters (2). For the capillary number, we have where it is supposed that σ = 1000 dyne/cm. Significant vertical displacement of a film shown in Fig. 1 begins in the point r m where the film is molten. Therefore, we suggest that the radius of appreciable deviation of a film from substrate is r m and we will regard below this radius as the separation radius r sep = r m . For gold (T m = 1.337 kK), the melting threshold on absorbed fluence according to (3) is compare with table composed for d f = 60 nm gold in [22]. Expression follows from the distribution (1) and formula (9). Let us calculate the thermal (2). From Eqs. (7) and (10) A beam radius (11); we take χ ≈ 0.5 cm 2 /s for molten gold to calculate v χ .
The map of numbers v 0σ , v 0χ is presented in Fig. 2. The curves 0.5 and 2 shows parametric dependence v 0σ (F c , d f = 60 [ nm] ) (8) and v 0χ F c , d f = 60 [nm] , R L (11) as function of the central fluence F c as a running parameter varying in the range F m < F c < 3F m for the radii R L equal 0.5 μm (curve 0.5) and 2 μm (curve 2). If F c = F m , then r sep = r m = 0. In this case, the thermal velocity v χ = χ/r m is infinitely large, while thermal number v 0χ = v c (0)/v χ is zero. Therefore, the curves 0.5 and 2 in Fig. 2 start from the horizontal axis. The hole in a film appears for the higher fluences F c > 3F m and regime of separation can change because internal spallation of a film begins.
Main new achievement used in the paper here (relative to the previous work [22]) is addition of the Monte-Carlo (MC) algorithm to the molecular dynamics (MD) algorithm [27]. The sense of the MC modeling is the following. The MC electron jumps from its host atom to the neighbor of the host thus transporting heat from hot to cold place along the electron subsystem. The quantity ν D gives frequency of these jumps per electron. The rate of electron-ion kinetic energy exchange is defined by two parameters: the second frequency ν ei and effective mass of electron m ef . The value ν ei is electron-ion collisional Fig. 2 The curves 0.5 (R L = 0.5 μm) and 2 (R L = 2 μm) present estimates of location of experimental parameters for a gold film d f = 60 nm on a glass substrate. Typically, they cover the range 2F m < F c < 3F m around the dashed grey curve separating the experimental non-ejecting (they are closer to the origin point) and ejecting regimes. The circles shows the simulation runs. The grey continuous curve separates non-ejecting (the blue circles) and ejecting cases (the red circles). We see that the simulations rather satisfactory correspond to the experiments frequency. Higher frequencies and mass m ef correspond to the stronger electron-ion thermal coupling. We achieve the necessary value for the thermal diffusion χ varying the triple of parameters ν D , ν ei , and m ef . The coefficients χ used in simulations are listed in Table 1. Number of millions of atoms is given in the second column. Monte Carlo (MC) algorithm allows to vary electron heat conduction κ and hence the heat diffusion coefficient χ. We adjust χ to correspond to the experimental values of v 0σ and v 0χ (see Fig. 2), while having smaller spatial scales in a simulation. The rectangle of the simulation box has a square cross-section L box × L box in the plane of a target. Values L box and d f are given in [nm] Decreasing diffusivity χ in comparison with an experimental value, we create smaller numbers of atoms N atom in computational system approximately equivalent to the real experimental system. This is equivalently relative to the capillary and thermal numbers (2). Mechanical and thermal dynamics of the real and the simulated objects should be similar if their numbers (2) are equal because the inertia/capillary and inertia/freezing processes govern dynamics. Parameters of the runs simulated by the molecular dynamics combined with Monte Carlo (MD-MC) program are listed in Table 1 and in Fig. 2. Numbers near the circles in Fig. 2 correspond to the numbers (#) of the runs in the first column of Table 1. We see that there is approximate correspondence between the experimental and simulation regions in Fig. 2 if we exclude the runs with the highest velocities v c (0).
As was said, the experimental and simulated objects are approximately equivalent. At the same time, the number of atoms N atom in simulation is hundred times less. The simulations shed light onto the internal processes running inside the illumined hot-spot invisible experimentally.

Solidification of bump, jet formation, separation of droplets, and destruction of bump
The idea of miniaturization of the MD-MC runs using v 0σ , v 0χ parameters (2) has been presented above. Experimental (the curves 0.5, 2) and simulation (the circles) situations on the map of these parameters are shown in Fig. 2. The curves 0.5, 2 have been obtained from the estimates (8) and (11). In experiments, the bumps appear at F c ∼ 2F m . The jets are formed at F c ∼ 2.5F m . Somewhere above 3F m , ejection and destruction begin. The thick dashed curve in Fig. 2 marks the approximate position of the line separating non-ejecting and ejecting regimes in experiments. At the present, in experimental papers, the values of energy of pulse are given. It is difficult to obtain precise values for the absorbed central fluence F c from these energies; F c is necessary for (8), (11). Using simulation results, let us consider how dynamics and the final shapes change with variation of position of the point on the map of the parameters in Fig. 2.
Blue symbols in Fig. 2 correspond to the non-ejecting cases. They are the low-velocity cases, when the rate of freezing is enough to crystallize the bump as whole before its destruction. Hydrodynamic velocity and hence influence of inertia increase as we go farther from the origin in Fig. 2. The jet begins to form in the run # 30 shown in Fig. 3. But kinetic energy is relatively low, and the capillary and freezing effects are stronger; thus, only the embryo of the jet had time to develop. The colors show a distribution of the symmetry parameter: the more green colors correspond to the more ordered structures, green is a crystalline state while red is liquid gold. The run # 23 marked in Fig. 2 behaves similar to the run # 30 in Fig. 3; again the underdeveloped jet is formed in the tip of the bump.
The longer (relative to the # 30 in Fig. 3) jet develops in the case # 20 presented in Fig. 4. This is the case with significantly larger initial kinetic energy (see Fig. 2) in comparison with the run in Fig. 3. The case # 20 is near the separation line bounding the non-ejecting cases: the long jet and the droplet ready for separation are grown enough. The freezing only slightly advances the detachment of the droplet. The approximate position of the separation or bounding line for the MD-MC simulations is given in Fig. 2 by the continuous grey curve separating the blue and reddish symbols. Snapshot of material states, including red-colored melt and green-colored solid, at time t = Fig. 4. Approximately 40 % of the droplet is frozen to this moment. The simulation # 20 is continued up to the final frame (not presented here) corresponding to the time t = 2.7 ns, when the droplet is solidified entirely.

ns is shown in
The mechanics of formation of the droplet at the end of a liquid jet or a thickened rim around the hole in a liquid membrane has been considered in [22] in the case of a bump and in [28] in the case of membranes in a foam; see also Figs. 8 and 9 below. If we cross the jet perpendicularly by two slightly separated planes, then there are two closed contours belonging to the crossed liquid surface of a jet. There are capillary forces acting on to the cutting of a jet between these two planes. These forces act on the first and the second contours and compensate each other if the lengths of the contours are equal; if the lengths are not equal, then this effect (together with a force acting perpendicular to a jet) favors the compression of the more narrow part of a jet during the Rayleigh instability. Let us now cut the neck of the droplet by the same plane. We see that the compensating forces of surface tension are absent. Therefore, the tension existing in the neck decelerates the droplet relative to the neighbor part of the jet. The droplet rakes in mass from the jet. This is the reason why a droplet or a rim is formed. In [22], calculation of corresponding velocities are given. It is surprising that these velocities are moderate in the sense that the shell of a bump has time to expand even if there is a hole in a shell. It is surprising because it seems that the whole shell or membrane is in strongly stressed state and thus the area of the membrane should contract immediately after the appearance of a hole. Indeed, it contracts but not immediately.
The run # 31 presented in Fig. 5 is slightly above the grey line in Fig. 2 isolating the non-ejecting runs (the blue symbols). Much longer jet develops in the run # 31 in Fig. 5 relative to the run # 20 in Fig. 4. As we wait for a long time, the large droplet is gradually formed at the end of the jet while the jet becomes thinner and thinner-compare the frames (b) t = 1.152 ns and (c) t = 4.536 ns. At the same time, the crystallization front moving along the jet is still far from the droplet, compare with Fig. 4. It seems that all  Table 1 where the sizes are given. Green is solid gold, while red is liquid. Liquid is strongly overcooled below the melting point. Thus, nanocrystallites appear not only at the boundary of continuous solid with the liquid part but also inside the liquid part. The cross-section of the bump is presented. Thickness of this cross-section is 3 nm in direction normal to the plane of picture. The symmetry parameter shown by the green/red colors (defining the solid/liquid states) is obtained by averaging along this 3 nm. a t = 216 ps, b 0.72 ns, c 1.37 ns, d 2.34 ns  Table 1. It is located near to the grey continuous line separating the simulated non-ejecting and ejecting cases in Fig. 2. The instant t = 2.3 ns is presented. Solid parts are painted by green color and liquid parts by red color. The figure is obtained by summing and averaging the symmetry parameter along a normal direction perpendicular to the plane of the figure. Therefore, the thin horizontal line (slightly above the bottom) appears; it marks the pedestal formed by the gold film remaining in contact on a surface of a substrate. The same pedestal is clearly seen in Fig. 7 below. Also, the thin curves near the upper boundary of the shell of the bump show the internal boundary of the shell from the side of the internal empty cavity. Let us mention that the solidified shell undergoes few damping oscillations up and down before it transits into a rest state is ready for separation of the first large and fast droplet. May be smaller and slower droplets will accompany the large one. If we suppose that the hardly probable case of full freezing of the jet in Fig. 5c will nevertheless takes place, then this will contradict to the experimental observations [1,7,9,17,23,24] where never such large ratios of a length of a jet to a bump separation radius r sep have been obtained.
Let us also compare the positions of the points corresponding to the runs # 30 and # 31 on the map in Fig. 2 and propagation of the green solidified area in Fig. 3c, d and in Fig. 5a. In the low-velocity case # 30 and in the higher velocity case # 31, the shapes of the bumps are similar (this is work of inertia and surface tension) while the solidified parts are very different. This is the relative delay of solidification in the faster case # 31.
Additional increase of velocity v c (0) above the value corresponding to the run # 31 in Fig. 5 causes faster decay of a jet into the train of successive droplets. This is the example # 25 shown in Fig. 6. They are spread in velocities, temperatures, and sizes of the droplets. The first droplets move faster and are hotter. The last droplet, which recently separates in the instant shown in Fig. 6, is semisolid. The liquid part of this droplet is overcooled below the melting temperature. Thus, this droplet will crystalize soon. The first droplets can fly in a liquid state for a long time because the rate of radiative cooling is extremely low. The angular distribution of the droplets in their cloud is very narrow if they are produced in the process of decay of a long jet. At even higher kinetic energies, the part of a shell can decay into cloud of droplets, then the angular distribution of droplets is wider.
The runs ## 22, 28, 29a, and 29 have the highest velocity values v c (0) investigated in our runs listed in Table 1 and in Fig. 2. Corresponding capillary numbers v 0σ are the same for these runs because a coefficient of surface tension is the same for all runs. Let us discuss the thermal numbers since their values are important for results. To cool a bump in the MD-MC runs, we use a Langevin thermostat. First, a molten film has been prepared by increase of temperature above the melting point. The temperature profile supported by the thermostat has the maximum temperature 2000 K in the center. Temperature decreases down to 1500 K in the periphery of the simulation box. We switch off the thermostat inside the circle r < R sep after creation of an equilibrium liquid. While in the periphery, we steadily use the thermostat (i) to keep a film at the periphery in contact with the substrate (therefore, the pedestals in Figs. 4 and 7 are in rest at the substrate) and (ii) to support temperature at the periphery below the melting temperature. This thermal condition corresponds to real situation where outside the molten hot-spot a film is solid at room temperature. The supported temperature is 500 K for all runs except the runs # 29a and # 29 where it is 300 K. The run # 22 has small thermal velocity v 0χ , see Table 1, relative to the cases ## 28, 29a, and 29. The runs # 29a and # 29 are below the run # 28 because the supported temperature is lower (300 K versus 500 K). The characteristic time of thermostat is 0.8 ps for all runs except the run # 29a where it is 1.2 ps. Therefore, cooling in the case # 29a is slower than in the case # 29. The circles in Fig. 2 are plotted according to their capillary and thermal numbers (2). Only the circles # 29a and # 29 are shifted down relative to the point # 28 in accordance with the estimates following from the differences in the temperature at a periphery and the characteristic time of thermostat.
Thus, the run # 29 has the largest kinetic energy and the fastest cooling. It is singled out by the blue square inside the red circle in Fig. 2. Figure 7 shows instant 792 ps of the run # 29. The run is significantly above the non-ejecting limit (the grey continuous line in Fig. 2). An extremely long jet is produced. Only a small part of the jet is presented in Fig. 7. The jet decays into multiple small droplets not visible in the truncated Fig. 7. The crystallization front moves fast enough to solidify the whole shell and the bottom part of the jet when there is still a significant stretching velocity gradient du/ds along a solid part; s here is length of an arc, and u is the tangential component of velocity. The inertial forces connected with this gradient breaks off the shell. Strong thinning of a shell does  Table 1 and in Fig. 2 is shown. The instants are (a) t = 0.72 ns, (b) 1.15 ns, (c) 4.536 ns. Solid is green; liquid is red. Significant redistribution of mass from the shell of a bump to the jet takes place. We see that the shell in its middle part and the part near the jet is much thinner than in the part near the substrate where the shell keeps its initial thickness. We also see the flattened vicinity of the tip of the shell of the bump, compare with previous Fig. 4. Such flattened tips are often observed in experiments, e.g., see Figure 2c in [17] and Figure 4c in [23] easier the break off. The thinning of a shell results from transfer of large amount of mass into the jet; see Fig. 5 above. By this means a large slowly moving piece of a solid shell may appear in the ejecta in addition to the droplet component.

Self-developed inhomogeneities and break through of shell
Initial data for a MD-MC simulation are smooth. A run starts from an ideal single fcc gold crystal with the exact plane boundaries. After that, the heating and melting by a thermostat take place and a smooth (analytical function) vertical (see Fig. 1) velocity profile is imposed. Therefore later, only the thermal fluctuations, possibly in combination with some kind of hydrodynamic instability, cause deviations from ideal smoothness-there are no any initial perturbations. Instability may be the modified Rayleigh-Taylor type one, but deceleration, necessary for this type, is caused by surface tension acting on the largest scale ∼ r sep . Thus, it is relatively weak to drive the much shorter than r sep wavelengths, because for them the surface tension plays stabilizing role overcoming destabilizing deceleration since wavelengths are much shorter than r sep .  But this conclusion is valid for the plane thin membrane when we compare the long and short wavelengths with the infinitesimally small amplitudes. It is unknown how the situation will change when the largest wavelength has a finite amplitude and/or the perturbation wavelength is less than r sep but is not much smaller. Theoretical study of a thermal excitation of a capillary wave spectrum and instabilities needs a separate work. Here, we present results exploring violation of smoothness in MD-MC simulations.
The top view onto the frozen shell and the flying up (relative to the plane of the bottom panel in Fig. 8) jet is shown in the bottom of the three panels in Fig. 8. The jet at this instant has a height 2.3L box and remains molten; L box is the length of the computational cross-section. To imagine the jet, see Figs. 5, 6, and 7 with the similar jet structures. We clearly see the net of domains (nanocrystallites) grown from the seed crystallization grains; see also Figs. 3, 4, and 7. The liquid is overcooled, thus the isolated grains appear in the transition solid-liquid zone; see Fig. 3. Figure 3 presents distribution of an averaged (along line of view) symmetry parameter in a thin vertical crossing while in Figs. 4 and 7 the averaging is made along a line of view intersecting thicker pieces of matter. They are especially thick in the lateral regions where the line of view passes matter tangentially to the shell. Thus in these regions, the net of the crystallites is smeared due to averaging, compare with Figs. 3 and 8 (bottom) where the smearing is absent. Figures 8 and 9 are designed to emphasize the pattern of compactions/depressions. The pattern is formed from capillary oscillations [29,30] during the stage when the shell was liquid. Rapid crystallization fixes the pattern. The reasons for appearance of the compactions in the map in Fig. 8 in the internal rectangle and in the rectangle A,C in Figs. 8 (bottom) and 9 are qualitatively different. In the internal rectangle, we see the random compactions; they are the remnants of the thermally induced capillary waviness [29,30] (may be in combination with a hydroinstability). While in the rectangle A,C the main compaction is regular (the rim), it is caused by raking of mass by the expanding hole. The hole expands when a shell is liquid. Freezing fixes the rim.
The deep blue atoms in Fig. 9 are the bulk atoms distantly located from the two surface monolayers. We see them in the cut made by the plane A-B. These atoms form the thickening called the rim. The plane A-B intersects the rim at an angle. The atoms at the edge surface of the rim are yellow-red colorized. We see also the spot of depression located at the left side in Fig. 9. Inside the depression, the visual lines pass through the empty gaps between atoms. There are only two monolayers remaining in the depression; therefore, the gaps are seen.
We pay attention to the chaotic net of the alternating compactions and depressions because namely progression of the depressions leads to perforation of a shell. The spots of depressions and the stripes of compaction are seen in Fig. 8 (bottom). The depressions correspond to the greyish spots. They look like perforated. The inspection of the random pattern is presented in the top and middle panels in Fig. 8. The inspection is made by the Linuxbased software AtomEye capable of visualizing atomistic configurations. This software allows us to see better the inhomogeneous structure.
In Fig. 8, the individual atoms are colorized in correspondence with potential energy. In the more blue regions, the shell is thicker (this is compaction or condensation); there are the bulk; most blue atoms are there. There are no gaps in the condensations. While we see the empty small gaps between the atoms in the depressions in The square is the cross-section of the computational box in a plane of a substrate. A distribution of a symmetry parameter defining an interatomic order is given. The parameter changes in the range from green (fcc lattice) color to red (molten gold, the central jet). The yellow rectangle marks the area of the shell for which the AtomEye (atomistic configuration viewer) middle and upper views in this figure are shown. An atomistic representation of the second rectangular (A,C) is given in Fig. 9. Middle. The top view onto the atomic ordering. In the middle and the upper panels, the individual atoms are colorized according to their potential energy-deep blue corresponds to the strongest binding, less blue colors mark weaker binding-these atoms are higher at the axis of potential energy. Top. The lateral AtomEye view onto the same piece of the shell. The lateral cross-section shown in the foreground is the upper horizontal side of the middle rectangle, compare correspondence of colors Fig. 8 (middle), where only two surface monoatomic layers remain. Figure 8 (top) clearly demonstrates that indeed there are only two monolayers in the thinnest spots of a shell. There are pairs of atoms in these spots (one above another) having light blue color indicating that potential energy is small. On the other hand, there are deep blue colors in the thicker spots where 3-4 monolayers are present. Fig. 9 Cumulation of mass in the rim around the hole according to the AtomEye view; again, atoms are colored by their potential energy. The cumulation is a result of raking of material particles forming a shell in a process of expansion of a hole, compare with cumulation of mass into the head droplet in Fig. 4. The particles having their instant positions outside relative to the rim do not "know" about the hole and about its expansion at this particular instant [22]. The rectangular box A-B-C has the basement A-C corresponding to the rectangle A,C in Fig. 8 (bottom). The thickening corresponds to the rim located around the boundary of the hole crossed by the rectangle A,C in Fig. 8

Conclusions
Production of solitary surface nanostructures by the ultrashort laser pulses is used in many applications listed in the "Background" section. Productive process is mainly governed by the interplay between inertia, capillarity, and solidification. Therefore, it can be described on the twodimensional plane of the capillary and thermal numbers. In the paper, the regimes corresponding to the experimental set (parameterized in experiments by energy of an incident pulse within a fixed focal radius) are studied. The regimes change with energy increases. We have described dynamics in these regimes starting from a simple bump, to a bump with a nipple (## 23, 30, Fig. 3), after that to a bump with a frozen jet (near the limit of non-ejecting cases, Fig. 4), and after that the ejecting regimes begin (Figs. 5, 6, and 7).
Description is based on a combination of numerical codes and physical approaches. Physical picture of inflation, stopping, and crystallization of bump includes the thermal and mechanical transport of heat and momentum from the hot shell to the external cold part of a film remaining in contact with substrate outside the bump. The combination of codes includes two-temperature hydrodynamics (2T-HD) code together with MD-MC code (molecular dynamics combined with Monte Carlo method). We run series of the 2T-HD simulations to explore dependence of separation velocity of gold film and its temperature on the absorbed fluence. This part of work is described in [22,31]. New results presented above are based on the usage of the Monte Carlo code for a heat transport problem for implementing the fast freezing of the bumps with material molten by an ultrashort laser heating. Formation of chaotic solidified inhomogeneities is described in Figs. 8 and 9.
Similarity on the capillary v 0σ and thermal v 0χ numbers allows us to use scaling in MD-MC simulation of a large number of involved atoms N atom in experimental system by a smaller atomic system. Typical experimental sizes vary from the beam radii R L of the order of few micrometer to the submicrometer radii in the cases where the diffraction-limited tight focusing UV lasers or higher harmonics of a Ti:sapp laser are used. With the scaling, we correctly describe formations of a bump, short and long frozen jets, and production of droplets by fragmentation of the long jets. Thresholds between these regimes agree satisfactorily with estimates of the experimental thresholds; see Fig. 2. Here it is worth to note that, unfortunately, accurate estimates of absorbed fluence F c are not available in the experiments.
Again, the scaling equivalence is exact for the listed processes: formation of bumps, jets, and droplets from jets. Indeed, even a complicated process of the droplet production is governed by the surface tension and the exponentially fast Rayleigh instability. Thus, it is enough for a jet radius in simulations to be larger than few interatomic distances. This demand is hold in our simulations. In this connection, it is necessary to say that a phenomenon of restoration of the surface tension coefficient with a number of interatomic layers is carefully studied; see Figs. 8 and 9 and Appendix in [22]. We are not sure only in accuracy of our description of the process of formation of a hole in a stretched film, Figs. 7, 8, and 9. It seems that our simulations here give the order-ofmagnitude estimates of timing and threshold. It is necessary to stretch a film down to 1-nm thickness to cause a rupture. In simulations, we use from few times to one order of magnitude initially thinner films. It is not clear if exponentially fast instability presented here determines the duration of rupture by its inverse increment (like the Rayleigh instability of cylindrical jets). If not, then the duration depends on initial thickness and will be faster for thin films. There are many cases where very strong thinning of films and formation of holes are observed in experiments.
In conclusion, we can say that our model successively describes several phenomena. Among them are