Identification and Space-Time Evolution of Vortex-Like Motion of Atoms in a Loaded Solid

The paper studies the redistribution of internal stresses and atomic displacements in a preloaded copper crystallite using the molecular dynamics method. It is shown that relaxation within the crystallite volume is accompanied by the formation of dynamic structures in which atomic displacements produce a coherent system of vortex lines. In so doing, the displacement of atoms in neighboring vortex structures has the opposite sign of the angular velocities. The evolution of the dynamic vortex structures is analyzed using an original technique for identifying the vortex motion in the space of a vector variable with a discrete step. It is shown that a system of dynamic vortices and antivortices can propagate inside the crystallite, ensuring the transfer of stresses from the bulk of the loaded material to its unloaded periphery in order to preserve continuity. The developed technique has revealed that the lifetime of such defects depends on their size and ranges from fractions to tens of picoseconds. The simulation results correlate well with the experimental electron microscopy data on the estimation of spatial parameters and lattice curvature during strain localization in the region of elastic distortions.


INTRODUCTION
The vortex-like and rotational motion of media is a topic of relevance across a vast number of applications, including dynamics of plasmas [1,2], modon structures in geophysics [3], astrophysical flows [4,5], the dynamics of quantum condensates [68], and the skyrmion magnetic structures [911], to name a few. At the same time, the appearance of circular motion in such fields of modern physics as physics of strength and plasticity still remains limited. Recent studies have shown that the vortex motion including the rotational mode of deformation is a very important object that in some cases can ensure the preservation of continuity and relaxation of a loaded body [12]. Indeed, various real and virtual experiments [1316] confirm the presence of grain boundary rotation along with other plastic flow mechanisms in nanocrystalline materials. Moreover, in the cases when the response of the loaded body becomes insufficient to provide necessary stress relaxation, the vortex-like motion becomes a precursor of the formation of a material discontinuity [17,18].
In addition to the conventional violation of the crystalline order, vortex-like displacements of atomic groups can be treated as dynamic defects because of their dynamic behavior. These defects correspond to a self-consistent circular motion of a large number of material points and exist only at some stages of loading, determining the consistency of displacements in different parts of the solid body. According to Refs. [12,18], the deformation of a heterogeneous material containing internal interfaces or free surfaces is accompanied by a collective circular motion near them. Such vortex-like patterns of elastic displacement are normally associated with the Rayleigh, Lamb, or Love surface waves, or Stoneley waves propagating along a PHYSICAL MESOMECHANICS Vol. 21 No. 5 2018 phase boundary. .or example, it was demonstrated by molecular dynamics simulations [19] that vortex structures in the velocity field are formed around grain boundaries of polycrystals under shear loading and can provide grain boundary migration. The nature of the vortex-like motion lies in the formation of local velocity gradients and stresses tangential to them near internal or external interfaces. It can therefore be expected that circular motion takes place at different scales: from atomic to macroscopic. Psakhie et al. [20] have shown that the vortex-like motion becomes a precursor of discontinuity formation in a loaded material on the macroscale. On the other hand, the theoretical results and experimental evidence of Zhang et al. [21] indicate that the contribution of the rotational mode of deformation can significantly increase under dynamic loading, especially in nanomaterials with a high concentration of grain boundaries. Nevertheless, such a fundamental factor as the elastic circular motion in a material under dynamic loading is still not fully explained. The above suggests that revealing the role of vortex displacement in elastic energy redistribution and, as a result, in the deformation and fracture of heterogeneous or homogeneous materials near free surfaces is a particularly important problem in materials science.
Although a direct observation of such defects is impossible because of their extremely small size and short lifetime, they strongly influence the entire deformation process at both the micro-and macroscale. In view of the principal significance of the free surface, internal or external interfaces, and dynamic nature of the considered vortex phenomena, the main approach to studying such defects is computer modeling [22]. There are difficulties in identifying such defects and especially in observing their time evolution. It is very difficult to visualize and identify the vortex motion of a discrete randomized system of a large number of particles. The most interesting task is to understand the role of vortices at different stages of the evolution. This would help us to find regularities in the motion and relate them to the general behavior of the material. We would better understand the relationship between different scale levels.
Jiang with coauthors [23] paid much attention to the question of identifying the vortex motion. They emphasize that it was a challenge to convert an intuitive description of the vortex into a formal definition. The authors reviewed existing detection methods and particularly discussed nine methods that are represen-tatives of the state of the art, each of which has its own advantages and drawbacks. The discreteness and randomized nature of a system of moving particles raises additional difficulties in the identification. However, there are a number of approaches that allow such analysis. .or example, Tordesillas et al. [24] analyzed the incompatibility of displacements of individual particles. Their approach revealed the existence of different vortices and allowed the estimation of the characteristic vortex sizes; some rules were established for the deformation and distortions in granular media.
Up to now it has been assumed that dynamic defects exist only in the active stage of loading, when forced displacement velocities are much higher than the typical time of conventional defect generation. At the same time, many fundamental studies consider the rotational mode of deformation as an important mechanism of internal stress relaxation [25]. The present paper focuses on the possibility of formation of dynamic vortex defects at the stage of internal structure relaxation in the material, i.e., after the stage of active loading. Molecular dynamics modeling is performed to study the deformation of a copper crystallite compressed to the deformation immediately preceding plastic flow. Possible vortex motion of a particle ensemble is identified using an original approach based on the vector field calculation in each point of the available space.
The paper is organized as follows. Initially, we describe the numerical model and the technique of vortex motion identification. Then, the numerical modeling results are combined with the application of the main approach, and the main features of the found vortex motion of particles are described. The final part presents conclusions of the study.

