Non-equilibrium Theory of Arrested Spinodal Decomposition
José Manuel Olais-Govea, Leticia López-Flores, Magdaleno Medina-Noyola
aa r X i v : . [ c ond - m a t . s o f t ] O c t Non-equilibrium Theory of Arrested Spinodal Decomposition
Jos´e Manuel Olais-Govea, Leticia L´opez-Flores and Magdaleno Medina-Noyola
Instituto de F´ısica “Manuel Sandoval Vallarta” ,Universidad Aut´onoma de San Luis Potos´ı,´Alvaro Obreg´on 64, 78000 San Luis Potos´ı, SLP, M´exico (Dated: July 3, 2018)
Abstract
The Non-equilibrium Self-consistent Generalized Langevin Equation theory of irreversible relax-ation [Phys. Rev. E (2010) , 061503; ibid. 061504] is applied to the description of the non-equilibrium processes involved in the spinodal decomposition of suddenly and deeply quenchedsimple liquids. For model liquids with hard-sphere plus attractive (Yukawa or square well) pairpotential, the theory predicts that the spinodal curve, besides being the threshold of the thermo-dynamic stability of homogeneous states, is also the borderline between the regions of ergodic andnon-ergodic homogeneous states. It also predicts that the high-density liquid-glass transition line,whose high-temperature limit corresponds to the well-known hard-sphere glass transition, at lowertemperature intersects the spinodal curve and continues inside the spinodal region as a glass-glasstransition line. Within the region bounded from below by this low-temperature glass-glass tran-sition and from above by the spinodal dynamic arrest line we can recognize two distinct domainswith qualitatively different temperature dependence of various physical properties. We interpretthese two domains as corresponding to full gas-liquid phase separation conditions and to the for-mation of physical gels by arrested spinodal decomposition. The resulting theoretical scenario isconsistent with the corresponding experimental observations in a specific colloidal model system. PACS numbers: 23.23.+x, 56.65.Dy . INTRODUCTION. Quenching a liquid from supercritical temperatures into the liquid/gas coexistence regionleads to the separation of the system into the equilibrium coexisting phases [1, 2]. This phaseseparation process is a paradigmatic illustration of non-equilibrium irreversible processes,whose description must be based on the general fundamental principles of molecular non-equilibrium statistical thermodynamics [3, 4]. When the final density and temperature ofthe quench correspond to a state point ( n, T ) inside the spinodal region, the separationprocess starts with the amplification of spatial density fluctuations, eventually leading tothe two-phase state in thermodynamic equilibrium [5–10]. However, under some conditions,involving mostly colloidal liquids [11–16], this process of phase separation is interruptedwhen the colloidal particles cluster in a percolating network of the denser phase, forming anamorphous sponge-like non-equilibrium bicontinuous structure, typical of physical gels [17].In understanding these experimental observations, computer simulation experiments inwell-defined model systems have also played an essential role. This includes Brownian dy-namics (BD) [18] and molecular dynamics (MD) [19] investigations of the density- andtemperature-dependence of the liquid-gas phase separation kinetics in suddenly-quenchedLennard-Jones and square-well [20] liquids, the BD simulation of gel formation in the col-loidal version of the primitive model of electrolytes (a mixture of equally-sized oppositelycharged colloids) [12], and MD simulations of bigel formation driven by the demixing ofthe two colloidal species of a binary mixture[15]. From the theoretical side, on the otherhand, a fundamental framework is missing that provides an unified description of both, theirreversible evolution of the structure and dynamics of a liquid suddenly quenched to an un-stable homogeneous state, and the possibility that along this irreversible process the systemencounters the conditions for its structural and dynamical arrest.Any theoretical attempt to model this complex phenomenon, however, must be guidedby the well-established experimental facts and simulation results above, according to which,arrested spinodal decomposition is indeed a possible mechanism for gelation in attractivecolloidal suspensions. Although there is still an ongoing discussion regarding specific details,such as the location of the arrest line [11, 14], one can expect that some of these detailsare actually system-dependent, and that they will eventually be fully understood. In themeanwhile, it is important to focus on the simplest experimental model systems, where2his phenomenology has been characterized as carefully and as thoroughly as possible. Thesimplest such model system should involve, for example, a perfectly monodisperse suspensionof spherical particles with attractive interactions that do not require a third species (ascolloid-polymer mixtures do), so that their interactions are state-independent, without theneed to translate the real temperature into any form of effective temperature. The solutionsof the globular protein lysozyme studied by Schurtenberger and collaborators [13, 14] are,to our knowledge, the closest experimental realization of such an ideal model system. Thus,the experimental scenario established by this group may provide an invaluable guide thatmight serve as a reference for the initial theoretical advances in this field. For this reason,we now quote the main features of such experimental scenario.According to Refs. [13, 14], the experimental equilibrium phase diagram of aqueouslysozyme solutions shows at high temperatures a stable gas-crystal equilibrium coexistence,and at lower temperatures, a metastable gas-liquid phase separation. Such equilibriumscenario is typical of systems of particles interacting through hard-sphere–like repulsionsplus very short-ranged attractions. The non-equilibrium experiments performed consist ofthe fast temperature quench at fixed overall protein concentration (or volume fraction φ ),starting with the system in equilibrium at an initial temperature T i above the gas-liquidcoexistence curve T gl ( φ ), cooling the system to a final temperature T lower than T gl ( φ ).The region below the coexistence curve “can be separated into three areas differing in theirkinetic behavior: a region of complete demixing (I), gel formation through an arrestedspinodal decomposition (II), and an homogeneous attractive glass (III)” [13, 14]. Morespecifically, for a shallow quench to a final temperature below the spinodal curve, T < ∼ T s ( φ ),the classical sequence of spinodal decomposition leading to complete phase separation isobserved (region I). However, for quenches below a threshold (or “tie-line”) temperature T ( φ ) ( < T s ( φ )), the spinodal domain structure initially coarsens but then completely arrests,forming gels via an arrested spinodal decomposition (region II). Microscopically, these gelscorrespond to a coexistence of the dilute fluid with a dense percolated glass phase. At lowertemperatures, region II is delimited by the boundary with the glass phase (region III), wherehomogeneous attractive glasses can be reached by quenches at sufficiently low T . Many otherobservations, which can be consulted in the original references [13, 14], enrich the observedexperimental scenario with more detail, but the existence of these three kinetically distinctregions constitute the most relevant feature calling for a fundamental explanation.3ome elements needed to construct the desired unifying theoretical framework may al-ready be contained in the classical theory of spinodal decomposition [5–9]. Unfortunately, itis not clear how to incorporate in these approaches the main non-equilibrium features of theprocess of dynamic arrest, including its history dependence and aging behavior. Similarly,it does not seem obvious how to incorporate the non-stationary evolution of the structureand dynamics of a liquid during the process of spinodal decomposition in existing theories ofglassy behavior [21]. For example, although many of the predictions of conventional modecoupling theory (MCT) [22–25] have found beautiful experimental verification, this theoryis unable to describe irreversible non-stationary processes. The reason for this is that, in itspresent form, MCT is in reality a theory of the dynamics of liquids in their thermodynamicequilibrium states. To overcome this fundamental limitation, in 2000 Latz [26] proposeda formal extension of MCT to situations far away from equilibrium which, however, hasnot yet found a specific quantitative application. In this context, let us also refer to themean-field theory of the aging dynamics of glassy spin systems developed by Cugliandoloand Kurchan [27]. This theory has also made relevant detailed predictions that have beenverified in experiments and simulations. Unfortunately, the models involved lack a geometricstructure and hence cannot describe the spatial evolution of real colloidal glass formers.MCT’s limitation above is shared by the self-consistent generalized Langevin equation(SCGLE) theory of dynamic arrest [28–30], a theory that predicts, in a manner completelyanalogous to MCT, the existence of dynamically arrested states. In contrast to MCT, how-ever, the SCGLE theory has been extended to describe the irreversible non-equilibriumevolution of glass-forming liquids, thus resulting in the “non-equilibrium self-consistent gen-eralized Langevin equation” (NE-SCGLE) theory [31–33]. This extended theory provides aconceptually simple picture of the crossover from ergodic equilibration to non-equilibriumaging [34, 35], offering a perspective of the glass transition in which the “waiting”-time t be-comes a relevant active variable and in which indeed the measured properties of the systemmay depend on the protocol of preparation and on the duration of the actual measurements.Thus, this non-equilibrium theory may provide the general framework in which to build themissing unified description of the processes of arrested spinodal decomposition.The main purpose of the present paper is to take the first steps in exploring this possibility.The first of these steps is to discuss the manner in which the NE-SCGLE theory may beadapted to address the description of arrested spinodal decomposition. The second refers to4btaining and presenting the concrete results of its application to a specific model systemthat permits at least qualitative contact with the observed behavior of an experimental modelsystem such as that described above. The physical interpretation of our theoretical results,however, is not as straightforward as, for example, calculating the equilibrium properties ofa system using century-old and well-established statistical thermodynamic methods (suchas calculating the partition function [2]). Instead, because of the non-equilibrium and non-linear nature of the phenomenon described, the interpretation of the specific predictions ismuch more subtle, involving features that might seem at first sight strongly counterintuitive.Thus, the third and most relevant step, is to properly interpret these results.Regarding the first of these steps, let us mention that the most general version of theNE-SCGLE theory of irreversible processes in colloidal liquids [31–33] consists of two fun-damental time-evolution equations, one for the mean value n ( r , t ), and another for thecovariance σ ( r , r ′ ; t ) ≡ δn ( r , t ) δn ( r ′ , t ) of the fluctuations δn ( r , t ) = n ( r , t ) − n ( r , t ) of thelocal concentration profile n ( r , t ) of a colloidal liquid. Thus, in the following section and inthe appendix we review this general description, to discuss the form adopted by these equa-tions in the absence of external fields, and the possibility that they allow the description of afingerprint of spinodal decomposition, namely, the early stage of the amplification of spatialheterogeneities. The resulting equations are thus restricted to allow the discussion of thesimplest explicit protocol of thermal and mechanical preparation, namely, an instantaneousisochoric quench, to describe the spontaneous relaxation of the system toward its equilibriumstate or toward predictable arrested states. To simplify the reading, Section II is essentiallya summary of the specific NE-SCGLE equations that we shall actually solve, while the ap-pendix contains a description of the main derivations and approximations leading to thesespecific equations.In Sect. III we apply the resulting approximate theory to the discussion of the possi-bility that the process of spinodal decomposition may be interrupted by the emergence ofdynamic arrest conditions. For the sake of concreteness, there we shall have in mind a veryspecific model system, namely, the “hard-sphere plus attractive Yukawa” (HSAY) potential,subjected to an instantaneous quench at t = 0 from an initial temperature T i to a finaltemperature T f = T ( < T i ), with the volume fraction φ held fixed. The main focus of thisapplication is the determination of the dynamic arrest diagram in the region of unstable ho-mogeneous states. According to the predicted scenario, the spinodal curve turns out to be5he borderline between the regions of ergodic and non-ergodic homogeneous states, and thehigh-density liquid-glass transition line, whose high-temperature limit corresponds to thewell-known hard-sphere glass transition, actually intersects the spinodal curve at lower tem-peratures and densities, and continues inside the spinodal region as a glass-glass transitionline.Some elements of this predicted scenario, however, may seem strongly counterintuitive,and not consistent with experimental observations. For example, this scenario predicts thatany quench inside the spinodal region will meet conditions for dynamic arrest, while we knowthat at least for shallow quenches the system will definitely phase-separate. A satisfactoryunderstanding of these features requires a detailed discussion of a number of issues whosesubtlety derives from the lack of familiarity with the solutions of the full NE-SCGLE non-linear equations. Thus, in Sect. IV we discuss the physical interpretation of this theoreticalscenario by analyzing in more detail the predicted long-time non-equilibrium arrested struc-ture factor S a ( k ; φ, T ), particularly the length scale ξ a ( φ, T ) = 2 π/k max associated with theposition k max of its small- k peak. This length represents the size of the spatial heterogeneitiestypical of the early stage of spinodal decomposition described by the Cahn-Hilliard-Cook(CHC) classical linear theory of spinodal decomposition [5, 6], shown to be a particularlimit of the NE-SCGLE theory. This analysis results in the recognition that indeed, shallowquenches lead to full phase separation whereas deeper quenches lead to the formation of gels.Thus, in Sect. V we finally establish a more detailed connection between our theoreticalscenario and the experimental observations in the lysozyme real model system, and in Sect.VI we summarize the main results of this work. II. REVIEW OF THE NE-SCGLE THEORY.
Let us start by considering a Brownian liquid formed by N particles that execute Brown-ian motion (characterized by a short-time self-diffusion coefficient D ) in a volume V whileinteracting between them through a generic pair potential u ( r ) = u HS ( r ) + u A ( r ) that is thesum of a hard-sphere (HS) term u HS ( r ) with HS diameter σ , plus an attractive tail u A ( r ).Let us imagine that we subject this system to a prescribed protocol of thermal and/or me-chanical treatment, described by the chosen temporal variation of the reservoir temperature T ( t ), of the total volume V ( t ), and of the applied external fields ψ ( r , t ). In principle, the6ltimate goal of the NE-SCGLE theory is to predict the response of our system to each pos-sible protocol of thermal and mechanical manipulation. For simplicity let us assume thatthe thermal diffusivity of the system is sufficiently large that any temperature gradient isdissipated almost instantly (compared with the relaxation of chemical potential gradients),so that the temperature field inside the system can be considered uniform, and instantlyequal to the temperature of the reservoir, T ( r , t ) = T ( t ).Let us now imagine that we use these time-dependent fields and constraints to drive thesystem to a prescribed (but arbitrary) initial state characterized by a given mean concen-tration profile n ( r ) and covariance σ ( k ; r ), and we then fix the field, the temperature, andthe volume afterward, ψ ( r , t ) = ψ ( r ), T ( t ) = T , and V ( t ) = V for t >
0. The NE-SCGLEtheory will then allow us to predict the most fundamental and primary information thatcharacterizes the intrinsic properties of the system, namely, its spontaneous behavior duringthe irreversible relaxation towards its new expected equilibrium state. This is described interms of the time-evolution of the non-equilibrium mean concentration profile n ( r , t ) andcovariance σ ( k ; r , t ) for t > instantaneous cooling or heating from an initial temperature T i to a final tem-perature T f , so that T ( t ) = T i for t ≤ T ( t ) = T f for t >
0, under isochoric conditions, V ( t ) = V (so that n ( t ) = n ) at all times, and in the absence of applied external fields, ψ ( r , t ) = 0. As explained in Ref. [34], the NE-SCGLE theory will then describe the spon-taneous response of the system in terms only of the non-stationary but uniform covariance σ ( k ; t ), now written as the non-equilibrium static structure factor S ( k ; t ) ≡ σ ( k ; t ) /n , whosetime-evolution equation for t > ∂S ( k ; t ) ∂t = − k D b ( t ) n E f ( k ) [ S ( k ; t ) − /n E f ( k )] . (2.1)In this equation E f ( k ) ≡ E h ( k ; n, T f ) is the Fourier transform of the thermodynamicfunctional derivative E [ r , r ′ ; n, T ] ≡ [ δβµ [ r ; n, T ] /δn ( r ′ )], evaluated at the uniform den-sity and temperature profiles n ( r ) = n and T ( r ) = T , in which case it can be writtenas E [ r , r ′ ; n, T ] = E h ( | r − r ′ | ; n, T ) = δ ( r − r ′ ) /n − c ( | r − r ′ | ; n, T ) or, in Fourier space, as E h ( k ; n, T ) = 1 /n − c ( k ; n, T ), where c ( r ; n, T ) is the so-called direct correlation function.7he non-stationary, non-equilibrium mobility function b ( t ) entering in Eq. (2.1) above,is the most important kinetic property. Within the NE-SCGLE theory this property isdetermined by (see Refs. [32, 33]) b ( t ) = [1 + Z ∞ dτ ∆ ζ ∗ ( τ ; t )] − , (2.2)with the t -evolving, τ -dependent friction coefficient ∆ ζ ∗ ( τ ; t ) given approximately by∆ ζ ∗ ( τ ; t ) = D π n Z d k k (cid:20) S ( k ; t ) − S ( k ; t ) (cid:21) × F ( k, τ ; t ) F S ( k, τ ; t ) (2.3)in terms of S ( k ; t ) itself and of the collective and self non-equilibrium intermediate scatteringfunctions F ( k, τ ; t ) and F S ( k, τ ; t ), whose memory-function equations are written, in Laplacespace, as F ( k, z ; t ) = S ( k ; t ) z + k D S − ( k ; t )1+ λ ( k ) ∆ ζ ∗ ( z ; t ) , (2.4)and F S ( k, z ; t ) = 1 z + k D λ ( k ) ∆ ζ ∗ ( z ; t ) , (2.5)with F ( k, z ; t ), F S ( k, z ; t ), and ∆ ζ ∗ ( z ; t ), being the corresponding Laplace transforms (LT)and with λ ( k ) being a phenomenological “interpolating function” [30], given by λ ( k ) = 1 / [1 + ( k/k c ) ] , (2.6)where the phenomenological cut-off wave-vector k c depends on the system considered. Sincein this paper we are only interested in semi-quantitative trends, here we adopt the value k c = 1 . π/σ ) = 8 . /σ , determined in a calibration procedure involving simulation dataof the hard-spheres system [41]. III. DYNAMIC ARREST DIAGRAM OF A SIMPLE MODEL SYSTEM.
In this section we apply the previous general theory to the discussion of the possibilitythat the process of spinodal decomposition may be interrupted by the emergence of dynamicarrest conditions. For this, let us consider a generic liquid formed by N particles in a volume8 interacting through a simple pair potential u ( r ) = u HS ( r ) + u A ( r ), which is the sum ofa hard-sphere term u HS ( r ) plus an attractive tail u A ( r ). Although all the results that weshall discuss here are independent of the specific functional form of u A ( r ), for the sake ofconcreteness, in what follows we shall have in mind a very specific model system, namely,the “hard-sphere plus attractive Yukawa” (HSAY) potential, defined as u ( r ) = ∞ , r < σ ; − ǫ exp[ − z ( r/σ − r/σ ) , r > σ, (3.1)where σ is the hard sphere diameter, ǫ the depth of the attractive Yukawa well at contact,and z − its decay length (in units of σ ). For given σ , ǫ , and z , the state space of this systemis spanned by two macroscopic variables, namely, the number density n = N/V and thetemperature T , which we express in dimensionless form, as n ∗ ≡ nσ and T ∗ ≡ k B T /ǫ . Inpractice, however, from now on we shall use σ as the unit of length, and ǫ/k B as the unitof temperature, so that n ∗ = n and T ∗ = T ; most frequently, however, we shall refer to thehard-sphere volume fraction φ ≡ πn/ A. Equilibrium thermodynamic properties.
The most fundamental thermodynamic property of this system, that we need to specifyin order to apply the NE-SCGLE theory described above, is the thermodynamic functionalderivative E [ r , r ′ ; n, T ] ≡ [ δβµ [ r ; n, T ] /δn ( r ′ )]. Evaluated at the uniform density and tem-perature profiles n ( r ) = n and T ( r ) = T , we have that E [ r , r ′ ; n, T ] = E h ( | r − r ′ | ; n, T ) = δ ( r − r ′ ) /n − c ( | r − r ′ | ; n, T ) or, in Fourier space, as E h ( k ; n, T ) = 1 /n − c ( k ; n, T ), with c ( r ; n, T ) being the direct correlation function. In the present paper we shall rely on theapproximate prescription proposed by Sharma and Sharma [42] to determine this thermo-dynamic property, which for the HSAY potential above reads c ( k ; φ, T ) = c HS ( k ; φ ) − βu A ( k ) , (3.2)in which c HS ( k ; φ ) is the FT of the direct correlation function of the hard-sphereliquid, approximated by the Percus-Yevick approximation with Verlet-Weis correction(PYVW) [43, 44] and βu A ( k ) is the FT of the attractive potential (in units of thethermal energy k B T = 1 /β ), i.e., βu A ( k ) ≡ π R ∞ σ drr [sin( kr ) / ( kr )] βu A ( r ). Thus,9 φ T Spinodal lineFreezing lineLiquid-glass dynamic arrest lineCoexistence lineCritical pointTriple point ¤ Bifurcation point ( φ b ,T b ) ★★ FIG. 1: Binodal and spinodal lines of the gas-liquid coexistence of the HSAY model fluid inEq. (3.1) with z = 2, calculated with the Sharma-Sharma approximation for E h ( k ; n, T ). Thecorresponding critical point (black star) is located at ( φ, T ) = (0 . , . S ( eq ) ( k ; φ, T ) =1 /n E h ( k ; φ, T ) reaches 2.85. The corresponding triple point (soft red circle) is located at ( φ, T ) =(0 . , . γ eq ( T, φ ) of the equilibrium “bifurcation equation”, Eq. (3.3). The dark circle denotesthe intersection ( φ b , T b ) = (0 . , .
45) of this arrest line with the spinodal curve. for the HSAY potential, βu σ ( k ) ≡ − (4 π/T ∗ ) R ∞ drr [sin( kr ) / ( kr )] exp[ − z ( r − /r = − (4 π/T ∗ ) [( k cos k + z sin k ) /k ( k + z )].For fixed z , the equilibrium phase diagram of the HSAY model system in the state space( φ, T ) contains the gas, liquid and (crystalline) solid phases. The previous approximationfor E h ( k ; n, T ) allows us to sketch the most prominent features of the fluid phases. Forexample, the gas-liquid transition involves a coexistence region, with its associated binodaland spinodal lines. The former can be determined by integration of β [ ∂p/∂n ] T = E h ( k =0; φ, T ) to get the equation of state, along with Maxwell’s construction. The latter can bedetermined from the condition for thermodynamic instability of uniform fluid states, namely, E h ( k = 0; n, T ) ≤
0. The resulting binodal and spinodal lines are illustrated in Fig. 1 forthe HSAY system with z = 2. On the other hand, the equilibrium static structure factor ofthe fluid state is determined by the equilibrium (Ornstein-Zernike) condition S ( eq ) ( k ; φ, T ) =1 /n E h ( k ; φ, T ). From S ( eq ) ( k ; φ, T ) the freezing line, which is another boundary of stability of10he uniform liquid state, can be readily sketched using the phenomenological Hansen-Verletcondition [45] that the height S ( eq ) max ( φ, T ) of the main peak of the equilibrium static structurefactor S ( eq ) ( k ; φ, T ) reaches the value S ( eq ) max ( φ, T ) = 2 .
85. Besides the spinodal, binodal andfreezing lines, Fig. 1 also exhibits the corresponding gas-liquid critical point, located at( φ, T ) = (0 . , . φ, T ) = (0 . , . S ( eq ) ( k ; φ, T ) = 1 /n E h ( k ; φ, T ). B. Dynamic arrest diagram from the equilibrium SCGLE theory.
The first relevant task in the description of dynamic arrest phenomena in a given sys-tem, such as that in our example, is the determination of the dynamic arrest diagram, i.e.,the region containing the state points ( φ, T ) where the system will be able to reach ther-modynamic equilibrium (ergodic region) and the region where the system, prevented fromcrystallizing, will be trapped in a dynamically arrested state (non-ergodic region). Undernormal circumstances this can be approached solving the so-called bifurcation equationsof mode coupling theory [22], or the corresponding equations of the equilibrium version ofthe SCGLE theory [30]. The latter consists of the following equation for the equilibriumlocalization length squared γ ( eq ) ( φ, T ),1 γ ( eq ) = 16 π n Z ∞ dkk (cid:2) S ( eq ) ( k ; φ, T ) − (cid:3) λ ( k )[ λ ( k ) S ( eq ) ( k ; φ, T ) + k γ ( eq ) ] [ λ ( k ) + k γ ( eq ) ] . (3.3)The only input of this equation is the equilibrium static structure factor S ( eq ) ( k ; φ, T ) =1 /n E h ( k ; φ, T ) at a given state point ( φ, T ). If the solution is γ ( eq ) ( φ, T ) = ∞ , we concludethat the state point ( φ, T ) lies in the ergodic region, whereas if γ ( eq ) ( φ, T ) is finite, the point( φ, T ) lies in the non-ergodic region of state space. In this manner we can determine thedynamic arrest transition line.This procedure has been performed in the past for several model systems, for which theequilibrium static structure factor S ( eq ) ( k ; φ, T ), given in terms of the equilibrium condition S ( eq ) ( k ; φ, T ) = 1 /n E h ( k ; φ, T ), is well-defined in the entire state space ( φ, T ). In the presentcase we have carried out this exercise for our HSAY model, and the result is the liquid-glass dynamic arrest transition line ( φ c , T c ) represented in Fig. 1 by the dark (blue) solidline, with the region to the left and above this curve corresponding to ergodic states. This11ynamic arrest line extends from its high temperature limit, corresponding to the hard-sphere glass transition at ( φ, T ) = (0 . , ∞ ), down to its intersection ( φ b , T b ) = (0 . , . φ, T ) inside the spinodal no equilibrium static structure factor S ( eq ) ( k ; φ, T ) exists that corresponds tospatially uniform states. It is at this point that the power of the non-equilibrium version ofthe SCGLE theory is manifested, since this general theory does provide a proper manner toovercome this limitation of the equilibrium theory thus offering an intriguing and unexpectedqualitative scenario of the interplay between dynamic arrest and gas-liquid phase separation.As we shall see below, the alternative method is NOT based on the use of the equilibriumbifurcation equation, Eq. (3.3), but on its non-equilibrium generalization, Eq. (3.7) below. C. Detecting dynamic arrest transitions using the NE -SCGLE theory. Let us now use the NE-SCGLE theory to detect dynamic arrest transitions. We startby noticing that from the very structure of Eq. (2.1) one can infer the existence of twofundamentally different kinds of stationary solutions, representing two fundamentally dif-ferent kinds of stationary states of matter. The first corresponds to the condition in whichstationarity is attained because the factor [ S ( k ; t ) − /n E f ( k )] on the right side of Eq.(2.1) vanishes, i.e., because S ( k ; t ) is able to reach its thermodynamic equilibrium value S ( eq ) ( k ; n, T f ) = 1 /n E h ( k ; n, T f ), while the mobility b ( t ) attains a finite positive long-timelimit b eqf . These stationary solutions correspond, of course, to the ordinary thermodynamicequilibrium states, in which the difference [ S ( k ; t ) − /n E f ( k )] decays (according to Ec. (2.1),and in a gross approximation) exponentially fast, [ S ( k ; t ) − /n E f ( k )] ∝ exp[ − t/t eq ], withthe equilibration time t eq estimated as t eq ( n, T f ) ∝ (cid:2) k max D n E f ( k max ) b eqf (cid:3) − .The second class of stationary solutions of Eq. (2.1) emerges from the possibility that thelong-time asymptotic limit of the kinetic factor b ( t ) vanishes, so that dS ( k ; t ) /dt vanishes atlong times without requiring the equilibrium condition [ S ( k ; t ) − /n E h ( k ; n, T f )] = 0 to befulfilled. This mathematical possibility has profound and general implications. For exam-ple, under these conditions S ( k ; t ) will now approach a distinct non-equilibrium stationarylimit, denoted by S a ( k ), which is definitely different from the expected equilibrium value12 ( eq ) ( k ; n, T f ) = 1 /n E h ( k ; n, T f ). Furthermore, the difference [ S ( k ; t ) − S a ( k )] is predictedto decay to zero not exponentially, but in a much slower fashion, namely, as t − x , with anexponent x near unity [33, 34]. Contrary to ordinary equilibrium states, whose propertiesare determined solely by the fundamental condition of maximum entropy, the properties ofthese stationary but intrinsically non-equilibrium states, such as S a ( k ), may depend on thepreparation protocol (in our example of the instantaneous isochoric quench, on T i and T f ).For a given instantaneous isochoric quench from an initial temperature T i to a lower finaltemperature T f , one can also predict if the system will equilibrate or will be trapped ina non-equilibrium arrested state by analyzing the stationary solutions of Eq. (2.1). Thisanalysis simplifies greatly with the change of variables from the actual evolution time t tothe so-called “material” time u ( t ) defined by [34] u ( t ) ≡ Z t b ( t ′ ) dt ′ . (3.4)This allows us to write, for example, the actual solution S ( k ; t ) of Eq. (2.1) as S ( k ; t ) = S ∗ ( k ; u ( t )), with the function S ∗ ( k ; u ) being the solution of ∂S ∗ ( k ; u ) ∂u = − k D n E f ( k ) [ S ∗ ( k ; u ) − /n E f ( k )] (3.5)i.e., S ( k ; t ) = S ∗ ( k ; u ( t )) ≡ [ n E f ( k )] − + (cid:8) S i ( k ) − [ n E f ( k )] − (cid:9) e − k D n E f ( k ) u ( t ) , (3.6)where S i ( k ) ≡ S ( k ; t = 0) = S ∗ ( k ; u = 0) represents the (arbitrary) initial condition.This expression interpolates S ∗ ( k ; u ) between its initial value S i ( k ) and the value of thethermodynamic property [ n E f ( k )] − . Thus, if nothing impedes the system from reachingequilibrium at the state point ( n, T f ), the non-stationary static structure factor S ( k ; t ) willeventually attain its equilibrium value S ( eq ) ( k ; n, T f ) = [ n E f ( k )] − . This, however, is onlyone of the two possibilities described above. In order to determine which of them willgovern the course of the spontaneous response of our system, we actually input S ∗ ( k ; u ) inthe “bifurcation” equation for the square localization length γ ( t ) at evolution time t (i.e.,the long- τ asymptotic value of the mean squared displacement). Denoting γ ( t ) = γ ∗ ( u ( t ))simply as γ ∗ ( u ), such an equation reads1 γ ∗ ( u ) = 16 π n Z ∞ dkk [ S ∗ ( k ; u ) − λ ( k ; u )[ λ ( k ; u ) S ∗ ( k ; u ) + k γ ∗ ( u )] [ λ ( k ; u ) + k γ ∗ ( u )] . (3.7)13s explained in Ref. [34], if we find that γ ∗ ( u ) = ∞ for 0 ≤ u ≤ ∞ , we conclude that thesystem will be able to equilibrate after this quench, and hence, that the point ( n, T f ) lies inthe ergodic region. If, instead, a finite value u a of the parameter u exists, such that γ ∗ ( u )remains infinite only within a finite interval 0 ≤ u < u a , the system will no longer equilibrate,but will become kinetically arrested, with the non-equilibrium asymptotic structure factor S a ( k ) given by S a ( k ) ≡ S ∗ ( k ; u a ) = [ n E f ( k )] − + (cid:8) S i ( k ) − [ n E f ( k )] − (cid:9) e − k D n E f ( k ) u a . (3.8)The parameter u a and the corresponding long-time asymptotic square localization length γ ∗ a ≡ γ ∗ ( u a ) are dynamic order parameters in the sense that if u a = γ ∗ a = ∞ , the system willreach equilibrium, but if u a and γ ∗ a are finite, the system will be dynamically arrested. Bothof them, as well as the non-equilibrium static structure factor S a ( k ), will depend on theprotocol of the quench (in our case, on φ , T i , and the final temperature T f , which from nowon will be denoted simply as T ). Thus, if one determines the functions u a = u a ( φ, T i , T ), γ ∗ a = γ ∗ a ( φ, T i , T ), and S a ( k ) = S a ( k ; φ, T i , T ), one can in principle draw the dynamic arrestdiagram.We have implemented this protocol and here we present the results of its application tothe concrete illustrative case involving the HSAY model. Let us first consider isochores with φ > φ b , which lie to the right of the intersection indicated by the dark circle in Fig. 1. Arepresentative isochoric quench of this class is schematically indicated in the inset of Fig.2(a) by the downward arrow along the isochore φ = 0 .
4. In the main panel of this figurewe plot the resulting u − a ( T ) as a function of T (solid line), together with the correspondingvalue of γ ∗− a ( T ) (dashed line).The first feature to notice is that indeed, for quenches to final temperature T above thecritical temperature T c = 0 . u a ( T ) is always infinite (so that u − a ( T ) = 0, and hence, donot appear in the logarithmic window of the figure), whereas for T ≤ T c we detect a finitevalue of u a ( T ). This value, however, diverges as T approaches T c from below, and decreasesmonotonically as the final temperature T decreases. According to the results in Fig. 2(a),for T ≤ T c also γ ∗ a ( T ) attains finite values, thus reaffirming the existence of an arrested phasein this temperature regime. In fact, these results for γ ∗ a ( T ) reveal a second relevant feature,namely, that this dynamic order parameter changes discontinuously at T = T c , from a value γ ∗− a ( T − c ) = 41 to a value γ ∗− a ( T + c ) = 0. We notice that this jump in the value of γ ∗ a ( T )14 φ T T u γ T i T f (a) T c ( φ=0.4) aa * -1 γ a - u * a; -1 - T -2 T i =0.8T i =1.0T i =1.3 (b) T c ( φ=0.4) u γ a-1 (T) a -1 * (T) γ a * - u a - FIG. 2: (a) Dependence of u a ( T ) (solid line) and of γ ∗ a ( T ) (dashed line) on the final temperature T ofan instantaneous quench with initial temperature T i = 1 . φ = 0 .
4. Theinset reproduces essentially Fig. 1, with the arrow indicating this illustrative isochoric quench. (b) u a ( T ) for the previous quench and for other two quenches that start at different initial temperatures,namely, T i = 1 . T i = 0 .
8. Notice that the three cases yield the same result for thetransition temperature, T c ( φ = 0 .
4) = 0 . corresponds to the jump in the solution γ eq ( T ) of the equilibrium bifurcation equation, Eq.(3.3), and it occurs at the same critical temperature T c = 0 . φ -dependence of thedynamic arrest transition temperature T c ( φ ), thus finding that indeed, outside the spinodalcurve (i.e., for volume fractions φ larger than the volume fraction of the intersection point( φ b , T b ) = (0 . , . u a ( T ) on the finaltemperature T do depend on the initial temperature T i of the quench. However, the resultingtransition temperature, T c ( φ = 0 .
4) = 0 . T i . This is illustrated in Fig.2(b), which compares the function u a ( T ) obtained with T i = 1 . u a ( T )obtained with T i = 1 . T i = 0 .
8. This independence of T c ( φ ) on the initial temperature T i is also a necessary condition for the present method to be a reliable approach to the15etermination of the liquid-glass dynamic arrest transition. In Fig. 2(b) we have alsoincluded the results for γ ∗− a ( T ) corresponding to the quench processes with different initialtemperatures, to show that γ ∗ a ( T ) is much less sensitive to the initial temperature of thequench. Let us also notice that for φ > φ b , the results for γ ∗ a ( T ) and u a ( T ) as a functionof the final temperature T , do not exhibit any discontinuous behavior when crossing thespinodal line, i.e., their dependence on T does not detect the spinodal. Thus, for φ > φ b thespinodal condition E h ( k = 0; φ, T ) = 0 does not seem to have any significant effect on thedynamics of the system. As discussed below, this rather astonishing result is a consequenceof the structure of the expression for S a ( k ) in Eq. (3.8), and is in dramatic contrast withwhat happens in the complementary regime, i.e., for volume fractions φ ≤ φ b , as we nowexplain. D. Dynamic arrest inside the spinodal region.
The main advantage of the NE-SCGLE methodology just explained, for determining thedynamic arrest transition line, is that its use can be extended to the interior of the spinodalregion. This subsection continues the description of the use of this methodology, but nowto unveil the structure of the dynamic arrest diagram of our model in the region inside thespinodal line. Thus, let us repeat the calculations illustrated in Fig. 2, but now for isochorescorresponding to volume fractions φ smaller than the volume fraction φ b = 0 .
33 of theintersection between the liquid-glass dynamic arrest transition and the spinodal curve. Thecorresponding illustrative example is contained in Fig. 3, which plots the functions u a ( T )and γ ∗ a ( T ) as a function of T for the HSAY model system subjected to an instantaneousquench from an initial temperature T i = 1 . T f < . φ = 0 . φ = 0 .
4, in Fig. 2, we notice that in the presentcase the functions u a ( T ) and γ ∗ a ( T ) remain infinite for all temperatures above a criticaltemperature T s = 0 .
76, and that this singular temperature coincides precisely with thetemperature of the spinodal curve for that isochore. The results in Fig. 3 also indicate thatboth parameters, u a ( T ) and γ ∗ a ( T ), are finite for T < T s , and both diverge as T approaches T s from below. This reveals a rather dramatic and unexpected conclusion, namely, thatthe spinodal curve, besides being the threshold of thermodynamic instability, now turns out16 φ T T -2 γ u T i T f T s ( φ=0.2) T c ( φ=0.2) * -1aa-1 γ - * ; u a - FIG. 3: Dependence of u a ( T ) (solid line) and of γ ∗ a ( T ) (dashed line) on the final temperature T ofan instantaneous quench with initial temperature T i = 1 . φ = 0 .
2. Theinset reproduces essentially Fig. 1, with the arrow indicating this illustrative isochoric quench.Notice that γ a and u a are discontinuous at a temperature T c , which implies the existence of atype-B glass-glass transition, whereas as T approaches to T s γ a and u a diverge, thus implying atype-A dynamic arrest transition at the spinodal. to be also the threshold of non-ergodicity. Furthermore, the fact that γ ∗ a ( T ) diverges as T approaches T s from below determines that this transition from ergodic to non-ergodic statesoccurs in a continuous fashion, i.e., that it is a “type A” dynamic arrest transition in MCTlanguage.Examining again the same results in Fig. 3, but now at temperatures well below thespinodal curve, we see that the parameters u a ( T ) and γ ∗ a ( T ) exhibit a discontinuity at a lowertemperature T c ( φ ). This discontinuity reveals the existence of still a second dynamic arresttransition, now corresponding to a glass-glass “type B” transition, in which the dynamicorder parameter γ ∗ a ( T ) changes discontinuously by about one order of magnitude, from avalue γ ∗ a ( T − c ) = 36 to another finite value γ ∗ a ( T + c ) = 4. Notice that the value γ ∗ a ( T − ) issimilar to the corresponding value at the isochore φ = 0 .
4. Repeating these calculations atother isochores, we determine both transition temperatures, T c ( φ ) and T s ( φ ), as a functionof volume fraction. The corresponding dynamic arrest transition lines are presented in Fig.4. This figure reveals a rather unexpected global scenario, in which (a) this low- T (typeB) glass-glass transition (dark blue dashed curve) turns out to be just the continuationto the interior of the spinodal region, of the liquid-glass transition previously determined17 φ T Spinodal lineGlass-glass transitionLiquid-glass transitionDynamic arrest line ( φ b ,T b ) FIG. 4: Dynamic arrest diagram exhibiting the low- T (type B) glass-glass transition (dark bluedashed curve), which continues the liquid-glass transition (dark blue solid curve) to the interior ofthe spinodal region. The (red) dotted line that superimposes on the spinodal curve is a continuous(type A) ergodic–nonergodic transition. These two transitions merge at the bifurcation point( φ b , T b ). The soft blue dashed curve to the right of the bifurcation point ( φ b , T b ) is the dynamicallyirrelevant part of the spinodal curve (see Subsection IV E). outside the spinodal (dark blue solid curve), (b) the spinodal line itself is a continuous (typeA) ergodic–nonergodic transition, (c) these two transitions merge at the bifurcation point( φ b , T b ) to continue outside the spinodal as the liquid-glass transition line that terminates atthe hard-sphere glass transition point ( φ, T ) = (0 . , ∞ ), and (d) for φ > φ b the spinodalline does not have any apparent dynamic significance.At first glance some of these predictions might appear counterintuitive. Let us remind,however, that they are based only on the analysis of the dependence of the dynamic or-der parameter γ ∗ a ( φ, T ) on the final temperature T of the instantaneous isochoric quenchesanalyzed. This order parameter is a functional of the long-time asymptotic value S a ( k ) ofthe non-stationary static structure factor S ( k ; t ). Thus, the analysis of the dependence of S a ( k ; φ, T ) on φ and T should provide additional pieces of structural information, while thewaiting-time dependence of S ( k ; t ) (and of all the other structural and dynamic properties)will provide valuable pieces of kinetic information. When these pieces are assembled togetherin a structural and kinetic jigsaw puzzle, one should expect that a sound and convincingscenario will emerge. For this, however, we must first identify and describe the most rele-vant of these pieces. The dynamic arrest diagram in Fig. 4 constitutes the initial and most18asic piece for this discussion, but in the rest of the paper we shall identify others, equallyfundamental. IV. PHYSICAL INTERPRETATION OF THE DYNAMIC ARREST DIAGRAM.
The scenario described by the dynamic arrest diagram in Fig. 4 is highly provoking, butalso at first sight difficult to reconcile with experience. First, we all know that simple fluids,when quenched to the inside of the spinodal region, are expected to fully phase-separateinto the gas phase and its corresponding condensed (liquid or crystal) phase. On the otherhand, as explained in the introduction, there are numerous experimental and numericalevidences that for sufficiently deep quenches, the process of spinodal decomposition maybe interrupted by its interference with the dynamic arrest of the condensed phase. Thus,the dynamic arrest diagram in Fig. 4 may contain the fundamental explanation of thisexperimental double scenario of full vs. arrested spinodal decomposition.To test this expectation, however, we must first understand what exactly is the detailedphysical meaning of these NE-SCGLE predictions and what are their most relevant andverifiable physical consequences. In this section we discuss several pieces of informationthat derive from the solution of the NE-SCGLE equations, and which will dissipate someof the most subtle puzzles of the interpretation of this dynamic arrest diagram. With thesedetails taken into account, we shall see that the resulting theoretical scenario is actuallyquite consistent with the main features of the experimental scenario of arrested spinodaldecomposition in lysozyme protein solutions, thus establishing a concrete link of our theorywith physical reality. A. S ( k ; t ) , S ( eq ) ( k ; φ, T ) , and S a ( k ) below the spinodal line. In this subsection, for example, we discuss the mechanism that allows the existenceof a well-behaved non-equilibrium S ( k ; t ) for quenches to state points ( φ, T ) inside thespinodal region, where the uniform fluid state has become thermodynamically unstable andthe equilibrium structure factor S ( eq ) ( k ; φ, T ) does not exist. According to Eqs. (3.6) and(3.8), the calculation of S ( k ; t ) and S a ( k ) requires the value of the thermodynamic property E h ( k ; φ, T ) for k ≥
0. As we know, the condition E h ( k = 0; φ, T ) > T abovethe spinodal temperature T s ( φ ). Thus, the spinodal condition E h ( k = 0; φ, T s ) = 0 definesthe threshold of thermodynamic instability of uniform states, so that for all states withtemperatures below T s ( φ ) we must have that E h ( k ; φ, T ) < k = 0, but also atleast within a finite interval 0 ≤ k ≤ k .This means that the equilibrium static structure factor S ( eq ) ( k ; φ, T ) = 1 /n E h ( k ; φ, T )will attain unphysical negative values in this interval and will exhibit a singularity at k = k . This unphysical behavior, which is a manifestation of the non-existence of spatiallyuniform equilibrium states inside the spinodal region, is illustrated in Fig. 5, which plots[ n E h ( k ; φ, T )] − (dashed curves) at the state points ( φ, T ) = (0 . , .
2) and (0 . , .
2) belowthe spinodal curve. One might then expect that the non-equilibrium static structure factor S ( k ; t ) and its asymptotic value S a ( k ), which are formally an exponential interpolationbetween S i ( k ) and [ n E f ( k ; φ, T )] − (see Eqs. (3.6) and (3.8)), will inherit these singularfeatures of the function [ n E f ( k ; φ, T )] − .In reality, this is not the case, since the appearance of [ n E f ( k )] in the exponent of theinterpolating function, allows us to show, for example, that near k = k we have that S ( k ; t )and S a ( k ) can be approximated, respectively, by S ( k ; t ) = S ∗ ( k ; u ( t )) ≈ S i ( k ) + 2 k D u ( t )and S a ( k ) ≈ S i ( k ) + 2 k D u a , which are always positive definite. This is illustrated in Fig. 5with the results for S a ( k ; φ, T ) (solid lines) for the quenches along the isochores φ = 0.2 and0.4, from the same initial temperature T i = 1 .
0, at which the initial static structure factoris given by S i ( k ) = S ( eq ) ( k ; φ, T i ) = 1 /n E h ( k ; φ, T i ) (dotted lines), to a final temperature T = 0 .
2. As we can see in Figs. 5(a) and (b), in both cases the non-equilibrium stationarystructure factor S a ( k ) does not exhibit any singular or non-physical feature at any wave-vector. Of course, for the same reason, the non-equilibrium static structure factor S ( k ; t )will evolve smoothly from S i ( k ) at t = 0, to S a ( k ) as t → ∞ , without exhibiting any hintof the singular behavior characteristic of the thermodynamic function [ n E f ( k )] − below thespinodal curve. This is illustrated in Fig. 5 by the light-gray solid lines, corresponding to S ( k ; t ) evaluated at representative finite values of the waiting time t . Notice in particularthat for the quench illustrated in Fig. 5(a) (corresponding to φ = 0 . k peak of S ( k ; t ) develops at k max ( t ) ≈
1. This is to be contrasted with the quenchalong the isochore φ = 0 . k peak is imperceptible.20 k S t r u c t u r e F a c t o r S eq (k;T f =0.2)S eq (k;T f =1.0)S(k;t)S a (k) a) φ=0.2 k S t r u c t u r e F a c t o r S eq (k;T f =0.2)S eq (k;T f =1.0)S(k;t)S a (k) b) φ=0.4 FIG. 5: Initial equilibrium static structure factor S ( eq ) ( k ; φ, T i ) = 1 /n E h ( k ; φ, T i ) (dotted red curve)of the isochoric quench from the ergodic state point ( φ, T i ) to the state point ( φ, T ) inside thespinodal region (panel (a) refers to the isochore φ = 0 . φ = 0 . n E h ( k ; φ, T )] − at the latter point, whose negativeand singular features indicate the nonexistence of uniform equilibrium states with static structurefactor S ( eq ) ( k ; φ, T ) = 1 /n E h ( k ; φ, T ). The dark (green) solid curves correspond to the (physicallyacceptable) asymptotic long-time static structure factor S a ( k ) ≡ lim t →∞ S ( k ; t ) = S ∗ ( k ; u a ) of thisquench. The soft gray curves represent S ( k ; t ) evaluated at an (arbitrary) sequence of values ofthe waiting time t . B. Relationship with the theory of Cahn, Hilliard, and Cook.
One of the main general features of the non-equilibrium evolution of S ( k ; t ) after anisochoric quench inside the spinodal region is precisely the development of this low- k peakof S ( k ; t ) located at the time-dependent wave-vector k max ( t ). This peak is associated withthe growing length scale ξ ( t ) ≈ π/k max ( t ) of the clusters formed in the early stage ofthe process of spinodal decomposition. According to the classical theory developed byCahn, Hilliard, and Cook (CHC) [5, 6], and reviewed in detail by Furukawa [7], this small- k peak evolves with waiting time t , moving to progressively smaller wave-vectors, so thatlim t →∞ ξ ( t ) = ∞ . The correct interpretation, however, is that this scenario only describesthe very early stages of the phase separation process, up to a point in which the size ξ ( t )reaches mesoscopic dimensions so that additional effects (such as surface tension, convection,21tc.), not contained in the theory, drive the system to full phase separation.As it happens, the NE-SCGLE theory contains this classical CHC theory as its linearand small- k limit. To see this more explicitly, let us notice that our fundamental timeevolution equation for S ( k ; t ) in Eq.(2.1) above, becomes Cook’s equation (Eq. (3.4) in Ref.[7]) in the linear regime and the long-wavelength limit. The linear regime is achieved inour theory by neglecting the non-linear dependence of the time-dependent mobility b ( t ) on S ( k ; t ) itself, i.e., by approximating b ( t ) = 1, whereas in the long-wavelength limit (small k ’s) we can expand the thermodynamic function n E h ( k ; φ, T ) in powers of k up to quadraticorder, n E h ( k ; φ, T ) ≈ ( a/k B T )+( K/k B T ) k , with a = a ( φ, T ) and K = K ( φ, T ) being wave-vector independent parameters. These two approximations, with M ≡ D /k B T , convertEq. (2.1) into ∂S ( k ; t ) ∂t = 2 k M (cid:2) k B T − ( a + Kk ) S ( k ; t ) (cid:3) , (4.1)which is Eq. (3.4) of Ref. [7]. Notice also that Cook’s approximate solution of this equationfor shallow quenches (Eq. (3.9) of Ref. [7]) is identical to the approximation S ( k ; t ) = S ∗ ( k ; u ( t )) ≈ S i ( k ) + 2 k D u ( t ) quoted above, provided that b ( t ) = 1.The fact that the NE-SCGLE theory contains the CHC theory as a particular limit alsodetermines the limitation of the NE-SCGLE theory to describe the full phase separationprocess. Thus, for shallow quenches, where we know that full phase separation occurs, weunderstand that the NE-SCGLE theory only describes the early stage of spinodal decom-position. Due to its non-linear nature, however, the NE-SCGLE theory complements thisclassical scenario (in which lim t →∞ ξ ( t ) = ∞ ) with the prediction that lim t →∞ ξ ( t ) = ξ a ,i.e., with the possibility that the process of spinodal decomposition is arrested when ξ ( t )reaches the length scale ξ a associated with the low- k peak of the non-equilibrium asymp-totic structure factor S a ( k ). Of course, if the predicted value of ξ a is sufficiently large forsurface tension and convection to take over, the system will phase-separate completely. Thepossibility exists, however, that the predicted value of ξ a is not that large, or to be even ofthe scale of a few particle diameters. Under such circumstances the dynamic arrest of thephase separation process will occur right at its early stage, freezing the current structure andpreventing the system from phase separating completely. These conditions are illustratedby the quench to T = 0 . φ = 0 .
2, whose S a ( k ) is represented by the solidline of Fig. 5(a), for which k max ≈ /σ , and hence, ξ a ≈ σ . Thus, we need a criterion tolocate the crossover temperature T ( φ ) below which the predicted dynamic arrest will in fact22ccur, and above which full phase separation will be driven by effects not contained in ourequations. In the rest of this section we shall search for such criterion in the temperaturedependence of the asymptotic structure represented by S a ( k ). C. Non-equilibrium small- k peak of S a ( k ; φ, T ) and diverging length scale ξ a ( φ, T ) . Let us thus analyze in more detail the dependence of S a ( k ; φ, T ) on the depth of thequench, i.e., on the final temperature T , as well as on the volume fraction φ . With thisintention, in Fig. 6(a) we present a set of results for S a ( k ; φ, T ) that illustrate the T -dependence of this asymptotic structural property for quenches along the isochore φ = 0.2.All of these quenches start with the system equilibrated at the same initial temperature T i = 1 .
0, whose corresponding equilibrium static structure factor S i ( k ) = S ( eq ) ( k ; φ =0 . , T = 1 .
0) is represented by the dotted curve. Each quench ends at a different finaltemperature T , and the solid lines represent the resulting S a ( k ; φ = 0 . , T ). Thus, thedotted line and the solid line corresponding to S a ( k ; φ = 0 . , T = 0 .
2) are the same as thedotted and solid lines of Fig. 5(a), but now in the vertical axis we use logarithmic scaleto visualize the strong increase in S a ( k max ; φ, T ) as the final temperature T approaches thespinodal temperature T s ( φ = 0 .
2) = 0 .
76 from below. This is clearly the most visible trendexhibited by this set of curves, together with the fact that the position k max ( φ, T ) of thisgrowing non-equilibrium peak of S a ( k ; φ, T ) moves monotonically to the left, so that thecorresponding length scale ξ a ( φ, T ) = 2 π/k max ( φ, T ) increases with increasing T .In fact, to further visualize these trends, in Fig. 6(b) we plot as a solid line the height S a ( k max ; φ, T ) of the low- k peak of S a ( k ; φ, T ) as a function of the final temperature T along the isochore φ = 0 .
2. This curve clearly illustrates that the behavior of S a ( k max ; φ, T )is distinctly different in the regime where the final temperature falls below the glass-glasstransition temperature T c ( φ = 0 .
2) = 0 .
22, and in the regime where T falls in the interval T c ≤ T ≤ T s . In the former, for the lowest temperatures S a ( k ; φ, T ) starts almost con-stant, with a value similar to S i ( k = 0) (recall that S a ( k max ) ≈ S i ( k max ) + 2 k max D u a ≈ S i ( k max ) ≈ S i (0) ≈ T approaches T c from below it increases sharply to a value S a ( k ; φ = 0 . , T − c ) = 55. In the interval T c ≤ T ≤ T s , on the other hand, S a ( k max ) startsfrom the value S a ( k ; φ = 0 . , T + c ) = 55 (i.e., S a ( k ; φ, T ) is continuous at T = T c ), and in-creases monotonically up to the spinodal temperature T s , at which it diverges as ( T s − T ) − µ ,23 k -1 S a ( k ) k -1 S a ( k ) S i (k)T=0.10T=0.20T=0.23T=0.40T=0.60T=0.75 φ=0.2 a) T -2 -1 S a ( k m a x ) k main k low k-peak φ=0.20 T c T s b) FIG. 6: (a) Non-equilibrium stationary structure factor S a ( k ; φ, T ) as a function of wave-vector k for a sequence of instantaneous isochoric cooling processes along the isochore φ = 0 . T i = 1 . T -dependence of the low- k peak ( k ≈ /σ ; also zoomed in the inset) and of the main peak ( k ≈ π/σ )of S a ( k ; φ, T ). (b) Height of the low- k peak (dark solid line) and of the main peak (dark dashedline) of S a ( k ; φ, T ), as a function of the final temperature T along the isochore φ = 0 .
2. The softgray solid line represents the diverging behavior (when T → T s from above) of the k = 0 value ofthe equilibrium static structure factor, S ( eq ) ( k = 0; φ, T ). with µ ≈ .
75. This divergence of S a ( k max ; φ, T ) is the non-equilibrium counterpart of thedivergence of S ( eq ) ( k = 0; φ, T ) = [ ∂ ( p/k B T ) /∂n ] − T when T approaches T s from above. Forcompleteness, although it has no direct relevance regarding the possibility of dynamic arrest,the divergence of S ( eq ) ( k = 0; φ, T ) is also shown in Fig. 6(b) (soft gray solid curve to theright of T s ).Let us also notice in Fig. 6(a) that in the regime T c ≤ T ≤ T s , while S a ( k max ) increasesand diverges as T → T − s , its position k max ( φ, T ) moves to the left and vanishes at T = T s , sothat the corresponding length scale ξ a ( φ, T ) ≡ π/k max ( φ, T ) also diverges as T approaches T s from below. This divergence is the non-equilibrium counterpart of the divergence of theequilibrium correlation length ξ eq ( φ, T ) as T → T + s , predicted in 1914 by Ornstein andZernike [36] to explain the phenomenon of critical opalescence along the critical isochore. InFig. 7 we present the temperature dependence of ξ a ( φ, T ) in the interval T c ≤ T ≤ T s alongthree illustrative isochores. We notice that in this interval ξ a ( φ, T ) exhibits two sub-regimes,24 .1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 T -2 -1 ξ − φ=0.05 a (a) T -2 -1 ξ − φ=0.25 a (b) T -2 -1 ξ − φ=0.30 a (c) FIG. 7: Characteristic length ξ a ( φ, T ) = 2 π/k max ( φ, T ) associated with the low- k peak of S a ( k max ; φ, T ), as a function of the final temperature T of the quench, for the isochores (a) φ = 0 . φ = 0 .
25, and (c) φ = 0 .
30. The dots represent ξ a ( φ, T ) determined directly from S a ( k max ; φ, T )at a discrete set of temperatures, while the composed (solid-dashed) curve is the fit of these numer-ical data obtained with the parametrization expression in Eq. (4.2). The vertical lines representthe temperatures T c ( φ ) < T ( φ ) < T s ( φ ), with T ( φ ) indicating the crossover of the T -dependenceof ξ a ( φ, T ) from exponential (solid magenta portion) to diverging power-law (dashed blue portion). characterized by a different dependence on the final temperature. In the first sub-regime ξ a ( φ, T ) increases exponentially with T , whereas in the second, which we shall refer to asthe “Ornstein-Zernike domain”, ξ a ( φ, T ) increases with T as an inverse power of the reducedtemperature T ∗ ≡ ( T s − T ) /T s , diverging at the spinodal line with an exponent close to theclassical (Ornstein-Zernike) exponent that describes the critical divergence of the equilibriumcorrelation length [36].In fact, we have noticed that the results for ξ a ( φ, T ) can be conveniently parametrizedwith an empiric expression that matches an exponential function below T with a power lawabove T , in a manner that both, the function and its derivative, are continuous at T . Suchan expression reads ξ a ( T ) = ξ a ( T + c ) e ν (cid:16) T − TcTs − T (cid:17) , for T c ≤ T ≤ T = ξ a ( T + c ) e ν (cid:16) T − TcTs − T (cid:17) (cid:20) T s − TT s − T (cid:21) − ν , for T ≤ T ≤ T s . (4.2)Since we already have the theoretical prediction for the spinodal and the glass-glass transitionlines T s ( φ ) and T c ( φ ), and we can determine the value ξ a ( T + c ) of ξ a ( T ) immediately above T c ,this functional form allows us to determine the exponent ν and the crossover temperature T ( φ ) by the best fit of the numerical data of ξ a ( T ). In fact, we actually fixed the exponent25 to its classical OZ value ν = 0 .
5, so that the only parameter to be determined by this fit isthe crossover temperature T ( φ ). As observed in Fig. 7, the quality of this parametrizationof the numerical results for ξ a ( φ, T ) along the three illustrative isochores is excellent. D. The threshold temperature T ( φ ) and the gelation line. The procedure just described, performed at a sequence of isochores, defines a curve T ( φ )in the dynamic arrest diagram, which divides the interval T c ( φ ) ≤ T ≤ T s ( φ ) in two well-defined sub-regimes: the subinterval T c ( φ ) ≤ T ≤ T ( φ ), where ξ a ( T ) increases exponentiallywith T , and the subinterval T ( φ ) ≤ T ≤ T s ( φ ), where ξ a ( T ) enters its diverging Ornstein-Zernike regime. Although we cannot jump immediately to the conclusion that these twosubintervals are to be identified, respectively, with the gel and the full spinodal decomposi-tion regions of state space, these findings provide important elements to consider in buildinga sound physical interpretation of the dynamic arrest diagram in Fig. 4. The expectation,however, is that indeed this crossover curve will eventually be identified with the gelationline.Further elements supporting the notion of the existence of this crossover are also revealedby the behavior of other properties. In fact, going back to Fig. 3, one can immediatelyappreciate that at least along the isochore φ = 0 .
2, also the asymptotic square localizationlength γ ∗ a ( T ) seems to exhibit a similar behavior as ξ a ( T ). To test this idea, let us proposea parametrization of its T -dependence similar to that employed for ξ a ( T ), namely, γ ∗ a ( T ) = γ ∗ a ( T c ) e α (cid:16) T − TcTs − T (cid:17) , for T c ≤ T ≤ T = γ ∗ a ( T c ) e α (cid:16) T − TcTs − T (cid:17) (cid:20) T s − TT s − T (cid:21) − α , for T ≤ T ≤ T s . (4.3)We found that a reasonable representation of the numerical results for γ ∗ a ( T ) in Fig. 3 couldbe obtained with the exponent α fixed as α = 1 .
5, and adjusting T to find the best overallfit, thus leading to the value T ( φ = 0 .
2) = 0 .
52. We tested the same procedure for otherisochores, also fixing the exponent α as α = 1 .
5, and determined T ( φ ) from the best overallfit. Fig. 8 illustrates the accuracy of this fit at other three illustrative isochores. Thisprocedure provides a second determination of the location of the crossover curve T ( φ ).Let us finally mention that not only the spatially-meaningful length scales ξ a ( φ, T ) and p γ ∗ a ( φ, T ) exhibit this dual temperature dependence. Looking back again at Fig. 3, we26 T -3 -2 -1 γ − φ=0.05 a (a) 0 0.2 0.4 0.6 0.8 1 T -3 -2 -1 γ − φ=0.25 a (b) 0 0.2 0.4 0.6 0.8 1 T -3 -2 -1 γ − φ=0.30 a (c) FIG. 8: Dynamic order parameter γ ∗ a ( φ, T ) as a function of the final temperature T , and itsparametrization according to Eq. (4.3) at the isochores (a) φ = 0 .
05, (b) φ = 0 .
25, and (c) φ = 0 .
30. As in Fig. 7, the vertical lines represent the temperatures T c ( φ ) < T ( φ ) < T s ( φ ), with T ( φ ) indicating the crossover of the T -dependence of γ a ( φ, T ) from exponential (solid magentaportion) to diverging power-law (dashed blue portion) see that the behavior of the “material time” u a ( φ, T ) also suggests a similar pattern. Al-though we do not provide further details, let us mention that we have carried out the sameparametrization exercise, with a similar composed functional form, leading to similar accu-racy of the overall fits. Although in this case the power-law exponent does depend slightlyon volume fraction, this parametrization with the exponent at the value of 2, does providestill a third route to determine the crossover temperature T ( φ ). Of course, there is no rea-son to expect full quantitative consistency between these three manners to determine T ( φ ).Nevertheless, they coincide in the general notion that the region bounded by the spinodalcurve T s ( φ ) from above and by the glass-glass transition line T s ( φ ) from below, is dividedin two physically distinct domains. E. High-density “ordinary” liquid-glass transition.
Let us also notice from Fig. 6(a) that besides the development of this low- k spinodal-decomposition peak, the asymptotic structure S a ( k ; φ, T ) also exhibits its main meak at k ≈ π/σ , associated with the short-range nearest-neighbor arrangement among the particles.This main peak also varies with the final temperature T , but in a far less dramatic manner.Thus, we see that the height of the main peak, plotted as the dashed line of Fig. 6(b),increases only slightly as the final temperature increases from 0 to T c , and then decreases27onotonically in the interval T c ≤ T ≤ T s . This behavior has some form of continuation inthe equilibrium regime (above the spinodal temperature) in the height of the main peak ofthe equilibrium static structure factor S ( eq ) ( k ; φ, T ) (soft gray solid line for T ≥ T s ). Thus,if the importance of these two competing modes is measured by the height of the respectivepeak of S a ( k ; φ, T ), then we conclude that for shallow quenches inside the spinodal region,but at volume fractions below φ b , the asymptotic structure is dramatically dominated bythe small- k peak, which diverges as we approach the spinodal line from below. Along thesame isochores, however, deeper quenches (below the glass-glass transition temperature T c )lead to an opposite scenario, in which the asymptotic structure is dominated by the mainpeak of S a ( k ; φ, T ). With the aim of understanding the nature of this glass-glass transition,let us first discuss its relationship with the “ordinary” liquid-glass transition predicted tooccur at T c ( φ ), but for volume fractions above the bifurcation point, φ ≥ φ b .Thus, Fig. 9 presents similar results as Fig. 6, but now at the isochore φ = 0.4, represen-tative of volume fractions well above the volume fraction φ b = 0.33 of the bifurcation point.Clearly, the scenario is now completely different from the previous one. As we can see inFig. 9(a), the main difference is that now the small- k peak is virtually absent. In realitythis peak is still there, but it is imperceptible in the scale of the main figure, as revealed bythe rather extreme zoom in the inset. Thus, in the virtual absence of this spinodal peak,the most visible feature of S a ( k ; φ, T ) is its main peak located at kσ ≈ π . As illustratedby the set of plots of S a ( k ; φ, T ) in Fig. 9(a), the height of this main peak exhibits a rathermild variation as a function of the final temperature T . This is also illustrated by the darkblue dashed line of Fig. 9(b).Let us notice that along the isochores with volume fraction φ < φ b , the dominatingfeature is the existence of the interval T c ≤ T ≤ T s , where the spinodal decomposition peakdetermines the nature of the process of dynamic arrest. This interval naturally defines theother two complementary regions of state space: the high-temperature equilibrium domainabove the spinodal curve, and the low-temperature glass region below T c ( φ ). With thisscenario in mind, it is natural to question if at isochores with volume fraction above φ b , wherethe order T c < T s has been inverted to T s < T c , we will again find three physically distincttemperature regimes: the high-temperature equilibrium domain above T c ( φ ), the window T s ≤ T ≤ T c , and the low-temperature domain below the spinodal curve T s ( φ ). In reality, theNE-SCGLE answer to this question was already given by the behavior of γ a in Fig. 2, which28 k S a ( k ) S i (k)T=0.30T=0.38T=0.46T=0.54T=0.62T=0.70 0 1 2 3 k S a ( k ) S i (k)T=0.30T=0.38T=0.46T=0.54T=0.62T=0.70 φ=0.4 a) T -2 -1 S a ( k m a x ) k main k low k-peak φ=0.40 T c T s b) FIG. 9: (a) Non-equilibrium stationary structure factor S a ( k ; φ, T ) as a function of wave-vector k for a sequence of instantaneous isochoric cooling processes along the isochore φ = 0 . T i = 1 . T -dependence of the low- k peak ( k ≈ /σ ; also zoomed in the inset) and of the main peak ( k ≈ π/σ )of S a ( k ; φ, T ). (b) Height of the low- k peak (dark solid line) and of the main peak (dark dashedline) of S a ( k ; φ, T ), as a function of the final temperature T along the isochore φ = 0 .
4. The softgray solid line represents the diverging behavior (when T → T s from above) of the k = 0 value ofthe equilibrium static structure factor, S ( eq ) ( k = 0; φ, T ). only identifies two physically distinct regions, namely, the high-temperature equilibriumliquid domain, T > T c ( φ ), and the region of dynamically arrested states, 0 ≤ T ≤ T c ,without the spinodal temperature T s ( φ ) having any relevance.In an attempt to better understand this scenario, let us notice that along the isochoreswith volume fraction larger than φ b , we have that the equilibrium static structure factor S ( eq ) ( k ; φ, T ) exists in the interval T ≥ T s ( φ ) ( S ( eq ) ( k = 0; φ, T ) is the soft gray solid linein Fig. 9(b)), and the dynamically-arrested static structure factor S a ( k ; φ, T ) exists in theinterval T ≤ T c ( φ ) (see the dark (blue) solid and dashed lines in Fig. 9(b)). This meansthat both stationary solutions of the time-evolution equation for S ( k ; t ), Eq. (2.1), existsimultaneously in the interval T s ( φ ) ≤ T ≤ T c ( φ ). However, it is not difficult to see [34]that under these circumstances, S a ( k ; φ, T ) is the only possible long-time asymptotic limit of S ( k ; t ). Thus, although S ( eq ) ( k ; φ, T ) exists and is well-defined in the interval T s ( φ ) ≤ T ≤ T c ( φ ), it does not have any influence in the structural relaxation of the system for quenches29n that interval. This is why the spinodal divergence of S ( eq ) ( k = 0; φ, T ) bears no relevanceon the long-time arrested structure of the system, and hence, the spinodal temperature T s ceases to have any dynamical significance in the regime φ > φ b .This means that if we revisit the predicted dynamic arrest diagram in Fig. 4, we mightas well remove completely the dashed portion of the spinodal curve to indicate its dynamicirrelevance and to make clear the predicted continuity of the glass phase throughout the re-gion of state space below the curve T c ( φ ). In the high-packing/high-temperature regime thisglass phase corresponds to the “hard-sphere” or “repulsive” high-packing glasses that resultfrom compressing hard spheres to φ HSg ≈ .
58, whereas in the low-packing/low-temperatureregime it corresponds to “attractive” glasses, in which the intense (compared with k B T ) at-tractive interactions allow the formation of macroscopically homogeneous spanning clusterswith large microscopic heterogeneities involving high local packing coexisting with voids ofsimilar dimensions, to allow mean packing fractions much lower than φ HSg ≈ .
58. Thispredicted scenario is fully consistent, at least qualitatively, with the existence of the glassregion of the experimental scenario reported in Refs. [13, 14], and referred to there as regionIII.
V. PREDICTED SCENARIO AND AGREEMENT WITH EXPERIMENT.
In this section we shall actually revisit and update the dynamic arrest diagram of Fig. 4.The resulting updated diagram will constitute the core of the scenario predicted by the NE-SCGLE theory, whose connection with the reported experimental scenario on the arrestedspinodal decomposition of lysozyme solutions will also be discussed. Let us thus incorporatethe two most relevant findings of the previous section, in the dynamic arrest diagram of Fig.4. The first is concerned with the dynamic irrelevance of the dashed portion of the spinodalcurve. This leads to the predicted continuity of the glass phase throughout the region ofstate space below (and to the right of) the curve T c ( φ ).The second refers to the existence of a crossover temperature T ( φ ) at which severalproperties change their temperature dependence from an exponential to a power-law thatdiverges at the spinodal line T s ( φ ). In Fig. 10 we present the updated dynamic arrestdiagram. There we have indicated the location of the crossover line T ( φ ) obtained in thethree manners described above (i.e., based on the T -dependence of ξ a ( φ, T ), of γ ∗ a ( φ, T ), and30 φ T Spinodal lineGlass-glass transitionLiquid-glass transitionDynamic arrest lineT ξο (φ) T γο (φ) T uo ( φ) I ( φ b ,T b ) II III
FIG. 10: Dynamic arrest diagram with the crossover line T ( φ ) determined by the parametrizationof the transition from exponential to power law of: the typical size ξ a ( φ, T ) of the heterogeneities( T ξ ), the square localization length γ a ( φ, T ) ( T γ ), and the the material time u a ( φ, T ) ( T u ). Eachof these three curves provide an approximate location of the crossover from full phase separation(region I) to gel formation (region II). Region III corresponds to the formation attractive glasses. of u a ( φ, T )). We do not observe, but we did not expect, any quantitative agreement amongthese three results for T ( φ ). However, it is clear that within their quantitative differences,the three of them suggest a common qualitative and semiquantitative scenario. In fact,their spreading may serve to remind us that we are talking about a crossover between twolimiting conditions, and not about a sharp transition between two well-defined states (likein the liquid-glass transition line).The referred two limiting conditions are, for shallow quenches ( T < ∼ T s ), the ineffec-tiveness of the dynamic arrest conditions to prevent full gas-liquid separation, whereas fordeeper quenches, slightly above the glass-glass transition ( T > ∼ T c ( φ )), the effectivenessof the condition of dynamic arrest to interrupt the process of spinodal decomposition atan early stage of its development. The result of this interruption cannot be other thanthe formation of glassy states whose asymptotic structure is strongly determined by theamplified heterogeneities represented by the frozen small- k peak of S a ( k ). These highly het-erogeneous structures contrast with the much more homogeneous structure of the attractiveglasses formed at deeper quenches, below the glass-glass transition ( T < T c ( φ )), in whichthe relevant spatial heterogeneities are determined by the main peak of S a ( k ), and not by31ts insignificant small- k peak. Besides the difference in the scale and degree of the spatialstructural heterogeneity of the two glass phases, we can highlight another relevant difference,this time in the localization length p γ ∗ a ( φ, T ). From the results of the previous sections wecan see that p γ ∗ a ( φ, T ) has the Lindemann-like value p γ ∗ a ( φ, T ) ≈ .
15 for the glass phaseright below T c ( φ ). In contrast, for the glass phase right above T c ( φ ), the correspondingvalue of p γ ∗ a ( φ, T ) is about ten times larger, thus indicating a sudden increase in the degreeof heterogeneity upon crossing the glass-glass transition line T c ( φ ). It is, of course, verydifficult not to associate the fluffier glass phase above T c ( φ ) with gel states, and the (ratherloosely defined) crossover line T ( φ ) with the corresponding gel line.In fact, if we now compare the topology of our updated theoretical diagram with itsexperimental counterpart in Fig. 4 of Ref. [13] (or Fig. 1(a) of Ref. [14]), we immediatelyrecognize that the qualitative mapping of one onto the other is impressive. To enhance thissimilarity with the experimental diagram, here we have also labeled, as in Refs. [13, 14], thethree kinetically distinct regions of the theoretical diagram as I, II and III, since the inter-pretation of each of these theoretical regions, as discussed in this section, are conceptuallyconsistent with the description in Refs. [13, 14] of the corresponding regions experimentallydetermined in the lysozyme system. The most relevant feature shared by the predicted andthe experimental scenarios, is the fact that the “ordinary” liquid-glass transition penetratesthe unstable region as a novel glass-glass transition. For φ > φ b , this glass phase borderswith the equilibrium liquid phase along the liquid-glass transition line T c ( φ ), so that weexpect that the predicted phenomenology of the corresponding liquid-glass transition willbe essentially identical to that of the purely repulsive soft-sphere system studied in Ref.[34]. Although there is still much to be understood about this “ordinary” glass transition,from the study in Ref. [34] we have learned that as soon as we enrich the scenario providedby the dynamic arrest diagram (obtained solely from long-time asymptotic properties like γ ∗ a ( φ, T )) with the kinetic dependence on waiting time, a number of important effects canbe predicted, such as the aging of the structure and the dynamics, and the experimentalimpossibility of observing the sharp and discontinuous liquid-glass transition implied by thedynamic arrest diagram, which can only be observed as a blurred and continuous transition[34].Let us finally mention other two additional points of comparison between our theoreticalscenario and the experimental results of Schurtenberger and collaborators [13, 14]. This32 .1 0.2 0.5 1 T * ξ a T o ( φ) T c ( φ) T *exp ξ a ( µ m ) (a) φ ξ a T=0.45T=0.30T=0.20 0.05 0.1 0.15 0.2 φ ξ a ( µ m ) (b) FIG. 11: Comparison between the theoretical scenario and the experimental results shown inFigs. 3(b) and 4(c) of Ref. [14] for the characteristic length ξ a ( φ, T ) = 2 π/k max ( φ, T ). In themain panel of (a) we show the theoretical results for ξ a as a function of T ∗ = ( T s − T ) /T s at φ = 0 .
15 for the HSAY model with z = 2, parametrized according to Eq. (4.2). In the insetthe dots are the experimental data in Fig. 4(c) of Ref. [14] for the fast quench of the lysozymesolution at the same isochore, and the composed solid-dashed curve is the fit of these data withthe parametrization expression in Eq. (4.2), yielding T = 12 . o C (vertical solid line). Fig. (b)illustrates the experimental observation reported in Fig. 3(b) of Ref. [14] that for fixed finaltemperature ξ a ( φ, T ) seems to be independent of volume fraction (inset). This feature is alsoqualitatively predicted by our theory, as illustrated by the results for ξ a ( φ, T ) along the isotherms T = 0 . , . , .
45 plotted in the main panel (circles, squares and diamonds, respectively). Theexperimental data in the insets were read from Figs. 3(b) and 4(c) of Ref. [14] with the permissionof IOP Publishing. group measured the temperature and density dependence of the heterogeneities length scale ξ a ( φ, T ) of the arrested structure of the gels. They reported that for fixed volume fraction( φ = 0 . ξ a ( φ, T ) could be adjusted by ξ a ( φ, T ) ∝ [ T s − T ] − ν with ν = 0 . ± . ξ a ( φ, T ), measured at fixed temperature, seemed to depend only weakly on volumefraction. Regarding the temperature dependence of ξ a ( φ, T ), in the main panel of Fig.11(a), we present the theoretical results for ξ a ( φ, T ) in the HSAY model with z = 2 alongthe isochore φ = 0 .
15 parametrized, as in Fig. 7, according to Eq. (4.2), but plotted hereas a function of the reduced temperature distance T ∗ ≡ ( T s − T ) /T s to the spinodal curve.33s in Fig. 7, the vertical (dotted and solid) lines represent, respectively, the glass-glasstransition temperature T c ( φ ) and the temperature T ( φ ) indicating the crossover of the T -dependence of ξ a ( φ, T ) from exponential (solid magenta portion) to diverging power-law(dashed blue portion). The prediction is that gels should be expected to form in the interval T c < T < T and to start to phase separate for temperatures above T . In the inset of Fig.11(a) we reproduce the log-log plot of the experimental data reported in Fig. 4(c) of Ref.[14] for the fast quench of the lysozyme solution at the isochore φ = 0 .
15. The composed(solid-dashed) curve is the parametrization expression in Eq. (4.2), in which we fixed theexponent ν at its Ornstein-Zernike value ν = 0 . T s (= 18 o C ), T c (= − . o C ) and ξ a ( φ, T c ) (= 1 . µm ), and varied T to obtain the best fit of these experimental data, yielding T = 12 . o C , a value 16 %below the experimental gel (or tie-line) temperature of 15 o C .The experimental observation that for fixed final temperature ξ a ( φ, T ) seems to be in-dependent of volume fraction is illustrated in the inset of Fig. 11(b), which plots the ex-perimental data of Fig. 3(b) of Ref. [14]. In the main panel of our figure 11(b) we plotthe theoretical ξ a ( φ, T ) as a function of φ at fixed temperature, in the same volume fractionwindow as the experimental report. We observe that indeed, in this particular window, thetheoretical prediction for ξ a ( φ, T ) along the three isotherms presented can be said to dependrather weakly on volume fraction. We must mention, however, that in a more extendedwindow that includes lower densities, the spinodal divergence of ξ a ( φ, T ) will lead to a morecomplex scenario than the simple independence of ξ a ( φ, T ) on volume fraction, as alreadyinsinuated by the curve corresponding to T = 0 .
45. Let us also stress that the comparison inFig. 11 between our theoretical predictions and the experimental results of Ref. [14] is onlymeant to illustrate the qualitative consistency between our theoretical predictions for a spe-cific model system and the results of an experimental study on the real system conceptuallyclosest to our idealized model. Although the difference in the regime theoretically studiedhere (relatively long-ranged attractions, z =2) and the experimental system (much shorter-ranged attractions) prevents a more detailed direct comparison, it is clear that our theoryis able to predict and describe the relevant measurable properties, and that the predictionsare consistent with the available experimental reports.34 I. SUMMARY.
In summary, in this paper we have introduced the Non-equilibrium Self-consistent Gener-alized Langevin Equation theory of irreversible relaxation in liquids, as a general theoreticalframework for the description of the non-equilibrium processes involved in the spinodal de-composition of suddenly and deeply quenched simple fluids. The main result of the paper isthe dynamic arrest diagram in Fig. 10 of the previous section, which is the non-equilibriumanalog of the equilibrium phase diagram in Fig. 1. This dynamic arrest diagram com-plements the equilibrium scenario in Fig. 1 with the kinetic perspective of the NE-SCGLEtheory, which predicts the dynamic arrest of a simple model liquid of attractive hard spheressuddenly and deeply quenched to its unstable region. Because of the non-linear and non-equilibrium nature of the phenomenon described, the initial interpretation of the specificpredictions turned out to be rather subtle. However, the detailed analysis of several proper-ties associated with the long-time asymptotic structure factor S a ( k ; φ, T ) provided a rathercoherent picture, summarized by the resulting dynamic arrest diagram of Fig. 10.The first main feature of the predicted dynamic arrest diagram of attractive hard spheresis the continuity of Region III, corresponding to the glass phase below and to the right ofthe curve T c ( φ ), which contains the high-packing/high-temperature “repulsive” glasses andthe low-packing/low-temperature “attractive” glasses. For φ > φ b , this glass phase borderswith the equilibrium liquid phase along the liquid-glass transition line T c ( φ ). For φ < φ b this liquid-glass transition penetrates the unstable region of the gas-liquid coexistence as aglass-gel transition, with the gel phase constituting region II of the diagram. Finally, weassociated the formation of gels with the early arrest of the growing spatial heterogeneities,whose size ξ a ( φ, T ) = 2 π/k max is determined by the position k max of the spinodal (small- k )peak of S a ( k ; φ, T ). The temperature dependence of this and other properties revealed theexistence of two regimes, separated by a crossover temperature T ( φ ). This crossover linewas conjectured to separate the region of gel formation from the region where the dynamicarrest is not able to prevent full phase separation.This theoretical scenario emerged only form the analysis of long-time asymptotic struc-tural properties, such as S a ( k ; φ, T ). The NE-SCGLE theory, however, provides abundantadditional information, including the detailed non-equilibrium evolution of the static anddynamic properties of the system at any instant after the quench. This kinetic and dynamic35nformation leads to a richer and more detailed scenario whose description, however, is outof the scope of the present paper. Nevertheless, even without this important information,we related the available theoretical scenario with the results of the study of dynamicallyarrested spinodal decomposition in the experimental model system conceptually closest tothe theoretical model analyzed here. As we demonstrated in the previous section, the mostrelevant predicted features are in satisfactory qualitative agreement with the experimentalobservations.Of course, rather than closing the theoretical discussion of arrested spinodal decomposi-tion, the results presented in this paper only open a number of question to further discussion.For example, we expect to expand this discussion to include the kinetic aspects not discussedhere, as well as the dependence of the resulting dynamic arrest scenario on the main fea-tures of the interaction potential, such as its spatial range and its specific form (attractiveYukawa vs. square well, hard-sphere vs. soft sphere, Lennard-Jones, etc.). Similarly, wehave not exhausted all the pertinent tests and comparisons of the theoretical predictionswith available relevant simulations, such as the thorough simulation study of Testard, Koband Berthier [19]. In future communications, however, we shall address these and otherissues that we could not properly treat in the present paper.ACKNOWLEDGMENTS: This work was supported by the Consejo Nacional de Cienciay Tecnolog´ıa (CONACYT, M´exico), through grants No. 182132 and 242364. Appendix. The NE-SCGLE theory: main approximations.
Let us have in mind a simple monocomponent liquid formed by N Brownian particlesin a volume V , which diffuse with a self-diffusion coefficient D between collisions whileinteracting through a pair potential u ( r ). The NE-SCGLE theory was originally derivedfrom a non-equilibrium extension of Onsager’s theory of thermal fluctuations [32]. This36bstract theory provides the general structure of the time-evolution equations for the meanvalue and the covariance of the thermal fluctuations of the macroscopic variables of a non-equilibrium system. Applied to the local number concentration of particles in our Brownianliquid, these two equations become ∂n ( r , t ) ∂t = D ∇ · b ( r , t ) n ( r , t ) ∇ βµ [ r ; n ( t )] (1)and ∂σ ( r , r ′ ; t ) ∂t = D ∇ · n ( r , t ) b ( r , t ) ∇ Z d r E [ r , r ; n ( t )] σ ( r , r ′ ; t )+ D ∇ ′ · n ( r ′ , t ) b ( r ′ , t ) ∇ ′ Z d r E [ r ′ , r ; n ( t )] σ ( r , r ; t ) − D ∇ · n ( r , t ) b ( r , t ) ∇ δ ( r − r ′ ) . (2)which describe the non-equilibrium evolution of the mean value n ( r , t ) and the covariance σ ( r , r ′ ; t ) ≡ δn ( r , t ) δn ( r ′ , t ) of the fluctuations δn ( r , t ) = n ( r , t ) − n ( r , t ) of the local concen-tration profile n ( r , t ) of a colloidal liquid.There are two fundamental elements entering in these equations. The first involves therelevant thermodynamic information, and is the electrochemical potential µ [ r ; n ( t )] of onecolloidal particle at position r , from which the thermodynamic stability function E [ r , r ′ ; n ] ≡ [ δβµ [ r ; n ] /δn ( r ′ )] derives [1, 3] ( k B T = 1 /β being the thermal energy per particle). The otheris the time-dependent local mobility function b ( r , t ), which in reality is also a functional of n ( r , t ) and σ ( r , r ′ ; t ), thus introducing a non-linear character to these equations. In whatfollows we provide a more detailed definition of these objects.
1. Thermodynamic elements.
The electrochemical potential βµ [ r ; n ( t )] is an ordinary function of the particle’s position r and a functional of the local concentration profile n ( r , t ). Although not explicitly stated, µ [ r ; n ( t )] is also a functional of the temperature profile T ( r , t ). Here, however, we shallonly have in mind uniform temperature profiles, T ( r , t ) = T ( t ), which will play the roleof the external control parameter describing the imposed thermal treatment of the system.The functional dependence of βµ [ r ; n ] on the arbitrary profile n ( r ) embodies the chemicalequation of state of the system, and is a fundamental thermodynamic input of the presenttheory. The functional derivative E [ r , r ′ ; n ] ≡ [ δβµ [ r ; n ] /δn ( r ′ )] is in general an ordinary37unction of both, r and r ′ , and a functional of n ( r , t ). Evaluated in the absence of externalfields and at an uniform profile n ( r ) = n , this two-point thermodynamic property onlydepends on the distance | r − r ′ | and its functional dependence on the profile n ( r ) becomesan ordinary dependence of the scalar argument n , i.e., E [ r , r ′ ; n ] = E h ( | r − r ′ | ; n ). The so-called direct correlation function c ( r ; n ) is defined in terms of this thermodynamic functionby the relationship E h ( r ; n ) = δ ( r ) /n − c ( r ; n ).In order to simplify the analysis of the equations above, at this point we introducethe “local-homogeneity” approximation, in which we approximate this two-point thermo-dynamic property by E [ r , r ′ ; n ( t )] = E h ( | r − r ′ | ; n ( r , t )), i.e., by the thermodynamic matrixevaluated at an uniform concentration profile, whose constant value is the local and instan-taneous concentration n ( r , t ). Under these conditions we expect that the covariance can alsobe approximated as σ ( r , r ′ ; t ) ≈ σ ( | r − r ′ | ; r , t ), so that Eq. (2) can also be written as ∂σ ( k ; r , t ) ∂t = − k D b ( r , t ) n ( r , t ) E h ( k ; n ( r , t )) σ ( k ; r , t )+ 2 k D b ( r , t ) n ( r , t ) , (3)where σ ( k ; r , t ) ≡ R d ke i k · x σ ( | x | ; r , t ) and E h ( k ; n ( r , t )) ≡ (2 π ) − R d ke − i k · x E h ( | x | ; n ( r , t )).To actually solve these equations, we must indicate the specific nature of the system, i.e.,the specific inter-particle interactions, embodied in the pair interaction potential u ( r ), andthe applied external force F ( r ) = −∇ ψ ( r ) acting on a particle at position r , whose effectsare introduced through the chemical equation of state. This is written in general as [37], βµ [ r ; n ] ≡ βµ ∗ ( β ) + ln n ( r ) + βψ ( r ) − c [ r ; n ] . (4)The first two terms on the right side are the ideal gas contribution to the chemical po-tential, whereas the term − c [ r ; n ] contains the deviations from ideal behavior due to inter-particle interactions. Thus, the functional dependence of c [ r ; n ] on the concentration pro-file n ( r ) must be specified. One possibility is to propose a theoretical approximation forthis dependence, which in the language of density functional theory [37] is actually equiv-alent to proposing an approximate free energy functional. We quote, for example, thesimplest such approximation, referred to as the Debye-H¨uckel or random phase approxima-tion, in which c [ r ; n ] is written as c ( RP A ) [ r ; n ] = − β R d r ′ u ( | r − r ′ | ) n ( r ′ ). An approximate c [ r ; n ] leads to a corresponding approximate thermodynamic matrix E [ r , r ′ ; n ]. For example, E ( RP A ) [ r , r ′ ; n ] = δ ( r − r ′ ) /n ( r ) + βu ( | r − r ′ | ). Other more elaborate approximations areavailable, however, in the literature of the equilibrium theory of inhomogeneous fluids [37].38e have assumed so far that the external fields are static and represented by ψ ( r ). Thereis, however, no fundamental reason why we have to restrict ourselves to these conditions.In fact, the general equations of the NE-SCGLE theory above can be used, within therange of validity of the underlying assumptions, to describe the response of the system toprescribed time-dependent external fields ψ ( r , t ) or to programmed thermal or mechanicalconstraints described by a time-dependent temperature T ( t ) of the reservoir and/or animposed time-dependent total volume V ( t ) of the system. Here we shall assume, however,that such time-dependent fields and constraints could have been used to drive the systemto a prescribed initial state with mean concentration profile n ( r ) and covariance σ ( k ; r ),but that afterward the external fields are set constant, ψ ( r , t ) = ψ ( r ).
2. Absence of external fields and spatial heterogeneities.
The NE-SCGLE Eqs. (1) and (3) are supposed to predict how the system relaxes to itsfinal equilibrium state whose mean profile and covariance are n eq ( r ) and σ eq ( k ; r ). Describingthis response at the level of the mean local concentration profile n ( r , t ) is precisely the aimof the recently-developed dynamic density functional theory (DDFT) [38, 39], whose centralequation is recovered from our theory by setting b ( r , t ) = 1 in Eq. (1). DDFT has beenapplied to a variety of systems, including the description of the irreversible sedimentation ofreal and simulated colloidal suspensions [40]. Similarly, the time-evolution equation for thecovariance, within the additional small-wave-vector approximation E h ( k ) ≈ E + E k + E k ,can be recognized as the kinetic equation that governs the irreversible evolution of the staticstructure factor of a suddenly quenched liquid, developed in the description of the earlystages of spinodal decomposition (see, for example, Eq. (2.11) of Ref. [8], in which E = 0,or Eq. (23) of Ref. [10]). In fact, the main motivation of the present work is to contributeto the development of the theory of spinodal decomposition starting from the general NE-SCGLE theory just introduced.The most essential aspect of the irreversible process of spinodal decomposition is thedevelopment of strong spatial heterogeneities. These originate in the thermodynamic spin-odal instability of thermal fluctuations of certain wavelengths, whose amplification generatesstrong spatial heterogeneities, leading eventually to the full separation of the system in twoseparate regions. At least in its initial regime, the main features of this irreversible phase39eparation process must be exhibited by the time-evolution of the two observable variablesinvolved in the present theory, namely, the mean value n ( r , t ) and the covariance σ ( k ; r , t ), asrecognized since the early work of Cahn and Hilliard [5] and Cook [6]. Our NE-SCGLE equa-tions actually contain Cook’s time-evolution equation for the non-stationary static structurefactor S ( k ; t ), but goes beyond this classical theory by incorporating the possibility of dy-namic arrest.To see this, let us imagine that we subject our system to a prescribed protocol of thermaland/or mechanical treatment, described by the chosen temporal variation T ( t ) and V ( t ) ofthe reservoir temperature and total volume, but in the absence of external fields, ψ ( r , t ) = 0.We also assume that the thermal diffusivity of the system is sufficiently large that any tem-perature gradient is dissipated almost instantly (compared with the relaxation of chemicalpotential gradients), so that the temperature field inside the system remains uniform andinstantly equal to the imposed temperature of the reservoir, T ( r , t ) = T ( t ).In contrast, since the essence of spinodal decomposition is the development of strongspatial heterogeneities, in this case we cannot in general assume that the mass diffusivityof the system is sufficiently large that any gradient of chemical potential can be dissipatedalmost instantly. If this were the case, the density field n ( r , t ) inside the system wouldalways remain uniform and instantly equal to the imposed bulk density, i.e., n ( r , t ) = n ( t ) ≡ N/V ( t ). Under such conditions, rather than solving Eqs. (1) and (3) for n ( r , t ) and σ ( k ; r , t ),we would only have to solve Eq. (3) for σ ( k ; r , t ), with n ( r , t ) substituted by its imposed bulkvalue n ( t ), which now becomes a control parameter. Although the resulting scheme mightseem too gross for the description of the spontaneous spatial heterogeneities represented bythe neglected deviations ∆ n ( r , t ) ≡ n ( r , t ) − n ( t ) (which are definitely non-zero), here weshall see that it provides a welcome mathematical and numerical simplification of Eq. (3)for σ ( k ; r , t ) ≈ σ ( k ; t ), which may now be written as ∂σ ( k ; t ) ∂t = − k D n ( t ) b ( t ) E h ( k ; n ( t )) σ ( k ; t ) + 2 k D n ( t ) b ( t ) . (5)40 . Determination of the t -dependent mobility b ( t ) . As explained in Refs. [32] and [33], the time-dependent mobility b ( t ) is in reality a functional of σ ( k ; t ), and is determined by the generalized Einstein relation b ( t ) = [1 + Z ∞ dτ ∆ ζ ∗ ( τ ; t )] − , (6)with the t -evolving, τ -dependent friction function ∆ ζ ∗ ( τ ; t ) representing the effects of theinterparticle interactions on the total friction coefficient of a representative tracer particle.The functional dependence of ∆ ζ ∗ ( τ ; t ) on σ ( k ; t ) introduces the essential non-linearities ofthe problem, and this functional dependence can only be determined in an approximatemanner.If we were dealing with an equilibrium system, we could simply adopt the correspond-ing approximation proposed by mode coupling theory or by the equilibrium version of theSCGLE theory. The non-equilibrium SCGLE theory extends the latter to non-equilibriumconditions, thus leading to an approximate expression for ∆ ζ ∗ ( τ ; t ), namely,∆ ζ ∗ ( τ ; t ) = D π n ( t ) Z d k k (cid:20) S ( k ; t ) − S ( k ; t ) (cid:21) × F ( k, τ ; t ) F S ( k, τ ; t ) (7)in terms of the non-stationary static structure factor S ( k ; t ) ≡ σ ( k ; t ) /n ( t ) and of the col-lective and self non-equilibrium intermediate scattering functions F ( k, τ ; t ) and F S ( k, τ ; t ).In Ref. [32] it was indicated how to extend the exact memory-function expressions for F eq ( k, τ ) and eq F S ( k, τ ) to non-equilibrium conditions, thus leading to the correspondingexact memory-function expressions for F ( k, τ ; t ) and F S ( k, τ ; t ). It is at the level of thecorresponding memory functions where one has to introduce the essential approximations,also extending the arguments of the SCGLE theory to non-equilibrium conditions. This con-verts the exact memory function expressions for F ( k, τ ; t ) and F S ( k, τ ; t ) into the following,approximate memory-function equations, written in terms of the Laplace transforms (LT) F ( k, τ ; t ) and F S ( k, τ ; t ), as F ( k, z ; t ) = S ( k ; t ) z + k D S − ( k ; t )1+ λ ( k ) ∆ ζ ∗ ( z ; t ) , (8)and F S ( k, z ; t ) = 1 z + k D λ ( k ) ∆ ζ ∗ ( z ; t ) , (9)41ith λ ( k ) being a phenomenological “interpolating function” [30], given by λ ( k ) = 1 / [1 + ( k/k c ) ] , (10)where the phenomenological cut-off wave-vector k c depends on the system considered. Inthe particular case of an instantaneous isochoric quench, which is the subject of the presentpaper, these equations become Eqs. (2.1)-(2.6) of Sect. II. [1] H. Callen, Thermodynamics , John Wiley, New York(1960).[2] D. A. McQuarrie
Statistical Mechanics , Harper & Row (New York, 1973).[3] J. Keizer,
Statistical Thermodynamics of Nonequilibrium Processes , Springer-Verlag (1987).[4] G. Lebon, D. Jou, and J. Casas-V´azquez,
Understanding Non-equilibrium ThermodynamicsFoundations, Applications, Frontiers , Springer-Verlag Berlin Heidelberg (2008).[5] J. W. Cahn and J. E. Hilliard, J. chem. Phys. , 688 (1959).[6] H. E. Cook, Acta Metall. , 297 (1970).[7] H. Furukawa, Adv. Phys., , 703 (1985).[8] J. S. Langer, M. Bar-on, and H. D. Miller, Phys. Rev. A , 1417 (1975).[9] J. K. G. Dhont, J. Chem. Phys. , 5112 (1996).[10] S. B. Goryachev, Phys. Rev. Lett. , 1850 (1994).[11] P. J. Lu, E. Zaccarelli, F. Ciulla, A. B. Schofield, F. Sciortino and D. Weitz, Nature , 499(2008).[12] E. Sanz, M. E. Leunissen, A. Fortini, A. van Blaaderen, and M. Dijkstra, J. Phys. Chem. B , 10861 (2008).[13] F. Cardinaux, T. Gibaud, A. Stradner, and P. Schurtenberger, Phys. Rev. Lett. , 118301(2007).[14] T. Gibaud and P. Schurtenberger, J. Phys.: Condens. Matter , 322201 (2009).[15] L. Di Michele, D. Fiocco, F. Varrato, S. Sastry, E. Eisera and G. Foffi, Soft Matter , 3633(2014).[16] Y. Gao, J. Kim and M. E. Helgeson. Microdynamics and arrest of coarsening dur-ing spinodal decomposition in thermoreversible colloidal gels. Soft matter (2015) DOI: , 323101 (2007).[18] J. F. M. Lodge and D. M. Heyes, J. Chem. Soc., Faraday Trans., , 437 (1997).[19] V. Testard, L. Berthier, and W. Kob, J. Chem. Phys. , 164502 (2014).[20] G. Foffi, C. De Michele, F. Sciortino, and P. Tartaglia, J. Chem. Phys. , 224903 (2005).[21] L. Berthier and G. Biroli, Rev. Mod. Phys. (2011).[22] W. G¨otze, in Liquids, Freezing and Glass Transition , edited by J. P. Hansen, D. Levesque,and J. Zinn-Justin (North-Holland, Amsterdam, 1991).[23] W. G¨otze and L. Sj¨ogren, Rep. Prog. Phys. , 241 (1992).[24] W. G¨otze and E. Leutheusser, Phys. Rev. A , 2173 (1975).[25] W. G¨otze, E. Leutheusser and S. Yip, Phys. Rev. A , 2634 (1981).[26] A. Latz, J. Phys.: Condens. Matter, (2000) 6353.[27] L. F. Cugliandolo and J. Kurchan, Phys. Rev. Lett. , 173 (1993).[28] P.E. Ram´ırez-Gonz´alez et al. , Rev. Mex. F´ısica , 327 (2007).[29] L. Yeomans-Reyna, M. A. Ch´avez-Rojo, P. E. Ram´ırez-Gonz´alez, R. Ju´arez-Maldonado, M.Ch´avez-P´aez, and M. Medina-Noyola, Phys. Rev. E , 041504 (2007)[30] R. Ju´arez-Maldonado et al. , Phys. Rev. E , 062502 (2007).[31] P. E. Ram´ırez-Gonz´alez and M. Medina-Noyola, J. Phys.: Cond. Matter : 504103 (2009).[32] P. E. Ram´ırez-Gonz´alez and M. Medina-Noyola, Phys. Rev. E , 061503 (2010).[33] P. E. Ram´ırez-Gonz´alez and M. Medina-Noyola, Phys. Rev. E , 061504 (2010).[34] L. E. S´anchez-D´ıaz, P. E. Ram´ırez-Gonz´alez, and M. Medina-Noyola, Phys. Rev. E , 052306(2013).[35] L. E. S´anchez-D´ıaz, E. L´azaro-L´azaro, J. M. Olais-Govea and M. Medina-Noyola, J. ChemPhys. , 234501 (2014).[36] L. S. Ornstein and F. Zernike, Proc. K. Ned. Akad. Wet. , 793 (1914).[37] R. Evans, Adv. Phys. : 143(1979).[38] U. Marini Bettolo Marconi and P. Tarazona, J. Chem. Phys. , 8032 (1999); ibid., J. Phys.:Condens. Matter , A413 (2000)[39] A. J. Archer and M. Rauscher, J. Phys. A , 9325 (2004).[40] C. P. Royall et al., Phys. Rev. Lett. , 188304 (2007).[41] L. L´opez-Flores et al., EPL, , 46001 (2012).
42] R. V. Sharma and K. C. Sharma, Physica A , 213 (1977).[43] J. K. Percus and G. J. Yevick, Phys. Rev. , 1 (1957).[44] L. Verlet and J. J. Weis Phys. Rev. A , 939 (1972).[45] J. P. Hansen and L. Verlet, Phys. Rev.
151 (1969).151 (1969).