Size effects in stress propagation and dynamics of dislocations: Fe–Ni–Cr steel

Movement of line dislocations in fcc steel 310S is found to depend on the size of nanometer sized structures, based on modeling within molecular dynamics (MD). The effect is attributed to time (and size) dependencies of pressure propagation into the medium interior. The observation is crucial in interpreting any MD studies of pressure effects since these are governed by time-dependent internal virial stresses. In particular, velocity of dislocations scales better with value of local internal shear component of virial stress S xy than with external shear pressure applied. Dynamics of stress penetration is described well within the model of damped harmonic oscillator, where characteristic oscillation frequency depends on number of crystallographic layers in direction along the wave propagation while the speed of stress propagation is the speed of sound. The minimal stress required for dislocation movement (Peierls stress) is determined to be 0.75 GPa. Pressure and temperature effects on dislocation movement are systematically investigated.


Continuum interpretation of virial stress in molecular simulations
The multiscale constitutive modeling numerical simulations are playing an important role in the understanding of the fundamental mechanisms governing microstructural evolution. The * Author to whom any correspondence should be addressed.
Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. problem of scale bridging between molecular dynamics (MD) and continuum mechanics analyses remains challenging, hindering the simultaneous study of physical processes at the atomic levels, especially in irradiated systems.
In the irradiated materials, all the processes take place during and soon after the interaction of energetic incident particles with lattice atoms. These processes are experimentally unobservable because a displacement phase of the collision cascade usually lasts picoseconds. The only approach that may address this issue is the use of MD simulations.
The efforts to develop validated methodologies capable of predicting microstructural evolution and mechanical property changes are hampered by problems associated with the stress definition and its dynamics, applicable to both continuum and discrete systems.
Stress is one of the most fundamental quantities in continuum mechanics, it is essential to introduce a stress definition applicable to both continuum and discrete systems. The most common one in discrete-particle systems is based on virial theorem. The virial stress consists of at least two components: a kinetic one depending on the atomic particle mass and velocity and a potential component depending on the interatomic forces and atom positions: where k and l indexes atoms in the volume V, m k is the mass of atom k, v k i , x l i is the ith component of the velocity and position of atom, while F kl j is the force along direction j on atom k due to atom l.
It has been reasoned that the virial stress represents an atomistic definition of stress that is equivalent to the continuum Cauchy stress [1], provided spatial and temporal averages are computed properly [2,3], while the problem attracted a broad discussion [4,5].