MODEL DESCRIPTION
The initial trial structure of the model crystallite without defects is shown in .ig. 1. It has the form of a parallelepiped measured 21.7 × 21.7 × 3.6 nm along the [100], [010], and [001] directions that respectively contain about 150 000 atoms. The crystallite was conditionally divided into three regions. Region II is a freely deformed part of the crystallite where new states of the particle ensemble were found from the Newtonian equations of motion. Atoms of regions I and III were exposed to an external influence using a two-stage loading scheme. This scheme includes an active and passive stage. The type of load at the active loading stage corresponded to uniaxial compression with the socalled string boundary conditions [26]. These conditions mean that in the [010] direction the atomic velocity projections in regions I and III are fixed, but in the directions different from [010] the velocities are determined by the corresponding atomic environment.
Such a type of boundary conditions is preferable for molecular dynamics modeling of a selected crystal lattice fragment. In contrast to conventional uniaxial compression where displacements of particles in the directions different from the loading direction are equal to zero, spring boundary conditions are used and lateral expansion under compression appears also in loaded parts of the system. The rate of deformation under compression reached up to 0.5% per picosecond. To achieve this, at each time step we assigned the values of 50 m/s and 50 m/s to the Y velocity projections of atoms from regions I and III, respectively. Loading continued to the deformation immediately preceding plastic flow. The total deformation reached in the active loading stage was varied between 5 and 10% in different calculations.
At the second loading stage (passive stage) the positions of atoms along the Y axis in the loaded layers were fixed to preserve the deformation achieved at the active stage. Here we modeled lattice relaxation based on the Newtonian equations, but without any additional sources of influence. In other words, a so-called microcanonical NVE ensemble with fixed total number of atoms, volume of the system, and its energy was applied.
Taking into account the realistic extension of the modeled fragment, we used periodic boundary conditions in the [001] direction and free boundaries in the [100] direction. To exclude the induced effect related to the symmetry of the ideal lattice, the copper crystallite was studied at the initial temperature 10 K. The final temperature was found a posteriori using the Maxwell distribution for the atomic velocities. Detail analysis of the relaxation process was carried out by studying the evolution of atomic configurations at different points of time and for different time intervals between them.
Molecular dynamics has been adopted as the most suitable method for simulating the acoustic response to indentation. Large-scale atomic/molecular massively parallel simulator software was used as a computational tool [27]. The interaction between atoms in the copper crystallite was described using an interparticle potential reconstructed in the framework of the embedded atom method [28]. The physical correctness of the software was verified by matching the computed and experimental values of the vacancy and stacking fault formation energies. The mismatch achieved was less than 1%. The visualization of the molecular dynamics simulation data and structure analysis were carried out using the open visualization tool OVITO [29].

