Classical Mobility of Highly Mobile Crystal Defects

Highly mobile crystal defects such as crowdions and prismatic dislocation loops exhibit an anomalous temperature independent mobility unexplained by phonon scattering analysis. Using a projection operator, without recourse to elasticity, we derive analytic expressions for the mobility of highly mobile defects and dislocations which may be efficiently evaluated in molecular dynamics simulation. The theory explains how a temperature-independent mobility arises because defect motion is not an eigenmode of the Hessian, an implicit assumption in all previous treatments.

Plastic deformation of crystals is effected by the motion of dislocations and point defects [1].Away from shock loading and the melting temperature this motion is usually modeled with the viscous damping law _ x ¼ γ −1 • f, employing a matrix of friction or drag coefficients γ, which set the time scales of defect dynamics [2].To reproduce the stochastic trajectories of highly mobile defects seen in experiment [3,4] this mobility law has been supplemented with a stochastic force to give the Langevin equation [5] [6,7].The stochastic force is usually more significant for small dislocation loops and point defects because the configurational force f λ is determined only by gradients in the stress field.For larger extended defects the configurational force usually dominates over the stochastic force.In both cases γ controls the rate of important microstructural processes such as swelling and post-irradiation annealing [8], but no universal theory for γ exists.
In this Letter we use the Zwanzig projection technique [18] to show that γ ¼ γ 0 þ k B Tγ w , in quantitative agreement with MD simulations of defects and dislocations.γ 0 arises because the defect displacement vector is not an eigenvector of the Hessian, so that thermal vibrations can induce a force on defects to linear order.This is missed in previous treatments [11,19] as by perturbing a quadratic integrable Hamiltonian one implicitly assumes that defect motion is an eigenmode, an assumption that we explicitly show to be false.Violation of this assumption is the origin of the anomalous mobility.
Defect coordinates.-Wedescribe a crystal using a 3Ndimensional vector of atomic positions X ∈ R 3N and velocities _ X ∈ R 3N .In this treatment crystal defects are not elastic singularities but localized deformations, which may be assigned a set of M ≪ N "position" labels x λ ∈ R 3M and "velocity" labels _ x λ ∈ R 3M to characterize the state of a defective crystal.Common methods for determining x λ , _ x λ include analysis of the atomic disregistry [20] or an energy filter [7], though in the following the only requirement is a repeatable protocol.By definition, the zero temperature configurations X ¼ Uðx λ Þ of the crystal potential energy VðXÞ may be entirely characterized by the parameters x λ , while variation of Uðx λ Þ with x λ can be determined through nudged elastic band calculations [21] or simply a finite difference derivative in the case of a defect with a wide core.To complete the discrete representation of a crystal at finite temperature, we must include displacements due to thermal vibrations Φ ∈ R 3N .The crystal configuration X at any given instant can now be expressed as . By introducing a defect position and velocity the coordinate set To rectify this we require the vibrational displacements Φ to be independent to the displacements caused by defect motion ∂ λ U, giving the 6M constraints [12] Published by the American Physical Society under the terms of the Creative Commons Attribution 3.0 License.Further distribution of this work must maintain attribution to the author(s) and the published article's title, journal citation, and DOI.
To obtain a dynamical equation for x λ , it now suffices [22] to project the exact equation of motion m ̈X ¼ −∇VðXÞ onto the direction ∂ λ U orthogonal to the crystal vibrations.Defining an effective mass tensor m ¼ m∂ λ U • ð∂ λ UÞ T , we exploit the time invariance of (2) to obtain . Similar equations of motion are standard in dynamical quasiparticle theories [12,22] and, in common with other authors, we will neglect the "hydrodynamic" term − _ . This is justified as we consider the motion of only subsonic defects, and it can be shown that these terms are of order j _ x λ j=c ≪ 1, where c is the speed of sound.As a result, the defect coordinates evolve according to where we have defined the instantaneous defect force f λ as the projection of the total force −∇V in the direction of defect motion ∂ λ U.The vibrational coordinates evolve in the subspace orthogonal to Removing the vibrational coordinates.-Fromthe form of the potential energy V½Uðx λ Þ þ Φ, it is clear that the evolution of the defect and vibrational coordinates are coupled, as they must be for a frictional force to exist.However, for highly mobile subsonic defects, which necessarily possess a wide defect core [23], the defect coordinates may be considered as slowly varying compared to the vibrational coordinates, a conclusion which will be explicitly demonstrated in molecular dynamics simulation below.Over a Debye period τ D ∼ a=c ∼ 0.1 ps, where a is the lattice parameter, the displacements of any atom due to thermal vibrations will approximately average to zero, with an oscillation amplitude of ∼τ D ffiffiffiffiffiffiffiffiffiffiffiffiffiffi k B T=m p .Since the defect speed will be approximately _ ≪ c, the displacement of any one atom due to defect motion in a time interval τ D will be at most , where jj∂ λ Ujj ∞ is the largest component of ∂ λ U.These calculations imply that if jj∂ λ Ujj ∞ ≪ j∂ λ Uj, then the displacement due to defect motion will be much less than the magnitude of displacements due to thermal motions, which implies that the Φ are effectively ergodic [24] over a time scale ∼τ D , where the defect coordinates are essentially stationary.But the condition ∥∂ λ U∥ ∞ ≪ j∂ λ Uj amounts to a requirement that the deformation associated with the defect is spread over many atomic sites, which is always satisfied by highly mobile defects with a wide core.We therefore assume that vibrational displacements average to zero over periods of ∼0.1 ps while the defect remains effectively stationary, an assumption that we will test explicitly when calculating the defect force autocorrelation.
We can exploit this separation of time scales to remove thermal vibrations from the defect equation of motion using the formalism of dimensional reduction by Zwanzig [18,25].In this formalism the solution of the "fast" equation of motion for Φ is substituted into the "slow" equation of motion for x λ .It may be shown, to order τ 3 D , that Φ, _ Φ are adiabatic with respect to x λ , _ x λ and ergodic over the partial Gibbs distribution where Zðx λ Þ ¼ exp½−βFðx λ Þ is the partial partition function and we integrate on the hypersurface defined by (2).
The defect coordinates now evolve on a coarse time scale τ D and follow the stochastic equation of motion It is usual in dislocation dynamics to neglect the inertial term m • ẍλ ðtÞ, which is valid when the potential energy landscape is slowly varying over the thermal length ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi k B T= mjγj p [5].In (5) we have introduced the expected force hf λ i ¼ −h∂ λ Vi ¼ −∂ λ F, the stochastic force ηðtÞ, where hηðtÞ ⊗ ηðt 0 Þi ¼ 2k B Tγδðt − t 0 Þ, and our central quantity, the friction matrix γ.In this timescale separated regime it is a standard result that γ is proportional to the time integral of the force autocorrelation CðτÞ, namely, where when x λ ∈ R, CðτÞ ≡ hf λ ðτÞf λ ð0Þi − hf λ ð0Þi 2 and may be expressed by ergodicity as We evaluate CðτÞ, and hence γ, in two ways: first by deriving in closed form the thermal averages (4) and second by numerical calculation of f λ ðtÞ in MD simulation.Analytic derivation.-Toderive an expression for γ we expand the potential energy V and the defect force f λ in powers of Φ.For the evaluation of the partition function the constraints (2) and the requirement that the Uðx λ Þ describes the zero temperature configurations results in an expansion where all inner products are with respect to Φ and all partial derivatives are evaluated at there is no restriction on the existence of mixed derivatives ∂ λ ∇ n Φ V ≠ 0. This is important as these mixed derivatives couple the defect and vibrational coordinates, as can be seen in the defect force expansion While we retain anharmonicity in the defect force, in order to perform analytical evaluation of expectation values we truncate V to quadratic order in Φ in the Gibbs distribution (4), allowing us to explicitly evaluate the expectation values in terms of the 3ðN − MÞ dimensional vibrational eigenset fω l ; v l g, where ∇ 2 Φ • v l ¼ mω 2 l v l .This truncation neglects any thermal expansion arising from the purely vibrational anharmonicities ∇ 3 Φ V and ∇ 4 Φ V.In the Supplemental Material [26] we systematically include these terms to produce an expression for γ up to linear order in temperature.It is shown that the anomalous temperature independent mobility γ 0 is unaffected by these additional terms.Using a quadratic Gibbs distribution, the expected force is found to be hf λ i ¼ −∂ λ ðV − TSÞ, where S is the harmonic entropy k B P l log ω l [27]; to evaluate CðτÞ we evolve the vibrational coordinates Φ from a given x λ .This is justified by the time scale separation and achieved by evaluating propagator terms of the form As appropriate for nonconservative dynamics, the propagator is evaluated using only the initial conditions hΦð0Þ ⊗ Φð0Þi ¼ P l k B T=mω 2 l v l ⊗ v l and, consequently, is closely related to the retarded Green's function GðtÞ ¼ ΘðtÞðk B TÞ −1 hΦðtÞ ⊗ Φð0Þi [28].All that now remains is to perform elementary Gaussian integrations to obtain our main result We see that the friction coefficient takes the form γ ¼ γ 0 þ k B Tγ w , with the new temperature independent γ 0 a function of the mixed quadratic derivative ∂ λ ∇ Φ V, and the temperature dependent k B Tγ w a function of the mixed cubic and quartic derivatives ∂ λ ∇ 2 Φ V and ∂ λ ∇ 3 Φ V.These terms may, in principle, be evaluated after diagonalizing ∇ 2 Φ V to obtain fω l ; v l g and computing the tensorial derivatives ∂ λ ∇ n Φ V.However, in common with modern methods to evaluate dispersion relations [29], we have found dynamical measurement of the thermal averages to be much more efficient than complete diagonalization of the vibrational Hessian ∇ 2 Φ V.
Numerical evaluation.-Wehave developed a method to calculate f λ ðtÞ by MD simulation, which yields CðτÞ and hence γ, yielding a numerical evaluation of the analytic expressions (11).In an ensemble of MD runs, with no stress applied, we time average the output for each run XðtÞ using a coarse-grained time step between τ D =4 and τ D to give hXi.To eliminate any errors, we find the zero temperature configuration U λ which minimizes j∂ λ hXi − ∂ λ Uj 2 .The calculated ∂ λ U is then used to project out the defect force f λ ðtÞ ¼ −∂ λ U • ∇V½XðtÞ over the same averaging time interval that produced hXi.We have found this method to be robust to variation in the averaging period and especially efficient for short line segments or nanoscale defects, where the zero temperature structures are typically related by rigid translation [30].An example of such calculations is shown in Fig. 1 for a 7 atom SIA cluster in tungsten, which exhibits the anomalous temperature independent mobility γ ¼ γ 0 [17], and in Fig. 2 for a highly mobile edge dislocation in iron, which exhibits a mixed temperature dependence γ ¼ γ 0 þ k B Tγ w [15].In both cases we see that CðτÞ loses all coherence after the first zero at ∼τ D =2, over which time the defect is observed to be essentially stationary.This validates our assumption of time scale separation between thermal vibrations and defect motion.We identify the subsequent force autocorrelation (FAC) signal as noise because it flattens with the system and ensemble size, limiting the integration CðτÞ only to the first zero.As shown in the figures, this method gives values in excellent agreement with conventional trajectory analysis.We also calculated the FAC for the 7-atom SIA  11) for a 7 atom SIA cluster in tungsten using LAMMPS [31] and an interatomic potential by Marinica et al. [32].We see a very similar peak in all methods which loses coherence after a time period ∼τ D =2, and we approximate the time integral in (11)  Φ V. We find excellent agreement with the dynamical method, as shown in Fig. 1.
Discussion.-Terms similar to (11) appear in phonon scattering predictions of γ, where they may be interpreted diagrammatically, with ∂ λ ∇ n Φ V approximately representing a vertex of one defect with n phonons [11,34].In this continuum picture, defects and phonons are separable to harmonic order, conserving energy and momentum in collisions.As a result, each term in (11) becomes dependent on the phase space available for the scattering process it represents.The anomalous term γ 0 is forbidden in such models as it represents the pure absorption or emission of a phonon, a process which has a vanishing phase space for subsonic defect speeds due to the linear phonon dispersion relation [34,35].It turns out that the second term in (11) dominates, describing a two-phonon elastic scattering process known as the phonon wind.With a cubic anharmonicity parameter A [36] this term has an approximate magnitude ∼k B TðA=μÞ 2 τ D , where μ is the shear modulus, in agreement with more detailed continuum treatments [11].However, the prediction γ 0 ¼ 0 from continuum analysis does not explain the observed simulation results.
To see how the present treatment allows an anomalous temperature independent mobility, we express γ 0 in the eigenbasis fv k g of the vibrational Hessian ∇ 2 Φ V. Using (10) and the expansion For this term to vanish, as in all continuum theories, we require ∂ λ ∇ k V ¼ 0. But this implies that the defect displacement operator ∂ λ U is an eigenvector of the total Hessian ∇ 2 V as the "off-diagonal" terms ∂ λ ∇ k V that mix ∂ λ U and the vibrational modes must vanish.We have explicitly demonstrated that this is not the case; it is precisely this effect, which relies on the weaker identification of a defect as a localized deformation that is not an eigenvector of the Hessian, in contrast to a canonical quantity separable from vibrations, that gives rise to γ 0 .Of course, anharmonic vibrations still affect the dynamics in a manner which becomes analogous to typical scattering theories in a continuum picture, giving the phonon wind term k B Tγ w in (11).These terms are appreciable for only extended line dislocations, which significantly deform the host lattice, while the anomalous γ 0 is the leading term for nanoscale defects which are typically elastically neutral in the far field.For some extended dislocations in close-packed crystals the defect translational operator is very nearly an eigenvector of the Hessian, implying that the anomalous mobility vanishes and γ ∼ k B Tγ w [13].But, in general, we have found this not to be the case, with the mixed dependence γ ¼ γ 0 þ k B Tγ w occurring across a wide range of crystal defects.
Concluding remarks.-Ourmain result is an explicit form (11) for the friction tensor γ of highly mobile crystal defects.We believe this is a new result.It may be used to parametrize accurately deterministic (_ ηðtÞg defect mobility laws.The result was obtained by identifying defects through a projection operator with no recourse to elasticity.An anomalous temperature independent mobility γ ∼ γ 0 arises because the displacement vector corresponding to defect motion is not an eigenvector of the Hessian, in violation of elasticity theory or solitonlike models, where vibrations are canonical.This finding highlights the importance of intrinsically discrete (i.e., atomistic) analysis to understand nanoscale crystal plasticity.We note that the form of γ 0 in (11) is closely analogous to the famous Kac-Zwanzig heat bath formula [18].But rather than a random array of harmonic oscillators with a constant coupling strength, we have here the vibrational modes of the entire crystal coupling to a localized deformation through ∂ λ ∇ Φ V.It is hoped that our explicit expression for γ and the method of evaluation may be used to provide further connections between analytic heat bath models and the thermal dynamics of real systems.

PRL 113 ,
215501 (2014) P H Y S I CA L R E V I E W L E T T E

FIG. 1 (
FIG.1(color online).Evaluation of the defect FAC in unbiased molecular dynamics simulation at three temperatures and the first analytic term in (11) for a 7 atom SIA cluster in tungsten using LAMMPS[31] and an interatomic potential by Marinica et al.[32].We see a very similar peak in all methods which loses coherence after a time period ∼τ D =2, and we approximate the time integral in(11) by the area under this first peak.Inset: Comparison of the predicted diffusivity D ¼ k B T=γ and the direct measurement D ¼ hx 2 i=2t.
FIG.1(color online).Evaluation of the defect FAC in unbiased molecular dynamics simulation at three temperatures and the first analytic term in (11) for a 7 atom SIA cluster in tungsten using LAMMPS[31] and an interatomic potential by Marinica et al.[32].We see a very similar peak in all methods which loses coherence after a time period ∼τ D =2, and we approximate the time integral in(11) by the area under this first peak.Inset: Comparison of the predicted diffusivity D ¼ k B T=γ and the direct measurement D ¼ hx 2 i=2t.
T. D. S. was supported through a studentship in the Centre for Doctoral Training on Theory and Simulation of Materials at Imperial College London, funded by EPSRC under Grant No. EP/G036888/1.This work was partially funded by the RCUK Energy Programme (Grant No. EP/I501045) and by the European Unions Horizon 2020 research and innovation programme under Grant Agreement No. 633053.To obtain further information on the data and models underlying this Letter please contact PublicationsManager@ccfe.ac.uk.The views and opinions expressed herein do not necessarily reflect those of the European Commission.This work was also partially funded by the United Kingdom Engineering and Physical Sciences Research Council via a programme Grant No. EP/ G050031.

FIG. 2 (
FIG. 2 (color online).Evaluation of CðτÞ for a 1=2h111ið10 1Þedge dislocation in Fe, using an interatomic potential by Gordon et al.[33], normalized to the unit length aj½1 21j ∼ 7 Å.The FAC increases with temperature such that γ ¼ γ 0 þ k B Tγ w , exhibiting both anomalous and phonon wind drag.Inset: Comparison with direct measurement of the diffusivity.The values are in quantitative agreement with finite stress simulations[15].