Stress propagation into a medium
There are two basic approaches towards modeling of stress effects in materials with the use of MD [6]. Either an external pressure is applied and its effects are studied as a function of time (so called dynamical simulations) or a continuous, changing linearly in time ramp force is applied (so called quasi-static simulations). In any of these cases we deal in fact with a dynamic situation and any interpretation of computer modeling results must be based on understanding the time scales involved and spacial distribution of stress penetration.
Longitudinal c L , and transversal c T speed of stress (sound) propagation is given by (see e.g. [7]): c L = E/ρ 1/2 , and c T = G/ρ 1/2 , where E and G is Young modulus and shear modulus, respectively, and ρ is material density. For finding elastic constants we used a package written by Aidan Thompson available in source code distribution of LAMMPS [8], 1 which is based on the modeling concept of Sprik et al [9]. At low temperatures for Bonny potential [10] elastic constants, C 11 , C 12 and C 44 , are, respectively, 327.06, 189.27 and 156.9 GPa. Poisson ratio obtained is 0.36. At T = 300 K density is 8.09 g cm −3 , which corrected for real atomic content of steel 310S by a factor 0.989 gives the density 7.999 g cm −3 at 300 K.
With computed elastic parameters we find the values of c L and c T at T = 0 K: c L is 6308 m s −1 and c T is 4368 m s −1 . Both, c L and c T almost do not depend on temperature and are of the order of 63 Å ps −1 and 44 Å ps −1 , in a reasonable agreement with known values for steel and iron [11].
We distinguish between microscopic pressure tensor components as reported by LAMMPS and pressure values applied at the surface, which is related to surface average of forces. Additionally, we ought to distinguish between the pressure reported by LAMMPS at the surface, and the external pressure applied to that surface. The last one will be denoted with superscript 0. Hence, P 0 xy (t), the applied shear pressure caused by force in X-direction on surface Y is in general not the same as the pressure exerted on the surface by sample interior, P {xy} (t), at any given moment of time.
When pressure is applied abruptly at the surface in a form of Heaviside function, P 0 xy = 0 for t < 0 and P 0 xy = const at t > 0), it propagates inside towards the opposite surface as a wave traveling with sound speed. It reflects from that opposite side with negative amplitude (that surface is fixed and it is movement is not allowed). The reflected wave interferes with the incoming wave forming a complex pattern of pressure inside of sample volume. A situation like that has a well known analytical solutions in case of a medium that is characterized by linear response. For instance, time dependence of amplitude, a(t) of sound wave is a solution of second-order differential equation and may be written as [12][13][14][15]: with phase ϕ given by cos ϕ = ζ, where ω = (k/m) 1/2 is called the undamped angular frequency, and ζ is called the damping ratio and it is related to energy losses. Quantities k and m in a simple classical model are spring constant and oscillating mass.
In case of a linear in time ramp pressure applied, P 0 where h is the pressure rate change in time, response of a linear medium is known analytically and, when to consider first order response only, it is given by: where P xy is component of pressure tensor as reported by LAMMPS. Equations like these above, (2) and (3), are often found in textbooks on automation and they describe a broad range of phenomena, mostly in engineering.
The above equations, written for xy component, have a more general nature, though parameters such as τ above may depend on components, due to material anisotropy and sample geometry. Hence, in case of linear ramp force applied, as it follows from equation (3), there is a time lag of pressure at the surface that additionally decays exponentially with time. equation (3) is however a solution of first order differential problem, with second order contributions (leading to oscillations) omitted.
We have performed a series of simulations for the applied pressure rate changes h from between 1 MPa ps −1 to 200 MPa ps −1 and found that surface pressure as reported by LAMMPS is well described by the function: For sample of 100 Å size in Y-direction the period of oscillations is around 10.7 ps, i.e. ω is around 0.6 ps −1 , while value of A is about 0.1 ps, when time step used in MD simulations is 0.001 ps. Similar oscillations are often found in MD studies when a linear in time load is applied [16].
These results show us that in case of a ramp pressure applied, P xy (t) is close enough (in a reasonable experimentally sense) to applied P 0 (t) for very low values of h only, of around just 1 MPa ps −1 or less, which is hardly computationally convenient for MD simulations that usually require a lot of testing yet in preliminary phase of any studies. On another hand ramp rates of the order of 100 MPa ps −1 are too fast for achieving quasi-stable physical conditions (dynamical stresses) inside the sample volume.
In this work we concentrate on using the first method of modeling, i.e. that one when pressure is applied abruptly at t = 0. Our aim is gaining a deeper understanding of stress dynamics and timescales involved for a proper design of the simulation setup and later analysis of results.

