Providing a physically sound explanation of aging phenomena in non-equilibrium amorphous materials is a challenging problem in modern statistical thermodynamics. The slow evolution of physical properties after quenches of control parameters is empirically well interpreted via the concept of material time (or internal clock) based on the Tool–Narayanaswamy–Moynihan model. Yet, the fundamental reasons of its striking success remain unclear. We propose a microscopic rationale behind the material time on the basis of the linear laws of irreversible thermodynamics and its extension that treats the corresponding kinetic coefficients as state functions of a slowly evolving material state. Our interpretation is based on the recognition that the same mathematical structure governs both the Tool model and the recently developed non-equilibrium extension of the self-consistent generalized Langevin equation theory, guided by the universal principles of Onsager’s theory of irreversible processes. This identification opens the way for a generalization of the material-time concept to aging systems where several relaxation modes with very different equilibration processes must be considered, and partially frozen glasses manifest the appearance of partial ergodicity breaking and, hence, materials with multiple very distinct inner clocks.
Amorphous states of matter that derive from the non-equilibrium solidification of supercooled liquids, such as most glasses and gels, are ubiquitous in nature and of significant technological and scientific relevance.1–4 One of the most essential fingerprints experimentally observed in these non-equilibrium amorphous states is the slow evolution of physical properties (e.g., volume and transport coefficients) with the “waiting” time t after preparation, a phenomenon commonly referred to as physical aging.5–16 The wealth of experimental data on the equilibration and aging processes of glass- and gel-forming liquids, thus, becomes a key to test the predictive power of available theories of glasses.17–20 Many such data are successfully fitted by phenomenological models, and it remains as a challenge for statistical-mechanics theories to reconcile the phenomenology with fundamental principles.21–24
Many of these phenomenological models were built out of the need to include essential effects such as the existence of several (not only one) relaxation modes or the presence of dynamical heterogeneities.14,16 For example, the need to describe measured thermal and viscoelastic properties during vitrification and devitrification under a wide variety of experimental conditions and thermal histories led to the development of powerful phenomenological tools, such as the Tool–Narayanaswamy–Moynihan (TNM5–7) and the Kovacs–Aklonis–Hutchinson–Ramos (KAHR12) models. These models are commonly used in industry to predict aging effects,25 and their development involved a rich discussion of many relevant issues.14 For our present purpose, however, we shall focus on the TNM model, and to simplify the discussion, we only highlight the essential feature of Tool’s original proposal,5 namely, to model the highly nonlinear aging processes observed in a given observable X(t) in terms of a simple mathematical representation: The model posits the generic nonlinear equation
where τ(X(t)) is a characteristic relaxation time that itself exhibits aging and that is postulated to depend on waiting time only through its dependence on X(t). A change of variable from the time t to the “material” (or “inner”) time ζ, defined as6
allows us to write the solution of Eq. (1) as X(t) = X*(ζ(t)) with the function X*(ζ) that obeys the ordinary linear differential equation
The solution of Eq. (1) with the arbitrary initial condition X(t = 0) = X0 is, thus, given by
Despite its general success, the applicability of this simple version of the TNM model is known to be limited, in practice, and this prompted the proposal of alternative independent empirical models.14 In fact, already in the original proposal, Narayanaswamy considered extended versions.6 The resulting description turns out to capture the essence of the behavior experimentally observed in many glass forming materials, although the fundamental reason for this is not fully understood.16,26 Determining the boundaries of validity of the TNM model and the origin of its corresponding limitations is still active matter of current research.27 Only recently, careful long-time experiments could establish the very existence of an internal “clock” as a manifestation of more fundamental physical principles in a real material.28
A general statistical mechanics theory of irreversible processes in liquids was proposed two decades ago.29,30 This approach, referred to as the non-equilibrium self-consistent generalized Langevin equation (NE-SCGLE) theory, predicts many relevant universal signatures of the glass and gel transitions—including aging phenomena—and more specific features reflecting the role of the particular molecular interactions acting in the specific system considered. For example, for simple liquids with purely repulsive interparticle interactions, this theoretical framework accurately describes the non-stationary and non-equilibrium process of the formation of high-density hard-sphere-like glasses,30–32 whereas for liquids with repulsive plus attractive interactions at low densities and temperatures, it predicts the formation of sponge-like gels and porous glasses by arrested spinodal decomposition.33–36
The NE-SCGLE was formulated in Ref. 29 and can be summarized in terms of a set of equations that describe the non-equilibrium evolution of the structural and dynamical properties of a simple glass-forming liquid [we refer here for simplicity to a generic system of N identical spherical particles in a volume V that interact through a radially symmetric pair potential U(r)]. The most central of these equations describes the time evolution of the non-equilibrium structure factor (SF) S(k; t) ≡⟨δn(k) δn(−k)⟩t, where the brackets indicate the average over the time-dependent non-equilibrium statistical ensemble and δn(k) is the Fourier-transformed (FT) fluctuating particle-number density . For a homogeneous system, instantaneously quenched at t = 0 from an initial equilibrated state corresponding to density and temperature , toward new values , the NE-SCGLE for the SF reads
The evolution of the non-equilibrium SF is driven by a thermodynamic force given by , the FT of the second functional derivative of a postulated Helmholtz free energy density-functional , i.e., ,37 evaluated at the uniform (bulk) density and temperature fields30, and T(r, t) = T, where N, V, and T are, respectively, the total particle number, volume, and temperature of the final state of the quench, and β−1 = kBT with the Boltzmann constant kB. The factor D0 in the NE-SCGLE relaxation [Eq. (5)] arises from the diffusive relaxation of density fluctuations and is related to the short-time friction coefficient ζ0 = kBT/D0 in colloidal suspensions (or a suitable kinetic or “Doppler” friction in the case of molecular liquids38,39).
More crucially, the factor b(t) in Eq. (5) is a non-equilibrium time-dependent mobility function29 that can be identified, in the NE-SCGLE framework, with the normalized instantaneous long-time self-diffusion coefficient b(t) ≡ DL(t)/D0. It is the crucial ingredient in the theory of aging, as it describes the kinetically slowed down evolution of structural quantities. It needs to be calculated by a suitable theory of (near-stationary) structural relaxation and kinetic arrest; in this context, b(t) is typically evaluated within the SCGLE framework, although a suitable extension of mode-coupling theory (MCT) would give qualitatively similar results.
Equations (1) and (5) bear an obvious resemblance: identifying S(k; t) as a length-scale dependent (k-dependent) observable X(t) within Tool’s model and as the asymptotic value X∞, we observe that NE-SCGLE encodes the mathematical structure of the generic relaxation postulated by Tool. In particular, a change of variable from time t to
allows us to write the solution of Eq. (5) as S(k; t) = S*(k; u(t)) in terms of the function S*(k; u) solving the linear differential equation
with a k-dependent relaxation factor and initial condition S*(k; u = 0) = S(k; t = 0) = S(i)(k); thus,
Hence, the NE-SCGLE proposes an unambiguous microscopic identification of the relaxation rate [1/τ(X(t))] corresponding to the variable X(t) = S(k; t). Remarkably, this factorizes, 1/τ(k; t) = α(k)b(t), and thus, the corresponding “material” time equally factorizes,
The quantity u(t) can be regarded as an “intrinsic” material time that governs the relaxation of density fluctuations on all length scales.
In this contribution, we aim to link the two different approaches to aging just highlighted—representative of a more practical-oriented materials-science and engineering methodology and of the statistical-physics method to predict macroscopic behavior from microscopic intermolecular forces—by exploiting the mathematical similitude between the phenomenological Tool model described by Eqs. (1)–(4) and the NE-SCGLE framework of Eqs. (5)–(9). We argue that this bears a deep physical significance as it reflects a remarkable connection in their fundamental roots, namely, the universal principles of Onsager’s theory of irreversible processes.
Relating classical and irreversible thermodynamics to glass science has, indeed, a long tradition, as illustrated since the 1950s by Davies and Jones,40,41 in the 1970s by Moynihan et al.7,8 and as discussed more recently by Bailey et al.42 In fact, the elegance and simplicity of the TNM model rest on its fundamental roots in the linear laws of irreversible thermodynamics, but it is notably augmented by Tool’s assumption that the relaxation times vary with time through their dependence on the state variables themselves. In irreversible-thermodynamic terms, this means that the kinetic coefficients of the transport equations (such as viscosity, conductivity, and diffusion coefficients) are treated as state functions. The empirical success of the TNM model clearly demands the theoretical investigation and prediction of the corresponding “kinetic equations of state,” which the NE-SCGLE theory incorporates in the very heart of its formulation.
As we will discuss below, the NE-SCGLE supports the definition of an intrinsic material time u(t) that is generic to a wide range of observables and, at the same time, suggests an extension of the TNM concept toward materials characterized by multiple intrinsic times, namely, those where fundamentally different modes of relaxation exist. We also hope to trigger by our discussion the reconsideration from a microscopic starting point of various phenomenological concepts of enormous practical and fundamental impact,14 whose discussion has been the subject of both historical arguments and enlightening reviews. Illustrative examples include the concepts of time–temperature superposition (TTS) and thermo-rheological simplicity (reviewed in Refs. 43 and 44) and the notions of fictive temperature and effective time, introduced decades ago by Struik.9 In close relationship with the fundamental description of Kovacs’ kinetic signatures in glass-forming systems and of the concept of fictive temperature4,10,11 (see Ref. 45 for a recent discussion), the present work is meant to advance along this program by focusing on the fundamental relationship between the empirically successful TNM model and the microscopic foundation of the NE-SCGLE formalism.
We start in Sec. II by discussing the Tool model as a manifestation of the linear laws of irreversible thermodynamics, modified by TNM’s recognition that the kinetic coefficients may be considered state functions and, hence, depend on the current value of the state variables. These linear laws, applied to several coupled state variables, naturally extend the original Tool’s one-relaxation time model5 to the multi-variable model proposed by Narayanaswamy.6 Section III addresses the problem of defining and predicting the kinetic equations of state, i.e., the dependence of the kinetic coefficients on the state variables. The route proposed by NE-SCLGE theory rests on the universal laws of irreversible thermodynamics, although with several important additions.29 The most important among them is the central role played by the kinetic equation of state L = L[X], which the resulting non-equilibrium microscopic theory predicts in an approximate, but systematic manner. For clarity, in Subsection III B, we provide an illustrative example on the use of the NE-SCGLE to calculate the intrinsic material time of a glass-forming model system. In Sec. IV, we show how the NE-SCGLE framework allows us to extend TNM’s concept of material time for the more sophisticated case of glass forming systems involving non-spherically interacting particles, leading to multiple relaxation modes originated in the dynamical coupling between the translational and orientational degrees of freedom. In Sec. IV, we also illustrate the use of this extended model for the description of aging phenomena in a model dipolar fluid. Finally, in Sec. V, we summarize our conclusions.
II. THE TOOL MODEL AND THE GENERAL LAWS OF IRREVERSIBLE THERMODYNAMICS
Let us start by reviewing the TNM framework from the perspective of the linear laws of irreversible thermodynamics, modified upon the recognition that the corresponding kinetic coefficients are state functions. We, then, explain how such linear laws applied to several state variables naturally extend Tool’s original model5 to the multi-variable model of Narayanaswamy.6
A. Linearity between “fluxes” and “forces” for one variable X(t)
Tool’s Eq. (1) can be understood in a simple manner as a direct extension of the universal laws of irreversible thermodynamics.46 For this, let the variable X(t) represent an extensive thermodynamic property. The most relevant of these laws states that the relaxation of this extensive variable toward its equilibrium value, Xeq, is governed close to equilibrium by a linear relationship,
between the “flux” [i.e., the time rate of change of X(t)] and the corresponding “thermodynamic force” ΔF ≡ F[X] − F[Xeq], where the function F[X] is the respective conjugate variable, F[X] = dS[X]/dX, with S[X] being the entropy47 and L being the corresponding Onsager kinetic coefficient. This fundamental law is one of the pillars of the TNM model.
The fact that “flux” and “force” are linearly related, however, does not imply that Eq. (10) is a linear equation. One obvious reason is that the thermodynamic equation of state F = F[X] does not involve a linear relationship between ΔF and ΔX ≡ X − Xeq. Such “thermodynamic” non-linearity is expected and well understood, and Boltzmann’s fundamental postulate S[X] = kB ln W[X]48 allows us, in principle, to determine the precise function F[X] = dS[X]/dX. A less obvious non-linearity of Eq. (10), however, arises if one allows for the proportionality “constant” L to be a state function L = L[X]. In contrast to the thermodynamic equation of state F = F[X], no general fundamental principle analogous to Boltzmann’s postulate is available to assure the existence of or to determine the “kinetic equation of state” L = L[X] from first principles.
B. Equilibrium vs kinetically arrested stationary states
The very postulate of existence of a kinetic equation of state L = L[X] already offers a glimpse of one of the most fascinating predictions of Eq. (10). Written as
this equation assures us that any equilibrium state Xeq is a stationary solution since in the limit t → ∞, the factor on the right-hand side vanishes whenever X(t → ∞) = Xeq. The inverse statement, however, is not true. In other words, not every stationary solution of this equation corresponds to an equilibrium state. The reason for this is the possibility that the kinetic coefficient L[X(t)] asymptotically vanishes at long times (i.e., L[X(t → ∞)] = 0), thus establishing conditions for stationarity without the equilibrium condition ΔF = 0 being required. Thermodynamic equilibrium states have long been understood as a consequence of the fundamental and general laws of classical47 and statistical48 thermodynamics. In contrast, the second class of stationary solutions, denoted by Xa, in which the transport mechanisms represented by L[X] become arrested (L[Xa] = 0), find their vivid physical realization in the formation of rather common non-equilibrium states of matter, such as glasses and gels. Demonstrating the existence of these states is one of the central aspects of the NE-SCGLE theory.
The notorious asymmetry in the understanding of equilibrium vs non-equilibrium stationary states of matter derives from the fact that equilibrium states can be understood solely in terms of a maximum entropy principle combined with Boltzmann’s postulate, without having to appeal to the fact that they are stationary solutions of any non-equilibrium process. This is in stark contrast with non-equilibrium stationary states, which can only be understood as stationary solutions of non-linear transport equations of the kind in Eq. (11).
C. A natural extension of the TNM model
The full significance of Onsager’s principle and its non-equilibrium extension comes to light in the multi-variable generalization of Eq. (10). Consider ν thermodynamic extensive variables X1(t), X2(t), …, Xν(t) grouped as the components of a vector X(t) ≡ [X1(t), X2(t), …, Xν(t)], and let S = S[X] be the fundamental thermodynamic relation.47 Then, Fi[X] ≡ (∂S[X]/∂Xi) are the corresponding intensive conjugate variables, and to first order,46
where ΔX(t) ≡ X(t) − Xeq, ΔF[X(t)] ≡ F[X(t)] − F[Xeq], and L is the matrix of Onsager kinetic coefficients.
For small enough derivation from equilibrium, the linearization of F[X] in [X(t) − Xeq] yields
where the elements of the matrix are the thermodynamic susceptibilities .
This linearization suggests a natural extension toward a multi-variable Tool model, based on the ad-hoc assumption at L and can be considered as kinetic state functions to describe the evolution of an observable X(t) toward is asymptotic value X∞,
with the relaxation-rate matrix .
Equation (14) is formally solved as
where the time-ordered exponential is the unique solution of the differential equation ∂tE(t, 0) = −τ−1(t)E(t, 0) to the initial condition E(0, 0) = 1. In principle, the matrix τ−1 can be diagonalized, and it will have non-negative eigenvalues because it is the product of a symmetric matrix L and a matrix that for reasons of thermodynamic stability is assumed to be positive (semi-)definite. The individual modes of relaxation in the set of observables Xi that correspond to the eigenbasis of the relaxation-rate matrix, then, follow
Note the subtle difference to Eq. (1): each of the corresponding relaxation rates now is a state function that, in general, depends on all the variables grouped in the vector X(t). In addition, the typical experimental setup prescribes a set of Xi that do not correspond to an eigenbasis of τ−1 so that, in general, the individual Xi(t) given by the matrix exponential in Eq. (15) do not decay as simple exponentials, but rather as superpositions of exponentials. This addresses Narayanaswamy’s concern6 that Tool’s original model5 in Eq. (1) “postulates a simple exponential relaxation mechanism governed by a single (time-dependent) relaxation time. In fact, most glasses exhibit a more complex nonexponential relaxation behavior that can, however, be characterized by a distribution or spectrum of relaxation times.” In addition, there need not exist a single scalar material time because different modes of τ−1 can have very different temperature dependence.
The formal solution, Eq. (15), suggests to define a “material-time matrix,”
As we exemplify below, for the discussion of the long-time fate of aging systems, in the limit t → ∞, we essentially distinguish between the following cases: if τ−1 approaches a matrix that is similar to a strictly positive definite matrix, asymptotically, ζ(t) ∼ ζ∞t, i.e., the material time becomes asymptotically linearly related to physical time, but with a non-trivial prefactor matrix ζ∞ that reflects the slow dynamics of the equilibrating system. In the second case, when the relaxation of all relevant modes of the system becomes kinetically arrested, τ−1 approaches the zero matrix, and in this case, the material is characterized by an asymptotically constant set of material times, which can be expressed as a set of eigenvalues ζi (or equivalently, a set of finite intrinsic material-time values ui). In this case, the physical aging of the system forever “gets stuck” at a finite internal clock-time.
The multi-variable generalization of the TNM model displays a third possibility: in systems with kinetic arrest but partial equilibration of some of the dynamical modes, the relaxation-rate matrix τ−1(t) may have blocks with different temperature dependence and, in particular, may develop zero-blocks for large t while remaining non-zero, in general. In this case, the matrix ζ(t) will have sub-blocks that grow asymptotically with t and ones that approach constants, i.e., a mixture between the two cases described above that are already present in the single-variable TNM model. As discussed below, this situation arises, for example, in the kinetic arrest of a dipolar hard-core fluid, where translational degrees of freedom display kinetic arrest, while the rotational degrees of freedom can remain ergodic and, thus, equilibrate. Other systems where such partial equilibration is known are mixtures of very disparate particles. The partial equilibration, then, manifests as a subset of material-time eigenvalues that remain stuck at finite values and a separate subset that grows indefinitely with physical time.
III. THE TNM MODEL AND THE NE-SCGLE THEORY FOR SIMPLE LIQUIDS
Section II dealt with the TNM model of aging, viewed as the phenomenological manifestation of the universal laws of irreversible thermodynamics. In this section, we now shift our perspective to that of the NE-SCGLE theory. We will focus, first of all, on the identification of the functions b(t) and u(t) appearing in Eqs. (5) and (6) as, respectively, the kinetic coefficient and intrinsic material time of the variable S(k; t), as determined from the solution of the NE-SCGLE system of equations.
A. The NE-SCGLEs
The main conceptual link between Tool’s Eq. (1) for a generic property X(t) and the relaxation equation [Eq. (5)] for the structure factor S(k; t) is the mobility function b(t) whose time integral is an intrinsic material time. Recall that, in essence, the mobility b(t) is the long-time self-diffusion coefficient and, hence, plays the role of the kinetic coefficient of this particular problem. Thus, in the context of the NE-SCGLE theory, it is a state function, which means that it depends on t through its functional dependence on S(k; t), a dependence that shall be denoted by b(t) = b[S(t)].
It has widely been recognized that the static structure factor is a “master property” of structural relaxation in the sense that—at least within the usual mean-field description of the glass transition as proposed by MCT or SCGLE—S(k) governs the long-time dynamics. Extending this notion to the non-equilibrium aging dynamics, we consider the kinetic coefficients in the TNM-like models to be functionals of S(k; t). In particular, within the NE-SCGLE, relaxation equations to determine b[S(t)] have been derived as a straightforward extension of the equilibrium SCGLE theory; a similar program could be carried out on the basis of MCT. In essence, the qualitative description given by both theories concerning the asymptotic structural-relaxation dynamics for long times will be identical.49
Within the NE-SCGLE, the mobility b(t) is expressed through a generalized friction kernel,
For the t-evolving and τ-dependent friction kernel Δζ*(τ; t), the following approximate expression is derived:29
where F(k, τ; t) ≡ N−1⟨δn(k, t + τ)δn(−k, t)⟩ is the non-equilibrium collective intermediate scattering function (ISF), whose self-part, FS(k, τ; t), is defined as , with being the displacement of a tracer particle.
These equations are complemented with the corresponding memory-function equations for the ISFs F(k, τ; t) and FS(k, τ; t), written approximately in terms of their Laplace transforms (LTs), F(k, z; t) and FS(k, z; t), as
where and kc is an empirically determined parameter.31
Equations (18)–(21) constitute a closed system of equations that summarizes the NE-SCGLE theory, together with Eq. (5). They determine the kinetic coefficient b(t) as a (rather involved) functional b(t) = b[S(t)] of the non-equilibrium static structure factor S(k; t) and, thus, constitute a kinetic equation of state associated with the variable S(k; t). Once the thermodynamic input has been determined, then the solution of these equations simultaneously yields predictions for S(k; t) and for the dynamic properties b(t), Δζ*(τ; t), F(k, τ; t), and FS(k, τ; t).
B. Illustrative concrete calculation of the intrinsic material time
The previous concepts are better illustrated by means of a concrete example involving the solution of the NE-SCGLEs. Despite their seemingly complex appearance and non-linear character, the mathematical features of Eq. (5)—and particularly Eq. (8)—render their numerical solution relatively manageable. The reader is referred to Ref. 30 for a thorough and detailed explanation of the numerical methods involved and to Ref. 50 for access to a code with the corresponding computational program. Here, we only illustrate the results for the glass-forming model liquid also studied in Ref. 30, which consists of N identical spherical particles interacting through short-ranged and purely repulsive interactions. We refer to the Weeks–Chandler–Andersen (WCA) model,51 whose pair potential U(r) vanishes for r larger than the distance σ of soft contact, but which for r ≤ σ is given by
The state space of this system is spanned by the reduced number density and the reduced temperature T* ≡ kBT/ϵ, here denoted simply as n and T, so that the volume fraction is ϕ = πn/6. For the determination of the thermodynamic property entering in Eq. (5), we employ the same approximation as in Ref. 30 based on the Percus–Yevick/Verlet–Weis (PY/VW) approximation52–54 for the hard-sphere liquid with an effective hard-sphere diameter σHS(T) determined by the blip function method.55
The solution of the NE-SCGLEs describes the evolution of the properties S(k; t), b(t), Δζ*(τ; t), F(k, τ; t), and FS(k, τ; t) for times t > 0 after the system, initially at equilibrium at temperature T(i) and volume fraction ϕ, is instantaneously and isochorically quenched (at t = 0) to a final temperature T(f). Here, however, we shall focus only on the behavior of the mobility function b(t) and on its integral u(t), the intrinsic material time of the variable S(k; t). As discussed in Ref. 30, the mobility function actually depends on both the initial and final temperatures, T(i) and T(f), and on the volume fraction ϕ of the quench, b(t) = b(t; ϕ, T(i), T(f)). Its long-time asymptotic limit ba(ϕ, T(f)) ≡ b(t → ∞; ϕ, T(i), T(f)), however, is independent of T(i) and serves to identify two qualitatively different scenarios predicted by the NE-SCGLE theory for the relaxation of the system after quenching.
The two scenarios are (i) equilibration, for which ba(ϕ, T(f)) is identical to its finite and positive equilibrium value, ba(ϕ, T(f)) = beq(ϕ, T(f)) > 0, thus defining a condition to predict whether the system would be able to reach a new equilibrium state, and (ii) aging, where ba(ϕ, T(f)) vanishes, and hence, the system is predicted to age endlessly. As illustrated in Fig. 1(a), using this criterion, the state space (ϕ, T) of the WCA system can be divided in two regions delimited by the dynamic arrest (solid) line T = Tc(ϕ), above which the system will equilibrate (ergodic region) and below which it will become arrested (non-ergodic or glass region). Let us mention in passing that full kinetic arrest at Tc is a mean-field simplification, as is well understood in the context of both MCT and SCGLE. For the present qualitative discussion, we do not consider additional relaxation modes that might be active in the idealized glass.
In the inset of Fig. 1(a), we plot the NE-SCGLE results for the asymptotic value ba(ϕ, T(f)) as a function of T(f) along the isochore ϕ = 0.66. These results determine that, at this isochore, the dynamic arrest temperature is Tc(ϕ = 0.66) ≈ 0.0915. At Tc(ϕ), the kinetic parameter ba(ϕ, T) passes continuously from its vanishing value below to its equilibrium non-zero value beq(ϕ, T) above. As also illustrated in the inset of Fig. 1(a), when T(f) approaches Tc from above, beq(ϕ, T) vanishes as with exponent μ = 2.2.
Figure 1(a) also shows the illustrative isochore ϕ = 0.66 and indicates with two downward arrows two representative isochoric quenches along this isochore. Both quenches start at the same initial temperature, T(i) = 0.04, but end at different final temperatures T(f). The shallower quench (T(f) = 0.2 > Tc) represents a typical equilibration process, whereas the deeper quench (T = 0.07 < Tc) is meant to illustrate a process of dynamic arrest. In Figs. 1(b) and 1(c), we plot the full non-equilibrium evolution of, respectively, the kinetic parameter b(t; ϕ, T) (in terms of its inverse) and the intrinsic material time u(t; ϕ, T(f)) as a function of the waiting time t, indicating with the thick solid and thick dashed lines these two illustrative quenches of Fig. 1(a). Figures 1(b) and 1(c) also include other similar quenches (thinner lines) to emphasize the crossover from the equilibration regime [solid lines, T(f) > Tc(ϕ)] to the aging regime [dashed lines, T(f) < Tc(ϕ)]. The thick dotted line corresponds to a quench to the critical temperature Tc(ϕ) ≈ 0.0915.
These results exhibit some relevant predictions of the NE-SCGLE theory whose detailed analysis, however, lies outside the scope of this article. Instead, we refer the reader to Ref. 30 for a detailed analysis of results analogous to those summarized in Fig. 1(b), but in which the various lines correspond to quenches with the same initial and final temperatures (Ti = 0.1 and T = 0), along different isochores. Here, we only highlight the predicted existence of the two regimes: equilibration for quenches with T(f) > Tc(ϕ) [solid curves in Figs. 1(b) and 1(c)] and aging for T(f) < Tc(ϕ) (dashed curves). In the equilibration regime, we also note that b−1(t; ϕ, T(f)) reaches its equilibrium value 1/beq(ϕ, T(f)) after an equilibration time-scale teq(T(f)), highlighted by the black circles in Fig. 1(b). The equilibration time teq grows asymptotically linearly with the relaxation time, teq ∼ 1/beq; both quantities, furthermore, increase strongly and eventually diverge as T(f) approaches the critical value Tc from above.
Figure 1(c) displays the corresponding evolution of the intrinsic material time . One observes that u(t) exhibits a short-time behavior given by u(t) ≈ bi × t, which is determined by the initial condition bi ≡ beq(ϕ, T(i)). At longer times, u(t) shows a crossover that can be easily understood in terms of the long-time behavior of b(t). For T(f) > Tc, for instance, the long-time limit of u(t) would be approximately given by u(t) ≈ A + beq(ϕ, T(f))t, where and teq(T(f)) is the corresponding equilibration time-scale. Hence, the material time u(t) continues to grow as a linear function of time when the system has reached equilibrium. For T(f) < Tc, instead, b(t) → 0 as t → ∞, so u(t) reaches an arrested plateau value , with ta being the time scale required by b(t) to reach its asymptotic arrested value ba = 0, which becomes shorter with decreasing T(f). The value ua depends on the specific quench and, hence, serves as a nonequilibrium state parameter, augmenting the set of thermodynamic control variables that are needed to uniquely specify a material.
C. Physically relevant functionals of S(k; t)
Let us mention in passing that the NE-SCGLE, in the spirit of a mean-field theory of the glass transition, expresses a number of physically relevant properties as functionals of S(k; t) so that these can, in principle, be easily evaluated on the basis of the results discussed above. Examples include the internal energy E(t) and the pressure p(t) of a simple liquid with only pair interactions whose configurational contributions ΔE(t) and Δp(t) are given, respectively, by the so-called “energy” and “pressure” equations,
In these equations, the non-equilibrium radial distribution function g(r; t) can be related with S(k; t) just as it is typically done in equilibrium conditions, i.e., [S(k; t) − 1]/n is the Fourier transform of [g(r; t) − 1]. Remarkably, these relationships, which are well known from equilibrium statistical physics, can be shown to also hold under non-equilibrium conditions by following their standard derivation,48,55 starting from the definition of g(r; t) as a statistical average, but avoiding the use of the canonical (or any other equilibrium) ensemble. The discussion of these aspects is beyond the scope of the present contribution and will be reported elsewhere.
Similarly, many other relevant dynamical, transport, and rheological properties may be written at least approximately in terms of S(k; t) or of the dynamical scattering functions, which through Eqs. (18)–(21) again depend on S(k; t). A simple example is provided by the mean squared displacement W(τ; t) ≡ ⟨(r(t + τ) −r(t))2⟩/6, which may be written within the NE-SCGLE in terms of the friction function Δζ*(τ; t) through the following exact equation:
An equally important example is provided by the non-equilibrium rheological or viscoelastic properties of the system, which may be represented by the complex-valued dynamic shear modulus G(ω; t) or by the dynamic shear viscosity η(ω; t). Either of these properties determines the macroscopic stress induced in a viscoelastic liquid upon the application of a low-amplitude oscillatory shear strain of frequency ω.56 These two properties are related by G(ω; t) = iωη(ω; t). An approximate determination of these properties may be provided by the non-equilibrium Stokes–Einstein relation,57 ζ(ω; t) ≈ 3πσ × η(ω; t), where ζ(ω; t) is the Fourier–Laplace transform of ζ(τ; t) ≡ ζ0[2δ(t) + Δζ*(τ; t)], with ζ0 being the short-time friction coefficient ζ0 ≡ kBT/D0. Thus, η(ω; t) may be approximated by
where η0 ≡ ζ0/3πσ is the short-time (or infinite frequency) viscosity. In its zero-frequency limit, this equation leads, using the definition of b(t) in Eq. (18), to the approximate identification of the non-equilibrium static viscosity η(t) ≡ η(ω = 0; t) with the inverse of the mobility function b(t), η(t) ≈ η0/b(t). This means that the non-equilibrium linear viscosity η(t) is also expected to exhibit the same qualitative behavior as the results for b−1(t; ϕ, T) in Fig. 1(b), which can also, then, be read as the predicted scaled viscosity η(t; ϕ, T)/η0.
IV. THE TNM MODEL FOR SYSTEMS WITH ROTATIONAL–TRANSLATIONAL DYNAMIC COUPLING
The simple glass-forming model liquid discussed in Sec. III B, formed by identical spherical particles, is, of course, an idealization. It allows us, however, to illustrate important features of the NE-SCGLE description of aging when a single relaxation mechanism is governed by one time-evolving relaxation time. A more interesting situation arises in systems where multiple relaxation mechanisms are active. This is the case, for instance, in dense and highly asymmetric colloidal mixtures of hard spheres58 and also in suspensions where the pair interactions between colloids depend on particle’s orientations.59–61 An interesting question that arises in this context, thus, is if the application of the NE-SCGLE theory to a simple model that incorporates, for example, rotational–translational dynamical coupling (or, alternatively, coupling in the dynamics of multiple species) could also serve to extend TNM’s model of aging to such more complex conditions.
With this idea in mind, in this section, we shall focus on the application of the NE-SCGLE framework to the description of aging in a model system whose interparticle interaction is not radially symmetric. We refer to a dipolar hard-sphere (DHS) fluid,55 formed by N axially symmetric identical particles in a volume V, interacting through a hard sphere repulsion of diameter σ, which shall be denoted as UHS(r12), plus a dipole–dipole interaction so that the total pair potential reads
where r1 and r2 are vectors describing the positions of the particles, r12 ≡ (r2 − r1), r12 = |r12|, and and and are unit vectors indicating the particles’ orientations so that the dipolar tensor contains all the orientational features of the interaction. A detailed characterization of this system has been provided recently in the context of the NE-SCGLE theory.62 Here, we only bring together some of the most important elements in the context of our present discussion.
A. Non-spherical NE-SCGLE theory
The non-spherical extension of the NE-SCGLE theory was formulated in Refs. 63 and 64 in terms of a set of equations describing the non-equilibrium evolution of the structural and dynamical properties of a system represented by Eq. (27) and has been systematically employed62 to understand the dynamically arrested states already detected in extensive molecular dynamics (MD) simulations of this model.65,66 We summarize the relevant ingredients of the non-spherical version of the NE-SCGLE theory in the Appendix.
Similar to the spherical case discussed previously, the most central equation of the extended NE-SCGLE describes the time-evolution of the non-equilibrium structure factor, which for liquids involving non-spherical interactions can be defined as . Here, is the Fourier transform of the fluctuation in the local number density of particles at the position r and with orientation at time t.67 In practice, however, one focuses on the coefficients of the spherical harmonic expansion of ,63,65–68 commonly referred to as tensorial density modes and denoted as [see Eq. (A1)]. These coefficients, in turn, allow us to define the corresponding projections of , namely,65,66 . From symmetry considerations and within well-defined approximations,62 only a finite subset of its isotropic and homogeneous diagonal components, denoted as Slm(k, t) ≡ Slm:lm(k, −k, t), are involved in the central NE-SCGLE, which reads64
Just as its spherically symmetric counterpart, Eq. (5), this equation describes the non-equilibrium isochoric response of the DHS system to an instantaneous quench at t = 0 from an arbitrary initial temperature T(i) to a final value T(f), while the system is constrained to remain spatially homogeneous and isotropic, with fixed bulk number density n = N/V. Alternatively, it also describes the isothermal response after an instantaneous compression from an initial volume fraction ϕ(i) to a final value ϕ(f).
In Eq. (28), and are, respectively, the translational and rotational free-diffusion coefficients of the constituent dipolar particles. The functions , on the other hand, are spherical harmonic projections of the second functional derivative of the free energy with respect to the local number density , evaluated at the final state point of a quench, i.e., (see Subsection 2 of the Appendix). As mentioned before, this intrinsically thermodynamic property is an externally provided input for the theory, which, for the DHS model, can be determined following the steps outlined in Refs. 62 and 66. The other relevant elements appearing in Eq. (28) are the time-dependent translational and rotational mobility functions bT(t) and bR(t). They are calculated with the non-spherical extension of the NE-SCGLE theory summarized in Eqs. (A4)–(A9). We highlight in the following the physical notions underlying the extended Tool–Narayanaswamy–Moynihan model suggested by the non-equilibrium evolution of the tensorial components Slm(k; t) in Eq. (28).
B. The rotational–translational Tool–Narayanaswamy–Moynihan model
One simple manner to identify an extended TNM model, inspired in the NE-SCGLE theory for fluids with coupled rotational and translational dynamics, is again based on the recognition of the mathematical structure of Eq. (28), which has exactly the same structure of the generic relaxation equation [Eq. (1)] postulated by Tool, namely,
with the relaxation time τlm(t) given by
which can also be written as
This immediately allows us to identify the corresponding TNM material time ζlm(k; t) for the relaxation of each tensorial component Slm(k; t) as
i.e., a linear combination of the translational intrinsic material time and a corresponding rotational intrinsic material time,
It is also straightforward to show that the solution of Eq. (29) can be written as with the function defined as
C. Equilibrium and non-equilibrium stationary solutions and arrested states diagram
Before considering specific results, let us highlight that Eq. (28) already reveals remarkable features regarding the possible non-equilibrium states of the DHS model system defined in Eq. (27) (and of any other that can be generically represented by this equation). The functional dependence of both bT(t) and bR(t) on the spherical harmonic projections Slm(k; t) [see Eqs. (A4)–(A7)] introduces strong non-linearities and, more crucially, implies the existence of two kinds of stationary solutions where ∂Slm(k, t)/∂t = 0, which are both fundamentally different from equilibrium states.
As mentioned before, we can distinguish ordinary thermodynamic equilibrium conditions, where stationarity is attained because the factor on the right-hand side of Eq. (28) vanishes, i.e., because all the spherical harmonic components Slm(k; t) are able to reach their corresponding thermodynamic equilibrium values, , at a finite t. This corresponds to quenches to states where both mobilities bT(t) and bR(t) attain finite long-time values and .
Other categories of stationary solutions of Eq. (28) are those for which the asymptotic value of the kinetic factor vanishes for t → ∞. Due to the structure of Eq. (28), this condition can be fulfilled if either both bT(t) and bR(t) vanish asymptotically or, for l = 0, only bT(t) vanishes since βlm(k) ∝ l vanishes in this case. These two possibilities, in fact, correspond to totally or partially arrested states of matter as predicted by the NE-SCGLE and observed in extensive molecular dynamics simulations:65 in the fully arrested state, both mobilities vanish, while in mixed-glass states, rotational degrees of freedom equilibrate, and thus, bR(t) never vanishes.
Thus, this simple analysis of the possible classes of stationary solutions of Eq. (28) already announces the possibility of only two qualitatively distinct non-equilibrium amorphous states for the DHS model. As just mentioned, the existence of these states has been confirmed recently with the assistance of extensive MD simulations.65 In fact, the use of the asymptotic long time solutions of the non-spherical NE-SCGLEs (see Subsection 4 of the Appendix) allows for the development of the corresponding arrested state diagram of the DHS fluid shown in Fig. 2(a). This diagram identifies the region in the (ϕ, T) plane (with T ≡ kBTσ3/μ2) corresponding to ergodic states, as well as two regions corresponding to the two aforementioned dynamically arrested states, and the boundaries between these regions. The different arrows in Fig. 1(a) illustrate two qualitatively distinct processes that could be used to bring the DHS system to each of the two arrested states, namely, isochoric quenches toward a fully arrested state (vertical arrows) and isothermal compressions toward partially arrested states (horizontal arrows). As in Fig. 1, solid arrows correspond to equilibration processes in which the final states also lie in the ergodic region (Q1 and C1), whereas the dotted arrows illustrate processes toward dynamically arrest states (Q2 and C2).
The full solution of the NE-SCGLE (outlined in Subsection 3 of the Appendix) provides the isochoric (or isothermal) evolution for t > 0 of the physical properties involved, for example, the projections Slm(k; t), the corresponding dynamic correlations Flm(k, τ; t), and , and also the friction memory kernels and . The detailed analysis of such of properties, however, is beyond the scope of this contribution and will be addressed elsewhere. Here, we only focus on the non-equilibrium behavior of the mobility functions bT(t) and bR(t), which, according to our previous discussion, define the kinetic coefficients , and also on the behavior of the corresponding intrinsic material times uT(t) and uR(t), which allow us to identify the corresponding TNM material time ζlm(k, t) of each projection Slm(k; t).
Figure 2(b) displays the evolution of the two relaxation times (upper panel) and (lower panel) for the isochoric quenches Q1 (thick solid lines) and Q2 (thick dashed lines) shown in Fig. 2(a), both starting from the same initial equilibrium state (ϕ = 0.4, T(i) = 1) and ending at distinct final state points (ϕ = 0.4, T(f) = 0.3) and (ϕ = 0.4, T(i) = 0.05), respectively. Thus, these processes consider state points that lie above and below the transition temperature Tc(ϕ = 0.4) ≈ 0.075. In Fig. 2(b), we also include results for other similar quenches in order to highlight the crossover from equilibration (solid lines, T(f) > Tc) to aging (dashed lines, T(f) ≤ Tc). One notes first that, qualitatively, both and display the same behavior, i.e., both become increasingly larger with the decreasing final temperature and both reach their corresponding equilibrium values, and , at a very similar time-scale, which increases with T(f). As emphasized by the black circles in the upper and lower panels, the equilibration times and follow a rather similar trend, increasing asymptotically linearly with the equilibration time, .
For quenches at-or-below Tc, instead, one observes an unbounded increase in both and , thus suggesting their simultaneous divergence as t → ∞. This indicates a gradual and simultaneous decrease in the translational and rotational diffusion (i.e., an increase in the viscosity of the system) and describes, thus, the formation of a fully arrested state in which both the orientational and translational dynamics slow down and undergo aging during the endless process of reaching a glassy state. Let us also mention that an analogous behavior is found for the two α-relaxation times, and , associated with the functions and , which describe the relaxation of, respectively, the translational and orientational dynamics.
These features are to be contrasted against the behavior that is found for the mobility functions upon isothermal compressions (such as C1 and C2) toward the region of partially arrested states. As shown in Fig. 2(c), for these processes, one observes a clear decoupling in both mobilities: again shows a crossover from equilibration (for compressions with final volume fraction ϕ(f) < ϕc = 0.582) to an aging regime (ϕ(f) ≥ ϕc = 0.582), whereas does not change appreciably in the whole volume fraction regime explored and, in fact, reaches an equilibrium value at a remarkably small time-scale. This scenario describes, thus, the formation of a partially arrested state, in which translational diffusion slows down and ages as in the process of formation of a hard-sphere glass, while the orientational dynamics of the dipolar particles remains ergodic, but in a disordered and only slowly evolving configuration of positions.
D. Time evolution of the intrinsic material times and physical trajectories on the plane (uT, uR)
where and . One requires, of course, that these two expressions yield the same value of the actual (i.e., real) waiting time, t, so that the differential form
must be satisfied along a (real) physical trajectory [uT(t), uR(t)] in the (uT, uR)-plane.
Thus, for each process considered (quench or compression), the system traces a corresponding physical trajectory of points [uT(t), uR(t)] on the material-time plane (uT, uR) as the waiting time evolves from t = 0 toward t → ∞. Some general features of these trajectories can be implied from the inspection of the possible long-time limits, and .
For instance, if a given protocol drives the system to an arbitrary equilibrium state (e.g., Q1 and C1), then both bT(t) and bR(t) eventually reach finite equilibrium values and within a finite waiting time window. Then, Eqs. (33) and (34) imply that both material times will evolve without bound as t runs from 0 to infinity so that the physical trajectory will start at the origin, (uT, uR) = (0, 0), and will tend to (uT, uR) = (∞, ∞). At asymptotically long times, the physical trajectory in the (uT, uR) plane will be a straight line with a prefactor that is determined by the specific quench. This is illustrated by the blue and purple solid lines in Figs. 3(c) and 3(f).
If, instead, the protocol drives the system to a fully arrested state (e.g., Q2), both asymptotic values and will be zero. In this scenario, Eqs. (33) and (34) imply that the real trajectory starts again at (uT, uR) = (0, 0) but now ends at a point , with and . This case is illustrated by the dashed line in Fig. 3(c).
Finally, a third possibility is a partially arrested state where vanishes, but attains an equilibrium value. In this case, uT(t) becomes finite, while uR(t) increases without bound; see the dashed line in Fig. 3(f).
We have discussed the concept of a material time in the context of aging in glasses, i.e., of the slow evolution of the physical properties of a system after a quench in control parameters. In particular, we have shown how the empirical Tool–Narayanaswamy–Moynihan model of aging can be seen as arising from the mathematical structure of equations that is suggested by the non-equilibrium extension of the self-consistent Langevin equation theory, NE-SCGLE. A key concept in this identification is that of a kinetic state function. Specifically, the theory suggests a structure where the relaxation rate of evolving observables during aging is expressed as a functional of the waiting-time dependent static structure function.
The generic mathematical structure of the TNM relaxation equation allows for two distinct stationary states: those corresponding to equilibration, where upon a quench within the ergodic-liquid state of a system, all relaxation rates remain non-vanishing, and those corresponding to quenches into the (ideal) glassy state, where the relaxation rates vanish. In the two cases, the material time evolves distinctly: upon equilibration, the inner clock of the material evolves asymptotically proportional to the physical time after the quench. For arrested states, the material time attains a finite asymptotic value.
More intriguingly, the extension of the NE-SCGLE to include multiple relaxation modes—exemplified in the present contribution through a model with distinct translational and rotational degrees of freedom—highlights the possibility of a third kind of stationary state, that of partial equilibration. Here, the aging process is characterized by two material times: one of which attains a finite long-time limit, while the other grows asymptotically proportionally with the real waiting time.
In principle, this highlights that a “material time” can only be defined with respect to a given observable, and not all observables need to follow the same internal clock; considering a set of observables that undergo aging, a corresponding set of material times would, in principle, be obtained. However, the NE-SCGLE theory encodes an interesting factorization property: for example, when considering the length-scale dependent relaxation of density fluctuations, the information on the length scale can be absorbed into a wave-number dependent prefactor. We expect, thus, that the aging process of a glass former, in general, is controlled by a (small) set of intrinsic material times ui corresponding to structurally distinct relaxation channels—such as those connected to translational and rotational degrees of freedom in a system with anisotropic interactions or those related to different particle species in mixtures of very disparately sized particles. The final material properties of a quenched glass are, thus, characterized by a vector ua of the asymptotic values of these intrinsic material times (some components possibly infinite in the case of partial equilibration), and the aging process is represented as a one-dimensional curve in the space spanned by the intrinsic material time components, parameterized by the physical waiting time.
We expect that the fundamental and universal laws of irreversible thermodynamics and the microscopic perspective provided by statistical mechanical theories of irreversible processes, such as the NE-SCGLE theory, may be instrumental in extending the descriptive power of the TNM model beyond its current limitations. Similarly, the use of the NE-SCGLE formalism should also enrich the discussion of other phenomenological concepts of great practical and fundamental relevance, such as the concepts of time–temperature superposition (TTS) and thermo-rheological simplicity,43,44 and the notions of fictive temperature and effective time.9 In the absence of fundamental principles that allow for the identification of their microscopic origin, it has been difficult to assess their validity and to establish a clear relationship between them, with the concept of material clock. In this regard, it is important to mention that the TTS scaling for equilibrium liquids is a distinct prediction of mode-coupling theory3,70–72 (and hence, also of the original equilibrium SCGLE theory49). These first-principles equilibrium theories break down in the very high viscosity regime. Future work will address the combination of the NE-SCGLE formalism with extensions of MCT-like theories that account for the avoidance of the ideal glass transition.73 This combination is expected to open a route to extend equilibrium glass-transition theories to the non-equilibrium regime and, hence, to provide a fundamental basis of a non-equilibrium extension of the notion of TTS and of its possible connection with the concept of a material clock.
At the same time, we should also expect that the success of the TNM model and of the other useful phenomenological notions referred to above will serve as a compelling motivation for the conventional statistical mechanical theoretical methodologies to surmount some of its own limitations, particularly those derived from the difficulty in identifying the universal physical principles that govern the macroscopic evolution of matter during the non-equilibrium process of amorphous solidification. In this regard, it will be interesting to reconcile the present description of aging based on the NE-SCGLE and linked to the TNM model, with more elaborate approaches rooted in the non-equilibrium statistical physics of spin glasses. The seminal work by Cugliandolo and Kurchan aside,74 most notably, the extension of MCT proposed by Latz75 offers a separate account of aging phenomena, which is less obviously recast in the form of empirical models. We also leave such discussion for future work.
L.F.E.-A. thanks Professor Dr. Andreas Meyer for helpful discussions and his crucial support to conclude this work and also acknowledges the Consejo Nacional de Ciencia y Teconología (CONACYT, México) for support through a Postdoctoral Fellowship (Grant No. I1200/224/2021). T.V. thanks the Glass and Time group of Roskilde University and, in particular, Jeppe Dyre for their kind hospitality during a research stay in summer 2021. M.M.-N. and R.P.-O. acknowledge the Consejo Nacional de Ciencia y Tecnología (CONACYT, México) for financial support through Grant Nos. CB A1-S-22362 and LANIMFE 314881.
Conflict of Interest
The authors have no conflicts to disclose.
The data that support the findings of this study are available from the corresponding author upon reasonable request.
APPENDIX: NON-SPHERICAL NE-SCGLE THEORY
In this appendix, we summarize the equations that constitute the extension of the NE-SCGLE theory for systems with non-spherical interactions and provide a numerical procedure that allows its solution. In addition, we briefly discuss the asymptotic long-time limit of the NE-SCGLEs, which allow for the determination of the dynamical arrest diagram of Fig. 2(a).
1. Fundamental equations
The extension of the NE-SCGLE to describe systems with non-spherical interactions was carried out in Refs. 63 and 64. As mentioned above, the core of the theory is Eq. (28), which describes the non-equilibrium time evolution of the structural correlations , defined in terms of the spherical harmonic projections, , of the Fourier transform of the local density, . More specifically, these coefficients are defined as63,65–68
Besides the projections Slm(k, t), the NE-SCGLE also deals with non-stationary two-time density correlation functions, also referred to as intermediate scattering functions (ISFs), defined as64
where τ denotes the correlation (delay) time. Note that the equal-correlation-time values of the collective ISF of Eq. (A2) are also the spherical harmonic components of the non-equilibrium SF at the evolution time t, i.e., Flm;l′m′(k, k′, τ = 0; t) = Slm;l′m′(k, k′; t). The dependence on the direction of the vectors k and k′ refers to the possible presence of external fields that would break conditions of global homogeneity and isotropy of space. Restricting to the case of disordered states (i.e., excluding the possibility of non-isotropic and non-uniform conditions, such as in crystalline solids), one might assume that the only relevant components of these functions depend on the magnitude of k and k′ with k′ = −k.
As discussed in Refs. 63–68, if all the orientations of particles, , and of the wave vector, , are described in the so-called k-frame,65,66 one is left with correlation functions of the form Fll′m(k, τ; t) ≡ Flm;l′m′(k, τ; t)δmm′, , and Sll′m(k; t) ≡ Slm;l′m′(k; t)δmm′. Furthermore, and on the basis of previous work for the specific case of the DHS fluid,65,66 one might restrict to consider only diagonal (l = l′) components Flm(k, τ; t) ≡ Fll′m(k, τ; t)δll′, , and Slm(k; t) ≡ Sll′m(k; t)δll′.
Within these considerations, the NE-SCGLE theory is summarized by a set of coupled evolution equations for the waiting-time dependent functions involved, whose solution for t > 0 describes the irreversible evolution of the structural and dynamical properties of an instantaneously quenched DHS liquid. As already mentioned, the most important of such equations is Eq. (28). The mobility functions bT(t) and bR(t) appearing in this equation are state functions, i.e., they are also functionals of the projections Slm(k; t). The exact functional dependence is, in general, unknown, but providing an approximate proposal for the functional dependence of bT(t) and bR(t) on Slm(k; t) is precisely part of the core of the NE-SCGLE approach. According to Ref. 64, these functions are defined as
that is, in terms of the non-stationary τ-dependent friction functions and , for which approximate expressions are provided,29,64,69 namely,
These equations are finally coupled with time-evolution equations that can be obtained for Flm(k, τ; t) and , which can be written in terms of their Laplace transforms Flm(k, z; t), , , and , as
The functions and are cutoff functions defined, respectively, as and , with kc(t) = α × kmax(t) and kmax(t) being the position of the main peak of S00(k; t) and α = 1.305. This ensures that, for radially symmetric interactions, one recovers the original NE-SCGLE theory for liquids of soft and hard spheres (contained as the case l = 0).
Provided that the fundamental thermodynamic inputs have been determined, Eq. (28) combined with Eqs. (A4)–(A9) constitute a closed set of equations for the non-equilibrium properties Slm(k; t), bT(t), bR(t), , , Flm(k, τ; t), and , whose self-consistent solution describes the non-stationary and non-equilibrium structural relaxation of glass-forming liquids formed by non-spherical particles.
2. Thermodynamic inputs
To solve the NE-SCGLEs, the thermodynamic functions must be externally determined. These functions derive from the spherical harmonic expansion of64
which is the second functional derivative of the Helmholtz free energy density functional F[n, T] (in units of the thermal energy kBT = 1/β) with respect to the local number density . In the absence of external fields, F[n, T] can be written, in general, as F[n, T] = Fid[n, T] + Fex[n, T], that is, as the sum of the ideal gas free energy (with Λ being the thermal wavelength48,55) plus the excess free energy, Fex[n, T], arising from the interparticle interactions, which reads
with being the two-particle direct correlation function.55,67 Therefore,
To determine , we evaluate at the uniform density and temperature of the final state. The Fourier transformed spherical harmonics projections of the resulting function read
Just as it happens in the case of the structural and dynamic correlation functions Slm(k; t), Flm(k, τ; t), and , the use of the k-frame simplifies the determination of the spherical harmonics projections in Eq. (A14), where one sticks only with projections of the form cll′m(k; n, T) ≡ clm;l′m′(k; n, T)δmm′. Hence, one can employ one of the available approximations generated in the context of the equilibrium theory of liquids for the determination of the direct correlation function . As shown in Refs. 53, 55, 63, and 65–68, an appealing choice is the use of the so-called mean spherical approximation (MSA) for the DHS fluid proposed by Wertheim.53 Within this approximation, it is straightforward to show that the functions cll′m(k; n, T) also become diagonal in l and l′, so one deals only with projections of the form clm(k; n, T), thus providing the diagonal elements of , namely,
3. Numerical solution of the NE-SCGLEs
To solve the NE-SCGLE system of differential equations in Eqs. (A4)–(A9), we appeal to standard integration schemes once the mobilities and have been determined. More specifically, we can use Eq. (5) to evaluate Slm(k; t) at a discrete sequence of waiting times t(j) (j = 0, 1, 2, …) or, equivalently, at the corresponding sequence of points in the (uT, uR) plane since for each waiting time t(j), Eqs. (33) and (34) assign values and to the material times uT and uR. In practice, we start by giving the initial value of Slm(k; t) as an input to Eqs. (A6)–(A9), whose solution determines and and, through Eqs. (A4) and (A5), also the initial values of the sequence . This is the j = 0 step of the evaluation of the sequence Slm(k; t(j)), bT(t(j)), and bR(t(j)) with (j = 1, 2, …).
In general, at any subsequent step j > 0 corresponding to a waiting time t(j) > 0, we first input the value Slm(k; t(j)) in Eqs. (A4)–(A9), whose solution determines bT(t(j)) and bR(t(j)), among all the other dynamic properties at time t(j) [or, equivalently, at the point ]. The next step starts with the arbitrary choice of a sufficiently small increment . Since Eq. (38) implies that , this choice also determines and Δt(j) ≡ t(j+1) − t(j), thus defining also t(j+1) and the new point along the (real) physical trajectory in the “material-time” plane (uT, uR). The determination of together with Eq. (35) allows us to evaluate , which we can also label as Slm(k; t(j+1)). This allows us to start the next step of the sequence by now giving Slm(k; t(j+1)) as the new input in Eqs. (A4)–(A9).
4. Asymptotic solutions of the NE-SCGLE: a non-equilibrium route for the determination of a dynamical arrest diagram
Under certain circumstances, one might be interested only in describing generic features summarizing the overall physical scenario for the distinct dynamical arrest transitions of the DHS fluid [or of any other model system whose potential of interaction can be described by means of Eq. (27)]. In fact, the similarity between Eqs. (A6)–(A9) and its equilibrium counterpart [cf. Eqs. (46)–(49) of Ref. 63] already suggests to consider the analysis of the solutions and of all the other quantities involved  in the asymptotic long-τ regime, a routinely exercise in both MCT and SCGLE.49 As it happens, in terms of such analysis, both the equilibrium SCGLE63 and the NE-SCGLE turn out to establish essentially the same physical scenario for the case of DHS fluid if one uses the MSA.
with α = T, R, which we may refer to as the non-equilibrium extension of the so-called non-ergodicity parameters (NEPs).49 Rather than labeling these t-dependent properties in terms of the (real) waiting time, we might relabel them in terms of the material times (uT, uR) defined above. Then, the time evolution of the NEP along the interval 0 < t < ∞ obviously will map onto the evolution of the point u(t) ≡ (uT(t), uR(t)) along the physical trajectories just introduced above in Sec. IV B.
involving the two “order parameters”
Equations (A19)–(A23) complemented with Eq. (35) yield a system of closed equations for the two dynamic order parameters and . In turn, the solution of these equations provides the values of and at any point u in the “material time” plane (uT, uR) and, in particular, along the “physical trajectories” discussed in Sec. IV D. As discussed there, the analysis of these real trajectories determines if the quench (or crush) considered, which ends in the (arbitrary) state point (η, Tf), will fall in any of the three scenarios (full equilibration, mixed states, or full arrest) described in Sec. IV C.
Following Refs. 30 and 63, we might proceed as follows: if for a given quench with parameters n, Ti, and Tf we find that both and for an arbitrary trajectory in the plane u = (uT, uR) (0 ≤ uT ≤ ∞, 0 ≤ uR ≤ ∞), we conclude that all the solutions , which implies that the system will be able to equilibrate after the specific quench considered. This condition also implies that the mobility functions and bR(t → ∞) = bR > 0, that is, they are free to reach their corresponding equilibrium finite values, so each particle is delocalized and mobile and free to rotate. In other words, the point (n, Tf) lies in the region of equilibrium (ergodic) states. If, instead, finite values of u exist such that one or both of and remain finite only within the intervals and , with and/or , then the system is expected to get dynamically arrested in a non-equilibrium asymptotic state, whose spherical harmonic components of the structure factor, Slm(k; ua) are, then, given by
which clearly differ from the corresponding equilibrium limit . In addition, one or both and will vanish in this limit stage.
As mentioned before, the parameters , , and depend, in general, on the fixed density, n, and on the initial and final temperatures of the quench (and similarly for the case of isothermal compressions), and such dependence contains the basic information from which the features of dynamical arrest emerge. As it happens, however, if our purpose is only the determination of the boundaries between the regions of ergodicity and dynamical arrest in the control parameters space, the influence of the initial state temperature becomes rather irrelevant since distinct values for Ti lead to the determination of the same boundary. The use of this protocol allows for the determination of the non-equilibrium phase diagram.