Unifying interatomic potential, g(r), elasticity, viscosity, and fragility of metallic glasses: analytical model, simulations, and experiments

An analytical framework is proposed to describe the elasticity, viscosity and fragility of metallic glasses in relation to their atomic-level structure and the effective interatomic interaction. The bottom-up approach starts with forming an effective Ashcroft-Born-Mayer interatomic potential based on Boltzmann inversion of the radial distribution function g(r) and on fitting the short-range part of $g(r)$ by means of a simple power-law approximation. The power exponent $\lambda$ represents a global repulsion steepness parameter. A scaling relation between atomic connectivity and packing fraction $Z \sim \phi^{1+\lambda}$ is derived. This relation is then implemented in a lattice-dynamical model for the high-frequency shear modulus where the attractive anharmonic part of the effective interaction is taken into account through the thermal expansion coefficient which maps the $\phi$-dependence into a $T$-dependence. The shear modulus as a function of temperature calculated in this way is then used within the cooperative shear model of the glass transition to yield the viscosity of the supercooled melt as a double-exponential function of $T$ across the entire Angell plot. The model, which has only one adjustable parameter (the characteristic atomic volume for high-frequency cage deformation) is tested against new experimental data of ZrCu alloys and provides an excellent one-parameter description of the viscosity down to the glass transition temperature.


I. INTRODUCTION
A. State of the art One of the most puzzling properties of glasses is the huge increase of viscosity, by many orders of magnitude, within a narrow range of temperature T upon approaching the glass transition temperature T g . As a consequence, considerable interest is being devoted to understanding this phenomenon in terms of the underlying atomic-level structure and dynamics of supercooled liquids. In recent years, much attention has been devoted to the roles of both short and medium-range order in promoting structural and dynamical arrest upon approaching the glass transition from the liquid side. Various order parameters have been proposed to quantify the local order in supercooled liquids, starting from the bond-orientational order parameter [1], until more recent proposals which showed evidence that glassy properties correlate strongly with the local breaking of inversionsymmetry at the atomic scale. As a matter of fact, the local inversion-symmetry breaking turns out to be the key microstructural aspect which controls both the nonaffine soft elasticity and the boson peak of glasses [2].
In terms of materials applications, there is little doubt that disordered glassy materials represent part of the future of materials science, due to their advanced applications, in particular in terms of their outstanding performance under mechanical loading [3]. Metallic glasses have emerged as the most important class of glassy ma-terials from this point of view, and it remains one of the major challenges of current research to understand the relationship between atomic dynamics and macroscopic mechanical response in these materials.
Due to the current limitations in terms of experimental and computational techniques, the most abundant and reliable information about the atomic-scale structure of liquid and glassy metals comes from two-point correlation functions such as the structure factor S(q), which can be extracted from neutron and X-ray scattering, and gives access to the radial distribution function, g(r). In the theory literature [4], attention has been devoted, over the last decades, to multi-point correlation functions such as the four-point correlation function χ 4 , which exhibits more significant changes upon crossing the T g , whereas the g(r) remains substantially unaltered upon going from the liquid into the glass. In reality, appreciable changes in g(r) can be detected upon vitrification, although the extent of these changes appears to vary from system to system, and this represents a possible way of linking structural evolution to dynamics [5].
Among the most popular pictures proposed to link the phenomenology of the Angell plot for the viscosity η versus T near T g , is the one which associates fragile glass formers (with the steepest dependence of η on T ) to an underlying steep interparticle repulsion at contact, whereas strong glasses (with Arrhenius dependence of η on T ) are associated with softer interparticle repulsion. This picture, which is largely based on the twopoint correlation dynamics and local structure, and on the Weeks-Chandler-Anderson [6] idea that the repulsive part of two-body interaction is what controls the overall structure of liquids, has been demonstrated convincingly for the case of soft colloidal glasses by the Weitz group [7]. B. "Soft atoms make strong glasses" Using a theoretical argument based on the highfrequency quasi-affine shear modulus and its relation to interatomic connectivity and thermal expansion, some of us [8] recently showed that this picture ("soft atoms make strong glasses") may be applicable to metallic glasses as well. In particular, we showed that a global interatomic repulsion parameter can be defined from theory and can be linked to the ascending part of g(r) to describe the full Angell plot (from strong to fragile) for metallic alloys of very different composition. In this theory, an important role is played also by thermal expansion: the fragile behaviour is linked with higher values of the thermal expansion coefficient α T . This is another "global" parameter, this time related to the longer-ranged anharmonic part of the interaction, which may also account, in a coarse-grained way, for effects of medium-range order. Both these global interaction parameters, λ and α T , are sensitive functions of the elemental composition and stoichiometry of the alloy, and may account for microalloying effects as well [9].