Creating oriented steel 310S with dislocations
For creation of samples of steel 310S (fcc structure) we used either LAMMPS [8] or Atomsk [17], as well custom-written Perl scripts. The weight contribution of Fe-Ni-Cr atoms is 0.55-0.20-0.25, as in specifications 2 . After testing available interatomic potentials ( [10,[18][19][20][21]) we choose the EAM potential of Bonny et al [10] as the one that reproduces correctly basic physical properties of the material and as the most efficient one computationally. Moreover, that potential has been developed to model defects in steel of similar composition and therefore is suitable for our planned future modeling. Some results were computed by using Artur potential [18] as well.
Samples have been imposed to potential energy minimization, checked for density and elastic constants at low temperatures and at room temperature (most of results we report on here were computed at T = 50 K).
It is important to include details of samples geometry and computational methodology since, as it will become evident, these have profound impact on results obtained and on their interpretation.
We use a typical [22] simulation setup applied in similar modeling, as shown in figure 1.
For sample sizes and positions of dislocations we use the following convention. Sizes in Y-direction are multiplicity of fcc unit cell in that direction. The sample is oriented, with X = [110] (which is direction of the Burgers vector b), Y = [111] (which is direction normal to the glide plane), and Z = [112]). That unit cell in Y-direction has size of 6.18 Å and it consists of atoms of hexagonal structure arranged in three layers in X-Z plane separated by 2.06 Å in Y-direction. Names of samples depend on sizes in Y-direction. M + N means that the position of dislocation is M unit cells from sample bottom (i.e. where Y = 0) followed by N unit cells to the surface, where pressure P 0 xy is applied. Hence, samples used here for dislocation dynamics studies have the geometries 2 + 2, 4 + 4, 8 + 8, 16 + 16. Additionally, there is a series of calculations on samples 8 + X, where X = n · 8 and n ranges from 1 to 7, and a sample named 48 (total size in Y direction) which is in three variants: 24 + 24, 16 + 32 and 32 + 16.
The shear surface stress P 0 xy is created by applying a force to atoms in the upper region while on the bottom surface of the lower region conditions are imposed of zero forces and velocities of atoms. The force F applied to atoms must be oriented in X-direction and have a value such that when summed for all atoms in the upper region and divided by the surface area in X-Z direction it will result in desired value of pressure P 0 xy . We find by testing that the width of upper and lower regions in figure 1 may be made as small as the distance between atomic layers in Y-direction, i.e. that these regions contain one only layer of atoms, not affecting quality of temperature control.
Temperature is stabilized with NVT thermostating in upper and lower sample regions only, while the interior of the sample remains under quasi-adiabatic conditions, with some only heat exchange with these two temperature-controlled regions.
Using NVT in entire sample volume would affect the dynamics of dislocations. Using NVE in whole sample volume would instead lead to temperature changes, due to energy losses, in particular when dislocations are present. For these reasons, in order to be able to simulate the effect of temperature on dynamics, we choose a method which is a mixture of NVT and NVE: with NVT at upper and lower regions and NVE where dislocation's cores are located between them.
Average sample temperature of entire sample volume remains constant, with RMS fluctuations not exceeding 0.7 K, for any sample size, at T = 50 K, provided that no dislocations are there and applied pressure is below a few GPa. Otherwise strong temperature instabilities may occur. External pressure is transferred to sample interior through NVE integration updating atoms positions and velocities.
For saving on data storage space and data analysis time it is sometime convenient to dump data from a narrow region only, named middle, which encompasses the core of dislocation. As a time step in MD simulations of dislocations movement we used mostly 0.001 ps, a time which is the largest possible to speed up the computation and is still not disturbing the accuracy of modeling. In some cases, as explained, another time step was used.
Dislocations of line type were created with the help of Atomsk or by displacing atoms with the method of Carpio and Bonilla [23] implemented ourselves. An overview of the methods of dislocations creation is found, e.g., in [24]. When using Atomsk, a perfect 1 2 [110] dislocation with a stacking fault is obtained by joining two crystals in Y direction that have length change by ±a/2 in oriented X direction, where a is lattice constant in that direction. That dislocation splits into two Shockley partials with Burgers vectors 1 6 [211] and 1 6 [121] [25] repelling each other and being attracted by stacking fault between them. In result a stable equilibrium position is reached when they are separated for about 66 Å in X-direction. Under stress they both move the same way, with some only fluctuations of the distance between them.
Dislocation visualization was done with Ovito [26], while for extraction of dislocation data Ovito were used, Ovitos (Python extension libraries to Ovito), Crystal Analysis tool 3 [27] and/or our own developed scripts in Perl.
Statistical averaging of physical quantities along X direction was done in Perl as well, with resolution of 2.06 Å in Y direction, which is the distance between neighboring crystallographic planes in that direction. Data analysis was often accompanied by creating movies with process automation utilizing bash, Perl, gnuplot scripts, as well ffmpeg.