DE.INITION O. VORTEX MOTION
Observations show that there are vortex-like displacements of atoms in the modeled crystallite. As an example, .ig. 2 illustrates the displacements of atoms of the modeled crystallite for different time intervals in the plane parallel to the (001) plane. .or better visibility, the line segments corresponding to the atomic displacements are enlarged 5 times. One can easily see from the framed region in .ig. 2a, where such displacements are shown for the interval between 8.4 and 8.9 ps, that 8 vortex lines are formed in the crystallite volume. The motion of atoms in all the lines takes place around the axes parallel to [001]. The diameter of the vortex lines is about 1520 distances between the crystal planes.
The rotation directions of the vortex lines are correlated as follows: the nearest pairs of the vortices rotate with opposite signs of the angular velocity and therefore the motion of atoms is correlated so that the continuity of the material is preserved. Our calculations have shown that such time intervals can be selected during relaxation in which the system of 8 vor-  is also a discrete array. .or further use, it is convenient to interpolate the array into a continu- The area of the continuous variables [x, y] must be chosen to completely cover the area of the discrete variables { , }.
k k x y The mutual correspondence of the surface J(x, y) and the discrete set can be controlled visually.
As an example, a particular implementation of the absolute value of the current J(x, y) = |J(x, y)| is depicted in .ig. 3a. .or further use, it is convenient to plot the value J(x, y) in the form of a planar density distribution shown as a grayscale density map in .ig. 3b, which clearly illustrates the fine structure of the current including the minima of the current J(x, y).
Mathematically, the points where J(x, y) turns to zero are indicators of probable positions of vortices. In some cases, these points correspond to saddle configurations of the current and hence they are not centers of the vortices. So, to find real vortices, we need to calculate the circulation of current around such points. However, since the field J(x, y) is defined as an interpolation of a randomized data array, exact zeros of the current J(x, y) cannot be reached. Taking this into account, it is convenient first to determine the circulation around every point of the regular array and only then to interpret the results obtained.
The circulation of displacements around all points of the array can be calculated by defining the radius of its contour R, which must belong to intermediate scales between the main characteristic distances of the system. The minimum distance is the mean distance between the nearest neighbors of the array { , }, k k x y hence mean(min | |). the cases, when the contour formally extends beyond the regular array [x, y], where data about J(x, y) are absent, the contour must be plotted along the corresponding boundary from one of the points of its intersection with the boundary to another. Numerically, the contour is a discrete data array with radius R and angle 0 2 , n ≤ α ≤ π where the subscript n = 1, 2, ..., , N α and the number of segments must be much larger than the unit 1 N α >> but is limited by the volume of available numerical information . with the respective value of the current ( , ). x y = J Now the study is reduced to the extraction of quantitative and visual information from consequently calculated distributions of the circulation C = C(t; x, y), which depend on the time t and regular coordinates [x, y].

Identification of Vortex Centers
.irst of all, we must formalize the calculation of the nominal positions of the centers of vortices and antivortices. .ormally, we must define the positions of (all) local extrema of the function C = C(t; x, y) for each time point. However, the direct use of this procedure in practical calculations can lead to artifacts related to some degree of randomization.
The above problem can be avoided if we make preliminary spatial smoothing of the function C = C(t; x, y) for the scales smaller or close to the characteristic scales of vortices. It is convenient to perform smoothing by a Gaussian convolution. Particularly, we can perform the fast .ourier transform C = C(t; x, y) → ( ; , ), x y C t q q multiply the obtained distribution by the Gaussian 2 2 2 ( , ) exp( ( ) ) x y x y G q q q q = − + ∆ with a given width ∆ in the .ourier space, define the new function ( ; , ) ( ; , ) ( , ) x y x y x y C t q q C t q q G q q = and perform the inverse transform ( ; , ) ( ; , ).
x y x y C t q q C t q q → Numerical check shows that for the absolute majority of the C(t; x, y) configurations such smoothing preserves all properties of the original distribution if ∆ is correctly chosen. It allows one to uniquely define the position of maxima and minima (i.e., centers of vortices and antivortices) correctly corresponding to the structure of the displacement current (see .ig. 4).
Using the described procedure, we can both visualize the positions of all vortices in the system and  A direct observation of the appearance and disappearance of vortices on numerical visualizations of spatially distributed circulations have shown that even in the cases when the total value of C Σ for the whole system is equal to zero, ( ; , ) 0, S C C t x y Σ = = ∫ its partial values for separate quadrants of the system can differ from zero. In other words, vortices and antivortices are generated separately in different parts of the system, but annihilate near its center.
This process can be quantitatively described by calculating the partial value of the total circulation C(t; x, y) for each quadrant: (c) (a) evolution of the system are represented in .ig. 5. The above described arguments, which relate to the nonbalanced appearance of vortices and antivortices in different parts of the system, their mutual compensation in different parts of the system, and finally annihilation in particular time points, are visualized by the curves ( ) C C t Σ Σ = (.igs. 5a and 5b) and igs. 5c and 5d).
The most valuable result of the described approach is perhaps the existence of the scalar but absolutely transparent language for the description of complex rearrangements in the compressed system, which is based on describing the position of the vortex centers, their total number ( ), C t Σ and the number of vortices in separate parts of the system 4 ( ). C t Σ The proposed approach provides an additional bonus. If the positions of the vortex centers are known, their motion can be tracked from the moment of their emergence to meeting and annihilation. In other words, the vortex motion can be described in a static manner as a map of their trajectories on the plane of the specimen cross section. It excludes long qualitative speculations and provides a quantitative language for further description of very different systems of this type.
As a result, we obtain a simple algorithm for processing experimental or numerical results for the system under external pressure, which includes: (i) interpolation of the distribution of discrete particles and their displacements to a regular data array, (ii) calculation of the circulation for every point of the obtained array using a fixed radius comparable with the characteristic scale of visually observed vortices, (iii) determination of all extrema (after preliminary smoothing of the circulation, if it is necessary) that define the positions of vortices and antivortices, (iv) calculation of the total number of the circulation extrema, their balance in the whole system, and their mutual redistribution in its parts, (v) determination of the spatial motion trajectories of vortices and antivortices from the moment of their emergence to annihilation, and finally derivation of a complete description of the system vortex state in a closed static form.