C. Medium-range atomic dynamics
Other very recent studies have appeared which substantially support this picture from different angles. Busch and co-workers [10] have shown, for various alloys, that fragile behaviour indeed correlates with higher thermal dilation, and with shallow changes in the first coordination shell. The first feature, is seen in the structural evolution on a scale of 3-4 atomic diameters, where longrange interaction effects and anharmonicity control the incorporation of volume upon changing T . The shallow structural change observed at the level of first coordination shell for fragile melts, instead, can be explained by the steepness of the potential: upon reducing T , systems with steeper repulsion experience an increased resistance towards further approach of two nearest-neighbours, thus leading to the observed shallow variation in the position of the first maximum of g(r). In recent work, this picture has also been connected to the Arrhenius crossover temperature, which marks the appearance of the high-T liquid regime where full Arrhenius behaviour of transport properties is recovered [11,12].

D. 5-fold symmetry picture
In a somewhat different picture proposed recently by Wang and co-workers [13], fragility was found to correlate with the rapidity at which 5-fold (icosahedral-like) symmetry develops upon lowering T [14,15]. Although apparently a different picture, this finding can be related to the interatomic softness parameter λ of Ref. [8]. Atoms tend to organize themselves in icosahedral clusters with their nearest-neighbours, as T decreases [16,17] . This process becomes faster with steeper interatomic repulsion (larger λ) because the re-organization energy, when two atoms are not too close to each other, is comparatively less than for softer potentials where at the same distance atoms still experience significant energy.

E. Analytical relations for η(T )
Finally, several of these recent works sought an analytical relationship for the increase of viscosity η as a function of T in the Angell plot. Kelton and co-workers [5] provided an extended Vogel-Tammann-Fulcher (VTF) relationship which involves a structural parameter related to the change in the first coordination shell. Wang and co-workers [13] proposed another extension of the VTF relationship which, instead, involves an order parameter for 5-fold symmetry. Both these relations have two free fitting parameters (in addition to the normalization constant η 0 ).
An altogether different relationship, with just one adjustable parameter, has been proposed by Johnson and co-workers [18] and further developed in Ref. [8]. This relation, furthermore, is not of the VTF type, but arises from a microscopic mechanistic picture given by the shoving model of the glass transition due to Dyre, made microscopic by using an atomic theory for the highfrequency shear modulus involving the interatomic repulsion parameter λ and the thermal expansion coefficient α T . The only adjustable parameter is the characteristic atomic volume V c for cage deformation, which was found to be related, on a larger scale, to the volume of shear-transformation zones (STZ) [19].
In this contribution, we apply theoretical ideas [8,18] to experimental data of metallic glasses to show that a bottom-up quantitative relationship can be built between the atomic-scale structure and interactions and the macroscopic T -dependent viscosity. The resulting framework, in which a major role is played by globallyaveraged interaction parameters related to short-ranged interatomic softness and longer-ranged anharmonicity, uses the short-range g(r) as input from simulation or experimental data, to arrive at the viscosity for which the theory provides an excellent one-parameter fitting. This good agreement shows that changes at the level of first-coordination shell, encoded in the short-range repulsion parameter λ, as well as changes at the mediumrange scale, encoded in the thermal expansion coefficient α T , are both important in determining the viscosity and fragility of metallic glasses.

II. INTERATOMIC POTENTIAL FOR THE ION-ION REPULSION IN METALLIC GLASSES
In recent work [8], we analysed several alloys in an attempt to extract an effective, averaged interatomic potential which describes the short-range repulsion between any two ions in a metallic alloy melt. Based on the systematic fitting of shear modulus and viscosity data for various three-and 5-component alloys we proposed the following interatomic potential which comprises two contributions: (i) the longer-ranged Thomas-Fermi screened-Coulomb repulsion modulated by the Ashcroft correction and (ii) the Born-Mayer closed-shell repulsion due essentially to Pauli repulsion. The Thomas-Fermi contribution is more long-ranged and is described by a Yukawa-potential type expression. The Born-Mayer contribution is a simple exponentially-decaying function of the core-core separation, motivated by the radial decay of electron wavefunctions for the closed shells. The effective interatomic potential reads as where is the Ashcroft factor [20], with R core being a typical value for the atom-specific core radius and Z ion the effective ionic charge number. Furthermore, a 0 is the Bohr radius andσ is the average ionic core diameter of the alloy, which corresponds to the average size of the ionized atoms constituting the alloy. The average ionic core diameter is obtained by averaging the respective ionic core diameter of the constituents with their contributing weights given by their volume ratios in the alloy. The values for the ionic core diameters of the atoms constituting the alloys are taken from Ref. [21]. The quantities A and B set the energy scales for the repulsive interaction from the Ashcroft and Born-Mayer term, respectively. The parameter q TF is the inverse of the Thomas-Fermi screening length given by Thomas-Fermi theory, and its value is known for different types of alloys [22]. We choose a representative value for q TF as 1.7Å −1 according to the values reported in Ref. [22]. The ionic core diameterσ is obtained by a weighted average of the core diameters of the atoms constituting the alloys taken from [21], where the weights correspond to the ratios of the respective atoms.
The characteristic range 1/C of the valence-shell overlap repulsion is not known a priori. However, its typical values are less sensitive to the atomic composition than the parametersσ, A and B. Different atoms have very similar values typically in the range C = 1.89 − 4.72Å −1 [23]. The latter cannot be easily estimated from firstprinciples or from literature. Similarly, the prefactor B of the Born-Mayer term, can be rigorously evaluated only from the exchange integrals of the various overlapping electrons belonging to the valence shells of the two interacting ions. This calculation, even in approximate form, is not feasible except for simple monoatomic crystals. Hence, both A and B were taken in our previous analysis as adjustable parameters in the mapping between our schematic logarithmic potential (to be introduced below in Sec.III) and the Ashcroft-Born-Mayer interatomic potential. We shall remark that the Born-Mayer prefactor B typically has non-trivial large variations from element to element across the periodic table, as shown in many ab initio simulation studies [23,24]. Consistent with this known fact, it turns out that B is the most sensitive parameter in our analysis, in the sense that small variations in B can lead to large deviations in the fitting of the experimental data. Conversely, the Ashcroft prefactor A is not a sensitive parameter, and its values are not crucial for the match with experiments.
In Ref. [8], it was found that, in order to fit shear modulus and viscosity data of various alloys, values of the Born-Mayer repulsion strength B are required which are between two and three orders of magnitude smaller than the Born-Mayer parameters tabulated for pure metals. This important difference has at least two reasons. One reason is that the Born-Mayer formula used for pure substances in the literature is written , which we use here. Evidently, exp [Cσ] > 1 partly contributes to explain this discrepancy. However, the fact that B fitted for multicomponent alloys is much smaller than B found for pure substances is also due to the so-called micro-alloying effect, whereby the addition of even a small amount of different elements with different ionic size and electronic structure induces a strongly nonlinear change in the inter-ionic potential. When there is an atomic size mismatch, this difference intuitively promotes softer repulsion due to the fact that short-ranged packing is more efficient. Of course there could be other important reasons related to the change in electronic structure, e.g. the change in anisotropy of closed electronic shell, which also make the effective short-range repulsion in alloys being effectively milder than in pure metals. Finally, the Born-Mayer parameters in Ref. [25] refer to pairs of atoms, in which the outer electronic structure is clearly much different from metals, where the atoms are significantly ionized.

