Determining Critical Temperature Tc in a Molecular Dynamics-Simulated Glass Forming Ni 0 . 8 Zr 0 . 2-System

In this paper molecular dynamics computer simulations is used to investigate a critical temperature Tc for a dynamical glass transition as proposed by the mode-coupling theory (MCT) of dense liquids in a glass forming Ni0.8Zr0.2-system. The critical temperature Tc are analyzed from different quantities and checked the consistency of the estimated values, i.e. from (i) the non-vanishing nonergodicity parameters as asymptotic solutions of the MCT equations in the arrested state, (ii) the gm-parameters describing the approach of the melt towards the arrested state on the ergodic side, (iii) the diffusion coefficients in the melt. The resulting Tc values are found to agree within about 10-15 %.


Introduction
The glass problem has a reputation of being controversial [1,2].It has even been said that "There are more theories of the glass transition than there are theorists who propose them" [3].Part of this controversy is fully justified.The problem is indeed extremely challenging, and it has long been difficult for theory to provide precise predictions that could be experimentally tested, in order to obtain conclusive answers [4].A lot of progress has nonetheless been achieved thanks to the impressive numerical and experimental developments of the last couple of decades, and theory is rapidly evolving in the aim to match this effort.
At the experimental level, the so-called glass transition is generally associated with a sharp increase in the characteristic relaxation times of the system, and a concomitant departure of laboratory measurements from equilibrium.At the theoretical level, it has been proposed that the transition from a liquid to a glassy state is triggered by an underlying thermodynamic (equilibrium) transition [5]; in that view, an "ideal" glass transition is believed to occur at the so-called Kauzmann temperature, T K .At T K , it is proposed that only one minimum-energy basin of attraction is accessible to the system.One of the first arguments of this type is due to Gibbs and diMarzio [6], but more recent studies using replica methods have yielded evidence in support of such a transition in Lennard-Jones glass formers [5,7,8].These observations have been called into question by experimental data and recent results of simulations of polydisperse hard-core disks, which have failed to detect any evidence of a thermodynamic transition up to extremely high packing fractions [9].One of the questions that arises is therefore whether the discrepancies between the reported simulated behavior of hard-disk and soft-sphere systems is due to fundamental differences in the models, or whether they are a consequence of inappropriate sampling at low temperatures and high densities.
Different, alternative theoretical considerations have attempted to establish a connection between glass transition phenomena and the rapid increase in relaxation times that arises in the vicinity of a theoretical critical temperature (the so-called "mode-coupling" critical temperature, T c ), thereby giving rise to a "kinetic" or "dynamic" transition [10].In recent years, both viewpoints have received some support from molecular simulations.Many of these simulations have been conducted in the context of models introduced by Stillinger and Weber and by Kob and Andersen [11]; such models have been employed in a number of studies that have helped shape our current views about the glass transition [7,12,13,14,15,16].
As in our previous paper [17] (it is referred as paper I), in this paper we investigate a critical temperature T c for a dynamical glass transition as proposed by the modecoupling theory (MCT) of dense liquids in a glass forming Ni 0.8 Zr 0.2 -system, while in paper I the system is a glass forming Ni 0.2 Zr 0.8 -system.The critical temperature T c are analyzed from different quantities and checked the consistency of the estimated values, i.e. from (i) the nonvanishing nonergodicity parameters as asymptotic solutions of the MCT equations in the arrested state, (ii) the g m -parameters describing the approach of the melt towards the arrested state on the ergodic side, (iii) the diffusion coefficients in the melt.
Thereforee it is tempting to consider how well the estimates of T c from different approaches fit together and whether the T c estimate from the nonergodicity parameters of the idealized MCT compares to the values from the full MCT.Regarding this, we here investigate a molecular dynamics (MD) simulation model adapted to the glass-forming Ni 0.8 Zr 0.2 transition metal system.The Ni x Zr 1-x -system is well studied by experiments [24,25] and by MD-simulations [26,27,28,29,30,31,32,33], as it is a rather interesting system whose components are important constituents of a number of multi-component 'massive' metallic glasses.In the present contribution we consider, in particular, the 0.8 = x composition and concentrate on the determination of T c from evaluating and analyzing the nonergodicity parameter, the g m (T)parameter in the ergodic regime, and the diffusion coefficients.
Our paper is organized as follows: In section II, we present the model and give some details of the computations.Section III.gives a brief discussion of some aspects of the mode coupling theory as used here.Results of our MD-simulations and their analysis are then presented and discussed in Section IV.

