Effectively one-dimensional phase diagram of CuZr liquids and glasses
EEffectively one-dimensional phase diagram of CuZr liquids and glasses
Laura Friedeheim, Jeppe C. Dyre, ∗ and Nicholas P. Bailey † “Glass and Time”, IMFUFA, Dept. of Science and Environment,Roskilde University, P. O. Box 260, DK-4000 Roskilde, Denmark (Dated: March 1, 2021)This paper presents computer simulations of Cu x Zr − x ( x = 36 , ,
64) in the liquid and glassphases. The simulations are based on the effective-medium theory (EMT) potentials. We find goodinvariance of both structure and dynamics in reduced units along the isomorphs of the systems. Thestate points studied involve a density variation of almost a factor of two and temperatures goingfrom 1500 K to above 4000 K for the liquids and from 500 K to above 1500 K for the glasses. Forcomparison, results are presented also for similar temperature variations along isochores, showinglittle invariance. In general for a binary system the phase diagram has three axes: composition,temperature and pressure (or density). When isomorphs are present, there are effectively only twoaxes, and for a fixed composition just one. We conclude that the thermodynamic phase diagramof this metallic glass former for a fixed composition is effectively one-dimensional in the sense thatmany physical properties are invariant along the same curves, implying that in order to investigatethe phase diagram, it is only necessary to go across these curves. ∗ [email protected] † [email protected] a r X i v : . [ c ond - m a t . m t r l - s c i ] F e b I. INTRODUCTION
Metallic systems constitute a very important category of glass formers due to their potential applications, as wellas their suitability as model systems for studies of the glass transition in computer simulations [1–5]. A well-studiedexample is the CuZr system, which is a fairly good glass former despite consisting of just two elements [3, 6]. Thispaper presents numerical evidence that both above and below the glass transition, CuZr systems are simpler thanhas hitherto been recognized. Specifically, for three different compositions of the CuZr system we show that curvesexist in the thermodynamic phase diagram along which the atomic structure and dynamics are invariant to a goodapproximation. The implication is that the two-dimensional thermodynamic phase diagram becomes effectively one-dimensional in regard to many material properties.The background of the investigation is the following. In liquid-state theory, a simple liquid is traditionally defined asa single-component system of particles described by classical Newtonian mechanics and interacting by pair-potentialforces [7–12]. It has been known for more than half a century that the hard-sphere (HS) model reproduces well thephysics of many simple liquids, both in regard to the radial distribution function (RDF) and to dynamic propertiessuch as the viscosity or the diffusion coefficient [12–18]. The traditional explanation of this is the “van der Waalspicture” according to which the repulsive forces dominate the physics of simple liquids [12, 15–17, 19].The HS model is a caricature simple liquid with pair forces that are zero except right at the particle collisions.In the HS model, temperature plays only the trivial role of determining particle velocities and thus the time scale;temperature is entirely unrelated to the geometry of particle positions. This implies that the thermodynamic phasediagram of the HS system is effectively one-dimensional with density as the only non-trivial variable: the dynamics oftwo different HS systems with the same packing fraction but different temperatures are identical, except for a trivialuniform scaling of the space and time coordinates. As a consequence, scaled RDFs are identical, scaled mean-squaredisplacements are identical, viscosities are trivially related, etc.The mapping of a simple liquid to a HS system presents the issue of identifying the effective HS packing fraction ata given thermodynamic state point of the liquid. Many suggestions have been made for how to calculate the relevanthard-sphere radius, yet no consensus was ever arrived at [20–28]. Already in 1977, Rosenfeld suggested a brilliantthermodynamics-based alternative by basically reasoning as follows [29]: Since the HS packing fraction determinesthe configurational part of the entropy, this quantity provides the required mapping between a simple liquid and itscorresponding HS system. Defining the excess entropy S ex as the system’s entropy minus that of an ideal gas atthe same density and temperature [12, 29], S ex quantifies the configurational entropy (note that S ex < T and theparticle number density ρ [31–36]. When density scaling is applied to experimental data, if γ is the so-called density-scaling exponent, plotting data for the dynamics as a function of ρ γ /T results in a collapse [32, 33, 35, 36]. This meansthat the dynamics depends on the two variables of the thermodynamic phase diagram only via the single number ρ γ /T . It should be emphasized that density scaling is not universally applicable; for instance, it works better for vander Waals liquids and metals than for hydrogen-bonded liquids [33, 36, 37]. An important extension of density scalingwas the discovery of isochronal superposition, according to which not only the average relaxation time is invariantalong the curves of constant ρ γ /T , so are the frequency-dependent response functions [38–40]. This indicates that theway atoms or molecules move about each other is identical at state points with same the same value of ρ γ /T . Onemay think of this as a “same-movie” property: Filming the atoms/molecules at two such state points results in thesame movie except for a uniform scaling of all particle positions and of the time.The above-mentioned findings can all be derived from the hidden-scale-invariance property stating that the orderingof a system’s configurations R ≡ ( r , ..., r N ) (in which N is the number of particles and r i is the position of particle i )according to their potential energy, U ( R ), at one density is maintained if the configurations are scaled uniformly to adifferent density [41]. The formal mathematical definition of hidden scale invariance is the following logical implication U ( R a ) < U ( R b ) ⇒ U ( λ R a ) < U ( λ R b ) . (1)Equation (1) implies that structure and dynamics, when given in proper reduced units, are invariant along the curvesof constant excess entropy, the system’s so-called isomorphs [30, 41, 42]. This result is rigorous if Eq. (1) applieswithout exception, but this is never the case for realistic models. However, isomorph invariance is still a goodapproximation if Eq. (1) applies for most of the physically important configurations and for scaling parameters λ relatively close to unity. This is believed to be the case for many metals and van der Waals bonded systems, whereassystems with strong directional bonds like hydrogen-bonded and covalently bonded systems are not expected to obeyisomorph-theory predictions [43]. For metals, however, the existence of isomorphs has only been validated in a fewcases [44, 45].In isomorph theory the density-scaling exponent γ is generally state-point dependent. This has recently beenconfirmed in high-pressure experimental data [46–48]. Systems with hidden scale invariance are referred to as R-simple in order to indicate the simplification of the physics that follows from this symmetry; most, though not all,pair-potential systems are R-simple and several molecular systems are also R-simple, necessitating a specific name forthis class of systems.The purpose of the present paper is to check for isomorphs in a typical metallic glass former. For this we havechosen to study three different CuZr mixtures. The systems have been computer simulated both in the liquid and glassphases, using the effective-medium theory (EMT) interaction potential [49–51]. We find good isomorph invarianceof structure and dynamics involving density changes up to a factor of two. This implies a significant simplificationin the description of the physics of this metallic glass former since the thermodynamic phase of diagram of CuZr iseffectively one-dimensional. II. THE EFFECTIVE-MEDIUM-THEORY POTENTIAL
The EMT potential [49–51] is one of several similar potentials aimed at describing metals with an accuracy compa-rable to that of a full density-functional-theory (DFT) treatment, but at a much lower computational cost. A widelyused class of potentials in this group of mean-field potentials is the embedded atom method (EAM) [52, 53]. EMTand EAM both write the total energy E as a pair-potential part plus a function of the local electron density at eachparticle. The EMT realizes this in a semi-empirical way, whereas the parameters of the EAM are determined byfitting to experimental properties of the bulk solid. For more on the relation between the EAM and EMT potentials,the reader is referred to Refs. 49 and 50, while Ref. 51 gives a detailed derivation of the EMT potential and itsparameters (see also the Appendix). A great advantage of the EMT is that the mathematical expression for theenergy is relatively simple. This made it straightforward to implement EMT in our GPU-code RUMD [54], whereasthe EAM typically involves tabulated data that are not easily implemented efficiently in GPU computing.The core of the EMT potential is a well chosen reference system defining the effective medium . The total energy ofthe system is the energy of the reference system plus the difference to the real system. Thus the EMT total energyis written E = (cid:88) i E c,i + (cid:32) E − (cid:88) i E c,i (cid:33) (2)where E c,i is the so-called cohesive energy, which is the energy of atom i in the reference system. The idea is now thatthe difference term should be small enough to be treated accurately by first-order perturbation theory. To obtain thisthe reference system must be as close as possible to the real system.The real and the reference systems are linked by a “tuning parameter”. In the first version of the EMT potential,the homogeneous electron gas was used as the reference system, with the electron density as the tuning parameter[55, 56]. The EMT version used in this paper is that of Ref. 51 for which a perfect fcc crystal is the reference system.Here, the lattice constant serves as the tuning parameter, i.e., is used to adjust the environment of an atom such thatthe average electron density surrounding the atom matches that of the real system.The EMT potential of a pure metal involves seven parameters: the negative cohesive energy E , a charge densityparameter n , the Wigner-Seitz radius s (defined in terms of the atomic, not the electronic density), a parameter η quantifying the influence of the density tail surrounding an neighboring atom, a parameter λ determined fromthe bulk modulus, and finally the product V ( βη − κ ) determined from the shear modulus. The density-relatedparameters n and η were originally calculated self-consistently by reference to the homogeneous electron gas, whilethe five other parameters are determined from experimental or ab initio data. More details on how the parameters aredetermined for single and compound systems can be found in the Appendix and in Refs. 55 and 57. The parametervalues for CuZr used in this work are those of Ref. 57, where parameters were adjusted to match DFT-determinedproperties—cohesive energies, lattice constants and elastic constants-of both the pure metals and a Cu Density [ Å -3 ] T e m p e r a t u r e [ K ] Cu Zr Cu Zr Cu Zr liquid Density [ Å -3 ] T e m p e r a t u r e [ K ] Cu Zr Cu Zr Cu Zr glass FIG. 1. Logarithmic density-temperature phase diagrams showing all state points simulated along isochores (vertical lines)and isomorphs (lines at an angle). The colors reflect the three different compositions studied. The left figure gives liquid statepoints, the right figure gives glass state points. Each isomorph is generated by means of the direct-isomorph-check method(see the text), proceeding in steps of 5% density changes starting from “reference” state points at temperature 1500 K for theliquids and 500 K for the glasses (full symbols). The reference state point densities were selected to have approximately thesame pressure ( ∼ III. REDUCED QUANTITIES
Structure and dynamics are invariant along isomorphs only when these are given in so-called reduced versionsbased on a “macroscopic” unit system that depends on the state point in question. The unit system defining reducedvariables reflects the system’s volume V and temperature T as follows. If the particle number density is ρ ≡ N/V ,the length, energy, and time units are, respectively, [42] l = ρ − / , e = k B T , t = ρ − / (cid:114) mk B T . (3)Here m is the average particle mass. Equation (3) refers to Newtonian dynamics; for Brownian dynamics one usesthe same length and energy units, but a different time unit [42]. All physical quantities can be made dimensionlessby reference to the above units. “Reduced” quantities are denoted by a tilde, for instance˜ R ≡ ρ / R . (4) IV. TRACING OUT ISOMORPHIC STATE POINTS
Three compositions were studied, Cu Zr , Cu Zr , and Cu Zr . Figure 1 presents the state points simulatedin a density-temperature thermodynamic phase diagram. Isomorph invariance is never perfect in realistic systems.In order to estimate to which degree this invariance holds, it is therefore useful to compare structure and dynamicsvariations along isomorphs to what happens along curves of similar temperature or density variation, which are notisomorphs. We have chosen to compare to isochores (lines of constant density) with the same temperature variationas the isomorphs. The isochores, of course, are the vertical straight lines in Fig. 1, the isomorphs are the lines witha slope. The high-temperature state points describe equilibrium liquids, the low-temperature points are glass-phasestate points.We now turn to the challenge of tracing out isomorphs. Recall that an isomorph is a curve of constant S ex for asystem that obeys the hidden-scale-invariance condition Eq. (1) at the relevant state points. To which degree thiscondition is obeyed may be difficult to judge because Eq. (1) always applies when λ is close to unity, but fortunatelya practical criterion exists: Eq. (1) applies to a good approximation if and only if the virial W and potential energy U are strongly correlated in their thermal-equilibrium constant-density ( N V T ) fluctuations [41]. These fluctuations arecharacterized by the Pearson correlation coefficient R defined by (in which sharp brackets denote canonical-ensembleaverages and ∆ is the deviation from the thermal average) R = (cid:104) ∆ U ∆ W (cid:105) (cid:112) (cid:104) (∆ U ) (cid:105)(cid:104) (∆ W ) (cid:105) . (5)As a pragmatic criterion, R > . R goes below 0.9 at high densities in the liquid phase, as wellas in most of the glass phase (Fig. 2), but at most state points studied R is above 0.8. Thus it makes good sense totest for isomorph invariance.Tracing out a curve of constants excess entropy is straightforward if one knows how S ex varies throughout the phasediagram. It is a bit challenging to evaluate entropy, however, because this involves thermodynamic integration (orthe Widom insertion method that is also tedious). In order to trace out an isomorph, one does not need to know thevalue of S ex , however, and one can therefore make use of the following general identity γ ≡ (cid:18) ∂ ln T∂ ln n (cid:19) S ex = (cid:104) ∆ U ∆ W (cid:105)(cid:104) (∆ U ) (cid:105) . (6)Here the quantity γ is the (state-point dependent) density-scaling exponent defined as the isomorph slope in alogarithmic density-temperature phase diagram like that of Fig. 1. The second equality sign is a statistical-mechanicalidentity that allows for calculating γ from N V T equilibrium fluctuations at the state point in question [42]. Figure 2shows how γ varies along the isomorphs of the CuZr systems studied below, plotted as a function of the densityrelative to that of the isomorph reference state point. All cases show similar behavior with γ decreasing significantlywith increasing density. This indicates a softening of the interactions at high densities.Equation (6) can be used to trace out an isomorph by numerical integration, for instance using the Euler algorithmfor small density changes of order one percent [42] or using the fourth-order Runge-Kutta algorithm that allows forsignificantly larger density changes [59]. Both methods are fool proof, but involve many simulations if one wishes tocover a significant density range. There are alternative “quick and dirty” methods which, depending on the systemin question, can be quite useful. For instance, isomorphs of the Lennard-Jones system are to a good approximationgiven by h ( ρ ) /T =Const. in which h ( ρ ) = ( γ / − ρ/ρ ) − ( γ / − ρ/ρ ) in which γ is the density-scalingexponent at a selected “reference state point” of density ρ [60, 61].A general and fairly efficient method for tracing out isomorphs is the “direct isomorph check” (DIC) [42], and weused this for generating the CuZr isomorphs. The DIC is justified as follows [41]. Hidden scale invariance (Eq. (1))implies that the microscopic excess entropy function S ex ( R ) is scale invariant, i.e., a function only of a configuration’sreduced coordinate ˜ R : S ex ( R ) = S ex ( ˜ R ) [41]. From the definition of S ex ( R ) it follows that U ( R ) = U ( ρ, S ex ( ˜ R )) inwhich the function U ( ρ, S ex ) is the average potential energy at the state point with density ρ and excess entropy S ex [41]. Considering configurations with the same density ρ and small deviations in the microscopic excess entropy fromthat of the givne state point S ex , an expansion to first order leads to U ( R ) ∼ = U ( ρ, S ex ) + T ( ρ, S ex ) (cid:16) S ex ( ˜ R ) − S ex (cid:17) . (7)Consider two state points ( ρ , T ) and ( ρ , T ) with the same excess entropy S ex and average potential energies U and U , respectively. If R and R are configurations of these state points with the same reduced coordinates, i.e.,obeying ρ / R = ρ / R ≡ ˜ R , (8)one gets by elimination of the common factor S ex ( ˜ R ) − S ex in Eq. (7) with T ≡ T ( ρ , S ex ) and T ≡ T ( ρ , S ex ) U ( R ) − U T ∼ = U ( R ) − U T . (9) γ Cu Zr Cu Zr Cu Zr liquid relative density ρ/ρ . . R γ Cu Zr Cu Zr Cu Zr glass relative density ρ/ρ . . R FIG. 2. Virial potential-energy Pearson correlation coefficient R (Eq. (5)) and density-scaling exponent γ (Eq. (6)), plottedalong the isomorphs as a function of the density relative to that of the reference state point (which has temperature 1500 Kfor liquids and 500 K for glasses). The left figure is the liquid, the right figure is the glass. We see similar pictures in thetwo cases, with γ decreasing significantly as density is increased, indicating an effective softening of the interactions. Thevirial potential-energy correlations are generally strong, with a maximum at densities close to the reference state point densitiesdenoted by ρ . The dashed lines mark R = 0 .
9, which is traditionally used for delimiting state points for which isomorph-theorypredictions are expected to apply [42, 58]. -5.09 -5.08 -5.07 -5.06
U(R) [eV] - . - . U ( ( ρ / ρ ) / R ) [ e V ] FIG. 3. Example of a direct isomorph check (DIC), here with state point 1 being the reference state point of the 64-36 mixture( T = 1500K) and state point 2 having a 5% higher density. At state point 1, a series of thermal-equilibrium configurations R are sampled. Each of these are scaled uniformly by the factor ( ρ /ρ ) − / = 0 . T /T , compareEq. (10); this determines the temperature T that makes the state point ( ρ , T ) isomorphic to state point ( ρ , T ). While not of direct relevance for the present paper, we note that Eq. (9) implies exp( − U ( R ) /k B T ) ∝ exp( − U ( R ) /k B T ),i.e., that the two configurations have the same canonical probability. This is a manifestation of the hidden scaleinvariance inherent in isomorph theory [42]. Equation (9) leads to U ( R ) ∼ = ( T /T ) U ( R ) + ( U − ( T /T ) U ). Forthe fluctuations about the respective mean values this implies∆ U ( R ) ∼ = T T ∆ U ( R ) . (10)Equation (10) implies that isomorphic state points may be identified as follows: First sample a set of equilibriumconfigurations at the state point ( ρ , T ). Then scale these configurations uniformly to density ρ . The temperature T of the state point with density ρ , which is isomorphic to state point ( ρ , T ), is now found from the slope of ascatter plot of the potential energies of scaled versus unscaled configurations. An example of how this works is shownin Fig. 3.Because the hidden-scale-invariance property is not exact, the DIC is less reliable for large density changes than forsmaller ones. We traced out the isomorphs studied below using step-by-step DICs involving density changes of 5%.The resulting isomorphs are shown in Fig. 1. The simulated isomorphic state points are listed in Table I (liquid) andTable II (glass). T [K] ρ [˚A − ] p [GPa] T [K] ρ [˚A − ] P [GPa] T [K] ρ [˚A − ] P [GPa]875 0.0508 1.7 870 0.0556 0.8 950 0.0620 1.51219 0.0551 9.6 1217 0.0602 9.0 1215 0.0659 8.4 Zr , Cu Zr , Cu Zr . The boldface row represents the reference state point for the isomorph of eachcomposition. T [K] ρ [˚A − ] p [GPa] T [K] ρ [˚A − ] P [GPa] T [K] ρ [˚A − ] P [GPa]339 0.0529 0.4 398 0.0602 3.1 400 0.0659 1.7432 0.0562 6.2 466 0.0627 7.7 467 0.0686 6.6
500 0.0585 11.0 500 0.0640 10.3 500 0.0700 9.5
587 0.0614 17.6 575 0.0672 17.3 578 0.0735 17.1683 0.0645 25.5 672 0.0706 25.9 667 0.0772 26.1781 0.0677 34.6 773 0.0741 35.8 758 0.0810 36.8880 0.0711 45.2 877 0.0778 47.5 858 0.0851 49.4980 0.0747 57.6 989 0.0817 60.7 968 0.0893 64.11088 0.0784 71.5 1103 0.0858 76.4 1083 0.0938 81.01200 0.0823 87.6 1218 0.0901 94.1 1203 0.0985 100.61304 0.0864 105.9 1340 0.0946 114.3 1325 0.1034 123.41442 0.0908 126.5 1440 0.0993 137.7 1453 0.1086 149.11558 0.0953 149.7 1575 0.1042 163.5 1596 0.1140 178.7TABLE II. Temperature, density, and pressure of glass isomorph state points. The vertical lines divide the compositions shownin the following order: Cu Zr , Cu Zr , Cu Zr . The boldface row represents the reference state point for the isomorphof each composition. The starting point for the glass isomorph is at the same density as the liquid isomorph starting point,but with temperature 500K. V. SIMULATION DETAILS
The three compositions studied in this work are Cu x Zr − x ( x = 36 , , ρ , but for each isomorph we also generated two isomorph state points at lowerdensities to ensure that samples close to zero pressure were included in the study (compare Tables I and II).For each state point on an isomorph, a state point was simulated at the same temperature at the reference-state-point density; these constitute the isochoric state points discussed below along with the isomorph state points. Forall three compositions, the glass-phase reference state points were obtained by cooling at a constant rate in 100000time steps from the liquid reference state point at 1500 K to the glass isomorph reference temperature 500 K. This isa high cooling rate, corresponding roughly to 1 K/ps.From these three glass reference state points, isomorphs were generated by the DIC method in the same way as forthe liquids.The simulations were carried out in RUMD [54], Roskilde university’s GPU Molecular Dynamics package that is FIG. 4. Snapshots of the glass configurations at the glass isomorph reference state points at 500 K. The Cu atoms are orange,the Zr atoms are grey. From left the figures are for the 36-64, 50-50, and 64-36 CuZr compositions. optimized for small systems. The
N V T ensemble was used to simulate cubic boxes containing 1000 particles. Theinitial configuration was a simple cubic crystal with particle types assigned randomly at the required ratios. At eachstate point of the liquid, 10 MD steps of equilibration were performed to melt and equilibrate the liquid. Followingthis we carried out 10 MD steps of the production run. For the glass-phase simulations, 10 MD steps of equilibrationwere performed before 10 MD steps of the production run. The time step in the simulations was 0.5 ˚A / (u / eV) / in which u is the atomic mass.Figure 4 shows the glasses prepared by cooling the liquids to the glass-isomorph reference state points. There areno signs of crystallization. VI. STRUCTURE AND DYNAMICS IN THE LIQUID PHASE
To investigate how the structure varies along isomorphs and isochores for the three CuZr compositions we probedthe radial distribution functions (RDFs), which in reduced units are predicted to be isomorph invariant. There arethree different RDFs, one for Cu-Cu, one for Cu-Zr, and one for Zr-Zr. Plotting a RDF in reduced units impliesscaling the distance variable according to the density (compare Eq. (4)). This results in peaks at roughly the sameplaces for all compositions because the scaling corresponds to taking the system to unit density.Figure 5 shows the reduced RDFs along the isomorphs and Fig. 6 shows the similar RDFs along isochores withsame temperature variation (compare Fig. 1 showing the similar state points). Comparing the two, we conclude thatthe structure is isomorph invariant to a good approximation, but varies significantly along the isochores. Deviationsare largest for the minority-minority RDFs for the non-equimolar compositions. Deviations from isomorph invarianceare also seen in some cases at the first maximum, where the maximum is generally lowered somewhat with increasingdensity. This is an effect that is well understood for pair-particle systems, for which it derives from the fact thata higher density-scaling exponent γ implies a steeper effective pair potential and therefore less likely particle closeencounters. This results in moving some of the low distance RDF to higher distances when γ is large, which is the caseat low densities. This explanation suggests that γ of the present non-pair-potential simulations can also be interpretedas an effective pair IPL exponent [44]. – For the isochores, there is a general “damping” of the RDFs at all distancesas temperature increases. This reflects the stronger thermal fluctuations at high temperatures.Focusing on the height of the first RDF peak, Fig. 7 shows these for all the data of Fig. 5 and Fig. 6. The fullsymbols are the peak heights along the three isomorphs, whereas the open symbols are the peak heights along thecorresponding isochores. Clearly, the variation is significantly larger along the isochores.Next we investigated the dynamics. Figure 8 shows results for the reduced-unit mean-square displacement (MSD)along isomorphs and isochores, respectively (two left columns). We focus on the majority atom MSD, but foundthat data for the minority atom are entirely similar (not shown). Clearly, the MSD is isomorph invariant and variessignificantly along the isochores. Note that the short-time ballistic-region collapse seen in all cases follows from thedefinition of reduced units, i.e., this collapse applies throughout the phase diagram of any system. The two right-handcolumns of Fig. 8 show similar data for the incoherent intermediate scattering function (evaluated at the wave vector2 πρ / that is constant in reduced units). Again, isomorph invariance is clearly demonstrated. C u - C u R D F Isomorph (liquid) r ~ Z r- Z r R D F C u - Z r R D F Cu Zr C u - C u R D F Isomorph (liquid) r ~ Z r- Z r R D F C u - Z r R D F Cu Zr C u - C u R D F Isomorph (liquid) r ~ Z r- Z r R D F C u - Z r R D F Cu Zr FIG. 5. Isomorph liquid-state radial distribution functions (RDFs), plotted as a function of the reduced pair distance ˜ r ≡ ρ / r .Each column gives reduced-unit RDFs along an isomorph of one composition. The color coding used is the xmgrace defaultordering with black for the data set corresponding to the reference state point, then at increasing density: red, green, blue,yellow, etc; orange is the color given to the 11th data set which is here that the highest temperature. The two state pointswith lower density than that of the reference state point are purple and brown. Generally, we see good approximate isomorphinvariance, with some deviation at the first peak maximum and the largest overall deviations at the lowest densities (at whichthe virial potential-energy correlation coefficient drops rapidly, compare Fig. 2). VII. STRUCTURE AND DYNAMICS IN THE GLASS PHASE
The above investigation was repeated in the glass phase of the three mixtures. Isomorph theory is traditionallyformulated by reference to thermal equilibrium [41–43], but we ignored this and proceeded pragmatically as if aglass were an equilibrium system. Each of the three glasses were prepared by cooling with a constant rate from aconfiguration at the reference state point ( T = 1500K) to temperature 500 K. For each composition, once a glassconfiguration was obtained at the reference state point, we generated an isomorph in the same way as the liquidisomorphs by repeated DICs involving 5% density changes. Again, for comparison we probed the RDF and thedynamics also at isochoric state points with the same temperature variation as that of the isomorphs (compareFig. 1).The RDFs are shown in Fig. 9 (isomorphs) and in Fig. 10 (isochores). The picture is similar to that of the liquidphase: overall, good invariance along the isomorphs is seen in contrast to a substantially larger variation along theisochores. This is also the conclusion from Fig. 11 showing all the first-peak heights as a function of the temperature.The dynamics of the glasses is investigated via the MSD and the incoherent intermediate structure factor inFig. 12. A glass consists mostly of atoms frozen at fixed positions, merely vibrating there. Thus the MSD is virtuallyconstant, though some particle motion is discernible at the longest times. This “glass flow” motion does not appear to0 C u - C u R D F Isochore (liquid) r ~ Z r- Z r R D F C u - Z r R D F Cu Zr C u - C u R D F Isochore (liquid) r ~ Z r- Z r R D F C u - Z r R D F Cu Zr C u - C u R D F Isochore (liquid) r ~ Z r- Z r R D F C u - Z r R D F Cu Zr FIG. 6. Isochore liquid-state RDFs plotted as in Fig. 5. There is a considerably larger variation of the structure, with thestructhre being more washed out, the higher the temperature is. This reflects the increasing thermal fluctuations, which in thecase of an isomorph are compensated by increasing forces (Fig. 5). be isomorph invariant, but we found no systematic variation with the density. This indicates that the non-invariancereflects statistical uncertainty. Along the isochores (second column), it is clear that the glass gradually melts astemperature is increased when moving from the black curve representing the 500 K reference state point. The factthat the particles in the glass virtually do not move except for vibrations is also visible in the incoherent intermediatescattering function (right two columns of Fig. 12), which along the isomorphs stabilize at a constant level at longtimes. In contrast, many of the isochore curves go to zero at long times, reflecting an ease of motion with increasingtemperature that is not achieved along the isomorphs.
VIII. SUMMARY
This paper has studied three different compositions of the CuZr system by computer simulations using the efficientEMT potentials. We have traced out isomorphs in the liquid and glass phases of the three systems. Good isomorphinvariance was observed of structure and dynamics in both phases, showing that the atoms move about each otherin much the same way at state points on the same isomorph. This means that the thermodynamic phase diagramof CuZr systems is effectively one-dimensional. Thus for many purposes, in order to get an overview of the CuZrsystem, it is enough to investigate state points belonging to different isomorphs. It should be noted, though, thatsome quantities like, e.g., the bulk modulus is not isomorph invariant even when given in reduced units [62]. Onthe other hand, most material quantities are isomorph invariant in reduced units, e.g., the shear modulus, the shearviscosity, the heat conductivity, etc [62].1
Temperature [K] m a x ( R D F ) liquid FIG. 7. Maximum values of all RDFs of Fig. 5 and Fig. 6, i.e., values of the RDFs at their first peak. The full symbolsrepresent isomorphs and the open symbols represent isochores. While the isomorph values are not invariant, in particular forthe Cu-Zr RDFs, the general picture is that of a visibly better invariance along the isomorphs than along isochores of the sametemperature variation.
It would be interesting to confirm the above results using the EAM potentials. Given that the EAM and EMTboth have been shown to nicely reproduce metal properties, we do not expect significantly different results. Indeed,the fact that isomorph theory describes metals well has been validated in DFT simulations of crystals [44].
ACKNOWLEDGMENTS
This work was supported by the VILLUM Foundation’s
Matter grant (16515).
APPENDIX: MORE ON THE EMT
Jacobsen, Nørskov, and Puska [51] and Puska et al. [49] discussed the ansatz used in the EAM with argumentsbased on the effective-medium approach. Jacobsen et al. demonstrated how the cohesive energy of a metallic systemcould be related to the embedding energies, with corrections accounting for the d-d hybridization in the transitionmetals. Their approach showed that with the neglect of the d-d hybridization (valid for simple metals and presumablyfor early and late transition metals), the EAM expression is recovered. For Al and Cu, for example, the practicalapplication of the effective-medium theory is mathematically equivalent to the EAM. The density of the effectivemedium was taken to be an unweighted average of the background electron number density over the Wigner-Seitz cellof the atom.As mentioned in Sec. II, the EMT potential is defined in terms of seven parameters: E the negative of the cohesiveenergy, an electron density parameter n , the Wigner-Seitz radius s defined from the atomic density, η quantifyingthe influence of the density tail surrounding an neighboring atom, λ determined from the bulk modulus, and theproduct of parameters V ( βη − κ ) determined from the shear modulus. The following summarizes the discussiongiven in Ref. 51 on how the values for these parameters are obtained for a system with an ideal fcc crystal as referencesystem.DFT shows that the correction term is a sum of atomic-sphere and one-electron corrections. The former termcorrects for the difference in electrostatic and exchange-correlation energies between real and reference system. Thisis approximately equal to the difference in pair interactions of the two systems. The one-electron correction is given bythe difference of sums over one-electron energies. This contribution does not require any fitting input from experimentsand can be calculated fully ab initio [63]. While this would yield an accuracy comparable to that of DFT, it is alsosimilarly computationally expensive. To circumvent this, the one-electron contribution is momentarily neglected.Including only the cohesive and the atomic-sphere correction parts, the total energy expression can be simplified to2 . Z r- M S D Isomorph (liquid)Cu Zr . Z r- M S D Isochore (liquid)Cu Zr . . . . . Z r-I SF Isomorph (liquid)Cu Zr . . . . . Z r-I SF Isochore (liquid)Cu Zr . C u - M S D Cu Zr . C u - M S D Cu Zr . . . . . Z r-I SF Cu Zr . . . . . Z r-I SF Cu Zr t ~ . C u - M S D Cu Zr t ~ . C u - M S D Cu Zr t ~ . . . . . C u -I SF Cu Zr t ~ . . . . . C u -I SF Cu Zr FIG. 8. Dynamics along the isomorphs and isochores of the liquids. The two left columns give data for the reduced MSDas a function of the reduced time. The two right columns give the analogous reduced-unit intermediate incoherent scatteringfunction, for simplicity evaluated at the same wave vector in all three cases 2 πρ / , which is constant in reduced units. Thereis good isomorph invariance of the dynamics, but a significant variation along the isochores (except in the short-time ballisticregion where the reduced MSD by definition is 3˜ t at all state points). E = (cid:88) i (cid:16) E c,i + ∆ E AS ( i ) (cid:17) = (cid:88) i E c,i ( n i ) + 12 (cid:88) i (cid:54) = j V ij ( r ij ) − ref (cid:88) i (cid:54) = j V ij ( r ij ) . (11)Here the atomic-sphere correction term ∆ E AS ( i ) is expressed as the difference between pair potentials V ij of the realand the reference system. The neglected one-electron contribution is not necessarily small enough to be neglected,but it is possible to include this contribution as a small adjustment to the potential V ij .The density argument of the cohesive energy E c,i is the parameter that connects the environment around a givenatom i in the real system with the reference system. This so-called embedding density, n i , is a sum of contributionsfrom the neighboring atoms, n i = (cid:88) i (cid:54) = j ∆ n j ( s i , r ij ) (12)in which the density “tail” ∆ n j ( s i , r ij ) is assumed to have the following exponential form3 C u - C u R D F Isomorph (glass) r ~ Z r- Z r R D F C u - Z r R D F Cu Zr C u - C u R D F Isomorph (glass) r ~ Z r- Z r R D F C u - Z r R D F Cu Zr C u - C u R D F Isomorph (glass) r ~ Z r- Z r R D F C u - Z r R D F Cu Zr FIG. 9. Isomorph glass-state RDFs, plotted as a function of the reduced pair distance ˜ r . The picture is pretty much as forthe liquid, with good isomorph invariance except for the somewhat more noisy data that are presumably due to the fact that,basically, only a single configuration and its vibrations is probed. Significant deviations from isomorph invariance are observed,though, for the 64% Cu mixture, both at the highest densities for the Cu-Cu RDF and at the lowest densities for the Zr-ZrRDF. ∆ n ( s i , r ij ) = ∆ n exp (cid:2) η ( s i − s ) − η ( r ij − βs ) (cid:3) . (13)Here s is the equilibrium Wigner-Seitz radius, s i is the neutral sphere radius at non-zero pressure, and β is a geometricfactor related to the nearest-neighbor distance d nn via βs = d nn (meaning that β = (16 π/ / / √ ≈ . E AS ( i ) = 12 (cid:88) j (cid:54) = i v ( r ij ) − v ( βs i ) (14)where the second term corresponds to the fcc lattice having 12 nearest neighbors at distance r = βs i . The potentialis given by v ( r ) = − V exp (cid:18) − κβ ( r ij − βs ) (cid:19) . (15)Input from experiments is required to determine the remaining parameters describing the elastic and structuralproperties of the metal in question. The cohesive energy and the parameter λ are related to the bulk modulus B via4 C u - C u R D F Isochore (glass) r ~ Z r- Z r R D F C u - Z r R D F Cu Zr C u - C u R D F Isochore (glass) r ~ Z r- Z r R D F C u - Z r R D F Cu Zr C u - C u R D F Isochore (glass) r ~ Z r- Z r R D F C u - Z r R D F Cu Zr FIG. 10. Isochore glass-state RDFs, plotted as in Fig. 9. There is here a considerably larger variation of the structure. B = − E λ πs (16)while the product V ( βη − κ ) = V δ is determined from the shear modulus C via the expression C = sV δκ πs . (17)Equation (16) and Eq. (17) are only correct in the nearest-neighbor approximation; Ref. 51 explains how to includemore neighbors using a cutoff implemented via a Fermi function. The standard choice, which is used also in this work,corresponds to the edge of the Fermi function being between the 3rd and 4th neighbor shell of the equilibrium fcccrystal.The above gives all parameters describing a pure metal. For two-component systems like CuZn mixtures, somegeneralization of the energy expression is required [51]. [1] A. L. Greer, “Metallic glasses,” Science , 1947–1953 (1995).[2] M. Zink, K. Samwer, W. L. Johnson, and S. G Mayr, “Plastic deformation of metallic glasses: Size of shear transformationzones from molecular dynamics simulations,” Phys. Rev. B , 172203 (2006).
500 1000 1500
Temperature [K] m a x ( R D F ) glass FIG. 11. Maximum values of all RDFs of Fig. 9 and Fig. 10, i.e., value of the RDF at the first peak. The full symbolsrepresent isomorphs and the open symbols represent isochores. The isomorph values are not invariant, in particular for theCu-Zr functions, but the general picture is that of significantly better invariance along the isomorphs than along isochores ofsame temperature variation.[3] M.I. Mendelev, M.J. Kramer, R.T. Ott, and D.J. Sordelet, “Molecular dynamics simulation of diffusion in supercooledCu–Zr alloys,” Phil. Mag. , 109–126 (2009).[4] W. H. Wang, “The elastic properties, elastic models and elastic perspectives of metallic glasses,” Prog. Mater. Sci. ,487–656 (2012).[5] W. H. Wang, “Dynamic relaxations and relaxation-property relationships in metallic glasses,” Prog. Mater. Sci. ,100561 (2019).[6] M.-B. Tang, D.-Q. Zhao, M.-X. Pan, and W.-H. Wang, “Binary Cu-Zr bulk metallic glasses,” Chinese Phys. Lett. ,901–903 (2004).[7] S. A. Rice and P. Gray, The Statistical Mechanics of Simple Liquids (Interscience, New York, 1965).[8] H. N. V. Temperley, J. S. Rowlinson, and G. S. Rushbrooke,
Physics of Simple Liquids (Wiley, New York, 1968).[9] S. M. Stishov, “The Thermodynamics of Melting of Simple Substances,” Sov. Phys. Usp. , 625–643 (1975).[10] J.-L. Barrat and J.-P. Hansen, Basic Concepts for Simple and Complex Liquids (Cambridge University Press, 2003).[11] T. S. Ingebrigtsen, T. B. Schrøder, and J. C. Dyre, “What is a simple liquid?” Phys. Rev. X , 011011 (2012).[12] J.-P. Hansen and I. R. McDonald, Theory of Simple Liquids: With Applications to Soft Matter , 4th ed. (Academic, NewYork, 2013).[13] J. O. Hirschfelder, C. F. Curtiss, and R. B. Bird,
Molecular Theory of Gases and Liquids (John Wiley & Sons (New York),1954).[14] J. D. Bernal, “The Bakerian lecture, 1962. The structure of liquids,” Proc. R. Soc. London Ser. A , 299–322 (1964).[15] B. Widom, “Intermolecular forces and the nature of the liquid state,” Science , 375–382 (1967).[16] J. A. Barker and D. Henderson, “What is ”liquid”? Understanding the states of matter,” Rev. Mod. Phys. , 587–671(1976).[17] D. Chandler, J. D. Weeks, and H. C. Andersen, “Van der Waals picture of liquids, solids, and phase transformations,”Science , 787–794 (1983).[18] J. C. Dyre, “Simple liquids’ quasiuniversality and the hard-sphere paradigm,” J. Phys. Condens. Matter , 323001 (2016).[19] J. D. van der Waals, Over de Continuiteit van den Gas- en Vloeistoftoestand (Dissertation, University of Leiden, 1873).[20] J. S. Rowlinson, “The statistical mechanics of systems with steep intermolecular potentials,” Mol. Phys. , 107–115 (1964).[21] D. Henderson and J. A. Barker, “Perturbation theory of fluids at high temperatures,” Phys. Rev. A , 1266–1267 (1970).[22] H. C. Andersen, J. D. Weeks, and D. Chandler, “Relationship between the hard-sphere fluid and fluids with realisticrepulsive forces,” Phys. Rev. A , 1597–1607 (1971).[23] H. S. Kang, S. C. Lee, T. Ree, and F. H. Ree, “A perturbation theory of classical equilibrium fluids,” J. Chem. Phys. ,414–423 (1985).[24] K. R. Harris, “The selfdiffusion coefficient and viscosity of the hard sphere fluid revisited: A comparison with experimentaldata for xenon, methane, ethene and trichloromethane,” Molec. Phys. , 1153–1167 (1992).[25] J. E. Straub, “Analysis of the role of attractive forces in self-diffusion of a simple fluid,” Molec. Phys. , 373–385 (1992).[26] D. Ben-Amotz and G. Stell, “Reformulaton of Weeks–Chandler–Andersen perturbation theory directly in terms of a hard-sphere reference system,” J. Phys. Chem. B , 6877–6882 (2004).[27] A. E. Nasrabad, “Thermodynamic and transport properties of the Weeks–Chandler–Andersen fluid: Theory and computersimulation,” J. Chem. Phys. , 244508 (2008).[28] T. Rodriguez-Lopez, J. Moreno-Razo, and F. del Rio, “Thermodynamic scaling and corresponding states for the self- . Z r- M S D Isomorph (glass)Cu Zr . Z r- M S D Isochore (glass)Cu Zr . . . . . Z r-I SF Isomorph (glass)Cu Zr . . . . . Z r-I SF Isochore (glass)Cu Zr . C u - M S D Cu Zr . C u - M S D Cu Zr . . . . . Z r-I SF Cu Zr . . . . . Z r-I SF Cu Zr t ~ . C u - M S D Cu Zr t ~ . C u - M S D Cu Zr t ~ . . . . . C u -I SF Cu Zr t ~ . . . . . C u -I SF Cu Zr FIG. 12. Dynamics along the glass isomorphs and isochores of the glasses. The two left columns give data for the reduced MSDas a function of the reduced time. The two right columns give the analogous reduced-unit intermediate incoherent scatteringfunction data, evaluated at the wave vector corresponding to the first peak of the corresponding RDF. Because the systemis a glass, along the isomorphs the MSD is constant over a very long time and, correspondingly, the incoherent intermediatescattering function does not decay to zero. For the isochore state points, raising the temperature takes the system closer to aliquid with a MSD that at long times is proportional to time and an intermediate incoherent scattering function that decaysto zero at long times. The overall picture is that, as for the liquid (Fig. 8), there is good isomorph invariance of the dynamics,but a significant variation along the isochores.diffusion coefficient of non-conformal soft-sphere fluids,” J. Chem. Phys. , 114502 (2013).[29] Y. Rosenfeld, “Relation between the transport coefficients and the internal entropy of simple systems,” Phys. Rev. A ,2545–2549 (1977).[30] J. C. Dyre, “Perspective: Excess-entropy scaling,” J. Chem. Phys. , 210901 (2018).[31] D. Kivelson, G. Tarjus, X. Zhao, and S. A. Kivelson, “Fitting of viscosity: Distinguishing the temperature dependencespredicted by various models of supercooled liquids,” Phys. Rev. E , 751–758 (1996).[32] C. Alba-Simionesco, A. Cailliaux, A. Alegria, and G. Tarjus, “Scaling out the density dependence of the alpha relaxationin glass-forming polymers,” Europhys. Lett. , 58–64 (2004).[33] C. M. Roland, S. Hensel-Bielowka, M. Paluch, and R. Casalini, “Supercooled dynamics of glass-forming liquids andpolymers under hydrostatic pressure,” Rep. Prog. Phys. , 1405–1478 (2005).[34] D. Gundermann, U. R. Pedersen, T. Hecksher, N. P. Bailey, B. Jakobsen, T. Christensen, N. B. Olsen, T. B. Schrøder,D. Fragiadakis, R. Casalini, C. M. Roland, J. C. Dyre, and K. Niss, “Predicting the density–scaling exponent of aglass–forming liquid from Prigogine–Defay ratio measurements,” Nat. Phys. , 816–821 (2011).[35] E. R. Lopez, A. S Pensado, J. Fernandez, and K. R. Harris, “On the Density Scaling of pVT Data and Transport Propertiesfor Molecular and Ionic Liquids,” J. Chem. Phys. , 214502 (2012).[36] K. Adrjanowicz, M. Paluch, and J. Pionteck, “Isochronal superposition and density scaling of the intermolecular dynamicsin glass-forming liquids with varying hydrogen bonding propensity,” RSC Adv. , 49370 (2016).[37] Y.-C. Hu, B.-S. Shang, P.-F. Guan, Y. Yang H.-Y. Bai, and W.-H. Wang, “Thermodynamic scaling of glassy dynamics and dynamic heterogeneities in metallic glass-forming liquid,” J. Chem. Phys. , 104503 (2016).[38] C. M. Roland, R. Casalini, and M. Paluch, “Isochronal temperature–pressure superpositioning of the alpha–relaxation intype-A glass formers,” Chem. Phys. Lett. , 259–264 (2003).[39] K. L. Ngai, R. Casalini, S. Capaccioli, M. Paluch, and C. M. Roland, “Do theories of the glass transition, in which thestructural relaxation time does not define the dispersion of the structural relaxation, need revision?” J. Phys. Chem. B , 17356–17360 (2005).[40] H. W. Hansen, A. Sanz, K. Adrjanowicz, B. Frick, and K. Niss, “Evidence of a one-dimensional thermodynamic phasediagram for simple glass-formers,” Nat. Commun. , 518 (2018).[41] T. B. Schrøder and J. C. Dyre, “Simplicity of condensed matter at its core: Generic definition of a Roskilde-simple system,”J. Chem. Phys. , 204502 (2014).[42] N. Gnan, T. B. Schrøder, U. R. Pedersen, N. P. Bailey, and J. C. Dyre, “Pressure-energy correlations in liquids. IV.“Isomorphs” in liquid phase diagrams,” J. Chem. Phys. , 234504 (2009).[43] J. C. Dyre, “Hidden scale envariance in condensed matter,” J. Phys. Chem. B , 10007–10024 (2014).[44] F. Hummel, G. Kresse, J. C. Dyre, and U. R. Pedersen, “Hidden scale invariance of metals,” Phys. Rev. B , 174116(2015).[45] L. Friedeheim, J. C. Dyre, and N. P. Bailey, “Hidden scale invariance at high pressures in gold and five other face-centered-cubic metal crystals,” Phys. Rev. E , 022142 (2019).[46] R. Casalini and T. C. Ransom, “On the experimental determination of the repulsive component of the potential from highpressure measurements: What is special about twelve?” J. Chem. Phys. , 194504 (2019).[47] A. Sanz, T. Hecksher, H. W. Hansen, J. C. Dyre, K. Niss, and U. R. Pedersen, “Experimental evidence for a state-point-dependent density-scaling exponent of liquid dynamics,” Phys. Rev. Lett. , 055501 (2019).[48] R. Casalini and T. C. Ransom, “On the pressure dependence of the thermodynamical scaling exponent γ ,” Soft Matter , 4625–4631 (2020).[49] M. J. Puska, R. M. Nieminen, and M. Manninen, “Atoms embedded in an electron gas: Immersion energies,” Phys. Rev.B , 3037–3047 (1981).[50] J. K. Nørskov, “Covalent effects in the effective-medium theory of chemical binding: Hydrogen heats of solution in the 3 d metals,” Phys. Rev. B , 2875–2885 (1982).[51] K. W. Jacobsen, P. Stoltze, and J. K. Nørskov, “A semi-empirical effective medium theory for metals and alloys,” Surf.Sci. , 394–402 (1996).[52] S. M. Foiles, M. I. Baskes, and M. S. Daw, “Embedded-atom-method functions for the fcc metals Cu, Ag, Au, Ni, Pd, Pt,and their alloys,” Phys. Rev. B , 7983–7991 (1986).[53] M. S. Daw, S. M. Foiles, and M. I. Baskes, “The embedded-atom method: a review of theory and applications,” Mater.Sci. Rep. , 251–310 (1993).[54] N. P. Bailey, T. S. Ingebrigtsen, J. S. Hansen, A. A. Veldhorst, L. Bøhling, C. A. Lemarchand, A. E. Olsen, A. K. Bacher,L. Costigliola, U. R. Pedersen, H. Larsen, J. C. Dyre, and T. B. Schrøder, “RUMD: A general purpose molecular dynamicspackage optimized to utilize GPU hardware down to a few thousand particles,” Scipost Phys. , 038 (2017).[55] K. W. Jacobsen, J. K. Nørskov, and M. J. Puska, “Interatomic interactions in the effective-medium theory,” Phys. Rev.B , 7423–7442 (1987).[56] J. K. Nørskov and N. D. Lang, “Effective-medium theory of chemical binding: Application to chemisorption,” Phys. Rev.B , 2131–2136 (1980).[57] A. Pˇaduraru, A. Kenoufi, N. P. Bailey, and J. Schiøtz, “An interatomic potential for studying CuZr bulk metallic glasses,”Adv. Eng. Mater. , 505–508 (2007).[58] N. P. Bailey, U. R. Pedersen, N. Gnan, T. B. Schrøder, and J. C. Dyre, “Pressure-energy correlations in liquids. I. Resultsfrom computer simulations,” J. Chem. Phys. , 184507 (2008).[59] E. Attia, J. C. Dyre, and U. R. Pedersen, “An extreme case of density scaling: The Weeks-Chandler-Andersen system atlow temperatures,” arXiv:2101.10663 (2021).[60] L. Bøhling, T. S. Ingebrigtsen, A. Grzybowski, M. Paluch, J. C. Dyre, and T. B. Schrøder, “Scaling of viscous dynamicsin simple liquids: Theory, simulation and experiment,” New J. Phys. , 113035 (2012).[61] T. S. Ingebrigtsen, L. Bøhling, T. B. Schrøder, and J. C. Dyre, “Thermodynamics of condensed matter with strongpressure-energy correlations,” J. Chem. Phys. , 061102 (2012).[62] D. M. Heyes, D. Dini, L. Costigliola, and J. C. Dyre, “Transport coefficients of the Lennard-Jones fluid close to thefreezing line,” J. Chem. Phys. , 204502 (2019).[63] N. Chetty, K. Stokbro, K. W. Jacobsen, and J. K. Nørskov, “ Ab initio potential for solids,” Phys. Rev. B46