Laser Shock Wave: The Plasticity and Thickness of the Residual Deformation Layer and the Transition from the Elastoplastic to Elastic Propagation Mode

Intense laser radiation leads to irreversible changes in the crystal structure of a target, which are used in laser shock peening technologies. Processes determining the thickness of the residual deformation layer and related residual stresses are studied in this work. It is known that the end of peening is caused by the decaying of the laser shock wave. New information on the transformation of the wave from the elastoplastic to elastic propagation mode under a picosecond impact is obtained. The elastic shock wave is inefficient for peening. The classical configuration with a plastic jump and an elastic precursor ahead of it turns out to disappear during transformation. In this case, the leading edge of the expanding plastic layer gradually decreases its velocity below the bulk velocity of sound, is smeared inside the rarefaction wave, and stops.


INTRODUCTION
Sufficiently intense laser radiation induces numerous hydrodynamic phenomena such as the expansion of a laser jet and the propagation of nonlinear acoustic waves from the impact area. In the case of an ultrashort pulse, the impact area is a layer with the transverse dimension determined by the diameter D L of the illuminated spot on the surface and with the thickness equal to the thickness of the heated layer d T . The supersonic heating for the time τ T t s = d T /c s creates the pressure p ini in the impact area; here, τ T is the time of formation of the layer with the thickness d T and c s is the speed of sound. Acoustic perturbation (shock wave, SW) carries reactive recoil p ini t s from the expansion of the jet (from the unloading of the pressure in the impact area) into the bulk of the target because the momentum of photons can be neglected. The strong SW with the amplitude exceeding the Hugoniot elastic limit (HEL) leaves a trace in the surface layer in the form of residual deformations, i.e., peens the surface layer [1][2][3][4][5][6].
Laser peening is widely used in modern industry to harden goods [1][2][3]7]; it is sufficient to google "laser shock peening." However, the corresponding physics is still unclear. Problems are due to processes resulting in the gradual decrease in the amplitude of the SW below the HEL. In addition, the HEL itself differs from that in shock wave experiments with millimeter strikers (see monograph [8] summarizing these experiments). The HEL in laser SW experiments is an order of magnitude higher [9][10][11][12]. This concerns not only femtosecond pulses [9][10][11] but also subnanosecond and nanosecond pulses [13][14][15]. An increase in the threshold significantly reduces the thickness of the residual deformation layer d plast .
Thus, in the cases interesting for hardening, the amplitude of the SW is first higher than the HEL. The thickness d plast (t) increases during the propagation of the SW. The propagating SW is weakened and its amplitude at the time t fin becomes below the HEL. After that, the expansion of the plastically deformed layer ends having the thickness d plast (∞) = d plast (t fin ), and the SW becomes elastic and disappears completely in the bulk of the target.
The weakening of the SW is due to the rarefaction wave. The shock wave with the rarefaction wave beginning immediately at the front of the SW is called the triangular SW according to its profile. These waves are just induced by ultrashort laser radiation.

OPTICS AND LASER PHYSICS
A system of the SW jump and rarefaction wave varies depending on the geometry. Indeed, the plane onedimensional (1D) geometry with one Cartesian coordinate x along the normal to the surface qualitatively differs [16] from the plane two-dimensional (2D p ) geometry with the coordinates x and y and the axisymmetric two-dimensional (2D a ) geometry with the coordinates z and r, where z is the coordinates along the cylindrical axis and r is the distance from this axis.
The shock wave decays in these geometries at different rates. It is the highest and lowest in the axisymmetric and one-dimensional cases, respectively [16]. However, this difference is manifested only at late times and large propagation distances L SW ≫ R L , where R L is the radius of the focal spot on the surface of the absorbing target. At such distances, lateral rarefaction waves begin to play a decisive role. At small distances L SW ≪ R L (near zone), the 1D approximation is applicable because the diameter of the spot D L = 2R L (1-1000 μm) is much larger than the thickness d T of the laser heated layer (about 0.1 μm for an ultrashort pulse).
In the 1D approximation, the difference between the elastic and plastic cases is masked. The longitudinal stress p xx really acts (elastic case, Hooke's law); it is similar to the dynamic action of the pressure p in the plastic case (Pascal's law, isotropic pressure). This allows one, using the plastic equation of state [17][18][19][20] and adding the shear modulus G, to approximate a laser-induced elastoplastic flow [21] in a hydrodynamic numerical code taking into account the twotemperature initial state (T e ≫ T i ) and the existence of the surface layer of the liquid phase [21].
Noticeable differences between the elastic, elastoplastic, and plastic cases appear in the non-onedimensional geometry. In particular, the trace on the surface of the target from the laser-induced shock wave flow in the bulk is different in the liquid and solid; the target is considered in the form of a halfspace with the irradiated front edge; i.e., waves reflected from the back edge are absent.