Simulations
The present simulations are carried out as state-of-the-art isothermal-isobaric calculations.The Newtonian equations of N= 648 atoms (518 Ni and 130 Zr) are numerically integrated by a fifth order predictorcorrector algorithm with time step ǻt = 2.5x10 -15 s in a cubic volume with periodic boundary conditions and variable box length L.
With regard to the electron theoretical description of the interatomic potentials in transition metal alloys by Hausleitner and Hafner [34], we model the interatomic couplings as in [27] by a volume dependent electron-gas term E vol (V) and pair potentials ‫(‬r) adapted to the equilibrium distance, depth, width, and zero of the Hausleitner-Hafner potentials [34] for Ni 0.8 Zr 0.2 [28].For this model, simulations were started through heating a starting configuration up to 2000 K which leads to a homogeneous liquid state.The system then is cooled continuously to various annealing temperatures with cooling ratet T = 1.5 x 10 12 K/s.Afterwards the obtained configurations at various annealing temperatures (here 1500-600 K) are relaxed by carrying out additional isothermal annealing runs.Finally the time evolution of these relaxed configurations is modelled and analyzed.More details of the simulations are given in [28].
In this section we provide some basic formulae that permit calculation of T c and the nonergodicity parameters f ij (q) for our system.A more detailed presentation may be found in Refs.[21,22,23,35,36].The central object of the MCT are the partial intermediate scattering functions which are defined for a binary system by [37] ,0) where The normalized partial-and incoherent intermediate scattering functions are given by , (5) where the S ij (q) = F ij (q, t = 0) are the partial static structure factors.
The basic equations of the MCT are the set of nonlinear matrix integrodifferential equations where F is the 2x2 matrix consisting of the partial intermediate scattering functions F ij (q, t), and the frequency matrix ȍ 2 is given by S(q) denotes the 2x2 matrix of the partial structure factors S ij (q), x i = N i /N and m i means the atomic mass of the species i.The MCT for the idealized glass transition predicts [41] that the memory kern M can be expressed at long times by where ȡ = N/Vol is the particle density and the vertex and the matrix of the direct correlation function is defined by The equation of motion for ) , ( t q F s i has a similar form as Eq.( 6), but the memory function for the incoherent intermediate scattering function is given by EPJ Web of Conferences 00011-p.2

Nonergodicity parameters
In order to characterize the long time behaviour of the intermediate scattering function, the nonergodicity parameters f(q) are introduced as .
These parameters are the solution of eqs.( 6)-( 10) at long times.The meaning of these parameters is the following: if f ij (q) = 0, then the system is in a liquid state with density fluctuation correlations decaying at long times.If f ij (q) > 0, the system is in an arrested, nonergodic state, where density fluctuation correlations are stable for all times.In order to compute f ij (q), one can use the following iterative procedure [23]: where the matrix N(q) is given by .(15) This iterative procedure, indeed, has two type of solutions, nontrivial ones with f(q) > 0 and trivial solutions f(q) = 0 .
The incoherent nonergodicity parameter ) (q f s i can be evaluated by the following iterative procedure: As indicated by Eq.( 16), computation of the incoherent nonergodicity parameter ) (q f s i demands that the coherent nonergodicity parameters are determined in advance.

m g -parameter
Beyond the details of the MCT, equations of motion like (6) can be derived for the correlation functions under rather general assumptions within the Lanczos recursion scheme [38] resp.the Mori-Zwanzig formalism [39].The approach demands that the time dependence of fluctuations A, B, ... is governed by a time evolution operator like the Liouvillian and that for two fluctuating quantitites a scalar products (B, A) with the meaning of a correlation function can be defined.In case of a tagged particle, this leads for ) , ( t q s i ) to the exact equation are hidden all the details of the time evolution of ) , ( t q s i ) . As proposed and applied in [19,20], instead of calculating ) , ( 0 t q M i from the time evolution operator as a continued fraction, it can be evaluated in closed forms once ) , ( t q s i ) is known, e.g., from experiments or MD-simulations.This can be demonstrated by introduction of ^`, ) with the Laplace transform of ) (t ) , and Eq.( 17) then leads to On the time axis, Adopting some arguments from the schematic MCT, Eq.( 17) allows asymptotically finite correlations , that means an arrested state, if remains finite where the relationship holds In order to characterize the undercooled melt and its transition into the glassy state, we introduced in [19 According to (23), ) , ( has the property that ) in the arrested, nonergodic state.On the other hand, if there is no arrested solution and the correlations ) , ( t q s i ) decay to zero for t ĺ , that means, the system is in the liquid state.From that we proposed [19] to use the value of g m as a relative measure how much the system has approached the arrested state and to use the temperature dependence of g m (T ) in the liquid state as an indication how the system approaches this state.

3.2
First we show the results of our simulations concerning the static properties of the system in terms of the partial structure factors S ij (q) and partial correlation functions g ij (r).
To compute the partial structure factors S ij (q) for a binary system we use the following definition [40] , 1) where are the partial pair correlation functions.
The MD simulations yield a periodic repetition of the atomic distributions with periodicity length L. Truncation of the Fourier integral in Eq.( 27) leads to an oscillatory behavior of the partial structure factors at small q .In order to reduce the effects of this truncation, we compute from Eq.( 28) the partial pair correlation functions for distance r up to For numerical evaluation of Eq.( 27), a Gaussian type damping term is included with Fig. 1-fig.3 shows the partial structure factors S ij (q) versus q for all temperatures investigated.The figure indicates that the shape of S ij (q) depends weakly on temperature only and that, in particular, the positions of the first maximum and the first minimum in S ij (q) are more or less temperature independent.To investigate the dynamical properties of the system, we have calculated the incoherent intermediare scattering function (ISF) ) , ( t q F s i and the coherent intermediare scattering function (ISF) ) , ( t q F ij as defined in equations ( 1) and (3).  of both species evaluated from our MD data for wave vector n q = L n/ 2S with n = 9, that means 24.4 = 9 q nm -1 .From the figure we see that ) , ( t q s i ) of both species shows at intermediate temperatures a structural relaxation in three succesive steps as predicted by the idealized schematic MCT [42].The first step is a fast initial decay on the time scale of the vibrations of atoms ( 0.2 < t ps).This step is characterized by the MCT only globaly.The second step is the E -relaxation regime.In the early E -regime the Between them a wide plateau is found near the critical temperature T c .In the melt, the D -relaxation takes place as the last decay step after the von Schweidler-law.It can be described by the Kohlrausch-Williams-Watts where the relaxation time D W near the glass transition shifts drastically to longer times.
The inverse power-law decay for the early E -regime is not seen in our data.This seems to be due to the fact that in our system the power-law decay is dressed by the atomic vibrations ( [19,20] and references therein).
According to our MD-results, ) , ( t q s i ) decays to zero for longer times at all temperatures investigated.This is in agreement with the full MCT.Including transversal currents as additional hydrodynamic variables, the full MCT [41] comes to the conclusion that all structural correlations decay in the final D -process, independent of temperature.Similar effects are expected from inclusion of thermally activated matter transport, that means diffusion in the arrested state.
At T = 900 K -800 K, the ) , ( t q s i ) drop rather sharply at large t .This reflects aging effects which take place, if a system is in a transient, non-steady state [43].Such a behaviour indicates relaxations of the system on the time scale of the 'measuring time' of the correlations.
The nonergodicity parameters are defined by Eq.( 13) as a non-vanishing asymptotic solution of the MCT Eq.( 6).Fig. 6 presents the estimated q -dependent nonergodicity parameters from the coherent scattering functions of Ni and Zr at T = 1005 K.In order to compute the nonergodicity parameters f ij (q) analytically, we followed for our binary system the self-consistent method as formulated by Nauroth and Kob [23] and as sketched in Section III.A. Input data for our iterative determination of f ij (q) = F ij (q,) are the temperature dependent partial structure factors S ij (q) from the previous subsection.The iteration is started by arbitrarily setting For T > 1100 K we always obtain the trivial solution f ij (q) = 0 while at T = 1000 K and below we get stable non-vanishing f ij (q) > 0. The stability of the nonvanishing solutions was tested for more than 3000 iteration steps.From this results we expect that T c for our system lies between 1000 and 1100 K. To estimate T c more precisely, we interpolated S ij (q) from our MD data for temperatures between 1000 and 1100 K by use of the algorithm of Press et.al.[44].We observe that at 1005 = T K a non-trivial solution of f ij (q) can be found, but not at 1010 = T K and above.It means that the critical temperature T c for our system is around 1005 K.The non-trivial solutions f ij (q) for this temperature shall be denoted the critical nonergodicty parameters f cij (q).They are included in Fig. 6.
By using the critical nonergodicity parameters f cij (q), the computational procedure was run to determine the critical nonergodicity parameters ) (q f s ci for the incoherent scattering functions at 1005 = T K .Here we present our results about the ) , ( function [19,20] described in section III.B.The memory functions ) , ( 0 t q M i are evaluated from the MD data for ) , ( t q s i ) by Fourier transformation along the positive time axis.For completeness, also 800 K data is included where the corresponding ) , ( t q s i ) are extrapolated to longer times by use of an KWW approximation.Fig. 7 and Fig. 8 show the thus deduced ) , ( 0 t q M i for 24.4 = q nm -1 .Regarding their qualitative features, the obtained ) , ( 0 t q M i are in full agreement with the results in [20] for the Ni 0.5 Zr 0.5 system.A particular interesting detail is the fact that there exists a minimum in ) , ( 0 t q M i for both species, Ni and Zr, at all investigated temperatures around a time of 0.1 ps.Below this time, ) , ( t q s i ) reflects the vibrational dynamics of the atoms.Above this value, the escape from the local cages takes place in the melt and the E -regime dynamics are developed.Apparently, the minimum is related to this crossover.In Fig. 9 and Fig. 10 we display )) , ( , ( ) . In this figure we again find the features already described for Ni 0.5 Zr 0.5 in [19,20].According to the plot, there exist (q-dependent) limiting values ) , ( is close to an universal behavior, while for ) , ( > ) , ( marked deviations are seen.) , ( 0 t q s i ) significantly decreases with increasing temperature.It is tempting to identify ) , ( with the polynomial form for ) , ( 0 t q M i assumed in the schematic version of the MCT [41].
By use of the calculated memory functions, we can evaluate the ) , ( , Eq.( 24).In Fig. 11 and Fig. has a maximum g m (q,T) at an intermediate value of ) .In the high temperature regime, the values of ) , ( T q g m move with EPJ Web of Conferences 00011-p.6

4.3.
decreasing temperature towards the limiting value 1.This is, in particular, visible in Fig. 13 where we present g m (q,T) as function of temperature for both species, Ni and Zr, and wave-vectors q 9 = 24.4nm -1 .At temperatures above 1000 K, the g m -values increase approximately linear towards 1 with decreasing temperatures.Below 1000 K, they remain close below the limiting value of 1, a behavior denoted in [19,20] as a balancing on the borderline between the arrested and the non-arrested state due to thermally induced matter transport by diffusion in the arrested state at the present high temperatures.) for q = 24.4nm -1 Linear fit of the g m -values for Ni above 950 K and for Zr above 1000 K predicts a crossover temperature T c * from liquid ( g m < 1) to the quasi-arrested ( g m = 1) behavior around 970 K from the Ni data and around 1020 K from the Zr data.We here identify this crossover temperature with the value of T c as visible in the ergodic, liquid regime and estimate it by the mean value from the Ni-and Zr-subsystems, that means by T c = 1000 K.
While in [19,20] for the Ni 0.5 Zr 0.5 melt a T c -value of 1120 K was estimated from g m (T), the value for the present composition is lower by about 120 K.A significant composition dependence of T c is expected according to the results of MD simulation for the closely related Co x Zr 1-x system [18].Over the whole x-range, T c was found to vary between 1170 and 650 K in Co x Zr 1-x , with T c ( x=0.2 ) ; 800 K. Regarding this, the present data for the Ni x Zr 1-x system reflect a rather weak T c variation.
Fig. 13.MD simulation results of the temperature dependence of g m (q,T) for q 9 = 24.4nm -1 (symbols).Linear fits to the g m (q,T) are included by full and dash lines ( for q 9 = 24.4nm -1 ); a) Zr-part with T c = 1020 K and b) Ni-part with T c = 970 K.
From the simulated atomic motions in the computer experiments, the diffusion coefficients of the Ni and Zr species can be determined as the slope of the atomic mean square displacements in the asymptotic long-time limit (31) with non-universal exponent Ȗ [45].
In order to estimate T c from this relationship, we have adapted the critical power law by a least mean squares fit to the simulated diffusion data for 1000 K and above.According to this fit, the system has a critical temperature of about 850-900 K.
Similar results for the temperature dependence of the diffusion coefficients have been found in MD simulations for other metallic glass forming systems, e.g., for Ni 0.5 Zr 0.5 [28], for Ni x Zr 1-x [18], Cu 0.33 Zr 0.67 [46], or Ni 0.81 B 0.19 [47].In all cases, like here, a break is observed ICASCE 2013 00011-p.7

Diffusion-Coefficients
in the Arrhenius slope.In the mentioned Zr-systems, this break is related to a change of the atomic dynamics around T c whereas for Ni 0.81 B 0.19 system it is ascribed to T G As in [47] T c and T G apparently fall together, there is no serious conflict between the obervations.There is close agreement between the T c values estimated from the dynamics in the undercooled melt when approaching T c from the high temperature side.The values are T c § 950 -1020 K from the m g -parameters, and T c § 950 K from the diffusion coefficients.As discussed in [18], the T c -estimates from the diffusion coefficients seem to depend on the upper limit of the temperature region taken into account in the fit procedure, where an increase in the upper limit increases the estimated T c .Accordingly, there is evidence that the present value of 950 K may underestimate the true T c by about 10 to 50 K, as it based on an upper limit of 2000 K only.Taking this into account, the present estimates from the melt seem to lead to a T c value around 1000 K.
The T c from the nonergodicity parameters describe the approach of the system towards T c from the low temperature side.They predict a T c value of 1005 K.This value is clearly outside the range of our T c estimates from the high temperature, ergodic melt.We consider this as a significant deviation which, however, is much smaller than the factor of two found in the modelling of a Lennard-Jones system [23].The here observed deviation between the T c estimates from the ergodic and the socalled nonergodic side reconfirm the finding from the soft spheres model [22] of an agreement within some 10 % between the different T c -estimates.
In the future to overcome these discrepancies, it will be included in computation the static triplet correlations, as suggested by D.Coslovich [48].From the theoretical side, prediction of three-point dynamic correlation functions has been possible within the so-called inhomogeneous mode-coupling theory (IMCT) [49], an extension of the standard mode-coupling theory (MCT) [2] of the glass transition, which accounts for the dynamic response of a fluid to an external field.The predictions of IMCT suffers from the same shortcomings of MCT, i.e., the dynamic correlation volume diverges at a finite temperature; furthermore, the predicted scaling behavior in the vicinity of the avoided, putative singularity is inconsistent with the available simulation data for four-point dynamic susceptibilities [50,51].

Fig. 4 and
Fig.4 and Fig.5 presents the normalized incoherent intermediate scattering functions ) , ( t q s i )of both species evaluated from our MD data for wave vector n q = L n/ 2S with n = 9, that means 24.4 = late E -relaxation regime, which appears only in the melt, according the von

Fig. 6 .
Fig. 6.Non-ergodicity parameter f ij for coherent intermediate scattering functions as solutions of eqs.(13)(red-line for Zr and black-line for Ni)

Fig. 11 .
Fig. 11.Ni-part: MD simulation results for the characteristic function g(ĭ s ) as a function ofs i

Fig. 14
Fig. 14 shows the thus calculated diffusion coefficients of our Ni 0.8 Zr 0.2 model for the temperature range between 600 and 2000 K.At temperatures above approximately 1000 K, the diffusion coefficients for both species run parallel to each other in the Arrhenius plot, indicating a fixed ratio D Ni /D Zr § 2.5 in this temperature regime.At lower temperatures, the Ni atoms have a lower mobility than the Zr atoms, yielding around 800 K a value of about 10 for D Ni /D Zr .That means, here the Zr atoms carry out a rather rapid motion within a relative immobile Ni matrix.According to the MCT, above T c the diffusion coefficients follow a critical power law c c i T T T T T D > for , ) ( ) ( J :(31) with non-universal exponent Ȗ[45].In order to estimate T c from this relationship, we have adapted the critical power law by a least mean squares fit to the simulated diffusion data for 1000 K and above.According to this fit, the system has a critical temperature of about 850-900 K.Similar results for the temperature dependence of the diffusion coefficients have been found in MD simulations for other metallic glass forming systems, e.g., for Ni 0.5 Zr 0.5[28], for Ni x Zr 1-x[18], Cu 0.33 Zr 0.67[46], or Ni 0.81 B 0.19[47].In all cases, like here, a break is observed

Fig. 14 .
Fig. 14.Diffusion coefficients D i as a function of 1000/T.Symbols are MD results for Zr (squares) and Ni (diamonds)