III. THE GLOBAL INTERATOMIC REPULSION PARAMETER λ
In the analysis of Ref. [8], it was found that different disordered alloys can be described by an effective interatomic potential given by Eq. (1)  mined by the elemental composition of the alloy, whereas q T F = 1.7Å −1 can be kept constant, independent of composition. Furthermore, an even simpler parametrization of the interatomic potential was obtained by mapping Eq.(1) onto a logarithmic expression with a single, global parameter λ, which contains all the information about the steepness (or its inverse, the softness) of the ion-ion repulsion. This simpler expression reads as The mapping between Eq.(1) and Eq.(3) is shown in Fig.1, for representative parameters of metallic glasses. In all instances examined thus far, Eq.(1) is very accurately represented by Eq.(3).
The main advantage of representing the repulsive part of the interatomic potential) with the compact Eq. (3) is that it allows us to pack various effects into a single, global interaction parameter λ which contains information about the overwhelmingly complex details of the ion-ion interaction in metallic alloys. The interaction between two ions in metallic glasses is in fact the result of the intricate underlying electronic structure as well as of many-body and non-local effects. It is a hopeless task to devise a theory of the interatomic interaction due to this complexity, and our aim here is to present an averaged non-local parameter which, similar to the effective mass concept in semiconductors, takes all these non-trivial effects into account while still allowing one to label different alloys and their properties in terms of the microscopic interaction.
A second important advantage of the compact form, Eq.(3), is that it allows us to relate the repulsion parameter λ directly to the short-range ascending slope of the radial distribution function g(r). Upon identifying the effective interatomic potential with the potential of mean force, Boltzmann inversion provides a link between the effective potential of mean force between two ions and the local structure Importantly, Boltzmann inversion [24] provides a definition of V (r) as a non-local interaction potential which contains important many-body contributions, since the potential of mean force defines the effective interaction of two ions in the field of all the other ions and manybody interactions thereof. Upon combining Eq.(4) with Eq. (3), we obtain a relationship between λ and g(r), in the form of the following simple power-law expression The link between the g(r) and the interatomic potential V (r) is schematically illustrated in Fig. 2.

