Local structural excitations in model glasses
LLocal structural excitations in model glasses
S. Swayamjyoti and P. M. Derlet ∗ Condensed Matter Theory Group, Paul Scherrer Institut, CH-5232 Villigen PSI, Switzerland
J. F. L¨offler
Laboratory of Metal Physics and Technology,Department of Materials, ETH Zurich, 8093 Zurich, Switzerland (Dated: November 8, 2018)
Abstract
Structural excitations of model Lennard-Jones glass systems are investigated using theActivation-Relaxation-Technique (ART), which explores the potential energy landscape of a localminimum energy configuration by converging to a nearby saddle-point configuration. PerformingART results in a distribution of barrier energies that is single-peaked for well relaxed samples. Thepresent work characterises such atomic scale excitations in terms of their local structure and envi-ronment. It is found that, at zero applied stress, many of the identified events consist of chain-likeexcitations that can either be extended or ring-like in their geometry. The location and activationenergy of these saddle-point structures are found to correlate with the type of atom involved, andwith spatial regions that have low shear moduli and are close to the excess free volume within theconfiguration. Such correlations are however weak and more generally the identified local struc-tural excitations are seen to exist throughout the model glass sample. The work concludes with adiscussion within the framework of α and β relaxation processes that are known to occur in theunder-cooled liquid regime. a r X i v : . [ c ond - m a t . m t r l - s c i ] N ov . INTRODUCTION Local structural excitations (LSEs) occurring at the atomic scale of Bulk Metallic Glasses(BMGs) are believed to mediate their unique plastic deformation properties. In the earlyseminal work of Spaepen , a free volume model was proposed in which single atoms wouldmigrate to a neighbouring free volume region. Subsequently Argon introduced the idea thatsmall groups of atoms undergo a structural transformation that could modify the local shear-stress state. The early atomistic simulations of Falk and Langer confirmed this picture andreferred to these LSEs as Shear Transformation Zones (STZs), leading to a series of staticand dynamical atomistic simulations to better understand how such STZs lead collectivelyto macroscopic plasticity . The structural transformations studied in these works wereobtained via application of strain rates that are many orders of magnitude higher than thatnormally seen in experiment, leading to localized material instabilities often identified asSTZs. In terms of the traditional thermal activation picture for plasticity, such simulationstherefore study the athermal limit of plasticity in bulk metallic glasses.A common question has been as to whether there exist local structural features withinthe glass configuration that facilitate the existence of LSEs. Indeed from the early athermalatomistic simulation work, those spatial regions which exhibited a localised material insta-bility and the subsequent observation of STZs have been referred to as pre-existing liquidlike regions - spatial domains that are not fully relaxed . In the work of Langer, a rela-tively low density of such regions was theoretically considered . From this perspective it isthen valid to ask whether there are structural features of the unstrained glass configurationthat allow identification of regions predisposed to LSEs - the so-called liquid like regions? Inthe work of L´eonforte et al , the reversible non-affine atomic displacements associated witha finite elastic distortion was considered. Such displacements give insight into the spatialextent of the low curvature of the unstrained PEL. This work found a percolative networkof large non-affine displacements, however when the applied strain was increased to promoteplastic activity, the location of the observed LSEs did not show a strong correlation withthe non-affine displacements suggesting the corresponding athermal energy barriers werenot directly related to the unstrained PEL curvature. In the work of Mayr , who analyzedthe local elastic constants (Born plus the kinetic fluctuation contributions) as a function ofapplied stress and temperature did find a strong connection between emerging elastic insta-2ilities and eventual plasticity as the temperature approached the glass transition. Whilstthese simulations were fully dynamical, the necessarily high strain rates place these resultsclose to the athermal regime.Atomistic simulation methods that avoid the high stress athermal limit, are the so-calledpotential energy landscape (PEL) exploration methods. The relevance of these methods tobulk metallic glasses relies on the idea that the atomic configuration of a structural glassspends much of its time in a local potential energy minimum, only occasionally transitingit via a saddle-point region to a new local minimum. Such structural changes are assumedto occur via thermal fluctuations and therefore do not involve a local material instabilitydriven by the applied stress. A well known example of this approach is the so-called Acti-vation Relaxation Technique nouveau (ART n ) which has been applied to both two andthree dimensional model glass systems . Starting from a well defined local minimum theART n method uses the local Hessian structure of the PEL to climb out of its current PELvalley and to eventually converge to a nearby saddle point region. The energy differencebetween the minimum and the saddle point gives the corresponding energy barrier of theassociated LSE. When applied to model glass systems a wide range of energy barriers areobtained producing a distribution of energy barriers that appear to converge when severalthousand such LSEs have been identified. For sufficiently relaxed glass configurations, thebarrier energy distribution is found to be single peaked and analysis of the associated LSEsrevealed their corresponding plastic strain increments to be uncorrelated with barrier en-ergy . Moreover, with respect to a particular loading geometry the plastic strains haveequal numbers of positive and negative sign. However, upon application of the loading ge-ometry strain, the distribution becomes slightly weighted towards positive strains with alarge number of negative strain LSEs still occurring. The situation is quite different for thecase of stress driven energy barriers which generally only yield LSEs that are positive instrain. Thus athermal stress driven simulations probe only a subset of possible LSEs —those compatible with the chosen driving mode.In this paper, the ART n method will be applied to three-dimensional model LJ structuralglass configurations to obtain a large number of LSEs and their associated barrier energies.The present work will only consider the case of zero external load. Using this data, the spatiallocation of the identified LSEs will be investigated in terms of local structural propertiesand a statistical analysis will be performed. The local quantities investigated will be atom3ype, local stress, local elastic constants, local pressure and Voronoi volume. In addition,the spatial nature of the structural changes and the number of atoms involved will bedetermined. All such information will also be correlated with the barrier energy. In sec. IIthe sample preparation procedure will be outlined and the local structural quantities to beinvestigated defined. Sec. III contains the major results of the present work encompassingthe ART n results and the ensuing statistical analysis of the identified LSEs. In sec. IV,the spatial nature of a number of representative LSEs will be atomistically visualized andcharacterised. Finally, sec. V will discuss the results in terms of contemporary pictures ofmicroscopic plastic deformation in BMGs. II. METHODOLOGYA. Sample preparation and ART n calculations Model glass samples have been prepared by molecular dynamics and statics using a 50/50binary mixture of a LJ system given by, V LJ ( r ) = 4 (cid:15) (cid:18)(cid:16) σ αβ r (cid:17) − (cid:16) σ αβ r (cid:17) (cid:19) , (1)where (cid:15) and σ set the microscopic energy and length scale of the model material. Theparameterization presently used is that of Wahnstr¨om parametrization , which is shown intable I. All simulation results are reported in LJ units where time is measured with respectto τ = (cid:112) mσ /ε and temperature with respect to ε/k b .Four samples, each with 1728 atoms, have been prepared using different quench-rates( η = 24 . / η = 24 . / η = 24 . / η = 24 . / × k b [ ε/k b ] and hydrostatic pressure of 8/160 [ ε/σ ],2) slow quenching of the sample from this well-equilibrated liquid state, which involves anincremental reduction in both temperature ( − . × k b [ ε/k b ]) and pressure (-0.158/160[ ε/σ ] ) by NPT molecular dynamics to form the disordered amorphous glass at a tem-perature of 100 × k b [ ε/k b ] and hydrostatic pressure of 0.1/160 [ ε/σ ] 3) Relaxation of theatomic coordinates to zero temperature and zero hydrostatic pressure by molecular staticsusing the Parrinello-Rahman method . For steps 1 and 2 the Parrinello-Rahman barostat4 arameter Parameter Description Value ε Depth of the potential well 1 σ Diameter of atom type 1 1 σ Cross-interaction diameter 11/12 σ Cross-interaction diameter 11/12 σ Diameter of atom type 2 5/6 r c Potential cut-off 2.5 σm Mass of atom type 1 2 m Mass of atom type 2 1TABLE I: Numerical values of the parameters of the LJ potential used in the present work , where σ ij represents the length-scale parameter between atoms of type i and j . has been used for pressure control and the Anderson-Hoover thermostat has been used fortemperature control.The ART n technique is then applied to the these samples. The ART n identification of anLSE involves randomly choosing one atom and displacing it by a small distance ( < . σ ) fromits equilibrium position — the starting condition for ART n . The lowest eigen-value of thecorresponding Hessian matrix is determined and the system is moved along the corresponding3 N dimensional eigen-vector until the eigenvalue of the Hessian becomes negative. Thispart of the algorithm is referred to as the Activation phase. At this point, the configurationhas passed an inflection region of the PEL and enters the new phase of Relaxation. Herethe configuration is moved in the direction of the eigenvector of the lowest eigen-valuewhich is now negative, with a new eigen-value and eigen-vector begin calculated at eachiteration. This is repeated until the dot-product of the total force of the configuration withthe eigen-vector is zero. When this occurs, the configuration has reached a saddle-pointwhich, following is referred to as the Activated state.To study the atomic scale environment of each identified LSE, an appropriate atomicweight for each atom is calculated via the displacement vectors between the initial state andeither the activated or final state configuration. Presently the normalised weights (which5um to unity) are defined to be w i = | ∆ R i | (cid:80) Ni =1 | ∆ R i | . (2)This form was motivated by the standard definition of participation number of aneigenvector , which gives information about the number of elements contributing to thenorm of the vector. Indeed, PN = 1 (cid:80) Ni =1 w i , (3)gives the effective number of atoms involved in the identified LSE.Given a Local Atomic Quantity (LAQ) for each atom, the weighted average of the quan-tity, according to eqn. 2, will give a representative value for the region occupied by a par-ticular LSE. That is, LAQ LSE = N (cid:88) i =1 w i × LAQ i . (4)With this quantity, an average over many LSEs can be made and compared to the un-weighted LAQ average where all atoms equally contribute to the average. Scatter diagramsfor the weighted LAQ with respect to barrier energy are also investigated to determine ifany correlation exists between local structural features and the activation energy of an LSE.To find out any possible linear correlation of such plots, the Pearson coefficient is used. For n data points, this is given byPC = (cid:80) ni =1 ( X i − ¯ X )( Y i − ¯ Y ) (cid:112)(cid:80) ni =1 ( X i − ¯ X ) (cid:112)(cid:80) ni =1 ( Y i − ¯ Y ) , (5)where X i and Y i are the data-sets, and ¯ X and ¯ Y are their respective arithmetic means.The value of PC ranges from − − linear correlation. B. Local Atomic Quantities
The local quantities (the LAQs) presently considered are • volume, calculated via an atomic scale Voronoi tessellation via the voro++ package • energy • pressure 6 dilation elastic modulus (three times the local bulk modulus) • the five linearly independent Kelvin eigen-shear elastic moduliSince a LJ pair potential is presently being used, the expressions for the local stress andelastic modulus are of a simple form, where respectively σ µνa = 12 V a (cid:88) ij V (cid:48) ( R ij ) R µij R νij R ij Λ a,ij (6)and C µναβa = 12 V a (cid:88) ij (cid:20) V (cid:48)(cid:48) ( R ij ) − V (cid:48) ( R ij ) R ij (cid:21) × R µij R νij R αij R βij R ij Λ a,ij + σ νβa δ µα + σ ναa δ µβ . (7)In the above expressions, Λ a,ij represents the proportion of the ij th bond within with thevolume of atom a . It is noted that bonds between two atoms, neither of which is atom a ,may also contribute to these two local quantities. Eqns. 6 and 7 properly partition volumeand therefore correctly take into account the contribution of each atomic bond . Thelocal pressure is obtained by taking one-third the trace of the local shear stress tensor. Toobtain the Kelvin elastic moduli from the fourth rank elastic stiffness tensor (eqn. 7), theusual Voigt elastic stiffness matrix is first constructed from which the Kelvin matrix (2ndrank tensor) is obtained via C µνK = A µν C µνV . For an explicit form of A µν see ref. . Thefive linearly independent eigen-shear moduli are obtained by first projecting out the puredilation distortions and then diagonalizing the resulting Kelvin matrix — for more details,see refs. .In contrast to the popular Voigt notation, the Kelvin notation preserves the norm of theactual elastic stiffness tensor and hence their eigen-values and eigen-vectors have geometricalsignificance . Since the Kelvin matrix is a tensor, the eigen-values of the stiffness matrixcan be computed to obtain the bulk modulus and the eigen-shear moduli. The invarianceof these eigen-shear moduli with respect to coordinate systems and thereby their role as anintrinsic material property has been highlighted in .7 . Natural Mode Analysis The natural modes of an N atom configuration can be obtained via the solution to (cid:88) jν (cid:0) m i [ ω n ] δ ij δ µν − ∆ µνij (cid:1) u νj,n = 0 (8)where ∆ µνij is the translationally invariant dynamical matrix obtained from∆ µνij = (cid:88) a,a (cid:54) = i H µνia δ ij − H µνij (1 − δ ij ) . (9)Here H µνij is the Hessian. In terms of the LJ interaction, V ( r ), the Hessian may be writtenas H µνij = (cid:20) V (cid:48)(cid:48) ( R ij ) − V (cid:48) ( R ij ) R ij (cid:21) R µij R νij R ij + V (cid:48) ( R ij ) R ij δ µν . (10)In eqn. 8, m i is the atomic mass of the i -th atom, and u νj,n is the eigen-vector associatedwith the eigen-frequency ω n of the n th natural mode. The number of atoms participatingin a particular eigenstate, u νj,n , may be obtained via the participation number PN n = (cid:34)(cid:88) i | (cid:126)u i,n | (cid:35) − , (11)where (cid:126)u i,n is the three-dimensional polarisation vector of atom i coming from the eigen-vector of eigen-frequency ω n . Assuming a normalised eigen-vector, PN n will range betweenunity (when the eigenstate is concentrated on just one atom) and N (when the eigenstateis distributed evenly over the entire sample). III. RESULTSA. ART n In total 1347 unique activated states where identified using the ART n method. To ver-ify that each activated state was directly connected to the initial atomic configuration, theactivated configuration was perturbed in a direction towards the initial atomic state config-uration and allowed to relax. If the resulting structure differed from the initial structure theLSE was discarded from the data-set, as was done in ref. . In a similar way, the final statecould be determined by perturbing the activated configuration in a direction away from the8 IG. 1: Red balls represent the centre-of-positions of all identified local structural excitations andthe green balls represent regions containing free volume within the simulation cell of the modelglass initial atomic state configuration and allowed to relax. In total 954 initial/activated/finalstate atomic configurations where obtained for subsequent analysis.To determine the spatial location of a particular LSE, the centre-of-position of thosesaddle-point atoms displaced relative to the initial configuration by more than 0 . σ , wascalculated. Fig. 1 displays these positions within a boundary box defining the three di-mensional periodic simulation cell. A general inspection of their spatial distribution revealssome heterogeneity and, in particular, regions where many LSEs are similarly located. Amore detailed inspection of such LSEs (as in a manner described in sec. IV) reveals themto be quite different in spatial extent and barrier energy, despite some LSEs having theircentre-of-position almost coincident. The goal of the proceeding section will be to see if anylocal structural feature correlates with this observed heterogeneity.Fig. 2a displays the distribution of barrier energies obtained from the 954 identifiedactivated states for the sample with the slowest quench rate. In agreement with refs. , thedistribution peaks at a non-zero barrier energy and appears to approach zero for small enough9arrier energies. The barrier energy scale is comparable to that seen in a previous fully threedimensional ART n simulations using the same LJ potential parametrisation . Inspectionof the participation number in fig. 2b reveals that the effective number of atoms involved inthe LSE is typically less than ten and that the LSEs with barrier energies in the lowest andhighest regime involve only a few atoms. Close inspection of the final state participationnumber (see inset of fig. 2b) reveals a clustering at integer values of the participation number.This does not occur for the activated state configuration demonstrating that whilst thenon-affine displacement field associated with the activated configuration is dispersed overmany atoms, the non-affine field of the final state can be localized on a discrete number ofatoms. Further inspection of all LSEs revealed that 81 initial/activated/final state atomicconfigurations had the feature that the final state configuration had an identical total energyto the initial state. Detailed inspection of these LSEs revealed that in these 81 cases, thefinal state involved a permutation of nearby atoms of the same type, where the activatedconfiguration involved a closed loop of displaced neighbouring atoms (see fig. 10a for anexample). Thus the final state is identical to the initial state and when not included infig. 2b the discrete participation numbers vanish. Since such LSEs cannot produce anystrain these were also removed from the data-set used in the LAQ analysis of proceedingsection.Fig. 2c now displays a histogram of the participation numbers indicating that the LSEsidentified by ART n generally involve somewhere between one and several atoms. Thus onaverage there is little difference in the number of participating atoms between a activatedand final state. On the other hand, fig. 2d displays the difference in participation numberbetween connected activated and final state configurations showing that the number caneither decrease or increase by several atoms. On average, however, the change in the partic-ipation number is close to zero indicating no strong bias to whether the final state containsmore or less participating atoms than the activated state, a result confirming fig. 2c. B. Statistical analysis of LAQs
Fig. 3 displays the normalised histograms of a number of LSE averaged LAQs with respectto the activated and final state configurations, also shown are the equivalent unweightedhistograms derived from the total sample. The two types of distributions will be referred to10
IG. 2: a) Normalised distribution of activation energies of LSEs in a 3D model glass system.b) Scatter plot of participation number (PN) of activated and final state configurations versusactivation energy where the inset is a blow up of the final state configuration participation numbers.c) Histogram of participation number and d) scatter plot of participation number difference betweenthe connected activated and final states as a function of activation energy. as LSE weighted and unweighted distributions of the LAQ. For the unweighted distributionsthe (atom-type resolved) partial histograms are also shown.Fig 3a represents the normalised distribution of potential energy of atoms. The doublepeak structure of the total unweighted distribution is clearly seen to arise from the singlepeaked distributions of each atomic type. The LSE weighted distributions, on the otherhand, do not exhibit a double peak structure, with the single observed peak coinciding withthe unweighted partial distribution of atoms of type 2. This result suggests that the atomsinvolved in an LSE are often of type 2. Fig. 3b now shows the corresponding normalised11
IG. 3: Normalized distribution of local a) energy b) volume c) pressure d) bulk modulus. Thered vertical lines represent the corresponding mean value derived from the total sample. distribution for the Voronoi volume. Inspection of the unweighted curves show that atomsof type 2 have lower volume compared to atoms of type 1, a feature that is expected giventhe nature of the Wahnstr¨om parametrization (see table I and ref. ). The distinct doublepeak feature is however absent for the LSE weighted distributions, both the activated andfinal state curves showing a single peaked structure approximately centered between theunweighted partial distributions. Closer inspection does however reveal a slight bias to thelower volumes of the type 2 atoms. Figs. 3c and d show respectively the distributions for bothlocal pressure and local bulk modulus. For both LAQs, the unweighted total distributionshow single peaked curves that are slightly asymmetric due to slightly different contributionsfrom the two atom types. That atoms of type 2 are involved in the LSEs is also reflectedhere because the LSE weighted curves are more biased toward the partial unweighted curves12f type 2 atoms.Fig 4a shows the distribution of Kelvin eigen-shear moduli for the activated and finalrelaxed states. Again, the unweighted total and partial distributions are shown for compar-ison. It is noted that for each atom, the local Voigt matrix is first constructed via eqn. 7,from which the local Kelvin matrix is built and is then diagonalized (after the dilation com-ponents are projected out) to obtain five eigen-shear moduli. These are then ordered andeach order is binned separately to produce the five panels of fig. 4. As seen in ref. , theleft tail of the distribution of lowest eigen-shear moduli extends into the negative modulidomain. That some atoms have a local distortion characterised by a negative modulus doesnot entail a local material instability since their calculation involves only those neighbouringatoms with a direct interaction and not the stabilising effect of the more distant surroundingmatrix. Such low or negative eigen-shear moduli do however indicate the presence of localshear distortions that are soft. Inspection of the unweighted total and partial LAQ singlepeak distributions reveal that atoms of type 2 are slightly biased towards regions of softermoduli. The weighted LAQ single peak distribution follows this bias, confirming that atomsof type 2 are often involved in an LSE.To gain more direct information on the type of atom involved, normalised atom type LSEweighted distributions are generated. Fig. 5a displays these for both the activated statesand the final states. Both distributions peak at the average atom-type of 1 and 2, andfor intermediate values are essentially a flat distribution with a slight bias towards LSEswith higher average atom-type. This suggests that the most probable specific chemicalcomposition will consist of only type 2 atoms, however on average an LSE will contain amixture of both atom types with a bias towards atoms of type 2. Inspection of those LSEscontaining only atoms of type 2 reveals them to consist of only one or two atoms. On theother hand, LSEs containing a mixture of both atoms types can consist of many more atoms.This is demonstrated in fig. 5b which shows the scatter plot of participation number withLSE weighted atom type. Here there is a clustering of small LSEs around atom type 2,but many larger LSEs consisting of both type 1 and type 2 atoms are evident. The lowerboundary in this figure reflects the fact that for an average atom type of 1.5 at least twoatoms must be involved.From the perspective of thermal activation the rate of occurrence for a particular LSEis given by ν exp( − E /k b T ) where ν is the so-called pre-factor (rate of attempt) and E IG. 4: Normalised distribution of five local eigenshear moduli where a) to e) represents thelowest to highest values. The red vertical lines represent the corresponding mean value derivedfrom the total sample. IG. 5: a) Normalised distribution of LSE weighted atom types for the activated and final stateconfigurations and b) Scatter plot of LSE weighted atom type against participation number. is the barrier energy of the LSE. Therefore it will be those LSEs with the lowest barrierenergy that will most likely occur. In the work of Koziatek et al this was demonstratedwithin harmonic transition state theory, where the corresponding pre-factor, ν , was alsocalculated, showing that despite a wide range of pre-factor values (spanning many orders ofmagnitude) those identified LSE with lowest barrier energy generally exhibited the highestrate of occurrences. It is therefore of interest to investigate the correlation between barrierenergy and the structural environment as presently defined by the LAQs of sec. II B.Fig. 6 displays scatter plots of the LSE weighted LAQs with respect to their correspondingbarrier energy. Data is shown for local cohesive energy, volume, pressure and bulk modulus.Inspection of these figures reveal strong scatter and no strong linear correlation. This isevidenced by the associated Pearson correlation coefficients shown in tab. II. The exceptionto this trend is that of local cohesive energy, which has a Pearson correlation coefficientequal to approximately -0.5, indicating a non-negligible linear correlation with the barrierenergy. Fig. 7 shows similar scatter plots for the five Kelvin eigen-shear moduli. Littlecorrelation is again evident, although the Pearson correlation coefficient for the lowest eigen-shear modulus is approximately 0.25 with the value progressively decreasing with increasingeigen-shear number — see tab. II which lists these coefficients for all five moduli. Whatcauses these weak correlations? A closer inspection demonstrates that the origin is againatom type, where those LSE with predominantly atoms of type 2 are not only more common,15 IG. 6: Scatter plots of LSE weighted local a) energy b) volume c) pressure and d) bulk modulusagainst activation energies. but they also appear to correspond to lower barrier energies. This is directly seen in fig. 7fwhich is a scatter plot of atom type and barrier energy. This plot shows again that there is abias towards LSEs with a majority of type two atoms with the corresponding barrier energybeing on average slightly lower than the barrier energies associated with LSEs dominated byatoms of type one. The scatter associated with this trend is however large, with both typesof LSEs having a spread in barrier energy comparable to the domain of the distributionshown in fig. 2a. 16 uantity Activated State PCC Final Relaxed State PCCBM 0.1593 0.1039E -0.4793 -0.4897P 0.1264 0.0084V 0.0821 0.0522ES1 0.2850 0.2473ES2 0.2170 0.1703ES3 0.1611 0.1343ES4 0.1051 0.0461ES5 0.0506 0.0154TABLE II: Pearson Correlation Coefficients (PCC) representing the correlation between localatomic quantity and activation energy
C. Correlation with Vibrational Modes
Fig. 8a shows the plot of the “vibrational participation numbers” (eqn. 11) against eigen-value number corresponding to each eigenstate, where the vibrational participation numberrepresents the effective number of atoms participating in a vibrational eigen-mode. Here theeigenvalues have been sorted from smallest to largest and therefore the horizontal axis is pro-portional to increasing vibrational frequency. Such data has been calculated before ,demonstrating that at low frequencies the eigen-modes are strongly heterogeneous indicat-ing quasi-localized mode behaviour that is believed to underly the well known Boson peakphenomenon of disordered systems . At the highest frequencies the participation numberagain drops reflecting a localization that is understood within the framework of Andersonlocalization . To investigate any correlation between the existence of such low-frequencyquasi-localized and high-frequency localized modes, the overlap between the vibrationaleigenstate and that of the LSE is determined. This is presently done by calculating thescalar product of the atomic weights (eqn. 2) with the eigenvector magnitude-squared, | (cid:126)u i,n | .Fig. 8b displays both the average overlap and the maximum overlap of all identified LSEswith each vibrational eigen-mode.The figure shows that there exists, on average, little overlap over the entire frequency17ange. In the low frequency regime, the average overlap is a well defined statistical quantityindicating that irrespective of the nature of the quasi-localized mode, the spatial locationand extension of the identified LSEs are similar for different modes. This is also reflectedin the maximum overlap which varies little with eigenstate, and is also a small quantity.At higher frequencies the situation is somewhat different in that there is much more scatterin the average value and that the maximum value. This however does not indicate anyimportant correlation since the small average and large maximum values are more likelyto indicate the scenario that statistically there will generally be one or some LSEs thatstrongly overlap with a one or some well-localized high frequency eigenstates. This doesnot occur at the low-frequency quasi-localized eigenstates since these are more extendedinvolving several tens to hundreds of atoms (see fig. 13 of ref. ). Thus the current analysisreveals little correlation between the spatial extent of the vibrational modes and the locationof the identified LSEs. D. Effect of Quench Rates
The ART n simulations have been also performed on glass samples prepared via fasterquench rates (sample 0a with η , sample 0b with η , sample 0c with η ), and a comparativestudy has been done to determine if the results of sec. III depend on the quench rate. Fig. 9adisplays the resulting activation energy distributions, showing that for the more rapidlyquenched systems the peak of the distribution shifts to lower activation energies and thatclose to the zero activation energy limit, the distribution does not limit to zero (a resultalso found in the work of Rodney and Schuh ). Fig. 9b now shows the participationnumber distribution (eqn. 3) revealing that with increasing quench rate there is a slightshift to a larger number of atoms being involved in the LSEs and that many of these involvea more mixed number of atom type — see fig. 9c which displays the average atom typedistribution of the LSEs. These results tend to suggest that the more rapidly the modelsystem is quenched the more shallow the local potential minimum is of each resulting 0Katomic configuration, and that the corresponding LSEs are somewhat larger involving bothtypes (sizes) of atom. In general however the weak correlation with local structural featuresseen in sec. III is insensitive to quench rate. For example, fig. 9d displays the LAQ weightedlowest Kelvin eigen-shear distribution as a function of quench rate, which displays only a18light shift towards lower elastic stiffness moduli. IV. ATOMIC VISUALISATION
In this section six identified LSEs are atomistically visualized. These examples have beenchosen since they represent typical features seen in all LSEs and are shown in fig 10. In all ofthe examples, the initial and final atomic positions are represented respectively by green andorange spheres, where the red arrows represent the displacement from the initial to activatedposition and the blue arrows the displacements from the activated to final position. Onlythose atoms are shown which are displaced by more than 0 . σ , either between the initialand activated or the activated, or final configuration. The large spheres represent atomsof type 1 and the small spheres atoms of type 2. Generally, the visualised atoms maybe classified into two groups, those central atoms that involve significant and irreversibledisplacement and those atoms that accommodate this activity either via reversible elasticor plastic activity. In all figures, the first class of initial atom positions are numbered, withthe dashed-corresponding-number labeling their final position.Fig. 10a represents an LSE with an activation energy of 11 . (cid:15) involving 8 atoms. Inthis LSE, the central atomic structure forms a symmetrical ring like (or closed chain-like)structure (1 → (cid:48) : 2 → (cid:48) : 3 → (cid:48) ) consisting mainly of smaller atoms (of type 2).Surrounding this plastic inner structure, there are mainly larger atoms (of type 2) whichmove back and forth during the initial to activated state and then from activated to final statetransition respectively. This is an example of elastic accommodation mechanism around theinner ring-like plastic rearrangement. This is an example of an LSE that results in a finalconfiguration identical to the initial configuration apart from a permutation of three labels.It is such LSEs, numbering 81, that have been removed from the data-set studied in sec. II B.Fig 10b represents an LSE with an activation energy of 6 . (cid:15) involving 17 atoms. It showsan extended chain-like atomic motion with the sequence being specified by (1 → (cid:48) : 2 → (cid:48) :3 → (cid:48) : 4 → (cid:48) ). Although smaller atoms (of type 2) are involved in the formation of thechain, there are a relatively large number of large atoms (of type 1) which are responsiblefor accommodating this structural excitation. There is evidence of both elastic and plasticaccommodation by the surrounding atoms, which is clearly seen by atoms moving backand forth as well as atoms moving irreversibly in the region surrounding the inner chain19ike formation. An observation that is inferred out of this LSE (and further confirmed inthe subsequent descriptions of LSEs) is that, excitations which involve a higher number ofatoms in the chain-like structure also involve a proportionally higher number of atomics inthe surrounding accommodation mechanism. Fig 10c represents an LSE with an activationenergy of 10 . (cid:15) involving 9 atoms. It also shows a chain-like atomic reconfiguration, now ofa strongly curved extension. The sequence is specified by (1 → (cid:48) : 2 → (cid:48) : 3 → (cid:48) : 4 → (cid:48) :5 → (cid:48) ). Here both sized atoms are involved in the re-configuration. There are two biggeratoms and one smaller atom surrounding this chain which undergo elastic displacement.One of the most spatially extended LSEs identified by ART n is shown in fig. 10d. ThisLSE has an activation energy of 16 . (cid:15) and involves 17 atoms, with the reconfigurationsequence being (1 → (cid:48) : 2 → (cid:48) : 3 → (cid:48) : 4 → (cid:48) : 5 → (cid:48) : 6 → (cid:48) ). It is noted thatsmaller atoms (of type 2) are involved at both ends of the chain sequence, and that bothtypes of atoms are involved in the elastic and plastic accommodation. Such mixed atom typechain-like activity is also seen in the smaller ring-like LSEs as shown in fig. 10e, which hasan activation energy of 10 . (cid:15) . Finally, fig 10f represents an LSE with an activation energyof 8 . (cid:15) involving 11 atoms. This LSE forms a chain (1 → (cid:48) : 2 → (cid:48) : 3 → (cid:48) : 4 → (cid:48) )which almost resembles a straight-line due to its low curvature.Upon inspection of these figures, the chain like sequence of an LSE generally involvesone atom replacing its neighbour (and so on) such that the chain or part of the chain isfully connected (with respect to the red and blue displacement arrows) or one atom canmove a previously unoccupied location with another atom doing the same with respect toanother atom (and so on) forming a disconnected chain (with respect to the red and bluedisplacement arrows). The smaller ring like structures of figs. 10a and e, fall into the firstcategory and the extended chains fall into both categories. Very low energy LSEs were alsovisualized (not shown) and these tended to involve just one atom changing its location withminor elastic and plastic accommodation in the surrounding regions. Such LSEs typicallyhave activation energies in the range of less than ∼ ε .Fig. 10 demonstrates that LSEs can be spatially extended. To better quantify this obser-vation, a non-dimensional quantity derived from the radius of gyration is used. The radiusof gyration (which is also commonly used in molecular applications) in the present context20s given by, R g = (cid:118)(cid:117)(cid:117)(cid:116) N LSE N LSE (cid:88) k =1 ( r k − r cop ) , (12)where N LSE is the number of atoms involved in the LSE, r k is the position of the concernedatom and r cop is the centre of position of the LSE. The dimensionless shape number is thendefined as R g /R max , where R max is the maximum distance from the centre-of-position ofan atom within the LSE. The shape number is computed using the initial positions of theatoms involved in the LSEs. A higher displacement cut-off of 0 . A was used to excludethe surrounding accommodating atoms (which have relatively smaller displacements), thusensuring only the central atoms within the LSE are considered.Fig. 11a displays the histogram of shape numbers derived from all of the identified LSEs.To understand this figure, some limiting cases of the shape number formalism are firstconsidered. When only one central atom is involved in an LSE, the shape number becomesindeterminate and when two atoms are involved the shape number becomes its largest valueof unity. An LSE consisting of an arbitrary number of atoms all situated on the surface(perimeter) of a sphere (circle) also gives a value of unity. Fig. 11b displays the shapenumber as a function of N LSE for the cases of a perfect linear chain and a spherical volume inwhich the positions of the atoms are chosen randomly. The data corresponding to the lattercase also includes the converged variance. In fig. 11a, the peak at a shape number of unity istherefore either due to a strong population of perfect rings containing an arbitrary numberof atoms or the many two atom LSEs discussed previously. Inspection of fig. 11c, whichplots the scatter diagram between shape number and participation number, demonstratesthe latter case, where there exists a dense line of LSEs with a participation number ofapproximately two at a shape number of one. For larger LSEs fig. 2c shows that typicalparticipation numbers are between 3 and 5. In this regime, the shape number for the limitingcases of fig. 11b admit shape numbers ranging between 0.7 and 0.85, which is precisely thelocation of the central peak in fig. 11a. Fig. 11c shows that LSEs consisting of greater thanfive atoms have shape numbers spanning the entire range of possible values indicating thatboth extended and more compact LSE structures are contributing to the peak, a conclusionthat is seen in fig. 10. The enhanced tail at low shape numbers does suggest larger extendedLSE chains. Importantly, in most cases an LSE may be characterized by an approximatesequence of atoms replacing atoms rather than a random rearrangment of atomic positions.21ig. 11d displays a scatter graph between shape number and activation energy showing thatthere is little correlation between the spatial extension of an LSE and its barrier energy.
V. DISCUSSION
The results of sec. III suggest that the location of an LSE is only weakly correlated withthe local structural features of those atoms involved. For the LJ system considered, the onlynon-negligible correlation is that the smaller atoms of type 2, are more often involved thanthe larger atoms of type 1, particularly when the LSEs consist of only a few atoms and areat the lower range of the activation energy spectrum. Despite the strong scatter, this latteraspect suggests a rather intuitive scenario where, because type two atoms generally involveless negative bond energies, the breaking of bonds that must occur in an LSE requires lessenergy and therefore the activation energy should be lower. Indeed this appears to be moreimportant than local Voronoi volume since in fig 3b, the volume LSE weighted distributionexhibits only a central peak structure not located at volumes typical of type 2 atoms, whereasthe local cohesive energy LSE distribution clearly correlates with the type 2 unweightedpeak (fig 3a). The remaining, somewhat weaker, correlation with activation energy is thata small or negative lowest local Kelvin eigen-shear tends to have a low activation energy.Again, an intuitive result since a small or negative Kelvin eigen-shear indicates a shallowpotential energy minimum and therefore a smaller activation barrier. In other words, theLSEs occurring in softer regions tend to have lower activation energies.The atomistic visualization shown in sec. IV generally demonstrates LSEs to be a sequenceof atoms that successively replace each others approximate location, with the surroundingatoms accommodating such movement through either elastic or plastic distortion. Thisappears to be a general result, although the spatial extension of the atomic sequence can bequite diverse, ranging from an almost linear extension, to strongly curved and closed ring-likestructures (for the smaller LSEs). Although those central atoms within the chain show noobvious decrease in their own local Voronoi volume, it is of interest to investigate whethernearby free volume is correlated with their existence. This is motivated by the originalassumption of Spaepen in his thermally activated free-volume theory . To determine thespatial location of free volume within the computer generated sample, the simulation cell isfilled with a fine regular cubic mesh of points, at a spacing much smaller than the typical22nter-atomic distance of ∼ σ . Those mesh points that have a distance to the nearest atomgreater than R max , and which are connected to each other, will then define the spatial extentof a region of free local volume. R max cannot be too small since then the normal interstitialregions, which span the entire simulation cell, will be identified. The parameter shouldalso not be too large since then no free volume will be identified. Such a method has beenused to identify free volume in grain boundaries . Using a value of R max = 0 . σ , fig. 1displays the identified regions as green balls. This figure shows that the computer generatedsample contains a few regions in which local free volume is above the normal background ofinterstitial regions.Fig. 12 now shows a histogram of, R fv , the nearest LSE centre-of-position (the redcoloured balls in fig. 1) to all identified free volume (the green coloured balls in fig. 1).Also shown are histograms of four random realizations of LSE center-of-positions derivedfrom a uniform distribution within the simulation cell. Inspection of this figure shows aslight bias of the ART n identified LSEs to be closer to free volume than that of an entirelyrandomly located LSEs. Thus there exists some correlation between the location of an LSEand nearby free-volume.The framework of fast β and slower α structural transformations, originally developedfor under-cooled liquids is now often applied to the regime of amorphous solids . Indeed,in the work of Harmon et al multiple microscopic β structural transformations (whichare assumed to be reversible) underly the emergence of irreversible α transformations inthe form of a less local release of elastic energy. In more recent work , which attempts toexplain dynamical-mechanical-spectroscopy data, some of these authors have postulated β structural transformations to consist of atomic chains of smaller atoms, rather like thoseencountered in the present work. On the other hand structural transformations that involveboth the small and large atoms, tend to reflect α transformations — a rather intuitive picturesince movement of the larger atoms will tend to involve more atoms due to accommodationissues and therefore be inherently less local. The usage of the terminology of α and β transformations in the regime of the amorphous solid is an interesting development giventhat from the under-cooled liquid perspective the α relaxations are assumed to be frozenout below the glass transition temperature . How this freezing occurs, and how far itextends to temperatures and affects plasticity below the glass transition has recently beenconsidered by one of the present authors from the perspective of a thermal activation theory23f deformation .Within the above framework a relevant question is, to which class of relaxation processes( α or β ) should the identified LSEs belong. Fig. 13a displays a histogram of the changein energy between the initial and final atomic configurations found by ART n . In mostcases this energy is positive, with a few LSEs leading to a decrease in energy and thereforea more stable atomic configuration than the initial configuration reached by dynamicalatomistic simulations. Given that only LSEs are considered which have a direct path betweenthe initial and activated states (that is, there exists no intermediate stable configuration),fig. 13a suggests that the initial configuration is in the basin of a much larger PEL valley andtherefore in the valley of the α landscape. From this context, the ART n method appears tobe probing primarily the β PEL involving the first LSE stage that would begin the atomicconfiguration’s journey out of its current α mega-basin. Fig. 13b shows the scatter plotbetween the change in energy between the initial and final atomic configuration and thecorresponding activation energy. The plot demonstrates the obvious fact that an activatedenergy cannot be less than the final state energy for LSEs that are directly connected tothe initial state. The figure also reveals that those final states that have an energy less thanthe initial state are separated by the full spectrum of possible activation energies with onlya very few having quite small activation energies. Generally, little correlation is seen apartfrom the observation that both energy scales are comparable demonstrating that, if theassumed surrounding α energy landscape does exist, the underlying “ripple” β energy scaleis that of the LSE energy scale. Thus it could be said that β processes are not so sensitiveto their local environment and that no such statement can be made for the α processes.Fig. 10 shows however that identifying LSEs as β processes has the consequence that bondsare broken for the latter — a result that is different from the view point that only α processesinvolve the breaking of bonds (see for example ref. and references therein). Clearly furtherwork is needed to confirm this picture.Whilst the results found here are not sensitive to a LJ sample realization, as well as thequench rate, the sample size is rather small containing only 1728 atoms — a cell size sidelength of only (cid:39) σ . The computational load of the ART n algorithm makes it difficult toapply the technique to much larger systems. Application of the ART n method to similar LJsystems containing 34,000 atoms (corresponding to side lengths of (cid:39) σ ) do reveal LSEs ofsimilar size to that encountered here suggesting that the LSEs accessible to ART n are at the24icroscopic length scale . Noting that ART n converges to nearby minima, this is entirelycompatible with the α − β potential energy landscape in which the larger and irreversible α transformation landscape is rippled with those of the underlying β transformation landscape.The usage of more realistic inter-atomic multi-component potentials is not expected tofundamentally change the current results, although some fine details of the LSE structureunique to the potential/system are expected to occur which are not present in the consideredmodel LJ studied here. Finally it is emphasized that it is assumed that ART n provides anunbiased probe to nearby saddle-point configurations within the potential energy landscape— an assumption that underlies all previous work applying ART n to structural glasses . VI. CONCLUSIONS