Pressure response of sample
The simplest way of studying stress dynamics is the measurement of stress at the sample surface after external pressure is applied, as shown in figure 2. Strong oscillations of P {xy} (t) are well  (5) or (6). The applied pressure at the surface, P 0 xy , is 1 GPa. In case of curves at the figure (b) the damping parameter ζ approaches (Y = 8) and exceeds the critical damping ratio (ζ = 1). The curve for Y = 12 is in overdamped regime. In this case we see that green lines, as computed by equation (6) start to depart from MD simulation results. That difference increases when ζ grows more above 1. described by the dumped harmonic oscillator response to a step time pressure change, as given by equation (2). In this case pressure P {xy} has a negative sign as a response to positive pressure applied at the surface, P 0 xy . Hence, to fulfill the boundary conditions let us rewrite equation (2) in the following form: Parameters ω and ζ are material specific. The damping parameter ζ depends also on details of the MD simulation process (in particular, it is found to be proportional to time-step). Angular frequency ω depends also on sample size in Y-direction. The all three green curves in figure 2 have been computed by using the same values of A and ω, and values of ζ scaled proportionally to time-step ts.
Parameter A in equation (5) has a dimension of pressure and it scales up linearly with applied pressure P 0 xy (at least for P 0 xy not exceeding a few GPa), but its value is lower than P 0 xy . It approaches P 0 xy only for ζ 1. This is seen better in figure 3, where P {xy} (t) dependencies are shown for samples of different sizes, from Y = 1 to Y = 12.
The damping parameter ζ in equation (5) must be lower than 1, since we have a 1 − ζ 2 factor there. In case of ζ > 1 the above equation is replaced by a one usually written as a sum of two exponential contributions. By using relations: sinh(x) = −i sin(ix) and  (5) and (6) on number of atomic layers in Y-direction, N, as deduced from analysis of MD simulations. The data points in green (full circles and curves) are obtained for simulations using Artur potential, while these in red and blue for Bonny potential. Data points and curves in blue are for direction [112] while all other for Y in direction [111].
, we may formulate equation (5) for ζ > 1 in a similar form: As seen however in figure 3, for large samples (for ζ > 1) the above equation does not reproduce exactly the results computed in MD simulations.
Equations (2), (5) and (6) are obtained as solutions of the second order differential equation. Therefore there must be two parameters for satisfying general boundary conditions. One of them is amplitude A of obvious physical meaning, while the second one is related to the initial velocity of wave propagation to the medium at the surface. This is why parameter ζ is found to depend on the time-step in MD simulations.
Detailed dependencies of ω, ζ and amplitude of oscillations A on number of atomic layers in Y-direction is shown in figure 4. The number N is simply multiplicity by 3 of the number of unit cells of [111] oriented fcc crystal.
Amplitude of oscillations does not depend on potential used or on crystal orientation and it has 1/N relationship, as shown in figure 4. Damping coefficient, shown there, is for MD timestep of 0.001 ps and it retains linear dependence on N, regardless of potential used or crystal orientation. It is slightly only different for potentials of Artur and Bonny, and it depends on crystal orientation. The largest change is found in ω: larger values are found for Bonny potential, as expected, since that potential produces a larger values of elastic parameters. Angular frequency ω depends on crystal orientation as well.
Since amplitude of oscillations is inversely proportional to N, it reaches 1 only in the limit of infinite size of sample, while for small sample sizes it is significantly lower than one. That means that volume average pressure value inside of a small sample is always lower than the pressure applied at sample surface.
Dependence of angular frequency on N, ω(N), may be interpreted within a simple classical model [13] in which all layers of atoms separated by about 2.06 Å in Y-direction are treated as masses connected by springs. This approach was found useful in case of explaining, for instance, Raman spectra in a few layers graphene [28][29][30]. A solution to this mathematical problem proposed by Tan et al [29] may be used as a qualitative description of the observed ω(N) dependence. In [29] they write angular frequencies of the possible set of solutions for N layers in the following form, for the ith vibrational mode: where i = 2 . . . N, and ω 0 is the angular frequency of the classical harmonic oscillator.
There is a significant difference however between their case and ours. In [29] the highest possible frequency mode is used (when i = N), as the one active in Raman spectroscopy, while we find the lowest possible frequency (when i = 2) as vibrational mode suitable here. Therefore we may write: ω = ω 0 · 1 − cos(π/N) 1/2 , which may be written as ω = 2 1/2 ω 0 · sin π/2N . As shown by solid lines in the upper part of figure 4, the angular frequency of damped oscillations are well fitted by the above relation when value of ω 0 of about 27/ps is used for Bonny potential and [111] direction.
For estimating the value of ω 0 in a model of harmonic oscillator we need to know k parameter in interatomic potential dependence on distance between atoms, when an atom embedded in a crystal lattice moving as a whole probes the potential well at various positions: U = kx 2 /2, where x is distance from the potential minimum. In order to find out this, a LAMMPS script was created with a frozen (not moving) fcc lattice and a set of 10 4 random potential values for iron atom were computed (for positions within a distance of 0.2 Å from potential minimum). It was found that k is about 8 eV Å −2 (potential is in fact not exactly harmonic and it is not even symmetric with respect to zero).
With Fe atomic mass of 55.845 we obtain ω 0 = k/m = 37.0/ps or in terms of frequency and energy it is 5.9 THz, or it is 24 meV. That value is in a reasonable agreement with computed [31] or experimentally determined phonon energies in Fe [32,33] or steel [34,35].
Detailed derivation and analysis of equations on ω(N), ζ(N) and A(N) will be published separately.