IV. ANALYTICAL EXPRESSION FOR THE HIGH-FREQUENCY SHEAR MODULUS
In this section we develop a link between λ and the high-frequency shear modulus. As is known at least since the seminal work of Zwanzig and Mountain [26], the highfrequency shear modulus of liquids can be efficiently described using affine or quasi-affine elasticity. In brief, the atoms forming the transient "cage" around a given atom are simply displaced proportionally to the imposed strain, with the high-frequency of the external driving not allowing for nonaffine rearrangements. The latter dominate instead the low-frequency response of liquids and glasses and are responsible for a dramatic softening of G, as they are associated with a negative contribution to the elastic free energy (in fact their contribution is internal work done by the system, hence negative by thermodynamic principles) [27][28][29][30][31][32].
Importantly, nonaffine response theory of Ref. [27,28,30] does correctly recover the Maxwell marginal rigidity criterion at the isostatic point at which the total number of constraints ZN/2 (where Z is the mean number of bonds per atom) is exactly equal to the total number of degrees of freedom per atom dN , with central-force interactions (leading to Z = 2d = 6 at the isostatic point). Conversely, previous approaches based uniquely on isostaticity [33] and ignoring the symmetry, are less general and of limited applicability. For example, unlike the more general nonaffine formalism, isostaticity-based approaches such as the one of Wyart [33] provide erroneous predictions of the moduli of centrosymmetric lattices; for instance they predict the shear modulus of perfect centrosymmetric lattices to scale as G ∼ (Z − 6) and not as G ∼ Z which is the correct scaling in view of the affine response ensured by local inversion-symmetry being enforced in centrosymmetric lattices.
Hence, upon taking the infinite-frequency limit of the frequency-dependent nonaffine shear modulus from nonaffine response theory [27,28], as shown in Ref. [8] and numerically confirmed in Ref. [31], the affine modulus is retrieved, which scales as G ∼ Z, where Z is the interatomic connectivity.
The link between the effective repulsion, encoded in λ, and the high-frequency G goes by the way of the atomic connectivity Z, which defines the average number of nearest-neighbours around a test atom and can be estimated by counting all neighbours contained within the first peak of g(r). This is a conservative way of counting nearest-neighbours, as opposed to e.g. also considering atoms beyond the maximum of g(r), since these effectively experience a much reduced restoring force. Furthermore, these farther apart neighbours are within the presumably attractive part of the potential of mean force and the physics of this region is already accounted for in our model by the thermal expansion coefficient α T . Therefore, at frequencies much larger than the inverse Maxwell relaxation time of the liquid, ωτ M 1, the shear modulus can be evaluated using affine elasticity theory [28], which gives Here, κ stands for the spring constant of a harmonic bond and R 0 for the average interatomic spacing at rest in the equilibrated glass. The coordination number Z refers to the average number of mechanically-active nearest-neighbours [28,29,32]. While coordination numbers are typically defined from an integral over the first peak of g(r), there is no consensus about the upper limit of the integral and on how this relates to the mechanical activity of the nearest-neighbours being counted. For example, it is well known that, if the first minimum of g(r) is taken as the upper limit of integration, then the integral yields Z ≈ 12 both in the liquid and in the glass, and this result is basically independent of the temperature of the system. Clearly, this definition is totally inadequate to estimate the average number of mechanical contacts which is required by the expression for G.
The value of g(r max ) evaluated at the first peak position r max , increases significantly with increasing the packing fraction φ (and with decreasing T ), which is typical of all dense liquids with an important repulsive component of the interaction [34]. This is related to the fact that nearest-neighbours located within the attractive minimum of V (r) are more likely to contribute to the rigidity. Metals definitely fall into this category, and their g(r) shares many features with the hard-sphere system, so that, for example, their g(r) can be expressed in terms of the hard-sphere g(r) using the celebrated Weeks-Chandler-Anderson method.
Most importantly, the increase of the value of g(r max ) with increasing atomic density is modulated by the steepness of the interatomic repulsion, since the slope of the ascending flank of the first peak of g(r) depends on the slope of the repulsive part of V (r) via the Boltzmann formula in Eq. (4). This fact can be intuitively understood by considering that the integral of g(r) up to the maximum g(r max ) must necessarily yield a larger number when the slope of the ascending part of g(r) is less steep compared to the case of a steeper slope (for the same maximum height). Hence, Z must be an increasing function of the atomic density modulated by λ.
Using the fact that the upper integration limit of r max increases with the packing fraction φ, integrating Eq.(5) up to a threshold which is proportional to φ, as done in Appendix A, yields the scaling law Z ∼ φ 1+λ . Although the upper limit of the integral could be perhaps identified with r max , since we are interested here in the qualitative behaviour we prefer to leave it as a generic threshold ∝ φ such that the limit Z → 0 is correctly recovered when φ → 0.
Moreover, the definition of the Debye-Gruneisen thermal expansion coefficient α T , in terms of the atomic packing fraction φ = vN/V (with v the characteristic atomic volume and N the total number of ions in the material) gives φ(T ) ∼ e −α T T , as discussed in Ref. [30]. According to this result, φ decreases with increasing temperature T , an effect mediated by the thermal expansion coefficient defined as α T = 1 V (∂V /∂T ) = − 1 φ (∂φ/∂T ). Replacing the latter relationship between φ and T in the equation for Z, see Eq. 6, we finally obtain a closedform equation which relates G to the two global interaction parameters, the short-range repulsion parameter λ and the attraction anharmonicity parameter α T , An exponential decay of the shear modulus with T is found also in approaches like Granato's intersticialcy theory which model the glass as a crystal with a high concentration of interstitials, see e.g. Ref. [35]. Equation (7) accounts, in compact form, for all the salient features of the interatomic interaction, and con-tains the effect of repulsion steepness (short-ranged part of V (r)) as expressed by λ, and of anharmonicity, expressed by α T . A schematic depiction of how the global parameters λ and α T are related to features of V (r), is presented in Fig. 2. It is important to note that the longer-ranged, anharmonic attractive part of the interaction in metals also stems from non-local, volumedependent terms in the interaction of the ions with the partly delocalized electron gas. Hence, a microscopic description in terms of pair-interactions alone is generally not valid, although for volume-preserving shear deformations, as considered here, it can still be used.
The above expression, Eq. (7), can be rewritten as where C G = ε 5π κ R0 e −α T Tg(2+λ) is defined as the shear modulus value at the glass transition temperature T g , i.e. C G ≡ G(T g ). The constant ε stems from the integration of α T and from the dimensional prefactor in the powerlaw ansatz for g(r). All the parameters in this expression are either fixed by the experimental/simulation protocol or can be found in the literature. The parameter λ has to be extracted from g(r) data, according to the protocol that we give in the Section VI.

