Atomic-scale expressions for viscosity and fragile-strong behavior in metal alloys based on the Zwanzig-Mountain formula

We combine the shoving model of $T$-dependent viscosity of supercooled liquids with the Zwanzig-Mountain formula for the high-frequency shear modulus, using the $g(r)$ of MD simulations of metal alloys as the input. This scheme leads to a semi-analytical expression for the viscosity as a function of temperature, which provides a three-parameter model fitting of experimental data of viscosity for the same alloy for which $g(r)$ was calculated. The model provides direct access to the influence of atomic-scale physical quantities such as the interatomic potential $\phi(r)$, on the viscosity and fragile-strong behavior. In particular, it is established that a steeper interatomic repulsion leads to fragile liquids, or, conversely, that"soft atoms make strong liquids".


I. INTRODUCTION
Different views of the glass transition have led to quite different descriptions of the viscosity of supercooled liquids.The kinetic view of the glass transition, which relies on a substantial continuity between liquid and solid glass, goes back to pioneering ideas of Y. Frenkel [1,2] and provides the basis for Dyre's shoving model and its ramifications.The entropic view of the glass transition, instead, based on the Adam-Gibbs scenario and later developed into a random first-order (ideal glass) transition, has led to a suitably modified Vogel-Fulcher-Tammann (VFT) equation with parameters that can be related to the entropy of cooperatively rearranging zones [3].Another approach based on the Adam-Gibbs scenario led to the Mauro equation for the viscosity [4].Yet a different type of approach based on Doremus' model [5] of viscosity where bonds are broken under shear flow leads to a two-exponential form for the viscosity as a function of temperature [6], which typically provides a better fitting to experimental data compared to single-exponential expressions [7].
In general, closed-form expressions for the viscosity of supercooled liquids contain three fitting parameters, which are typically related to microscopically poorly defined quantities such as free volume or entropy.The shoving model developed by Dyre provides a different approach in this sense, as it links viscosity to thermallyactivated jumps of atoms out of the nearest-neighbor cage, as in Frenkel's and Eyring's early approaches, with an activation energy which is described rigorously by means of continuum mechanics.In particular, the activation energy is expressed via the product of the highfrequency shear modulus of the liquid and an activation volume, which follows from the analysis of the work done by a particle to shove around the surrounding atoms to escape from the cage.A similar relation between energy barrier and shear modulus is provided by the Cooperative Shear Model (CSM) where, however, the characteristic volume is larger and can be connected to the concept of shear transformation zones (STZs) [8].
The shoving model provides the starting point for a more microscopic description of the viscosity and relaxation time of supercooled liquids.In particular, upon approximating the shear modulus G ∞ with Born-Huang (affine) lattice dynamics (as appropriate for the high-frequency modulus), G ∞ can be directly related to the short-range part of the radial distribution function (RDF) g(r), and hence to the interatomic potential.This led to the Krausser-Samwer-Zaccone (KSZ) equation [9], which expresses the T -dependent viscosity in closed-form in terms of the thermal expansion coefficient α T , the interatomic repulsion steepness parameter λ (obtained from a power-law fitting of the RDF up to the maximum of the first peak) and the activation volme V c mentioned above.The KSZ equation reads as where C G is the value of the G ∞ at T g , again evaluated analytically with the Born-Huang formula.This equation provides a two-parameter fitting of viscosity data since λ is determined by fitting of the g(r) data, α T is an experimentally determined quantity, which leaves η 0 , T g and V c as the only parameters, with the important constraint that V c ∼ 10 −30 m 3 .The glass transition temperature T g , in the above formula, defines to the temperature at which the glass state loses its rigidity [10].Furthermore, it has been shown that the repulsion parameter λ can be mapped accurately on realistic microscopic parameters of Thomas-Fermi screening length and Born-Mayer repulsion for the electron-gas screened repulsion between two ions in metals.
In spite of its success in providing, for the first time, a direct connection between viscosity (and fragility) and the microscopic physics of atomic-scale structure and interactions [11], the KSZ equation suffers from an intrinsic ambiguity in determining the λ parameter from experimental or simulated g(r) data.While the fitting protocol of Ref. [9] provides a consistent assessment of the effect of interatomic repulsion on T -dependent viscosity and the fragility m (defined as slope of η(T ) near T g ), other protocols [12] have provided ambiguous results.In particular, the original protocol of Ref. [9] has shown that a steeper interatomic repulsion results in a fragile behavior of the glass-forming liquid, whereas a softer repulsion is associated with strong liquids.A different protocol for extracting λ with less prescriptive constraints on the fitting was used in Ref. [12].In that approach, λ was taken to be a free parameter which also depends on temperature, and the opposite scenario, i.e. softer repulsion leads to fragile behavior, was found.However, it has been later demonstrated that λ does not depend on temperature [13], which invalidates this fitting protocol.
Here we develop a different, perhaps more sophisticated, approach, which combines the shoving model with the microscopic Zwanzig-Mountain formula for the G ∞ of liquids.This leads to semi-analytical expressions for η(T ) and for m, which directly link these quantities to the g(r) and to the interatomic potential φ(r).Upon successfully calibrating these expressions for the case of Cu 50 Zr 50 , a new interatomic repulsion paramter l is identified which is unambiguously linked with the repulsive part of φ(r).Upon letting this parameter vary, fictive materials with different interatomic repulsion softness are generated.The model analysis demonstrates that steeper interatomic repulsion leads to fragile behavior, thus reaffirming the correctness of the fitting protocol of Ref. [9] for the identification of λ in the KSZ equation above.It also confirms the qualitative increasing trend of fragility m increasing with potential repulsion steepness l or λ, and recovers the linear trend already seen for m(λ) in Ref. [9].