Dynamics of internal virial stress
Pressure P {xy} discussed so far is a physical quantity which is observed at the surface of sample and it is value is a result of summation of contributions to stress within entire sample interior. Dynamics of dislocations depends however mostly on stress distribution at dislocation core and around it. Therefore what we need to know for studying it is the internal stress distribution, which additionally changes in time. Introducing dislocations into sample volume changes distribution of stress there and changes also propagation dynamics of externally applied pressure. We cannot assume that contributions to local stresses is a simple sum of stress field caused by dislocation and a one due to penetrating the material interior external wave of pressure, since the stress at the core of dislocation is beyond the linearity limit. At the same time stress field changes due to dislocation movement.

Z Kozioł
Throughout this article we analyze the virial contribution from interatomic potential as given by the second term in equation (1), since in case of EAM potential of Bonny et al [10] no other terms there are present, while the kinetic energy term contribution to stress is of the order of ∼10 −3 only of that of potential energy, at temperature T = 50 K, and therefore it can be neglected.
The LAMMPS implementation of virial stress by Thompson et al [36], of a microscopic quantity that may be interpreted as a local, atomistic stress tensor, equation (1), passes the basic test of validity, in our case for steel 310S, which is the condition that the sums of all diagonal virial stress tensor components for all atoms in sample volume and at any time must be equal to total pressure with minus sign: We find that a similar relation is observed for the surface pressure shear component P xy as well, i.e.: In equation (9), S {xy} is value of virial stress tensor component xy averaged over all sample atoms.
One needs to take care of the conversion between S {i j} and P {i j} to avoid confusion: virial stress, as reported by LAMMPS, must be divided by an average volume of atoms (in units of Å 3 , when metal units are used), and hence it is material specific, while in LAMMPS it is reported in units of pressure. In our case, for steel 310S that conversion factor is 11.21, i.e. S {xy} (GPa) = −11.21 · P {xy} (GPa).
The relation expressed by equation (9) is fulfilled with high accuracy during all dynamic simulations and it can be used as a test of validity of any computation: any deviations from it must be treated as an almost certain signature of computational instabilities or methodological errors in the modeling process itself. Figure 5 shows typical results on virial stress distribution S {xy} in a 16 + 16 sample containing a dislocation line in its geometrical center in Y-direction.
The profiles of S {xy} averaged over X-direction (which is direction of dislocation movement, i.e. it is perpendicular to dislocation alignment along Z-axis) are shown in figure equation (5) as the pressure penetrates successively inside of the sample interior. Values of S {xy} saturate for times larger than about 30 ps for values of Y-coordinate that are below the dislocation line position (which in this case is at around Y = 100 Å). A large peak of S {xy} near the surface initially increases with time and starts to decrease slowly at around the same time when saturation of S {xy} is observed below the dislocation line position.
Results of figure 5 clearly indicate the existence of a profile of pressure spreading into the interior of sample. Therefore we performed simulations for samples of various sizes in Y-direction in order to determine accurately the speed of pressure wave penetration. Similar analysis must always be done on a large number of data, as fluctuations of internal stresses are large: probability Distribution function for S {xy} , as computed by us for all atoms, resembles a Gaussian function and at T = 50 K and at zero applied pressure has a half-width of about 3.0 GPa.
Results of figure 6 indicate that the penetration depth, as marked approximately by arrows, is linear with time and the penetration speed is about 75 Å ps −1 , which is close to the speed of sound waves, as estimated by us.
It is interesting to see how the stress penetrates in case of samples without dislocations. This is shown in figure 7(a), where profiles of S {xy} across 8 + 8 sample with no dislocations,  at time 50 ps after the external pressure P 0 xy has been applied to the sample surface located at Y of around 100 Å. At certain range of applied pressure (2.0 to 2.5 GPa in that figure) well observed spatial oscillations of S {xy} are found. At higher P 0 xy values the stress profile become nearly linear with y-coordinate.
These results suggest the existence of spatial oscillation in sample volumes at certain values of applied pressure. Therefore we performed a search for time oscillations as well. After performing computation on a sample of 8 + 8 with a time resolution of 0.005 ps, as shown in figure 7(b), we observed with a good accuracy the wave propagation at speed of 72 Å ps −1 . It was checked that this speed does not depend on sample size in Y-direction. A few waves that coincide are plotted in figure 7(b) with thicker lines. The period of oscillations and distance between wave maxima is well determined, with angular frequency ω = (2π/Δt)/ps = 29.2/ps. Hence, the observed oscillation frequency is in a good agreement with that one expected, based on classical oscillator equation, and in agreement with angular frequencies found out from measured ω(N), as shown in figure 4.