V. ANALYTICAL EXPRESSION FOR THE VISCOSITY
We can now use our model for G(T ) to evaluate the activation energy E(T ) involved in restructuring the glassy cage and, hence, the viscosity η(T ) of the melts. Within the framework of the cooperative shear or elastic model of the glass transition [36][37][38][39], the activation energy for local cooperative rearrangements is E(T ) = GV c . The characteristic atomic volume V c appearing here is accessible through the theoretical fitting of the viscosity data, although its value cannot be arbitrary and it must be representative of the atomic composition of the alloy and of the atomic sizes of its constituents.
Replacing the expression for E(T ) in the Arrhenius relation given by the cooperative shear model of the glass transition, and using Eq. (8) for G(T ) inside E(T ), we obtain the following analytical expression for the viscosity, where η 0 is a normalisation constant set by the high-T limit of η.
It is important to consider how the double-exponential dependence of the viscosity on the temperature arises. The first exponential stems from the elastic activation described in the framework of the cooperative shear model, whereas the second exponential is due to the Debye-Grüneisen thermal expansion rooted in lattice-dynamical considerations of anharmonicity. This formula accounts for both anharmonicity, through α T , and for the repulsion steepness λ (or softness 1/λ).

VI. ESTIMATING λ FROM THE RADIAL DISTRIBUTION FUNCTION OF BINARY ZrCu ALLOYS
We can now apply the tools introduced above to find an analytical connection between the interatomic interaction parameters and the η(T ) and G(T ) of the melt up to T g . To this end, we will apply our model to the binary system of ZrCu alloys. The atomic-level structure of this system is studied by means of numerical simulations, and we present here also ad-hoc experimental data for η(T ) and g(r) for the special case of Cu 50 Zr 50 . Using our analytical model and the g(r) data from simulations, we can build an analytical connection between g(r) and η(T ), which makes use of the parameter λ introduced above. The model will be tested against ad-hoc experimental data of viscosity as a function of temperature for the Cu 50 Zr 50 melt.
We also should notice that the λ value is sensitive to the elemental composition and stoichiometry, but is less sensitive to temperature changes. This fact can be readily understood by considering that the interatomic potential results from the electronic structure of both valence and conduction electrons. In the temperature regime that we consider here and in the comparison with viscosity data below, it is quite unlikely that the electronic and ionic structure change with temperature, and this is reflected in the fact that the slope of the ascending flank of the the first peak of g(r) does not change significantly in the range of T under consideration. What changes with T , is the position of the first peak, of course, because the average distance between ions increases due to thermal expansion. For these reasons, in our comparison for the viscosity, below, we will use the value of λ determined near T g .