MODELING RESULTS
The proposed approach to identifying dynamic vortices has been tested using a time-space analysis of defect motion in an originally ideal copper crystallite at the active stage of uniaxial loading. Some examples of vortex and antivortex trajectories from the moment of their appearance to annihilation are presented in .ig. 6. In .ig. 6a the vortex center trajectories are shown in the interval from 0.5 to 5.5 ps which starts exactly from the beginning of relaxation. It is easily seen how the diagonal vortexantivortex pairs are generated in the regions away from the center of the crystallite. Later they move towards the center and annihilate. The defects from the left side of the crystallite annihilate with the opposite-sign defects from the right side. The trajectories of the dynamic defects in the interval between 5.5 and 9.5 ps are shown in .ig. 6b. There is one dynamic defect (type 1 defect) in each quadrant which is generated in the regions maximally distant from the center and moves gradually to the central region. In contrast to the previous time interval, the annihilation of such defects takes place in the vertical, not horizontal, direction. The dynamic defects from the bottom part of the figure annihilate with the defects from the upper half. Another type of defect (type 2 defect) is generated in the bulk of the crystallite and moves to  It should be noted that in some cases the vortices and antivortices moving in the bulk of the crystal lattice can elongate and split into two ones. Later, one of the centers disappears and the whole vortex jumps into a new point in space. This effect looks like quantum tunneling. The corresponding parts of the trajectories are shown by the dashed lines (see intervals A and B in .ig. 6b). The fine structure of such dynamic transformations is presented in .ig. 7. The moments of deformation and fast spatial motion of the vortex cores are marked by the arrows. Now let us analyze the evolution of the total number of vortices and antivortices. Generally, we find that there is often only one vortex in each quadrant, or two of them with different signs of circulation. The corresponding data are represented in .ig. 8 showing the time dependences of the number of vortices and antivortices (.ig. 8a) and their difference (.ig. 8b). It is seen that the difference between their numbers is normally equal to 0 or 1, and rarely reaches 2 or more. The number of vortices is obviously limited by the  small size of the system. Their deformation and core splitting can influence the formally calculated number of vortices (probably with higher original circulation). However, taking into account these limitations, it is possible to make general conclusions about the vortex motion of atoms at the active loading stage.
Vortex motions promote the transfer of internal stresses/displacements from the crystallite bulk to the open surfaces. The stresses redistribute either because of the formation of solitary defects or due to neighboring defects of different signs which rotate like bearings and transfer the stresses/displacements in a zigzag manner. This process preserves continuity of the system and is schematically shown in .ig. 9a.
Comparison of atomic trajectories in the bulk and near the open surface also illustrates the redistribution of stresses/displacements. In particular, the atomic trajectories in the bulk of the crystallite for two consecutive time points are shown in .ig. 9b. One can see that the atom moves almost cyclically along a curve around its equilibrium position. The initial position of the atom is marked by the circle. In contrast to the internal atoms, the same trajectories for the atoms close to the open surface (see .ig. 9b) are more extended and preferably oriented in the X direction of the main transfer of stresses/displacements from the center to the periphery of the material.
A detailed analysis revealed a definite correlation between the vortex and antivortex sizes and their life cycle duration. It has been established for a particular system that the sizes of the observed vortices vary from 56 to 2025 interplanar spacings. In this case, the vortex size is used to mean a half of the distance between two neighboring defects. Close defects with the distance between their centers smaller than 15 inter-planar spacings are treated as a solitary defect with a split.
It has been found that larger size defects live longer. Their maximum life duration is about 10 ps. The position of such a defect is not fixed and changes in space during this time interval. Small size defects usually exist for fractions of picoseconds and form a pair with an opposite-sign defect, as it is shown in .ig. 10.
Thus, short time intervals and small defect sizes allow the defects to be treated as a support tool that helps the main vortices and antivortices to perform the redistribution of stresses and displacements in order to preserve the material continuity.