PHYSICAL MODEL AND NUMERICAL
SIMPLIFICATION (i) In this work, we performed the molecular dynamics calculations for aluminum, which can be described by the existing well-approved broadband interatomic interaction potential [22]. This allowed us to depart from the usually applied (to describe SW flows [8]) semiempirical elastoplastic models of solid mechanics based on the Mie-Grüneisen equation of state. The point is that, first, the parameters of semiempirical models are poorly defined at spatiotemporal scales of interest in laser problems (e.g., the HEL differs by an order of magnitude [9,10]). Second, it is necessary to join flow parts with evaporation, melting, and dynamically deformed solid in a single flow. It is difficult to perform such joining in solid mechanics. To describe phase transitions, hydrodynamics with multiphase equations of state is used [17][18][19][20], but it is necessary in this case to solve problems with the elastoplastic mode [21] and with non-onedimensional geometry.
(ii) To reduce the load of a multiprocessor computer (involving about 100 cores), the 2D p geometry is considered instead of the 2D a geometry because we are interested in the transition between the near, L SW ≪ R L , and far, L SW ≫ R L , zones. The difference between the plane and axisymmetric geometries is absent in the near zone because the heating disk is thin, d T ≪ R L . Differences between the geometries in the transient layer L SW ~ R L are insignificant [16]. Asymptotic decay laws in the far zone certainly depend on the geometry [16,23]. We do not consider this problem because it will be seen that the transition from the elastoplastic to elastic propagation mode of the SW occurs in the near zone. Furthermore, we emphasize that the problem with the 2D p geometry is of a particular interest for the irradiation conditions with a cylindrical lens [24].
(iii) In the hydrodynamic representation of the flow, the heated layer was a rectangle with the dimensions D L = 100 nm and d T = 20 nm in the y and x directions, respectively. The x × y × z parallelepiped heated in the molecular dynamics code had a dimension of 20 nm in the z direction. Calculations were performed with the heated Gaussian spot and with the heated rectangular spot; the y axis was directed along the surface of the target, and the point (x = 0, y = 0) was in the center of the spot. To distinctly separate longitudinal acoustic waves and transversely polarized waves, we present the results of calculations with the heated rectangular spot. Transversely polarized waves are associated with side effects and propagate from the edges nm of the rectangle. The film sample had the dimensions of 200, 500, and 20 nm in the x, y, and z directions, respectively. Free boundary conditions were imposed at the front (x = 0 at t = 0) and back (x = 200 nm) edges. Periodic boundary conditions were imposed at the lateral edges nm. The number of atoms was 120 × 10 6 . The fcc aluminum single crystal target was oriented so that its [111], , and directions coincided with the x, y, and z axes of the Cartesian laboratory coordinate system, respectively.
The D L × d T heated rectangle was in the center and was adjacent to the front edge. The laser beam axis coincided with the x axis and divided the rectangle into two equal halves. The initial temperature of the rectangle was 17500 K, and the initial pressure was 49 GPa.
(iv) Logarithmic scaling. Since the computational resources were limited, the length of the rectangle D L = 100 nm was approximately an order of magni- [ 1 12] tude smaller than the minimum diameter D L > 1 μm in the experiment. To ensure the relation D L /d T ≫ 1 (as in the experiments), the depth d T was reduced to 20 nm compared to the typical experimental value of about 0.1 μm. Nevertheless, these parameters approximately correspond to the experiment because the HEL in the limit of ultrafast deformations is determined primarily by thermal fluctuations. In this case, the dependence of the HEL on spatiotemporal scales is logarithmically weak. Thus, a change in the threshold, as well as the thickness of the residual deformation layer d plast , under the variation of the scale by an order of magnitude is small.
(v) The description of the distribution of the thermal energy over the depth measured from the surface of the target is simplified. Instead of the solution of the real two-temperature hydrodynamics problem, we suggest that a uniform stepwise distribution of the absorbed laser energy with the thickness d T is instantaneously established in the heated rectangle; i.e., this distribution is an initial condition in the Cauchy problem.
Two-temperature hydrodynamics calculations were performed in our previous works [25][26][27]. The two-temperature relaxation time in aluminum is about several picoseconds, whereas the acoustic scale in the 100-nm-long heated layer is about 20 ps. Consequently, the replacement of the two-temperature stage by the initial condition is acceptable. Moreover, the thickness of the dislocation layer d plast under significant laser loads used in peening is much larger than the thickness of the heated layer d T . Correspondingly, the initial distribution of the absorbed energy over the depth is insignificant and it can be accepted in the stepwise form with the conservation of the absorbed fluence F abs . Figure 1 shows the state soon after the decay of the initial state. The SW passed a path equal to the thickness of the heated layer. All presented plots correspond to the experimental conditions with the absorbed energy of approximately three or four thermomechanical ablation thresholds. This threshold for the bulk aluminum target is F abs | abl ≈ 70 mJ/cm 2 [26]. The initial temperature of the target is equal to room temperature.