A. Numerical simulations of g(r) for ZrCu alloys
In order to extract the λ parameter from sensible data, we carried out Molecular Dynamics (MD) simulations of the Zr 100−x system (where x = 20, 35, 46, 50, 60, 65, 80), by employing an interatomic potential proposed by Duan et al. [40]. This potential is a semi-empirical manybody potential developed in analogy to the tight-binding scheme in the second-moment approximation [41,42]. The equations of motion were integrated by using the Verlet algorithm with a time step of 5 fs. The configurations were prepared starting from a cubic cell box of 1.28 × 10 5 atoms in the B2 structure and periodic boundary conditions in all directions. In each case, the posi-tions of all atoms were redistributed randomly within the simulation cell and the resulting systems were first equilibrated at 300 K in the isothermal-isobaric ensemble (NPT) for 100 ps and subsequently heated up to 2000 K for melting. After sufficient equilibration in the liquid state, the configurations were cooled down to 300 K (always in the NPT ensemble) with a cooling rate of 10 K/ps, and they were finally equilibrated for 100 ps (always in NPT). The structural changes of the system were studied by calculating the total g(r) together with the Faber-Ziman partials every 100 K upon cooling. In Fig.3 we show plots of the radial distribution function g(r) for different compositions of the binary alloy.
B. Fitting protocol of g(r) data Also shown in Fig. 3, are the fittings obtained using the power-law ansatz g(r) ∼ (r − σ) λ . In order to make quantitative fittings, the numerical coefficients in this ansatz need to be specified, so that we write The numerical coefficients g 0 and a take care of two important facts. Firstly, the g(r) is dimensionless by its definition, hence the prefactor g 0 takes care of dimensionality and has dimensions ofÅ −λ , if we express lengths in units ofÅ in the above formula. Secondly, it is not realistic that g(r) = 0, exactly, at r = σ. The probability density of two ions at that distance will be extremely small, but not identically zero. Hence, the parameter b > 0 takes this fact into account. Finally, both parameters g 0 and b control the value of g(r) at the hard-core distance σ = 2a 0 , since g(σ) = g 0 b λ . Hence, the two parameters must satisfy the condition g(σ) = g 0 b λ ≪ 1. In Ref. [8], we found that, in order to simultaneously fit shear modulus and viscosity data of different alloys, values of the exponent λ in the order of magnitude λ ∼ 100, are required. In turn, this constraint on the order of magnitude of λ, puts constraints on the orders of magnitude of the coefficients g 0 and b. Hence, for example, the only acceptable fitting of the g(r) of the Cu 50 Zr 50 alloy with λ ∼ 100, can be achieved with λ = 80. This value, in turn, constrains the values of the other coefficients to be g 0 = 1 × 10 −28 and b = 2.16Å.
Using the same values of the g 0 and b coefficients also for the other compositions for which g(r) was obtained from simulations, we produce the fittings reported in Fig.  3 for different values of stoichiometry.  0 1 . 2 1 . 4 1 . 6 1 . 8 2 . 0 1 . 0 1 . 2 1 . 4 1 . 6 1 . 8  Cu 50 Zr 50 . As the alloy becomes richer in Cu, a trend becomes visible, where the λ value increases as Cu becomes the dominant component. This effect might be explained in terms of electronic structure of the ions, along the lines of [8]. In particular, Zr-Zr interactions may be effectively softer because of the pronounced d-wave character of the outer electron shells of Zr which is a transition element, compared to the effectively more steeply repulsive Cu-Cu interaction, dominated by the more pronounced s-wave character of the Cu ions. Based on this picture, the Zr-Cu interaction may be comparatively the most steeply repulsive, as the smaller Cu atoms tend to nest in the corners, in between the lobes of the d-wave structure of Zr.
We also see from Fig. 3, that for the alloy richest in Zr, Cu 20 Zr 80 , two slopes in the repulsive ascending flank of the first peak are visible. In order to obtain more insight into this phenomenon, we have studied also the Faber-Ziman partials g αβ (r), where α, β = Zr, Cu, for this particular composition. The results are shown in Fig. 4. It is seen that the left-most slope is due to the Cu-Cu contribution, which is steeper and shorter-ranged and thus contributes a higher value to the effective λ which is larger than the partial contribution due to the Zr-Zr contribution which is softer, as argued above because of the more pronounced d-wave character of the Zr ion. The lower slope of the second flank is in fact due to the Zr-Zr contribution. From this analysis, it appears that the short-range repulsive part of g(r) is dominated by the Zr-Cu contribution which is responsible for the λ = 80 overall value, which varies only weakly with composition. This can be explained if one considers that smaller Cu atoms tend to surround larger Zr atoms and to nest into the corners of the d-wave outer shell of Zr. We also note here that similar electronic interactions affect the slowing down of the dynamics in Zr-(Cu/Ni/Co) melts with Al additions via an enhanced short-range packing between the Al and late-transition metal species [43][44][45].
D. The interatomic repulsion parameter λ is independent of T In order to explore a possible dependence of the interatomic repulsion steepness parameter λ on the temperature T , we performed numerical simulations on the same alloy melts at different temperatures. To this end, we have determined g(r) from simulations of the Cu 50 Zr 50 melt at T = 600 K and T = 800 K, in addition to the T = 700 K data reported in Fig. 3.
The comparison is shown in Fig. 5. As is evident, the only difference that comes from varying T is the decrease in the height of the first peak of g(r) upon increasing T . This effect is trivial and is attributed to the less pronounced smearing of atomic correlations as T is increased [46]. Moreover, this observation clearly supports our main assumption that the connectivity changes with T , a manifestation of which is the fact that the parameter Z in our model decreases upon increasing T . Thermal expansion greatly enhances this effect in experimental systems, as volume is not constant but is allowed to change following a change of temperature. It is also evident, that T has no effect whatsoever on the slope of the repulsive left-hand-side flank of g(r). The conclusion of this analysis is, therefore, that our repulsion steepness parameter λ, which is directly related to the slope of the left flank of g(r), is independent of T . In the following steps of our model we will thus take λ as independent of T .