II. THEORY A. The shoving model
We base our derivation of the viscosity of liquid metals here upon the so-called shoving model [14].The assumption at the basis of this model is that, within the transition state theory [1], the activation energy of the average relaxation time is determined by the work done in shoving aside the surrounding liquid to allow "flow events".The main result of this model is the temperature dependence of the viscosity, which is related to the temperature dependence of the high frequency limit of the shear modulus G ∞ (T ): where η 0 is a constant prefactor and V c is the characteristic volume of the group of atoms involved in the shoving event (on the same order of magnitude of the nearestneighbour cage), and is a weakly dependent function of T such that this dependence is typically neglected.The T -dependence of viscosity η is thus directly controlled by the T -dependence of G ∞ , while V c has a dependence on temperature which is more difficult to assess and for simplicity is normally taken as T -independent [14].
The values of V c typically found are on the order of the atomic size, which is much smaller than the values of the characteristic volume normally found within the CSM model.In Frenkel's original derivation, activation energy in the Arrhenius exponential of Eq. ( 1) is given by E = 8πG ∞ r∆r 2 , where r is the cage radius and ∆r is the increase of the cage radius due to thermal fluctuations which enables the atom to escape the cage (see page 193 in Ref. [1]).

B. High-frequency shear modulus from the Zwanzig-Mountain formula
Zwanzig and Mountain derived a popular expression for the high-frequency shear modulus of liquids, which reads as [15]: where g(r) denotes the radial distribution function (RDF) of the atoms in the liquid (average of the atomic species in an alloy), ρ(T ) is the total atomic density, and φ(r) is the average interatomic potential between any two atoms.
Two crucial input to evaluate G ∞ are, therefore, the RDF g(r) and the interatomic potential φ(r).In lieu of a theoretical expression for g(r) valid for real metal alloys, which is obviously beyond reach, we base our analysis on MD simulations data of the average g(r) in Cu 50 Zr 50 alloys, from Ref. [11].
In order to provide an analytical handle on the various features of the g(r), which in turn are connected to features of the interatomic potential φ(r) as shown below, we proceed to the following parameterization of the g(r).