Displacement as a function of time
Once we have gained the basic understanding of the dynamics of S {xy} field penetration, we may start analysis of dislocation displacements.
Dislocations have been detected and analyzed by using Ovito [26] and Crystal Analysis tool [27]. These tools allow to extract the data about position of dislocation nodes for later analysis by Ovitos or by our own scripts in Perl. In most cases dislocations are detected only after pressure is applied. Figure 8 shows typical frames of visualization performed with Ovito, when pressure of 5 GPa is applied at t = 0. At that time no dislocations are found. Only after very short time of exposure of the sample to pressure ( figure 8(a), at t = 0.1 ps) fragments of dislocation lines are formed in entire sample volume. At t = 0.3 ps, figure 8(b), two very irregular lines are formed at positions as expected, in geometrical center of sample in Y-direction. At t = 0.4 ps, as in figure 8(c), a well developed pair of line dislocations is formed. Their shape changes all the time during simulations, regardless of pressure values, but no particular dependence of that shape on time or applied pressure was noticed. At some time intervals additional fragments of line or loop dislocations are found, mostly near the top Y surface where pressure is applied.
As the position of dislocations in X-direction we assumed an average of coordinates of their nodes, computed in a semi-automatic way by using Perl scripts and a set of CA data files [27], typically around a 1000 of them for every run as a function of time. Figure 9(a) shows typical dependencies of dislocation displacement as a function of time, after pressure P 0 xy is applied, for pressure values as indicated there. For the most bottom data points, at pressure of 0.6 GPa, dislocations move at some random moments only and it is ambiguous determining its velocity. The data points at 0.8 G Pa move almost all the time. At the highest speeds (for 2.4 GPa in this case) displacements are monotonic in time and form nearly ideal straight lines. At higher pressures however there is an abrupt change to curves with a strong curvature, and again there is some ambiguity in determining the speed of dislocation movement (curve for 5.0 GPa). We used for that an approximate slope of curves at t = 60 ps. At yet higher pressures again an irregular movement is found and dislocations gradually dissociate.
Hence, we observe that there exists a maximum of dislocation velocity at pressures around 2.4 GPa in case of the sample 8 + 8 of figure 9(a). A question arises if sample geometry has any influence on computational results. To check that modeling was performed on a series of samples named 8 + X, where X = n · 8 and n ranges from 1 to 7, as shown in figure 9(b). In these samples dislocation core is located eight lengths of unit cell from sample bottom while pressure is applied at the top surface which is X = n · 8 lengths of unit cell from dislocation core.
We find that distance of surface from the dislocation line not only influence dislocation velocity but also the time when it starts to move. Figure 10 shows S {xy} profiles for the same samples as in figure 9(b), at time moments indicated there by green arrows, when dislocations start to move. Each of these profiles intersects with a straight horizontal line at S {xy} = 0.75 GPa, indicating that this is the critical value of stress needed for dislocation movement. Moreover, position of these points on Y-axis is the same as position of the core of dislocations.  Figure 11(a) shows a collection of velocity data points from computations performed on samples of different sizes and pressures P 0 xy as well various distances between the surface where pressure is applied and dislocations positions. Clearly three ranges of applied P 0 xy can be distinguished (these ranges change with sample geometry).