VII. COMPARISON WITH EXPERIMENTAL DATA OF VISCOSITY VERSUS T
We are now able to evaluate our analytical model linking the interatomic potential and the g(r) with the macroscopic viscosity of the alloys upon approaching T g from the high-T end. Since, unfortunately, experimental viscosity data over the wide range of CuZr compositions modeled here are not readily available, we use experimental data for the Cu 50 Zr 50 composition taken from ad-hoc measurements as detailed in the Appendix.
The value of the interatomic repulsion steepness λ = 80.0 for this system was obtained from the fitting of g(r) in Fig. 3. Using Eq. (9) it is therefore possible to obtain a one-parameter fitting of the experimental viscosity data for the Cu 50 Zr 50 alloy as a function of T , as shown in Fig. 7. Furthermore, we recall from Eq. (8) that C G is defined as the value of the high-frequency shear modulus G at T = T g . For the Cu 50 Zr 50 alloy, the value of C G = 31.3 GPa was obtained from experimental measurements in the literature [47]. The only adjustable parameter in Eq.(9) is thus the characteristic atomic volume V c . From the fitting we obtain V c = 0.0039 nm 3 .
In Fig. 7 (left panel) we also reported a fitting done with the VFT-type model of Ref. [4], which uses the concept of kinetic fragility D * to provide a more physical interpretation to the VFT relation (which normally has three adjustable parameters). The VFT-type relationship of Ref. [4] is given by η(T ) = a exp[D * T 0 /(T − T 0 )], where a and T 0 are free parameters, while D * was ad- justed to the experimental Cu 50 Zr 50 data in Ref. [4] and interpreted as the kinetic fragility. In the fitting shown in Fig. 7, T 0 = 578.6 K was used to obtain the best fitting, which is very far from the T g value of this system, T g = 770 K as determined from our simulations and experiments for the Cu 50 Zr 50 system. Therefore, the VFT-type fitting of Ref. [4] uses two free parameters, D * and T 0 (without considering the parameter a which plays the same role as η 0 in our Eq. (9) and is set by the high-T behaviour of viscosity). In contrast, our Eq. (9) has only one free parameter, V c . In spite of having one adjustable parameter more than our model, the fitting of Ref. [4], provides a less accurate fitting in comparison, especially in the regime T > 1200 K. The characteristic volume is defined within the framework of the cooperative shear model as V c = (2/3)(∆V ) 2 /V and was derived by Dyre [39] using linear elasticity for an expanding sphere. The quantity ∆V is the activation volume or the local volume change associated with a cooperative shear event and can be expressed as ∆V ≈ √ 480V c . Furthermore, V in these relations can be identified with the volume of a shear-transformation zone (STZ), and in Ref. [8] it was found that, for various alloys, the approximate relationship V ≈ 320V c holds.
The V c = 0.0039 nm 3 value obtained here for Cu 50 Zr 50 is on the same order of magnitude but smaller than the value obtained in our previous analysis [8] for Zr 41.2 T i 13.8 N i 10 Cu 12.5 Be 22.5 , for which a characteristic volume of 0.0085 nm 3 was found, which compares well with independent determinations on similar alloys [8]. Moreover, for Cu 50 Zr 50 , we obtain an activation volume of ∆V ∼ 100Å 3 , which is somewhat larger than ∆V = 84.3Å 3 , which was determined in MD simulations on a similar binary Cu 56 Zr 44 alloy [48].
Finally, in Fig. 7 (right panel) we have also shown the underlying effective interatomic potential corresponding to the repulsion steepness value λ = 80 obtained from the analysis of the g(r), and its mapping to the Ashcroft-Born-Mayer potential in Eq.1. Some parameters used in the mapping are fixed and independent of the alloy composition, such as q T F = 1.7Å −1 , from Ref. [22], and C = 2.8Å −1 , perfectly within the range 1.89−4.72Å −1 as discussed in Ref. [8]. The parameters A and B, instead, set the characteristic energy scales of the Ashcroft and Born-Mayer terms, respectively, in Eq. (1). The values that we found here for Cu 50 Zr 50 are quite close to the values found in Ref. [8] for other metallic glasses, which reflects the robustness of this approach.