C. Analytical parameterization of g(r)
The parameterization of g(r) is given by a sum of three terms, and it depends on 6 parameters.The first term describes, at the same time, the short-range repulsion and impenetrability of the atoms as well as behavior for large interatomic separations.This is mathematically expressed by the following asymptotic limits: g(r) goes to zero for small r and tends to one for r → ∞.These asymptotic behaviors are encoded in the hyperbolic tangent.The second term represents the first nearest-neighbor shell, that corresponds to the first peak of the RDF, while the third term describes the decrease of structure with the distance, that corresponds to the decreasing height of the second peak and of the following Friedel oscillations.The first peak is given by a Gaussian, while the second peak and the Friedel oscillations are described by a decaying exponential multiplied by an oscillating function.Here we report the parameterization: This expression contains six parameters for which we now give a qualitative description.Parameter a is approximately the position of the first peak (i.e. the center of the Gaussian).Parameter b is linked to the width of the first peak since it's proportional to the variance of the Gaussian.Parameter l, which is the most important for our subsequent analysis, gives the asymptotic power-law trend of the first term, hence it represents the steepness of the ascending part of the first peak of g(r).Parameter k is the height of the Gaussian distribution, so it is linked to the height of the peak, which is also influenced by the other two terms.Finally, w is related to the frequency of the Friedel oscillations while h controls the decay of the envelop of the oscillations.
The fitting of the MD simulations data of g(r) for the Cu 50 Zr 50 system is shown in Fig. 1 1 (a)) that the ascending part of g(r), before the first peak, is perfectly described by a power-law ∼ (r/a) l , with l = 20, over more than two decades in r.

D. Determination of the interatomic potential φ(r)
The RDF is related to the potential of mean force, φ(r), via the reversible work theorem: the proof of which can be found in the textbooks [16].
The mean interatomic potential between two atoms φ(r) determined in this way thus accounts for many-body effects from the surrounding electronic and atomic environment.Using the fitting of the MD g(r) data for the Cu 50 Zr 50 alloy we obtain the interatomic potential profile shown in Fig. 1 (b).The potential features a rather steep interatomic repulsion due closed electron shell repulsion followed an attractive bonding minimum mediated by the nearly-free electrons.After the minimum, the oscillations represent the Friedel oscillations in the electronic density.

III. RESULTS AND DISCUSSION
A. T −dependent viscosity Using Eqs. ( 3)-( 4) calibrated on the g(r) data for the Cu 50 Zr 50 alloy inside the Zwanzig-Mountain formula Eq.
(2), we are now able to evaluate the T -dependent viscosity η by means of the shoving model, Eq. (1).
Upon denoting the k B T -normalized integral in the Zwanzig-Mountain formula Eq. ( 2) as I, we therefore arrive at the following expression for the viscosity: Note the cancellation of a k B T factor contained in G ∞ (and recall that φ(r) = −k B T ln g(r)) with the k B T factor in the denominator of the argument of the exponential in Eq. ( 1).Here the T -dependent density ρ is expressed in terms of the thermal expansion coefficient α T , using the definition of the latter α T = V −1 (∂V /∂T ) T , leading to: where ρ 0 = 7 • 10 3 [kg/m 3 ] is a known value of density at the reference temperature T 0 = 298K [17].The above relation is normally linearized due to the small value of α T ,  Viscosity (mPas) Figure 2: Model fitting of experimental viscosity data from Ref. [11] using Eq. ( 5) with V c = 5.6 • 10 −29 m 3 , η 0 = 14 mPa•s, I given by the integral in Eq. ( 2) with φ(r) given by Eq. ( 4) and g(r) from the analytical parameterization Eq. ( 3) of MD simulations data of g(r) from Ref. [11].An effective value of thermal expansion coefficient has been taken, α T = 0.00377K −1  for the comparison.
]. Also, one should note that the integral I in Eq. ( 2) is also T -dependent due to the k B T factor in the definition of φ(r) in Eq. ( 4).In Fig. 2 we report the fitting of the viscosity data of Cu 50 Zr 50 as a function of temperature ( measured with a levitating drop method in Ref. [11]) by means of Eq. ( 5).There are two fitting parameters in the comparison, one is the effective thermal expansion coefficient α T = 0.00377K −1 , which is significantly larger than the typical values for metallic melts (∼ 10 −4 K −1 ) and effectively compensates for the neglected T -dependence of the activation volume V c and of the atomic structure given by g(r), since their T -dependence is not known.In particular, the T -dependence of V c may play an important role, since it was shown to be a rapidly decreasing function of T upon approaching T g from below [18].However, the dependence of V c on T in the high temperature liquid phase is not known, and this may explain the larger value of the fitted α T coefficient, which makes up for neglecting the decrease of V c upon increasing T .
The other fitting parameter is V c which is found to be equal to 5.6 • 10 −29 m 3 , i.e. in close agreement with typical values of the shoving volume found in previous works [9], and corresponds to the characteristic size of the nearest-neighbor cage in disordered metals.However, it is significantly smaller than the typical size of a cooperative flow event [8].
All in all, Eq. ( 5) provides a three-parameter fit of viscosity data over a broad range of T , and, unlike other popular three-parameter models such as VFT, the Avramov-Milchev (AM) equation [19] and the Mauro equation [4], all the parameters can be traced back to atomic-scale structure and interactions, in that being similar to the KSZ equation [9].The latter still retains the favorable advantage, over Eq. ( 5) of being in simple, fully analytical form.