CONCLUSIONS
Although recent data indicate that rotational deformation occurs at different scales (from atomic to mac- roscopic), this study focused on some aspects of a collective rotational motion of atoms. Molecular dynamics was used to model the evolution of a preloaded copper crystallite at the stage of relaxation. The results demonstrate the possibility to generate a coordinated vortex-like motion of atoms with the opposite signs of circulation in conjugated regions and alternating directions of rotation over time. In view of their dynamic behavior, such structures can be treated as dynamic defects. Despite the fact that these structures are within the elastic range, they can lead to non-elastic consequences and induce the formation of conventional static defects in the crystal lattice. Such defects are formed, in particular, by increasing the duration of the active loading stage (not considered in this paper), mainly in the zone of conjugation of loaded layers I and III with a free surface.
The experimental detection of vortex dynamic defects is objectively difficult because of their extremely small sizes and very short lifetime. However, recent electron microscopy studies of deformation localization in the region of elastic distortions during torsion in Bridgman anvils of nickel have proved the validity of numerical modeling results [30]. This suggests that numerical simulation methods are an effective way to study the role of circular movement in a loaded solid. There are, however, technical difficulties in the practical identification of vortex motion due to a discrete distribution of the vector field on randomized arrays in space. In this study, we proposed an original approach to solving the problem. It is based on the interpolation of the field to a regular array, calculation of the circulation in every point of the array, and then the use of a Gaussian convolution to smooth the results. This technique allows both visualizing the vortex motion of particles in space and tracking its evolution in time.
The proposed approach was used to determine the characteristic sizes of vortices and antivortices (circulations of atomic displacements of opposite signs) and their lifetime duration. Possible vector field transformations in space and preferable directions of motion of vortex cores were determined. It was found that the vortexantivortex pairs can annihilate if they approach close to each other. Some effects of vortex core splitting and instant vortex jumps to new positions at large vector field deformations were observed. In the case when a number of relatively close vortices with opposite circulations appear in the proximity of each other, an additional effect of a roller transfer mechanism of atomic displacements can arise. In this case, particle displacements along boundaries between vortices become correlated, and the self-consistent flow of particles transfers internal stresses from the bulk of the crystal to the open boundaries on its periphery.
Limitations of the proposed approach are the following. The accuracy of identifying the vortex motion is limited by the natural step between the points of the array corresponding to the particles for which the vector variable must be applied. The Gaussian convolution is also limited by itself and makes it impossible to identify small-sized and short-living dynamic defects. Increasing the algorithm accuracy in this limit sometimes leads to phantom centers of collective rotational motion with a nonzero but extremely low angular velocity. The last problem implies further modification of the proposed approach to parameterize it correctly depending on the angular velocity in a wide range of its values. The specific choice of boundary conditions also influences the value of the spatially distributed circulation obtained in the close vicinity of the crystallite boundaries.

ACKNOWLEDGMENTS
The work was financially supported by the .undamental Research Program of the State Academies of Science for 20132020 (Project III.23.2.4). Molecular dynamic modeling was carried out with the financial support of RS. Grant 17-19-01374. The supercomputer of Tomsk State University was used in the framework of the TSU competitiveness improvement program. The development of the algorithm to visualize the circular motion of particles in space was supported by the German Academic Exchange Service (DAAD). We express our gratitude to Prof. S.G. Psakhie for a fruitful discussion of the results.