VIII. CONCLUSION
In summary, we developed a protocol which takes the short-range repulsive part of g(r) from simulation (or experimental) data as the input to qualitatively infer the viscosity of metallic alloys as a function of temperature upon vitrification, and their fragility. Using a combination of simulation and experimental data for the case of the binary Cu-Zr alloys, it has been possible to produce a one-parameter theoretical fit of the viscosity as a function of temperature for the Cu 50 Zr 50 system. Importantly, the steepness of the viscosity rise upon approaching the glass transition, is controlled by two averaged interaction parameters: the λ interatomic-repulsion parameter introduced in Ref. [8], and the thermal expansion coefficient α T . This result confirms the picture according to which high thermal expansion coefficients favour fragile behaviour, alongside with steep short-ranged interatomic repulsion. The repulsion steepness is related to shallow changes at the level of nearest-neighbours upon decreasing T , whereas thermal expansion, associated with the long-range part of the interatomic potential, in metallic alloys is related to medium-range order effects and to the ability of the system to expel free volume at the level of the 3rd-4th coordination shells [10].
In contrast to VFT-type relations for the viscosity as a function of T , our analytical theory is derived from the underlying atomic dynamics as shown in Ref. [8], and gives the viscosity as a double-exponential decreasing function of T . The outer exponential comes from the shoving model of the viscosity derived by Dyre, while the inner exponential is due to the exponential decrease of interatomic connectivity with increasing T due to thermal expansion [30]. By estimating the repulsion parameter λ from numerical simulations of g(r) (validated at higher T against experimental data), we have shown here how a quantitative fitting of viscosity as a function of T can be obtained with just one fitting parameter, which is the characteristic atomic volume in the shoving model. Also, all parameters involved in our viscosity expression have a clear microscopic physical meaning.  9), with the calculation parameters λ = 80 determined from the g(r) fitting and CG = 31.3 GPa from Ref. [47]. Tg = 770 K is the glass transition temperature for the Cu50Zr50 system and is therefore not a fitting parameter. The only adjustable parameter is, hence, the characteristic atomic volume of the cooperative shear model Vc = 0.0046 nm 3 , which turns out to be a meaningful value in comparison with previous estimates for similar alloys [8]. The dashed line is a fitting made with the VFT model η(T ) = a exp[D * T0/(T − T0)] of Ref. [4]. Here three fitting parameters are used: a = 0.223, T0 = 578.6 and the kinetic fragility D * = 4.8. (b) The logarithmic potential of mean force from Boltzmann inversion of the repulsive flank of the g(r) first peak, dashed line, and its fitting (solid line) using the Ashcroft-Born-Mayer interatomic potential, Eqs. (1)-(2), with parameters A = 0.121 eV, qT F = 1.7Å −1 , as in Ref. [8], B = 5.61 eV, C = 2.8Å −1 , as in Ref. [8], andσ = 1.35Å from the average of ionic core sizes for Cu50Zr50.
Appendix A: Scaling between the connectivity Z and the global interatomic repulsion steepness λ When increasing T , the average spacing between atoms in the coordination shell becomes larger, and the probability of nearest neighbours leaving the connectivity shell increases. It is then possible to use the radial distribution function g(r) to relate the change in atomic packing fraction φ, due to an externally imposed change in temperature T , to the change in connectivity Z. Following along the lines of Ref. [30], the connectivity can be written as usual as an integral over the g(r), typically written as Z ≈ rmax 0 r 2 g(r)dr up to a threshold which could be perhaps set conservatively as the maximum of g(r). In amorphous systems, the peak of g(r) increases upon increasing the packing fraction (and upon decreasing T as is visible in Fig.5) due to many-body excluded-volume interactions, and this effect translates into a larger value of Z. A way to describe this effect, is to shift the upper limit of the integral which gives Z in proportion to φ, as r max = cφ, where c is some proportionality constant, and one can write Z ≈ cφ 0 r 2 g(r)dr, where r represents the separation distance between two ions in the system. This relation correctly recovers the limit Z → 0 in the limit φ → 0. From Fig.5 it is also seen that the ascending part of the first peak of g(r) happens within a narrow r-interval, which is quite close to the metallic diameter of the system, ≈ 2.2Å. At this very short separation between the ions, one can approximate the local geometry as being effectively Cartesian, instead of spherical, which removes the metric factor from Eq. (A1). This is a customary simplification in dealing with small gaps between two spheres, and is widely in colloidal systems, where it is known as the Derjaguin approximation [49]. Upon replacing our approximation for the ascending part of g(r) Eq. (9) we thus obtain Z ≈ cφ 0 (r − σ + b) λ dr. (A2) Next, we make the standard change of variable in the integral x = r − σ + b, and the integral becomes where q is some coefficient. We then realize that the following approximation can be used Hence we obtain the key asymptotic scaling used in our model to link the g(r) with the high-frequency shear modulus G and then with the viscosity η.
Appendix B: Experimental measurement of g(r) and viscosity for the Cu50Zr50 system The g(r) of the Cu 50 Zr 50 melt was determined previously in combined neutron and synchrotron X-ray diffraction experiments, described in detail in Ref. [50] . From the total neutron and X-ray static structure factors S(q), two of the three Bhatia-Thornton partial structure factors [51], S nn (q) and S nc (q), were able to be determined with good precision. The corresponding g nn (r) and g nc (r) were calculated by Fourier transformation of S nn (q) and S nc (q), respectively.
The viscosity values were determined in both the equilibrium and undercooled states of the Cu 50 Zr 50 melt at the Institute for Materials Physics in Space of the German Aerospace Center (DLR) in Cologne, Germany. An electrostatic levitation (ESL) apparatus under ultra-high vacuum conditions of ≈ 10 −8 mbar was used to levitate an electrically charged sample with a mass of 50 mg. Heating was achieved via two 25 W IR lasers and the temperature measured contact-free by a pyrometer directed at the sample side. A high-speed camera allowed determination of the vertical sample radius R z (t) as a function of time t. The levitated sample was brought in to the liquid state and heated to an initial temperature of 1308 K, which is 100 K higher than the liquids temperature T liq = 1208 K, taken from the published phase diagram [52] . Measurements of the viscosity were carried out isothermally during cooling using the oscillating droplet method [53] . The viscosity was determined at each temperature by measuring the decay of surface oscillations induced by a sinusoidal electric field described by R z (t) = R 0 + A exp(−t/τ 0 ) sin(ωt + δ 0 ), where R 0 is the quiescent sample radius, A the oscillation amplitude, τ 0 the decay time constant, ω the frequency and δ 0 the phase shift. Using Lamb's law [54] , the viscosity was calculated as where ρ is the macroscopic density of the droplet, also determined in these ESL experiments combined with video diagnostic techniques [53,55].