B. Effect of interatomic potential on viscosity and fragile-strong behavior
Thanks to the direct connection that the above model provides between η(T ) and microscopic atomic-scale parameters, it is possible to analyze the effect of the interatomic potential φ(r) on the viscosity and on the fragilestrong behavior.We use the interatomic repulsion steepness parameter l in Eq. ( 3) as a proxy to design fictive materials of varying interatomic repulsion, with the aim of studying the effect of the interatomic repulsion on viscosity and fragility.
We start from the level of the RDF g(r) and vary the repulsion steepness parameter l around the value (l = 20) that we found in the fitting of MD simulation data (Fig. 1 (b)).We thus obtain the fictive RDF's shown in Fig. 3 (a).It is clear that large values of l correspond to a sharp rise of the first peak of the RDF, whereas low values of l correspond to a less steep rise of the peak.Using Eq. ( 4) we then obtain the corresponding potentials of mean force describing the interatomic potential, shown in Fig. 3 (b) for the same values of l shown in Fig. 3 (a).It is evident that large l values (l 20) result in steep interatomic repulsion, whereas low values of l (i.e.l < 10 − 20) result in softer repulsive potentials.Based on Eq. ( 3), when r/a < 1 (i.e. at short range before the first peak of the RDF) the tanh can be approximated as a linear function, hence ∼ (r/a) l .Before the peak, the parameter l is thus closely related to the λ parameter of the KSZ equation, defined as g(r) ∼ (r − σ) λ where σ is a hard-core atomic size.Even though the power-law trend to describe the interatomic repulsion is a common feature of KSZ and of the present approach, the presence of a hard-core cut-off size σ in the KSZ model is a significant difference, which quantitatively may lead to different values of l and λ, although it should not affect the qualitative trends.
It is interesting to note that the parameter l affects only the repulsive part of the potential, whereas, remarkably, the attractive part of the potential, from the minimum on, is almost unaffected.This allows us to single out the effect of the interatomic repulsion, which strongly depends on the atomic composition of the metal.
The corresponding viscosities as a function of T calculated based on the RDFs and φ(r) profiles of Fig. 3 are shown in Fig. 4 (a).Upon varying the interatomic steepness l, it is clear that the slope of the η(T ) changes.In particular, larger values of steepness l correspond to larger (steeper) slopes of the viscosity curves.The microscopic explanation for this behavior resides in the k B Tnormalized integral I in the ZM formula Eq. ( 2).The value of I in Eq. ( 4) increases systematically upon increasing l, which results in a G ∞ more steeply rising with T , hence in a larger slope of η(T ).This is mainly due to the factor d dr r 4 dφ dr inside the integral, which clearly gives larger contributions as the repulsive decay of φ(r) becomes steeper.This fact can be easily checked on the simple example of a power-law repulsive decay φ(r) ∼ r −n : the larger the exponent n the larger the contribution of this factor to the integral.
The fragility of a glass-forming liquid is given as the slope of the viscosity at glass transition temperature, that is: m = ∂ log 10 (η/η0) ∂(Tg/T ) T =Tg [20].Upon applying the definition to our viscosity formula Eq. ( 5), we obtain the following expression: which can be evaluated by computing I for different values of l, by keeping all the other model parameters the same as in the fitting of Fig. 2.
The fragility evaluated according to Eq. ( 7) is plotted in Fig. 4 (b) as a function of the interatomic repulsion steepness parameter l.It is evident that the fragility m increases monotonically with the repulsion steepness l, or, in other words, the fragility m is lower (the liquid is stronger) with softer repulsion steepness.This fact reaffirms the conclusions of Ref. [9] that "softer atoms make stronger liquids", with a surprising robustness of this law across different materials, from colloids [21] to metals.It is important to note that the estimate via Eq.( 7) is not quantitatively predictive because we used a viscosity fitting calibrated at significantly higher T than T g where, by definition, the fragility should be evaluated (recall T g = 700K for Cu 50 Zr 50 ).However, Eq.( 7) is valid and reliable as a prediction of the qualitative trend for m as a function of l.This is because the only parameter in Eq. ( 7) which depends on l is the factor I, while all other parameters are independent of the interatomic repulsion steepness and act mainly as scale factors.This includes α T which depends on the attractive part of the potential but not so much on the repulsive part [22].
It is to be noted that the fragility m in Fig. 4 is predicted to go through a crossover from a rapidly varying increasing trend for l < 20 to a linear trend for l > 20.The latter linear trend perfectly recovers what has been found with the KSZ equation in [9], indeed in a range of larger m values.Instead, the crossover from a quickly rising initial trend of m with repulsion l into the (more slowly growing) linear m vs l regime is a new prediction of the present work.Here the data point with l ≈ 20 corresponding to Cu 50 Zr 50 is highlighted by the arrow.