TRANSITION TO THE ELASTIC MODE
The amplitude of the SW is insufficient for melting of aluminum (see Figs. 1-3). According to [28], a pressure of about 1 Mbar is required for melting in the SW. For this reason, a fairly sharp transition from aluminum in the supercritical state to plastically deformed solid aluminum occurs near x = 20 nm.
The SW in the target maintains the pressure in the heated layer. After a significant expansion of heated aluminum, the pressure in it decreases strongly and the profile of the SW gradually becomes triangular. Thus, the triangular profile is formed because the action of the piston (the pressure in the hot layer), which supported the SW in the solid, vanishes.
The emission of the elastic wave (ESW in Fig. 4) from the traveling elastoplastic structure is illustrated in Figs. 1-5. In the process of emission, the increase in the thickness of the plastic layer d plast (t) slows down and stops (see Fig. 4). The front of the plastic SW is associated with the beginning of the region of a decrease in the symmetry parameter s. This parameter is the number of symmetrically located pairs of neighbors of each atom. In the ideal fcc single crystal, s = 6.
The deceleration and broadening of the plastic front are presented in Fig. 5. The transverse stress near the front of the plastic SW begins to approach the longitudinal stress (see Figs. 2 and 3). As seen, the front of  the plastic SW is decelerated, broadened, and stops. Such a situation was observed in [29,30], but beyond laser peening. Figures 2-5 illustrate the evolution of defects of the lattice. In Figs. 2 and 3, this evolution is shown by blue lines S. The evolution of the plastic region is shown on the (x, t) diagram in Fig. 4. All these plots refer to the x axis at y = 0, i.e., the continuation of the laser beam axis in the target. Figure 6 shows the spatial picture of defects formed around the impact rectangle. Arrows indicate slip planes. The upper arrow marks a set of planes parallel to the edge of the target. Such an orientation is related to the chosen orientation of the crystal, in which the flat edge of the target coincides with the (111) plane of the crystal. Only slip planes that are perpendicular to the plane of Fig. 6 according to the chosen symmetry of the crystal are seen in this figure. The contribution from the plane that is not perpendicular to this plane vanishes after averaging over z. An (x, y) pixel map for the symmetry parameter was obtained by atomic averaging in the z direction perpendicular to the plane of Fig. 6. Gliding occurs along close-packed planes. As seen, gliding begins in lateral flows, whereas slip planes are not observed in the central (front) region.