Velocity scaling with pressure and sample sizes
I region: at lowest pressures, below around 1 GPa, velocity is small (and determined with large uncertainty) due to non-monotonic dislocation movement, marked with periods of no movement and abrupt jumps in speed. The value of 1 GPa coincides well with the value of internal stress of around 0.75 GPa, needed for displacing dislocations, as deduced from figure 10.
II region is between around 1 GPa and 2.4 GPa or 2.2 GPa for samples 8 + 8 and 16 + 16, respectively. Dislocation movement in this region is monotonic in time, and displacement is characterized by a dependence well approximated by a linear function of time.
The III region of pressures starts at the point where an abrupt drop in dislocation velocity is observed. For instance, for sample 2 + 2 that drop is from the maximum velocity to nearly zero within pressure range of 0.1 GPa. Initially, for larger samples, at pressures just after that drop, displacement as a function of time changes abruptly from linear to concave, with slope decreasing with time, and at certain moment dislocations start to dissociate.
Straight lines in figure 11(a) illustrate an attempt to approximate velocity in region II by a linear dependence on pressure, V = a · P 0 xy + b. It is found (not shown here) that coefficient a in this function scales well with the distance d between the surface at which pressure is applied and dislocations: a = 26.9 − 3.6 √ d where units of V are in Å ps −1 while distance d is the number of unit cells in Y-direction.
We find that a better scaling of velocity in regions I and II is obtained with S {xy} instead, as shown in figure 11(b). For that an average stress S {xy} in the sample space region between the sample bottom and location of dislocations is used. Profiles of S {xy} in that region become quasi-stable for all studied samples for times greater than about 30 ps since the time when pressure P 0 xy is applied. For a better accuracy the volume average S {xy} has been computed in that entire region for times in interval 40-60 ps.