IV. CONCLUSIONS
In summary, a microscopic model of the viscosity and fragile-strong behavior of liquid metals in the supercooled regime has been developed, based on combining the shoving model with the Zwanzig-Mountain (ZM) formula for the high-frequency shear modulus, G ∞ of liquids.The model has been evaluated using MD simulations data for the g(r) of the Cu 50 Zr 50 alloy, from which the interatomic potential of mean force, φ(r) has been determined upon analytically parameterizing the g(r) data.These input allowed us to evaluate the G ∞ with the ZM formula, and in turn the viscosity as a function of T .This led to a semi-analytical formula for η(T ) which provided an excellent three-parameter fit of experimental viscosity data for the Cu 50 Zr 50 alloy from the literature [11].Compared to other popular three-parameter models, such as VFT and the Mauro equation [4], this expression, Eq. ( 5), features only microscopic parameters and provides direct access to atomic-scale structural and interaction quantities, such as g(r) and φ(r).
The analytical parameterization of the g(r) led to the possibility of studying the slope of viscosity as a function of T and hence the fragile-strong behavior of the liquid, in terms of microscopic interaction potential parameters.It has been shown that a crucial parameter which controls the fragility is the interatomic potential repulsion steepness l (related to Born-Mayer repulsion and to the λ parameter of the KSZ equation [9]).It has been demonstrated that the fragility m of an atomic liquid is a monotonically increasing function of the potential repulsion steepness l.This result is independent of the values of the other parameters entering the fragility formula Eq. ( 7), and depends exclusively on the integral I in the ZM formula, which gives a direct insight into the microscopic explanation for this phenomenon.This analysis thus reaffirms that "soft atoms make strong liquids", in full agreement with previous claims on vastly different materials such as metal alloys [9] and colloids [21,23].
Finally, we also note that the above double-exponential form for η(T ) Eqs.( 5)-( 6) has the potential to effectively describe the non-Arrhenius to Arrhenius crossover that has been observed in supercooled liquids [24][25][26][27][28][29][30], since a crossover to a single-exponential Arrhenius form (Eq. ( 1)) is predicted to occur when the thermal expansion coefficient α T is very low.In future work, this consideration may be the starting point to rationalize the extreme variability of the non-Arrhenius to Arrhenius crossover with material chemistries and bonding in terms of the bonding anharmonicity.

Figure 1 :
Figure 1: (a): Analytical parameterization of the g(r) simulations data for the binary Cu 50 Zr 50 alloy.The values of the parameters are reported in the main text.The inset shows the analytical fit of the main panel together with a power-law trend line ∼ (r/a) l , with l = 20.(b):The average interatomic potential (potential of mean force) obtained from the analytical fitting of g(r).

Figure 3 :
Figure 3: (a): Evolution of the g(r) upon varying the interatomic steepness parameter l around the value l ≈ 20 used in the fitting of the g(r) in Fig. 1 (a).From left to right l = 4, 18, 30, 60.(b):The average interatomic potential (potential of mean force) for the same values of the l parameter shown in panel (a).

Figure 4 :
Figure 4: (a): Viscosity as a function of temperature calculated upon varying the interatomic steepness parameter l around the value l = 18.77(red symbol) used in the fitting of Fig. 1 (a).From left to right l = 4, 18, 30, 60. (b): Fragility m calculated upon varying the interatomic repulsion steepness parameter l.Here the data point with l ≈ 20 corresponding to Cu 50 Zr 50 is highlighted by the arrow.