SPATIAL DISTRIBUTION OF DISLOCATIONS
The symmetry parameter s = 5.5 deviates from the ideal value s = 6 by approximately 10%. This deviation in the cold crystal is primarily due to dislocations rather than to thermal fluctuations. The shock wave weakly heats the crystal (see Fig. 7). The elastic SW increases the temperature by 50 K. The front of the elastic SW has a width of approximately 20 interatomic layers. In this width, atomic layers vibrate and the longitudinal and transverse temperatures are gradually equalized (see Fig. 7). The temperature in the layer with dislocations increases with the concentration of dislocations approximately by 200 K, but it remains much lower than a melting temperature of 934 K.
Deviation of the symmetry parameter by 10% corresponds to a very high concentration of dislocations of about a percent of an extreme concentration of about 10 14 cm -2 . Figure 8 shows the spatial structure of the wave field in the region between the near and far zones. The    The analysis of works on the laser generation of SWs shows that results for the spatial configuration of the complex wave field (Fig. 8) are new. Most of the previous studies are devoted to the 1D geometry. A. Chemin reported experimental results with the picture of spatial waves at the ANGEL conference in summer 2021 (http://angel-conference.org/). However, these are waves in the far zone, where the asymptotic structure consisting of the compression (с) wave followed by a pair of the shear (s) wave and Rayleigh (R) wave is already formed (see [31]).

LATERAL WAVES IN THE BULK
The compression wave is represented the outer rounded elastic shock jump in Fig. 8 and is shown by the dashed line in Fig. 6. The reflection of the central region of the SW (c) from the back edge begins at the time t = 25.2 ps (Fig. 8). The back edge begins to move to the right in Fig. 8. The reflected wave propagates inside the film. However, the situation near the point of interest and near the Rayleigh wave does not change: the reflected wave is still far.
The results obtained by Chemin develop the results obtained in [32]; see Fig. 4a in [32]. Figure 4a in [32] demonstrates a nucleated bubble, the SW in water, and ρ / = B Π* an acoustic wave in water generated by the fast SW propagating in the target at the velocity much higher than that of the SW in water. The authors of [32] paid no attention to the SW in the target.
The scheme of the c-s-R wave separation is described in [31], where the asymptotic propagation and decay behaviors of c, s, and R waves in the 2D a case are presented. This scheme is typical of an earthquake localized on the surface of a homogeneous halfspace. The c wave propagates first at the longitudinal . The decay of the rectangle into the half-space of the target generates a shock wave and a dislocation field. Symmetry parameter profiles shown by blue lines S in Figs. 2 and 3 were recorded along the laser beam axis corresponding to the vertical passing through the middle of the rectangle. The symmetry parameter in the range from 5 to 6 is given in green to red. Above the last red contour, this parameter is identically equal to 6; i.e., the crystal is not destroyed (see blue lines S in Figs. 2 and 3), although the longitudinal stress on the axis is as high as 15 GPa (see Fig. 3). velocity of sound and a pair of s and R waves follow it. The R wave is slightly behind the s wave because the velocity of the R wave c R is 5-13% lower than the transverse velocity of sound c t . This difference in aluminum is 6%; see the next section.
We consider situation up to the transitional stage, the s and R waves are not yet separated, and they form a joined complex (see Fig. 8). In this work, we con-sider regimes generating intense shock waves, which cannot be described within linear acoustics.
The authors of recent work [33] considered laser impact in the 2D p geometry in the linear acoustics approximation. As in [31], the reported results refer only to the far zone, i.e., after the separation of waves. The authors of [23] applied linear acoustics to describe the mechanics of a deformable solid in the far zone in the 2D a geometry.
A large cycle of works on laser generation of sound were performed during more than two decades at Hokkaido University (Sapporo, Japan), where optoacoustic devices, phonon crystals, and generation of coherent photons were considered [34]; see also longterm studies by Kozhushko, Lomonosov, and Hess [24]. The Rayleigh wave is the most important element of the asymptotic picture of waves [31] because the amplitudes of other waves decrease with time much more rapidly than the amplitude of the Rayleigh wave. This property stimulates applications of the Rayleigh wave in opto-acoustics (see [24,33]). The quasihemispherical elastic SW dominates in the near zone (Figs. 6,8), whereas the Rayleigh wave dominates in the far zone.
The elastic constants in a cold (T = 0) Al crystal for the interatomic potential [22] used in our molecular dynamics calculations (in comparison with reference values at room temperature given in parentheses) are c 11 = 118 (108.2) GPa, c 12 = 62.4 (61.3) GPa, and c 44 = 32.5 (28.5) GPa. The densities at room temperature and at T = 0 are ρ 0 = 2.7 and 2.734 g/cm 3 , respectively. The difference from the isotropic elastic body model is determined by the difference of the parameter = 1.17(1.22) from unity. As seen, aluminum can be described in the first approximation within the isotropic model.
In the case of our crystal, where the [111] direction coincides with the x axis, the (x, y) plane (z = const) coincides with the plane of the fcc lattice of Al. An anisotropic crystal is characterized by two low quasi-transverse velocities of sound c t1 > c t2 and one high quasi-longitudinal velocity of sound c L , see [35]. c t is approximately half the velocity c L . Waves are polarized in the general case of an arbitrary direction of the wave vector with respect to the principal axes of the crystal. For the elastic constants of Al presented above, the velocities calculated by formulas presented in [35] are c t1 ≈ 3.2 km/s, c t2 ≈ 3 km/s, and c L ≈ 6.3-6.5 km/s. The point of crossing П* of the SW with the surface in Fig. 8 moves on the surface with deceleration. The velocity of the point П* at a time of 25.2 ps shown in Fig. 8 is u cross = (7.3 ± 0.1) km/s. Using Huygens' principle and Mach cone/wedge and neglecting the velocity of the flow compared to the velocity of the point П*, we determine the angles of the wedges θ s ≈ arcsin(c t /u cross ) ≈ 22° and θ L ≈ arcsin(c L /u cross ) ≈ 60°, which form with the surface rarefaction waves if they propagate at the velocities of longitudinal and transverse sounds, respectively. We compare these angles with the angles of wedges indicated by two arrows in Fig. 8a.
The angle θ s is in satisfactory agreement with an angle of 22°-24° indicated by the lower arrow in Fig. 8a. The angle θ L is significantly different from a wedge angle of 35°-36° marked by the upper arrow in Fig. 8a. The difference between the angles associated with waves c t1 and c t2 is approximately 1° because the velocities of sounds c t1 and c t2 are close to each other; i.e., the parameter of anisotropy is small. In our calculation, it is difficult to reveal this small difference (related to small anisotropy of Al) if it really exists.
According to the velocity jump behind the narrower wedge θ s (the lower arrow in Fig. 8a), this wedge is related to the softest mode c t2 . This wedge likely does not include the addition of the mode c t1 . The reason is that the mode c t2 is completely polarized in the z direction; i.e., change in the velocity behind the jump occurs only along the z axis, whereas the velocities of perturbations of the modes c L and c t1 are in the (x, y) plane. The molecular dynamics calculation (Fig. 8) shows that the velocity jump behind the wedge θ s has z polarization.
The main SW makes an angle of 97° with the surface at a time of 25.2 ps shown in Fig. 8a. This angle is measured in the same way as for wedge rarefaction waves indicated by arrows in Fig. 8a. Wedge rarefaction waves are generated by the main SW. 7. CONCLUSIONS Processes in the near zone that precede the emission of a complex system of linear elastic waves in the bulk and on the surface in the far zone have been described within molecular dynamics simulation. The case of intense irradiation where the absorbed energy is three or four times larger than the ablation threshold has been studied. In this case, the hot material is emitted in the form of a plasma jet (see Fig. 8) and a com-plex shock wave motion is formed inside the target owing to the reactive recoil from the jet. The molecular dynamics approach has allowed us to remove linear acoustics constraints, which prevent the consideration of processes in the hot area, and to take into account nonlinearities (strong shock waves, dispersion of triangular waves).
The region of crystal modification by a plastic shock wave has been studied. This region is a dense cloud of dislocations and its dimensions are determined by the dynamics of deceleration, broadening, and stop of the plastic shock wave. Deceleration occurs against the background of the rarefaction wave following behind the elastic shock wave.