Velocity scaling with temperature
For a more complete description, in this section we report on temperature scaling of dislocation velocity by a sample of size 8 + 8. In this case we obtain results very close to these already presented for T = 50 K. The difference is that dissociation of dislocations is seen only at these applied pressures when dislocation velocity becomes close to the speed of sound (the region III is narrow). In figure 12(a) velocity is shown as a function of applied pressure for several values of pressure and for a broad range of temperatures. We see that a critical stress needed for dislocation movement, τ p of about 0.75 GPa does not depend on temperature in a noticeable way. This value is in accordance with estimations for line dislocations in fcc lattice, as it is of the order of 10 −3 -10 −2 of shear modulus [37]. An exact value of τ p is however difficult for computation and a broad range of theoretical predictions is known [38,39]. In the classical Peierls-Nabarro model of stress flow [40], it is given by: where G is shear modulus, ν is Poisson ratio, b is magnitude of Burgers vector and d is distance between glide planes. For Burgers vector 1 6 [211], ν = 0.36, and d = 2.06 Å (which is a/ √ 3), we obtain τ p /G ≈ 0.03. Dislocation velocity V(T ) is determined by the speed the system can dissipate energy with phonon mediation. This can be seen well in MD simulations, which are of the realm of classical physics: reducing temperature to zero leads to no dislocations movement. In phenomenological models V(T ) is written in the form [41,42]: where τ is stress, b is value of Burgers vector, and B(τ , T) is approximated by B 0 + B 1 T. Curves of velocities of dislocations in figure 12(a) as a function of pressure remain similar, change of their slope is small (if any) and it is of about 8Å (GPa ps) −1 , regardless of temperature, albeit these curves shift upward when T increases. This is different from what is expected based on commonly accepted models [43,44], and different than for instance MD results observed for bcc Fe [16,45]. It is possible that large stress fluctuations for atoms in crystal lattice, due to atomic inhomogeneity, hinder in our case the effect of thermal fluctuations on dislocations mobility.
The figure 12(b) shows the same data as on figure 12(a), with the difference that velocities and temperature are in log-log scales. The solid straight lines illustrate a fitting of velocity to a power-law temperature dependence, for a constant value of external pressure applied. We find that velocities at constant pressure may be approximated by dependencies of the kind: V(T) ∼ T α , where α is around 0.13 for applied pressure of 2 GPa and it is around 0.09 for 6 GPa, in agreement with results reported in the literature.
Many features of dislocation movement as described here have been reported in the past, but they were not interpreted consistently in terms of dynamic internal shear pressure propagation.
For instance, Rodary et al [44] propose four-region approximation for velocity vs pressure diagram. We however do not observe a well defined saturation of velocity for large pressures. Instead dislocations dissociate when velocity is close to sound speed. Also, we do not find well defined velocities in region I. Instead we observe in some conditions that movement of dislocations occurs at certain moments only. In region II, which is usually defined as a region of monotonic movement with constant velocity we sometime observe periodic oscillations of velocity with time. These effects are well explained by time-oscillations of pressure profile entering material volume. We find also that instead of linear V(P 0 xy ) relation in region II a somewhat better approximation of the data would be provided by a power-law dependence of the kind: V = V 0 (T) · (P 0 xy − τ p ) β , where τ p , a Peierls threshold stress needed for the movement is in our case around 0.75 GPa, and exponent β is around 0.7.
Despite theoretical efforts and simulation results, still no good understanding exists of physical mechanisms behind dislocations mobility and dissociation at velocities comparable to the speed of sound [46]. Frank [47] found that the elastic energy of a moving dislocation diverges at the transverse speed of sound, in an analogy to divergence of energy of masses approaching the speed of light in relativistic theory (see also [40]). That suggested the existence of a limiting speed for dislocations. There is however evidence based on MD simulations that in some cases dislocations may achieve supersonic speeds [46], but no experimental observations confirmed that so far. It is thought that dislocation dissociation and supersonic speeds are accompanied by potential energy instabilities during the travel of pressure waves [48]. We observed energy (and temperature) instabilities in analysis of high-speed dislocations and treat these as a possible aim of the future studies.

Summary
In order to model in MD line dislocations evolution under the local stress field in fcc steel 310S, an analysis of pressure penetration into the sample interior has been carried out. It was shown that pressure (stress field distribution) enters the sample volume with sound speed (of about 70 Å ps −1 in our case) in a form of waves that can be described well by the damped harmonic oscillator equations, while oscillations frequency of atoms displacements are related to interatomic potential. Angular frequency of these oscillations, of about 29.2/ps, is in agreement with that determined for classical linear harmonic oscillator (for an atom of Fe in a fcc lattice), which we estimate at about 37.0/ps. Hence, stress distribution depends on time since the pressure has been applied at the surface and therefore it depends on size of samples. That results on sensitivity of any MD modeling on time scales involved (time step used) and on size of samples studied. Large oscillations of shear surface pressure are shown to be related to intrinsic material properties and sample size (number of atomic crystallographic planes) along the direction of sound propagation.
Velocity of dislocations scales well with spatially averaged local shear component of virial stress tensor value around dislocation core and not with the external shear pressure.
The maximum dislocation velocity approaches sound speed for high enough stress field. The existence of minimal value of stress (Peierls stress) of about 0.75 GPa is necessary for causing line dislocations movement. In most cases dislocations move with constant speed after the critical stress is exceeded and that speed is linearly proportional to applied pressure. It is confirmed that dislocation speed dependence on temperature may be described by a power-law formula.