Formation of solitary microstructure and ablation into transparent dielectric by a subnanosecond laser pulse
Laser ablation in liquid (LAL) is important technique used for formation of nanoparticles (NP). The LAL processes cover logarithmically wide range of spatiotemporal scales and is not fully understood. The NP produced by LAL are rather expensive, thus optimization of involved processes is valuable. As the first step to such optimizations more deep understanding is necessary. We employ physical models and computer simulations by thermodynamic, hydrodynamic, and molecular dynamics codes in this direction. Absorbing light metal expanding into transparent solid or liquid dielectrics is considered. We analyze an interplay between diffusion, hydrodynamic instability, and decrease of surface tension down to zero value caused by strong heating and compression transferring matter into state of overcritical fluids. The primary NPs appear during expansion and cooling of diffusion zone when pressure in this zone drops below critical pressure for a metal. Long evolution from the overcritical states to states below a critical point for a metal and down to critical point of liquid and deeply down to surrounding pressure of 1 bar is followed. Conductive heating of liquid from hot metal is significant.
keywords: laser ablation, instability, condensation
Paper presents theoretical and simulation results on laser ablation in liquid (LAL) or into initially solid transparent dielectrics; gold-water and gold-silica pairs are investigated. We consider and compare ultrashort fs and subnanosecond ps pulses and vary absorbed fluence in a wide range up to values limited by optical breakdown of dielectrics; threshold values for breakdown are W/cm for the subnanosecond case with eV [1, 2]. LAL is used to prepare colloidal solutions of nanoparticles (NP). These processes have many technological applications, see detailed reviews [3, 4]. Feeling about the today’s level of achievements and unsolved problems may be obtained reading approximately hundred one-page abstracts from the ANGEL 2018 conference .
In spite of numerous applications and many peoples involved in working in this direction, a lot of problems remain poorly understood. This is because the ultrafast early stages still are unobservable experimentally. Only rather late stages starting from the microsecond s) time scale are carefully documented . These late stages relates to formation and development of a cavitation bubble. While pulse durations from 0.1 ps to a few ns are from 7 to 3 orders of magnitude shorter.
Recently computer simulations shed light onto the beginning of the ”hidden era” lasting up to s before the bubble begins to form [7, 8, 9, 10, 11, 12]. Ultrashort pulses with duration ps were considered, and the processes have been followed up to few nanoseconds in [7, 8, 9, 10] and up to 200 ns in [11, 12]. In all these papers the absorbed energies were comparable and of the order of J/cm In the present paper we analyze new case with the three order of magnitudes longer pulse.
At first glance, it seems surprising that dynamics at the stage beginning after instant ps (when the 50 ps pulse is ended) is similar in the cases with so different durations (0.1 ps versus 50 ps) if energies are the same (in more details this question is discussed below). This is what we find comparing the cases with ps and ps, see Fig. 1. Thermal diffusivity of gold is large at the 2T stage  relative to its 1T values in solid and liquid states. But duration of the ultrashort pulse and duration of a 2T stage are shorter than duration 50 ps of the subnanosecond pulse. Therefore thicknesses of the heat affected zones are approximately the same in both cases, see Fig. 1.
Thicknesses of the liquid layers are approximately the same in Fig. 1. The thin layer separating liquid and solid has approximately the position nm in both cases at the instant shown. Our model of heat conduction includes the change of thermal conductivity caused by melting. Therefore the gradient is 2-3 times smaller in the solid. For the same reason the gradient has a sharp kink at the Au-water contact boundary (CB), see Fig. 1.
As was said, the thermal distributions are similar in the cases with the durations 0.1 ps and 50 ps shown in Fig. 1. But there are differences also. Before saying about them let’s describe a global structure of the 1D flow. In the both cases the foam region and atmosphere are formed [7, 8, 9, 10, 11, 12]. The atmosphere is a layer contacting with water and being in quasi-hydrostatic equilibrium with CB as it is described in [11, 12]. The atmosphere separates foam and water. In the case of a rather long subnanosecond pulse the atmosphere forms when CB begins to decelerate. As it will be seen below, the CB accelerates during the rise of the pulse ps and some time after; we begin our simulations at the negative time instant equal to
In the ultrashort case shown in Fig. 1 the nucleation process and formation of foam begins earlier than in the subnanosecond case. This is the difference between the cases in Fig. 1. In the ultrashort case nucleation begins at the stage ps while in the subnanosecond case the starting stage of nucleation is ps.
Nevertheless the amount of the ejected gold mass g/cm is approximately the same in the both cases. The ejected mass forms a layer separated by a gap from the rest of a bulk gold target up to the stages ns. The layer contacts with water. Velocity of the layer drops down to very small values m/s at the stage ns [11, 12]. Fate of the layer is unknown at the late stages s. The gap is filled by saturation vapor of gold and has pressures bars until gold has temperatures above kK [11, 12]. While pressure in liquid and later in gaseous water is bars at the stage of few microseconds. Thus the gold layer may moves back under water pressure and closes the gap if the layer is cooled down below few kK. The layer may stick back with the bulk gold if gold in the layer remains liquid during the shrinkage of the gap. Velocity change of the layer under pressure difference acting onto the layer from its internal (vapor of Au in the gap) and external (CB, water) boundaries is large km/s for bars, g/cm and s. Such velocities easily may close the gap during s because the gap thickness is m only. This problem obviously needs separate consideration.
Below we will limit ourself with comparative analysis of the early 10s ps) and middle ns) stages of the ultrashort Au-water and subnanosecond Au-water and Au-glass cases. There are two main conclusions following from this analysis.
First. All these three cases are similar in the sense of (i) energy distribution in the surface layer, (ii) the value of ejected mass, and (iii) the trajectory of the CB where the stands for kinematic shift, velocity, acceleration) or for thermodynamics quantities.
Second. The states near the CB during the middle stage ns are hot and compressed, thus both the gold and, of course, the water (with its low critical parameters) are in overcritical states. Therefore Au-water surface tension coefficient is zero and diffusion is strong. Rather thick diffusively mixed Au-water layer forms around the CB thanks to that. Later in time the pressure drops due to expansion of the pressurized layer near the CB. Condensation into clusters and liquid droplets begins as pressure decreases below the critical values for gold.
From the other hand, strong diffusion suppresses development of the Rayleigh-Taylor instability (RTI), see Section 7 below. There is a limited in time interval when the RTI develops [11, 12], see Section 3 and Fig. 11. Thus if the development is delayed during this temporal interval then contribution of RTI remains insignificant. Significant contribution from RTI into the nanoparticles (NP) production was observed for ultrashort actions and smaller (than in this paper) absorbed energies in [9, 10, 11, 12].
The paper is organized as follows. First in Section 2 we consider separation of large distances covered by shocks in gold and water from small distance passed by contact boundary (CB). Velocity of a CB decays in time quickly while shocks propagate at velocities above speed of sound. In Section 3 the acceleration stage is considered. This stage differs the case of ultrashort pulse from the case of subnanosecond pulse. Later in time (that is after the acceleration stage) these cases behave similarly. In Section 4 duration of the stage when parameters of a metal at the CB are in overcritical states (overcritical stage) is analyzed. It is found that this duration is a function of absorbed energy.
We use two-temperature (2T) and one-temperature (1T) hydrodynamic codes to obtain results listed in the Sections 2-4. Embedded atom method (EAM) potentials of interatomic interactions and molecular dynamics (MD) simulations are applied to reveal another part of valuable information in the next Sections. Saturation pressure, curves of equilibrium between phase states, surface tension coefficient, and description of diffusion on the base of the EAM/MD approach is presented in Section 5. Condensation of gold and formation of nanoparticles is described in Section 6. Section 7 demonstrates results of large scale MD simulation of ablation of gold in water.
2 Global and internal structures
Global structure of flow is clear from Fig. 2. Before action of a pulse we have two motionless half-spaces in contact. The left one is filled with water, while the right one corresponds to the bulk gold target; here is Lagrangian coordinate; it corresponds to initial (before laser action) positions of material particles. Laser beam penetrates through transparent water and heats absorbing gold in its skin-layer. Heating and pressure rise cause motion. After finishing of a pulse the spatial expansion gradually decreases pressure. Two shock waves (SW) in Fig. 2 run into water and into gold target. The pulse is rather long ps. Therefore it takes some time to form the shocks. The SWs form thanks to nonlinear effects after focusing of characteristics and wave breaking.
Velocity ”jump” 0.23 km/s in Fig. 2 covers the layer where the atmosphere and foam locate. The ”jump” characterizes expansion of the CB relative to the bottom of the foam at the instant 0.5 ns. This velocity is small (just few %) in comparison with velocity km/s giving the lowest limit of velocity of increase of a distance between the SWs in Fig. 2. Expansion proceeds as a result of easy stretchability of a vapor component in foam. Structure of the thin layer near the CB (the ”jump” in Fig. 2) is revealed in Fig. 3.
Comparison of the Au-water simulations for ultrashort pulse fs [11, 12] and subnanosecond pulse with ps is demonstrated in Fig. 4. Absorbed energies are the same. We cut out the slowly moving part of bulk gold in the ultrashort case [11, 12]. Only really significant for mixing with water and nanoparticles production part of a gold target was kept in the simulation of ultrashort action. Comparing two simulations in Fig. 4 we conclude that:
First, at the nanosecond time scale and later the duration varied 500 times from 0.1 to 50 ps is insignificant.
Second, the passive part of the gold target kept in one simulation (pay attention to the SW running into bulk) and removed in the another one is, indeed, insignificant for dynamics of foam and water.
But internal structures of density distributions inside the foamy layer differ, see Fig. 5. This is because of random origin of nucleations and separation into vapor and liquid phases and (may be) due to different acceleration/deceleration behavior at the stage of the order of a few tens of ps.
Expansions in water (liquid dielectrics) and in glass (initially solid dielectrics) qualitatively are similar up to the stage when in the case of water the bubble begins to form. Corresponding situations are shown in Fig. 6. In Figures 6 and 7 the cases with a bulk gold target and a thick film target (initially 370 nm thick) are compared. We call a film thick if its thickness is above thickness of a heat affected zone (HAZ) [13, 14]. Thickness 370 nm is larger than thickness of HAZ, see Fig. 1. The case with a film 370 nm and a glass substrate is considered in connection with experiments [15, 16]. Discussion in connection with these experiments needs separate description. Figures 6 and 7 demonstrate dramatic difference between the cases with a bulk and a film targets if we consider motion of gold; dynamics of dielectrics are similar as it follows from Fig. 6. Initially the system from dielectrics and a target was motionless, thus its momentum is zero. Laser beam introduces energy into the system, but not momentum because momentum of a photon cloud is negligibly small. Therefore total momentum of a system is zero after laser impact.
The momentum directed into the side of dielectrics produces motions in dielectrics which are more or less similar in both cases, see Fig. 6. But the motions in the gold targets triggered by the right side momentum are very different depending on thickness of a film. In a bulk target this momentum is carried into the bulk side and out from the near CB layer by approximately triangular shock wave SW in Figures 2, 4, and 6. In the case of a film this carryover is impossible - thickness is limited. Thus all amount of the powerful push to the right remains in a film. This causes decay of a film into a sequence of quasi-spallative-plates separated by foamy layers in case of a film, see Fig. 7. Velocity of the rear-side plate is large, see Fig. 6.
3 Active dynamics: acceleration, deceleration, and comparisons
Dynamics of expansion of CB in the case of Au-wt pair is shown in Fig. 8. The cases with ultrashort and subnanosecond pulses are compared. Absorbed fluences are the same 400 mJ/cm In both cases intensity is given by the Gaussian law But durations differ 500 times. Thus if we impose these Gaussian laws onto the plane in Fig. 8 then the subnanosecond pulse will have visible left and right wings. And the rise of its left wing is proportional to the growing shift of CB in the subnanosecond case. While the Gaussian for the ultrashort pulse looks like delta-function at the temporal scale used in Fig. 8.
The speed up interval is very short for the pulse with ps. We don’t see it because the time scale in Fig. 8 is too large. At this scale, it seems that in the ultrashort case the CB starts to move immediately at the instant with large velocity. Later in time it only decelerates thanks to resistance of water to expansion of gold. But in the subnanosecond case the definite acceleration stage exists. Of course, it is related to the increase in time of absorbed intensity at the left wing of Gaussian pump.
Influence of absorbed energy, resistive ability of dielectrics, and thickness of a gold target on the trajectory of the contact boundary (CB) is illustrated in Fig. 9. Obviously, the shift increases if we increase absorbed energy fixing two other parameters: density of dielectrics and thickness of a gold layer. If we increase density of surrounding medium then deceleration is faster; initial density of the pyrex silica used in our Au-glass simulations is 2.2 g/cm Compare the cases wt/400 mJ/cm and gl/560 mJ/cm in Fig. 9. This is also obvious.
Less obvious are consequences of variation of thickness of a layer of Au. It seems, at first sight, that a thick target will decelerate more slower. But there is foam formation and appearance of a contact layer of gold (atmosphere) [11, 12]. Speed of sound in foam is low. The first portions (atmosphere) of fragmented Au are decelerated by surrounding water or glass, while the second portions move faster thus adding their mass to the first portions located at the CB. This is said to underline that (i) there is a decelerated layer (atmosphere) and (ii) its mass may change with time.
Analysis of simulations show that asymptotically velocity and pressure at the CB are acoustically connected where is acoustic impedance in surrounding medium. Writing we find that
where is CB position, and are density and speed of sound of surrounding matter. The estimate (1) gives for the deceleration time constant a value ns for density ratio height of atmosphere nm, and speed of sound in water km/s.
The difference between ultrashort and subnanosecond laser drives is distinctly clear from Fig. 10. Whereas there is a relatively prolonged acceleration interval in the subnanosecond case, in the ultrashort case the velocity trajectory starts from high value and after that only decelerates with time. Let’s mention that acceleration continues after the instant when the maximum of a laser pump is achieved. Thus intensity begins to decrease while velocity still increases. Increase of velocity continues approximately ps after maximum of a pump in
Velocity drop in the ultrashort case shown in Fig. 10 is well described by a function km/s where time should be taken in ps. The subnanosecond velocity trajectories obey the approximately same law but only after few tens of ps after their maximum of velocity Deceleration dependence on time is plotted in Fig. 11. There is relatively slow acceleration in the subnanosecond case. The acceleration interval is lost for the Rayleigh-Taylor (RT) amplification of perturbations of CB. Thus development of the RT instability is significantly weaker in the subnanosecond case, compare the curves ”gain USP” and ”gain subNS” in Fig. 11. The temporal interval of a few tens of ps is important for development of the RT disturbances because during this interval the deceleration is huge, see the curve ”decel USP” in Fig. 11. These decelerations overcome values cm/s typical for surface of a neutron star. Atomic decelerations in interatomic interactions are of the order of cm/s The curves ”gain USP” and ”gain subNS” in Fig. 11 corresponds to the wavelength 60 nm of a harmonic perturbation. They were calculated using linear RT theory including surface tension and viscosity according to the method described in [11, 12].
4 Strong heating, inertial geometrical confinement, overcritical states
Gold is strongly heated at the end stage of absorption of a subnanosecond pulse, see Figures 1 and 12. Simulations are done using equation of state for gold, water and glass taken from [17, 18, 19, 20, 21]. We use heat conduction description for gold from [12, 22]. The description includes (a) calculated contribution from electron-electron collisions important at high temperatures and (b) decrease of conductivity during melting. Critical parameters of gold are 7.8 kK, 0.53 GPa, and 5.3 g/cm [17, 18, 19, 20, 21].
Heated gold is cooled due to heat conduction into deeper layers and thanks to adiabatic cooling during expansion. Inertial confinement by rather dense dielectric limits a rate of expansion and thus the rate of cooling. It takes approximately 0.2 ns to cool the CB below critical temperature in the case with mJ/cm In the case with mJ/cm the time interval for such cooling is longer, see Figures 13-15.
The temporal interval, when gold near a CB is in overcritical state, plays a special role in the process of nanoparticles (NP) formation. See discussion in the next Chapter. Here we have to discuss strong increase of temperature near the CB (see Figs. 12, 13) in the case of the most powerful action mJ/cm from the sequence of the considered cases. It is obvious that temperature at the CB after absorption of a subnanosecond pulse increases as energy increases while density drops down. How high is temperature at the CB? From energy balance we have an estimate where is heat capacity of gold per unit of volume. Our model of heat conduction [12, 22] describes dependence of thermal diffusivity on density. Values of and decreases as density decreases. Thus temperature isn’t simply proportional to as it is in the case with permanent values of and
Namely nonlinear dependence of a function on causes appearance of the overheated contact layer in Figures 12 and 13, see also Fig. 15. In this layer temperature increases 5-7 times above the average level. While density drops 5-7 times below the average level, see Fig. 15. Thus pressure is approximately the same across the matter of glass near the CB, in the overheated layer, and in the dense gold near this layer, see Fig. 14.
Transport properties of gold in the states with g/cm and few eV are poorly known. Thus we cannot pretend that our conduction model [12, 22] is very accurate at these hot states of intermediate density between condensed and gaseous values. But we are sure that at least at the semi-quantitative level the model is true in these states. Therefore phenomenon with the nonlinear overheating have to take place but may be at slightly different energies If we suppose absorption of gold at the level than incident intensity is W/cm that is somewhere not far from the values corresponding to optical breakdown of dielectrics. This understanding seems limits achievable level of the highest temperatures
5 Disappearance of capillary barrier and strong enhancement of diffusion
As follows from Figures 1, 12-15 the CB is above the critical point for gold for some temporal interval. Critical point of water (220 bar, 647 K, 0.3068 g/cm is much lower than the critical point of gold. Thus both liquids are in the overcritical states. Then they behave like the non-ideal gases. Surface tension at the contact boundary (CB) between them becomes zero. While diffusion coefficient increases. In condensed matter the diffusion process is slow because an atom loses time doing many oscillations inside its potential well before it tranships to the neighbor well.
We use molecular dynamics (MD) simulation with EAM (embedded atom method) interatomic potential for gold to calculate saturation pressure (Fig. 16), the boiling and condensation curves (Fig. 17), and surface tension at the boiling curve (Fig. 18). We modify our previous EAM potential for Au-Au interaction described in . The modified potential provides the correct dependence of surface tension coefficient on temperature shown in Fig. 18 and dimer Au-Au binding energy. Mechanically both potentials are close to each other (equilibrium density, compressibility, Hugoniot adiabat). Melting temperatures are very near to the experimental value for both potentials. We refer to these EAM potentials as old and new EAM potentials in Figures 16 and 17 and below in the text.
Figure 16 presents the binodal curve at the phase plane. We see that both EAM potentials underestimate saturation pressure at high temperatures near the critical point relative to the data from [17, 18, 19, 20, 21]. While critical temperature is or higher for the new EAM potential, or lower for the old one. It should be said that the exact value of critical pressure and temperature of gold still remains unknown. Thus in reality we cannot definitely conclude what set of data is more close to the real gold.
Achievement of the thermodynamic equilibrium between liquid and its vapor in the MD run becomes more difficult as temperature and thus pressure decrease (because the rate of evaporation drops down sharply thus simulation time increases). But the equilibrium is achieved in all cases presented. The lowest pressures are of the order of 1 bar, see Fig. 16. The difference between gas and liquid pressures in Fig. 16 for the same EAM potential follows from large errors in definition of pressure in liquid when pressure is low. Indeed, in liquid pressure is a sum of two large values having opposite signs: thermal and virial pressures. Unfortunately these errors leave rather large systematic shift. Therefore pressures in liquid in Fig. 16 are below the gaseous pressures. We have to use the gaseous pressures as saturation pressure for a given EAM potential. It is interesting that the pair interatomic potentials are less sensitive to this problem than the EAM (embedded atom method) potentials (they include many-body interactions). Liquid pressures are presented to familiarize readers with quality of results and real problems existing in MD approach.
Fig. 17 shows boiling and condensation curves at the density-temperature plane. Equation-od-state (EoS) [17, 18, 19, 20, 21], and two EAM potentials are compared. We see that at the plane the new and old potentials limit the EoS data from above and from below.
Molecular dynamics (MD) simulation of atom diffusion in the hot binary mixture Au-water was performed with using three different potentials describing interactions between atoms Au-Au, water-water, and Au-water. The new EAM potential was used for describing of Au-Au interactions. Water is considered in our simulation as a monatomic system, for which the EAM potential was developed in , see the data files on the project https://www.researchgate.net/project/Development-of-interatomic-EAM-potentials. It provides the correct mechanical properties of liquid water including the density and sound speed in ambient conditions, and shock Hugoniot as well.
Due to lack of information, the interaction of Au atom and water molecule is greatly simplified to a pairwise potential between Au and “water” atoms using the van der Waals radii 0.166 nm for Au and 0.152 nm for oxygen atom, and a strong repulsion is assumed. The usage of larger van der Waals radius 0.19 nm of water for our “water” atom reduces the diffusion coefficient of Au atoms in the simulated 50/50 binary mixture by about 10-20%; 50/50 is the ratio of the numbers of atoms per unit of volume; partial densities are 4.1 and 0.4 g/cm for Au and water, respectively. Results of MD simulations relating to calculations of diffusion coefficient in the overcritical states of gold and gold-water mixture are presented in Table 1. Corresponding MD simulations were done using new EAM potential. Positions of corresponding points (where coefficient was calculated) are marked as magenta rhombus in Fig. 17.
|of Au [g/cm]||4.1||8.2||16.4|
6 Cooling, condensation, and formation of nanoparticles
There is significant gold-water atomic inter-diffusion mixing process when matter near the CB exist in the overcritical states; see Figures 12-15 about duration of residence in the overcritical states and Fig. 17 with Table 1 about the enhanced value of diffusion coefficient in these states. Thus the jump like CB diffuses into a smearing layer. Thickness of this layer is Taking durations from few to few tens of nanoseconds we obtain values of from few tens to nm.
Hydrodynamic expansion causes cooling of mixture thanks to work done for expansion. Gradually temperature drops below the critical value for gold. Then condensation into clusters of Au atoms begins. These clusters are in Brownian movement, sometimes meet each other and stick together. Atoms surrounding a cluster condensate to a cluster and evaporate from it, thus mass and energy exchange between atoms and a cluster takes place. Our MD simulations presented below show that temperatures of atoms and clusters are approximately the same. Thus collisional removal of the positive energy excess from growing in size (and hence heating) cluster is fast enough to support equality of temperatures.
The hydrodynamic expansion rate is defined by the gradient of velocity near the CB. At the nanosecond time scale this gradient is of the order of 1/s, see Figures 2, 4, and 6. Therefore expansion velocity of the CB diffuse layer is a few meters per second.
Following these conditions we run several MD simulations to estimate temporal range of nucleation and growth of nanoparticles (NP) as pure gold or the 50/50 mixture (in concentrations of Au and water) intersect a condensation curve shown in Fig. 17. Calculations are done using the old EAM potential. Condensation starts with some delay in time after the intersection. Of course, this delay depends on the rate of expansion - in the case of very slow expansion the condensation begins very close to the intersection point.
Figures 19 and 20 show how the adiabatic curves kink as they pass into the two-phase region. The kink corresponds to the transition from the fast temperature drop with expansion to the slow temperature drop. The fast drop (rather large effective adiabatic index relates to the one-pase behavior (Au is in a gas state). The slow drop (decrease of begins when liquid component appears in a two-phase vapor-liquid mixture. Release of evaporation energy through condensation during volume expansion reduces the rate of temperature decrease.
Fig. 19 presents global view of the EoS and old EAM boiling curves and the MD simulated adiabatic curve crossing the boiling curve. The MD simulation is intended for definition of the transition piece from the one-phase to the two-phase behavior along the adiabatic expansion process for our particular conditions: velocity gradient 1/s and the diffusion thickness of the order of few tens of nm. We are interested to know how long this transition lasts and how large nanoparticles are at the end of our temporal range. The MD simulation is rather difficult because duration is long (tens of ns) and system should consists from large number of atoms Therefore we limit our self to a narrow time interval just before the intersection (omitting trivial one-phase piece of the gaseous adiabatic curve) and we follow the process during a time scale typical for our problem. This explains why the piece of the adiabatic curve is short in Fig. 19 relative to the field shown.
Fig. 20 demonstrates shape of the transient kinks for pure Au and for a mixture Au-wt. This is the local vicinity of the points where the adiabatic curves for pure and mixed Au cross the curves of thermodynamic equilibrium for pure and mixed Au, compare with Fig. 19. MD simulations with old EAM interatomic potential were used to obtain these two adiabatic curves. The simulation boxes for the both cases have fixed in time periodic boundary conditions at the lateral walls, that is in the directions. Lateral dimensions are nm. The moving walls are imposed in the direction. There are reflecting potentials connected with these -walls. They return the Au or Au and water atoms colliding the -wall back into box. Reflection from the expanding walls extracts small (velocity of expansion is much less than thermal velocity) part of kinetic energy of colliding atoms thus adiabatically cooling a whole system in the box.
Initial densities and temperatures are given in Fig. 20 and Tables 2 and 3. Initial distance between the -walls is nm. This starting distance is the same in pure Au and the mixture cases. Expansion velocity is 2 m/s. It is also the same in the both case. Distance between the -walls grows with time and density decreases This is a model of spatial expansion of diffusion mixture in real conditions during expansion and cooling.
Figures 21 and 22 present pictures of evolutions of pure Au (Fig. 21) and mixed systems, respectively. Detailed description of simulation approach and obtained results is given in the paper  presented at the ANGEL-2018 conference. Movies of evolution are attached to this paper. Figures 21 and 22 show how nanoparticles nucleate and grow. They corresponds to the view along the -axis. The simulation box is 100 nm long in the -direction. Figures 21 and 22 are composed from the set of the individual pixels. Every pixel contains information about potential energy averaged along the -direction.
More light green pixel marks position where potential energy is deeper. The light green spots are the nanoparticles (NP). Their sizes are of the order of ten nanometers at the final stages in both cases of pure Au (Fig. 21) or Au-water mixture 22). Atoms sticking together into NP have profit in their interaction energy. Sticking of part of Au atoms in the two-phase system does deeper potential energy of the whole system (averaging along -directions). This is shown in Fig. 20 by the dependence relating to potential energy. Figures 21 and 22 show the same but with averaging only along the -direction. Temperatures in Figures 21 and 22 are higher than melting temperature of gold. Thus NP remains in liquid states in the temporal range considered.
7 Molecular dynamics simulation of ablation in water
From the sequence of three shots compared in Figures 9, 9 and 12-15 we see that increase of absorbed energy causes serious consequences. The main feature is that with this increase a metal near a CB transfers from undercritical to overcritical conditions and remains their for rather long time. Undercritical (condensed) and overcritical (uncondensed, gaseous) conditions differ qualitatively. In the undercritical case (relatively small the liquid gold contacts with liquid; our fluences are above melting threshold for Au. Thus to penetrate into liquid and to form NPs the Au atoms have to overcome surface barrier through evaporation/dissolution (the first step) and after that they diffusively drift out from a CB (the second step). At the last step the expansion and cooling lead to condensation and formation of NPs, see Figures 20 and 22.
In the overcritical states the first step disappears together with surface tension. This is very important. Indeed, low values of saturated pressures of evaporated/dissolved Au atoms outside the boundary of a continuous condensed phase (outside the CB) strongly limits amount of Au atoms in a diffusion layer and thus limits amount of NPs. In the overcritical states only the second and third steps are present. This enhances diffusion - now diffusion starts not from evaporated/dissolved atomic number density of Au near the CB but from much larger value of atomic number density directly in the continuous phase.
We perform a large scale MD simulation to follow ablation of gold into water at absorbed energy larger than in the case previously considered in paper . Comparison of these two cases together with detailed description of evolution of the flow is given in the paper  presented at the ANGEL 2018 conference. Here in Fig. 23 rather late instant from MD simulation with large is presented. We see that Au-water interpenetration significantly overcomes the estimate nm for cms and ns. May be enhancement of interpenetration is due to presence of some embryos of Rayleigh-Taylor instability. Indeed, it is known that diffusion is a laminar like process which produces very flat atomically mixed layer; diffusion smooths away lateral roughness. But in Fig. 23 we see significant inhomogeneities in the direction transverse to the direction of longitudinal expansion.
Above ablation of gold into transparent media is analyzed. Initial stages of formation of nanoparticles (NP) are considered. These stages relate to tens and hundreds of nanoseconds (see Section 6) well before bubble dynamics comes into play. All thermodynamical (with first order phase transitions) and hydrodynamical aspects have to be included to come to understanding of appearance of primary NPs. It is found that
* Deceleration of the contact boundary (CB) at the nanosecond time scale is approximately exponential with the e-folding time of the order of 0.5-1 ns. This law describes deceleration in the 1-0.1 km/s range of velocities after ultrashort and subnanosecond pulses.
* Maximum shift of the CB from its initial position is rather limited and is of the order of micron.
* Temperature of gold at the CB nonlinearly depends on absorbed energy Temperature strongly increases for energies J/cm up to values K for subnanosecond pulses. These values are much above critical temperature. While density of gold decreases down to values g/cm
* Ablation motion triggered by a subnanosecond pulse is order of magnitude more stable against Rayleigh-Taylor instability (RTI) in comparison with an ultrashort action, see Fig. 11. This is simply caused by the shortening of deceleration stage because more than a half of the heating pulse is spent to accelerate the CB.
* Another important reason to limit the RTI perturbations growth is connected with high temperatures, decrease of the density ratio at the CB, and enhancement of diffusion. Enhanced diffusion is a powerful stabilizer of RTI. Diffusive smearing grows as while for RTI in the case with permanent deceleration the law of growth is Thus diffusion always dominates at the early stages while the RTI wins this competition later in time. But in laser ablation the deceleration quickly decays in time. Thus there is only limited interval of times when perturbations can develop thanks to RTI. In this situation the diffusion process suppresses RTI. Then RTI is removed from the process of fabrication of NPs at the capillary scale. Appearance of NPs of the capillary scale was observed in simulations [9, 10, 11, 12].
-  B. C. Stuart, M. D. Feit, A. M. Rubenchik, B. W. Shore, M. D. Perry, Laser-induced damage in dielectrics with nanosecond to subpicosecond pulses, Phys. Rev. Lett. 74 (1995) 2248–2251. doi:10.1103/PhysRevLett.74.2248.
-  A. V. Smith, B. T. Do, Bulk and surface laser damage of silica by picosecond and nanosecond pulses at 1064 nm, Appl. Opt. 47 (26) (2008) 4812–4832. doi:10.1364/AO.47.004812.
-  D. Zhang, B. Gökce, S. Barcikowski, Laser synthesis and processing of colloids: Fundamentals and applications., Chem Rev. 117(5) (2017) 3990–4103. doi:10.1021/acs.chemrev.6b00468.
-  J. Xiao, P. Liu, C. Wang, G. Yang, External field-assisted laser ablation in liquid: An efficient strategy for nanocrystal synthesis and nanostructure assembly, Progr. Mater. Scienc 87 (2017) 140 – 220. doi:https://doi.org/10.1016/j.pmatsci.2017.02.004.
-  J. Lama, J. Lombard, C. Dujardin, G. Ledoux, S. Merabia, D. Amans, Dynamical study of bubble expansion following laser ablation in liquids, Appl. Phys. Lett. 108 (2016) 074104. doi:10.1063/1.4942389.
-  M. Povarnitsyn, T. Itina, P. Levashov, K. Khishchenko, Mechnisms of nanoparticle formation by ulra-short laser ablation metals in liquid environment, Phys. Chem. Chem. Phys. 15 (2013) 3108–3114. doi:10.1039/c2cp42650a.
-  M. Povarnitsyn, T. Itina, Hydrodynamic modeling of femtosecond laser ablation of metals in vacuum and in liquid, Appl. Phys. A 117 (1) (2014) 175–178. doi:10.1007/s00339-014-8319-1.
-  C.-Y. Shih, M. V. Shugaev, C. Wu, L. V. Zhigilei, Generation of subsurface voids, incubation effect, and formation of nanoparticles in short pulse laser interactions with bulk metal targets in liquid: Molecular dynamics study, J. Phys. Chem. C 121 (30) (2017) 16549–16567. doi:10.1021/acs.jpcc.7b02301.
-  C.-Y. Shih, R. Streubel, J. Heberle, A. Letzel, M. V. Shugaev, C. Wu, M. Schmidt, B. Gökce, S. Barcikowski, L. V. Zhigilei, Two mechanisms of nanoparticle generation in picosecond laser ablation in liquids: the origin of the bimodal size distribution, Nanoscale 10 (2018) 6900–6910. doi:10.1039/C7NR08614H.
N. Inogamov, V. Zhakhovsky, V. Khokhlov,
Laser ablation of gold into water:
near critical point phenomena and hydrodynamic instability, arXiv (2018)
-  N. A. Inogamov, V. V. Zhakhovskii, V. A. Khokhlov, Dynamics of gold ablation into water, J. Exper. Theor. Phys. 127 (1) (2018) 79–106. doi:10.1134/S1063776118070075.
-  N. A. Inogamov, V. V. Zhakhovskii, V. A. Khokhlov, Jet formation in spallation of metal film from substrate under action of femtosecond laser pulse, J. Exper. Theor. Phys. 120 (1) (2015) 1548. doi:10.1134/S1063776115010136.
-  X. W. Wang, A. A. Kuchmizhak, X. Li, S. Juodkazis, O. B. Vitrik, Y. N. Kulchin, V. V. Zhakhovsky, P. A. Danilov, A. A. Ionin, S. I. Kudryashov, A. A. Rudenko, N. A. Inogamov, Laser-induced translative hydrodynamic mass snapshots: non-invasive characterization and predictive modeling via mapping at nanoscale, Phys. Rev. Applied 8 (4) (2017) 044016. doi:10.1103/PhysRevApplied.8.044016.
-  Q. Li, A. P. Alloncle, D. Grojo, P. Delaporte, Generating liquid nanojets from copper by dual laser irradiation for ultra-high resolution printing., Opt Express 25(20) (2017) 24164–24172. doi:10.1364/OE.25.024164.
-  Q. Li, A. P. Alloncle, D. Grojo, P. Delaporte, Laser-induced nano-jetting behaviors of liquid metals, Appl. Phys. A 123(11) (2017) 718. doi:10.1007/s00339-017-1308-4.
-  A. V. Bushman, G. I. Kanel’, A. L. Ni, V. E. Fortov, Intense dynamic loading of condensed matter, Taylor & Francis, 1993.
-  K. V. Khishchenko, S. I. Tkachenko, P. R. Levashov, I. V. Lomonosov, V. S. Vorobev, Metastable states of liquid tungsten under subsecond wire explosion, Int. J. Thermophys. 23 (5) (2002) 1359–1367. doi:10.1023/A:1019821126883.
-  I. V. Lomonosov, Multi-phase equation of state for aluminum, Laser and Particle Beams 25 (2007) 567–584. doi:10.1017/S0263034607000687.
-  Y. V. Petrov, N. A. Inogamov, K. P. Migdal, Thermal conductivity and the electronion heat transfer coefficient in condensed media with a strongly excited electron subsystem, JETP Lett. 97 (1) (2013) 20–27. doi:10.1134/S0021364013010098.
-  S. I. Anisimov, D. O. Dunikov, V. V. Zhakhovskii, S. P. Malyshenko, Properties of a liquidgas interface at high-rate evaporation, J. Chem. Phys. 110(17) (1999) 8722–872. doi:10.1063/1.478779.
-  V. K. Semenchenko, Surface Phenomena in Metals and Alloys, Pergamon, New York, 1961.
-  I. Egry, G. Lohoefer, G. Jacobs, Surface tension of liquid metals: Results from measurements on ground and in space, Phys. Rev. Lett. 75 (1995) 4043–4046. doi:10.1103/PhysRevLett.75.4043.
-  E. B. Webb, G. S. Grest, Liquid/vapor surface tension of metals: Embedded atom method with charge gradient corrections, Phys. Rev. Lett. 86 (2001) 2066–2069. doi:10.1103/PhysRevLett.86.2066.
-  V. V. Zhakhovskii, N. A. Inogamov, Y. V. Petrov, S. I. Ashitkov, K. Nishihara, Molecular dynamics simulation of femtosecond ablation and spallation with different interatomic potentials, Appl. Surf. Sci. 255 (24) (2009) 9592–9596. doi:10.1016/j.apsusc.2009.04.082.
-  V. V. Zhakhovsky, N. A. Inogamov, Molecular dynamics approach to laser ablation in liquid, Appl. Surf. Sci. (2018) sent to journa.