Intermittent dynamics and logarithmic domain growth during the spinodal decomposition of a glass-forming liquid
IIntermittent dynamics and logarithmic domain growth during the spinodaldecomposition of a glass-forming liquid
Vincent Testard, Ludovic Berthier, and Walter Kob
Laboratoire Charles Coulomb, UMR 5221 CNRS and Universit´e Montpellier 2, 34095 Montpellier, France (Dated: September 25, 2018)We use large-scale molecular dynamics simulations of a simple glass-forming system to investigatehow its liquid-gas phase separation kinetics depends on temperature. A shallow quench leads to afully demixed liquid-gas system whereas a deep quench makes the dense phase undergo a glass tran-sition and become an amorphous solid. This glass has a gel-like bicontinuous structure that evolvesvery slowly with time and becomes fully arrested in the limit where thermal fluctuations becomenegligible. We show that the phase separation kinetics changes qualitatively with temperature, themicroscopic dynamics evolving from a surface tension-driven diffusive motion at high temperatureto a strongly intermittent, heterogeneous and thermally activated dynamics at low temperature,with a logarithmically slow growth of the typical domain size. These results shed light on recentexperimental observations of various porous materials produced by arrested spinodal decomposition,such as nonequilibrium colloidal gels and bicontinuous polymeric structures, and they elucidate themicroscopic mechanisms underlying a specific class of viscoelastic phase separation.
I. INTRODUCTION
The behavior found in nonequilibrium kinetic phenom-ena such as self-assembly [1], pattern formation [2], andphase ordering kinetics [3] is typically much richer thanthe one encountered in equilibrium since a broader rangeof morphologies and a more complex relaxation dynamicscan be observed in systems that are far from equilibrium.While this has been an active research area since manydecades [2, 4], the field has seen a surge of interest in thelast few years, since progress in the synthesis of colloidalparticles with complex shapes and tunable pairwise in-teractions permits the self-assembly of materials of evergrowing complexity, i.e. systems that exist only if onemasters their intricate formation process [1].One of the most important phenomenon occurring innonequilibrium systems is spinodal decomposition [3, 5].The case of shallow temperature quenches has been stud-ied in great detail using theory as well as experiments.Various regimes for the time evolution of the average do-main size have been predicted, and observed in controlledexperiments and simulations of simple fluids using a mul-titude of computational approaches [6–10]. However, amuch greater complexity can be expected if the densephase is not a simple liquid [11] but is itself a “com-plex” fluid [12]. Of particular importance is the so-called“viscoelastic phase separation” [13, 14] which is charac-terized by a strong physical (mainly, rheological) asym-metry between the two coexisting phases. For instance,one can study the phase separation between two fluidsof unequal viscosities, or the coexistence between a solidand a fluid phase. In the following we shall be concernedwith a situation intermediate between these two, whereone component will evolve from being a simple fluid tobecome a highly viscous liquid or an amorphous glassphase, thus giving rise to the phenomenon of a glass-gasphase separation [15–18].The phase ordering process between a fluid and anamorphous solid can be expected to display a complex phenomenology, since even the equilibrium bulk behav-ior of (homogeneous) amorphous glasses is not well un-derstood [19, 20]. When suddenly quenched from hightemperature to below the glass temperature, a glass-forming material evolves slowly with time, undergoing anonstationary aging dynamics characterized by intermit-tent, heterogeneous dynamics occurring far from equilib-rium [20–24]. It is not clear how this aging glass statewill evolve if it is given by the dense phase in the complexbicontinuous structure formed during spinodal decompo-sition. Since a glass behaves mechanically like a solid, onecan expect that the bicontinuous structure formed afterthe quench into the coexistence region becomes rigid, andas a result will become kinetically arrested into a bicon-tinuous porous structure [16]. In that case, the glass-gas phase separation would thus be a conceptually verysimple way of producing porous media. However, onecan also expect that the aging of the glass phase enablesintermittent, thermally activated microscopic rearrange-ments, which could potentially make the porous materialvery fragile. To advance our understanding regardingthese questions we present here the results of our studyconcerning the interplay between phase separation kinet-ics and aging behavior resulting from deep quenches atlow temperatures in a simple glass-forming model.Although our main motivation is to obtain a funda-mental understanding of the glass-gas phase separationkinetics, there are also several experimental considera-tions that motivate such a study. Firstly, a number ofcolloidal systems and protein solutions can be modelledas spherical particles with nearly hard-core repulsion andlonger-ranged attraction, whose range and strength canbe controlled. Thus, they will undergo a phase sepa-ration in some part of the phase diagram, which mightpossibly interfere with the colloidal glass formation oc-curring at higher densities. Therefore, the idea that, atleast in some materials, gelation results from a kineti-cally arrested glass-gas spinodal decomposition has beenexplored in several experimental works [25–33]. This pro- a r X i v : . [ c ond - m a t . d i s - nn ] S e p cess has also been the subject of a number of numeri-cal studies, mostly aimed at reproducing realistic coarse-grained pair interactions characterizing colloidal systemsthat were specifically studied experimentally [27, 34–39].More recent studies, in line with our own work [15], haveconsidered more diverse systems, such as the phase sep-aration kinetics of a coarse-grained model for buckyballC carbon molecules [40]. Moreover, bicontinuous dis-ordered structures reminiscent of the ones obtained inspinodal decompositions may also be found in colorfulbird feathers, and were recently interpreted as incom-pletely phase separated polymeric glasses [41, 42], thusdemonstrating the broad relevance of the problem of theglass-gas phase separation. Finally, even more complexmixtures that are relevant to food processing [43] or solarcell technology [44–46] also exhibit kinetically arrestedspinodal decompositions, that might result from the factthat one component becomes mechanically rigid.There are several possible ways to study coarseningprocesses by means of theory [3]. One efficient approachis to study a coarse-grained model of a biphasic materialusing a Ginzburg-Landau free energy functional of thetwo-phase system complemented with phenomenologicaldynamical equations to incorporate relevant dynamicaland mechanical properties of each phase [13, 47]. For sim-ple fluids, this amounts to studying model H in the classi-fication scheme proposed by Hohenberg and Halperin [4].Possible extensions to viscoelastic materials have indeedbeen considered in the past [13], and numerically inte-grated to obtain insights into some specific viscoelasticphase separation processes. Still, it remains difficult tofaithfully incorporate in this approach the complex (typ-ically highly nonlinear and history dependent) physicalproperties of real glass-forming materials in their agingregime. Moreover, such coarse-grained equations cannotprovide direct information on the microscopic dynamicsresponsible for the evolution of such bicontinuous mate-rials.To avoid this drawback we use here a microscopicallyrealistic description of the homogeneous glass using anatomistic model combined with molecular dynamics tech-niques. We study its behavior during phase separationover a broad range of control parameters, mainly densityand temperature. Indeed, numerous successful atomic-scale simulations of the simpler situation of a liquid-gasspinodal decomposition have been reported [6, 8, 48, 49].In Ref. [50], a Lennard-Jones system was quenched to lowtemperature in the coexistence region, and the resultingcrystal-gas phase separation was studied, but no kineticarrest was reported (see Ref. [51] for a recent related ex-perimental study). Simulations of realistic colloidal in-teractions have been reported in Refs. [34, 35, 37, 38],but the quenches have been performed to very lowtemperatures and hence particles aggregate nearly irre-versibly and thermal fluctuations play little role, whichcorresponds effectively to zero-temperature quenches inour approach. Thus, a careful study of the crossoverregime between ordinary spinodal decomposition and ir- reversible aggregation is so far not available.In this work, we fill this gap and provide a detailed nu-merical study of the phase separation kinetics between agas and a glass-forming material at various temperaturesencompassing the glass transition of the bulk material.We describe how the phase separation kinetics changesqualitatively with temperature, the microscopic dynam-ics evolving from the well-known diffusive motion drivenby surface tension for shallow quenches, to a qualitativelydifferent coarsening regime in which dynamics becomesstrongly intermittent, spatially heterogeneous and ther-mally activated at low temperature, leading to logarith-mically slow growth of the typical domain size. A shortaccount of our results has been published [15].Our paper is organized as follows. In Sec. II we definethe model, provide technical details about our numericalsimulations, and discuss the phase diagram of the sys-tem and the relevant parameters to be explored in thiswork. In Sec. III we provide a qualitative descriptionof the temperature influence on the spinodal decompo-sition kinetics. In Sec. IV we present several structuralcharacterizations of the bicontinuous structures and inSec. V we discuss the time evolution of these structures,characterizing in detail the influence of temperature onthe growth law. In Sec. VI we provide insights into themicroscopic mechanisms responsible for the coarseningstructures at high and low temperatures. In Sec. VII wediscuss the nature of the coexistence line below the glasstransition temperature, that was the subject of a recentcontroversy in the literature. In Sec. VIII, we close thepaper with some perspectives for future work. II. MODEL AND PHASE DIAGRAM
In this section, we describe details of the Lennard-Jones model used in this study and provide some tech-nical informations about the numerical simulations. Wethen describe the relevant features of the phase diagramof the model. Finally we introduce a coarse-grained den-sity field that will be useful to the analysis of the simu-lations.
A. Model and technical details
To study the interplay between liquid-gas phase sepa-ration and the liquid-glass transition we consider a sim-ple Lennard-Jones model for a liquid that was first de-vised to study the dynamics of glass-forming materialsin the bulk [52]. The model is a 80:20 binary mixtureof Lennard-Jones particles with asymmetric interactionparameters chosen such that the minority componentfrustrates, and therefore efficiently prevents, the crys-tallization of the majority component. The details ofthe interaction parameters are as in the original pub-lication [52]. In the following we use Lennard-Jonesunits corresponding to the majority component, express-ing length in units of the particle diameter, σ , and timein units of τ = (cid:112) mσ /(cid:15) , where m is the particle massand (cid:15) the energy scale of the Lennard-Jones interactionbetween particles of the majority component.We integrate Newton’s equation of motion usingLAMMPS [53] for N particles enclosed in a volume V ,working with periodic boundary conditions. We work atconstant number density, ρ = N/V , and adjust the tem-perature T using periodic velocity rescaling as a simplethermostatting procedure. The equations of motion areintegrated using a standard velocity Verlet scheme, us-ing a discretization time step of 0 .
01 in reduced Lennard-Jones units.To study the kinetics of phase separation, we first pre-pare homogeneous samples at the desired density ρ , work-ing at high temperature, T = 3 .
0, well above the criticalpoint T c ≈ .
2, until thermal equilibrium is reached. Wethen instantaneously quench the temperature T to thedesired final value, where the phase separation dynam-ics starts. In the following we will denote as the “age”of the system the time elapsed since the quench to thefinal temperature. To improve the statistics of our quan-titative measurements and the robustness of our findings,we have repeated simulations at each state point using 5to 10 samples, using independently prepared samples athigh temperatures.Specific attention has been paid to system size. Whilenumerical studies of the glass transition in the homoge-neous liquid typically require simulating about N = 10 particles, we have found that up to N = 10 particles areneeded to obtain results devoid of finite size effects dur-ing the phase separation. We have obtained most of ourquantitative results using N = 3 · particles. Whereappropriate, we will discuss the N -dependence of our nu-merical results. B. Coarse-graining the density
Since the bicontinuous structures produced during aspinodal decomposition are characterized by a typicallength scale that is often much larger than the typicalparticle size, the local structure of the fluid is largely irrel-evant. Therefore, it will prove useful to first coarse-grainthe density field before quantifying the spatial fluctua-tions of the obtained field [54, 55]. This coarse-graineddensity will also be used to facilitate the visualizations ofthe particle configurations.We start from the microscopic density field, which isdefined as ρ micro ( r ) = N (cid:88) i =1 δ ( r − r i ) , (1)where r i is the position of particle i . Our first step isto discretize space by dividing it into boxes of linear size ξ b , so that continuous space is now replaced by a discrete ρ _ p ( ρ _ ) T=1.75 T=0.1 ρ =0.6 FIG. 1: Probability distribution function of the coarse-grained density field for ρ = 0 . T = 1 .
75 and for a phase separating systemat T = 0 .
1. The dashed lines are a fit with a Gaussian. lattice containing
V /ξ sites. We consider a discrete den-sity, ρ ( r ), defined for discrete positions r located at thecenter of the boxes described above, as ρ ( r ) = 34 πξ N (cid:88) i =1 θ ( ξ s − | r − r i | ) , (2)where θ ( x ) is the Heaviside function and a second coarse-graining length, ξ s , is introduced. Thus, the density atposition r takes into account all particles located in asphere of radius ξ s centered at r .Finally, we obtain the desired coarse-grained density¯ ρ ( r ) by using the following weighted average of ρ ( r ) overthe boxes surrounding r :¯ ρ ( r ) = 18 (cid:2) ρ ( r ) + ρ ( r + ξ b e x ) + ρ ( r + ξ b e y ) + ρ ( r + ξ b e z )+ ρ ( r − ξ b e x ) + ρ ( r − ξ b e y ) + ρ ( r − ξ b e z ) (cid:3) . (3)Here e α is a unit vector in direction α .The procedure described by Eqs. (2, 3) is easy to im-plement. It returns for each lattice site a coarse-graineddensity field which is a much smoother function than themicroscopic density field ρ micro ( r ) in Eq. (1). The twocoarse-graining length scales ξ b and ξ s can be adjustedby seeking a compromise between having a smooth fieldwithout losing too much spatial resolution. After havingtried several combinations [56], we have settled to thevalues ξ b = σ/ ξ s = σ , so that the lattice spac-ing is equal to the particle radius, while the density fieldis coarse-grained by taking into account the immediateneighbourood of each particle.Having defined the coarse-grained density field ¯ ρ ( r ),we can now easily transform a set of particle coordinatesfrom a single particle configuration into the probabilitydistribution of the coarse-grained density, p (¯ ρ ). In Fig. 1,we show this distribution for two configurations obtainedat density ρ = 0 .
6. The first example is measured for T = 1 .
75, where the system is in the homogeneous fluidphase. As shown by the dashed line, the distribution iswell described by a Gaussian functional form with a max-imum located at ¯ ρ = 0 .
6, as expected. More interestingis the second case at T = 0 . p (¯ ρ ) directly reflects this phase coexistence, since itis characterized by two peaks. One peak is located atvery low density, representative of the gas phase (whosevery small density is not resolved in the scale chosen inFig. 1). A second peak corresponds to the dense phaseand has a maximum near ρ ≈ . ρ and p (¯ ρ ), between cells that belong toeither of the two phases, and those belonging to the in-terfaces. To this end, we need to determine a thresholddensity delimiting the dense phase from the gas phase.By careful visual inspections of several configurations atvarious state points, we have chosen ¯ ρ = 0 .
42 as givingthe most faithful representation of the particle configura-tions. Thus, we define from now on cells with ¯ ρ > .
42 asbelonging to the fluid phase, and cells with ¯ ρ < .
42 asthose forming the gas phase. Interfaces are representedby cells that are in the dense phase and that have at leastone neighboring cell that is not in the dense phase. Wehave used for instance these definitions to produce theimages shown in Figs. 4 and 5, where only cells belong-ing to the interfaces were shown, using different colorsfor gas and fluid phases.
C. Phase diagram
Using a combination of direct visual inspection ofequilibrium configurations coupled to more quantitativemethods to analyze the morphology of biphasic atomisticconfigurations (as described above), we have determinedthe phase diagram shown in Fig. 2, which we comple-ment with some relevant literature data. In this phasediagram, the control parameters are the temperature, T ,and the number density, ρ .When density is high enough, roughly ρ ≥ .
2, the ρ T gas+liquid liquid glassgas+glass b i noda l s p i noda l g l a ss t r an s i t i on li ne gas FIG. 2: Temperature-density phase diagram of the binaryLennard-Jones mixture, showing the fluid, liquid and glasshomogeneous phases. The “glass transition line” has been ob-tained from a mode-coupling analysis of the glassy dynamicsin the liquid phase [57]. The “binodal” line separates homoge-neous from biphasic states. It is obtained by direct inspectionof the state points marked by symbols at high temperature,and from the evaluation of the density in the dense phase ofbiphasic states at low temperature. The “spinodal” line istaken from Ref. [58]. system is always homogeneous. Similarly, the system isa homogeneous fluid at high temperature, T ≥ T c ≈ . T c ≈ . ρ ≤ .
2, the homogeneous system is unstable andphase coexistence is observed, see Fig. 2. To determinethe shown binodal line we have performed quenches to anumber of state points (symbols in the figure) and an-alyze whether the system remains homogeneous at verylong times, i.e. we determined whether or not the distri-bution of coarse-grained density discussed in the previoussubsection has only one peak. Very close to the binodal,we payed attention to metastability and hysteretic ef-fects and performed additional numerical tests to assessthe location of the binodal [56]. Unfortunately, this di-rect method becomes inefficient if temperature becomessmall, T ≤ .
4, because complete phase separation doesnot occur in the limited time window of our simulations,and hysteresis effects become more pronounced at lowertemperatures. In this region, therefore, the coexistenceline has been determined by performing quenches wellwithin the coexistence region, and by measuring the av-erage density of the dense phase of the biphasic config-urations, using the method described above in Sec. II B.We have made sure that both methods yield consistentresults in the vicinity of T = 0 . III. QUALITATIVE OVERVIEW OF RESULTS
In this section we present a general overview of the dis-tinct types of morphologies obtained in the course of ournumerical studies. We then describe qualitative aspectsof the kinetics of the phase separation process at variousstate points.
A. Biphasic morphologies at long times
We start by discussing the morphologies observed whenwe quench the system to various state points in the co-existence region of the phase diagram in Fig. 2. Fig-ures 3 and 4 show that for a broad range of densi-ties, 0 . ≤ ρ ≤ .
0, and for low enough temperatures, T ≤ .
1, the morphologies obtained at long times af-ter the quench to the final temperature, t = 10 , arebicontinuous gel-like structures. These particle configu-rations are strongly reminiscent of the nonequilibriumcolloidal gels observed for instance using confocal mi-croscopy [26, 31]. By changing the density and the tem-perature, the typical length scales characterizing these FIG. 3: Snapshots of representative bicontinuous gel-likeconfigurations obtained at long time, t = 10 , for T = 0 . ρ = 0 . ρ = 0 . porous structures change and a central goal of the presentpaper is to quantify these changes.As can be noticed from the snapshots shown in Fig. 3,it is not trivial to visualize the complex morphologiesof porous, bicontinuous materials in three dimensions.Therefore, to ease the visualization, we have implementeda numerical method to localize the two phases, and no-tably the interfaces between them. The procedure relieson the definition of a coarse-grained density field, as ex-plained above in Sec. II B. While we primarily developedthis procedure to quantitatively characterize the numeri-cally obtained bicontinuous structures (described below),we find that it is also useful for visualization purposes, asdemonstrated in Fig. 4 where now the entire simulationbox is shown and the geometries are more easily visual-ized than by showing particle configurations directly.In the left column of the figure we show representa-tive configurations obtained at long times for a low tem-perature, T = 0 . T c ≈ .
2) and various densities. For low density, ρ = 0 .
15, one obtains disconnected droplets of densefluids that slowly coarsen with time. For ρ ≥ .
2, abicontinuous structure is obtained, with a dense phasewhich occupies an increasing fraction of the total vol-ume as density increases from ρ = 0 .
2, to ρ = 0 . ρ = 0 .
8. For ρ ≥ . t = 10 for quenches at fixed intermediate density, ρ = 0 .
6, and different temperatures. While the phase sep-aration proceeds rapidly for shallow quenches, T ≥ . T = 0 .
4, seeFig. 4. (Note that the fact that in this panel the finalconfiguration is not the expected spherical object mightalso be a finite size effect.) Decreasing further the tem-perature, T ≤ .
2, one observes that even at the endof our simulations an intricate bicontinuous structure re-mains apparent, indicating that phase separation is far
FIG. 4: Representative configurations obtained at longtimes, t = 10 , and various state points indicated in the fig-ure. Left: constant low temperature, T = 0 .
1, and variousdensities. Right: constant intermediate density, ρ = 0 .
6, andvarious temperatures. For all snapshots, periodic boundaryconditions are used, as is obvious for instance in the bottomright configuration. from being reached. Even smaller domains are obtainedat large times for very low temperatures, T = 0 . T the phase separation kineticsis slowed down dramatically and thus the typical domainsize remains relatively modest even at very long times. Inthe subsequent sections of the paper we will characterize these observations in a quantitative manner. B. Kinetics of phase separation
We now make a qualitative description of the kineticsof the phase separation process at various state points,as illustrated in the time series shown in Fig. 5. Eachsnapshot in a given row is separated from the previousone by a factor of 10 in time.The first column in Fig. 5 reproduces the typical timeevolution observed in molecular dynamics studies of theliquid-gas spinodal decomposition [6, 10, 48]. For thisdensity, ρ = 0 .
4, the gas and liquid phases occupy sim-ilar volumes. The temperature is T = 0 .
5, about halfthe critical temperature T c , and thus the quench is notvery deep. Using the phase diagram in Fig. 2, we seethat the liquid phase is still well above the glass transi-tion line. It thus behaves as a simple liquid, and so thespinodal decomposition takes place with no interferencefrom the physics related to the glass transition. Indeed,one observes that a bicontinuous structure forms over avery short time, t ≈
10, and then coarsens slowly withtime. For a finite system as the one used in our numer-ical simulations, the phase separation proceeds until asimple geometry is reached with a flat interface separat-ing the two phases. For an infinite system, of course,spinodal decomposition would proceed indefinitely in ascale-invariant manner [3].When temperature is decreased to T = 0 . T = 0 . t with thediffusion constant. We shall see later that such a simplerescaling is insufficient.At even lower temperature, T = 0 . t ≤ , proceedsas before, the slowing down is now dramatic, as demon-strated by the fact that the two snapshots for t = 10 and t = 10 are virtually identical (at least to the eye).This implies that at this temperature phase separation is FIG. 5: Time evolution of the phase separating system after a temperature quench at various state points indicated in eachcolumn. The selected timescales are identical in the four column. From top to bottom, t = 1 .
45, 13.2, 120, 1096, and 10 . nearly arrested at intermediate times, and well before thetypical domain size has reached the system size. Finallythe fourth column in Fig. 5 illustrates that a similar slow-ing down of the phase separation is observed for a broadrange of densities, the example shown being ρ = 0 . T = 0 . IV. STRUCTURAL ANALYSIS
To quantify the above qualitative observations we mustfirst characterizate the observed bicontinuous structuresin terms of quantitative observables and then determinehow these depend on time and temperature. In this sec-tion we introduce and compare several structural indi-cators, and show that the so-called “chord length distri-bution” represents an efficient choice to describe phaseseparating systems in our particle-based numerical sim-ulations.
A. Pair correlation function
Since we know at each timestep of the simulation theposition of all the particles, an obvious choice of a mi-croscopic function to quantity the large scale structuresis to record the pair correlation function, defined as [11] g ( r , t ) = 1 ρN N (cid:88) j =1 N (cid:88) k =1 (cid:104) δ ( r − ( r j ( t ) − r k ( t ))) (cid:105) , (4)where the brackets represent an average over indepen-dent initial conditions and trajectories, and r j ( t ) is theposition of particle j at time t . Since our configura-tions are isotropic we further perform a spherical aver-age and divide by the phase space factor 4 πr to obtain g ( r = | r | , t ). We have considered also partial pair corre-lation functions, specializing the sums in Eq. (4) to eitherone of the two species of the binary Lennard-Jones mix-ture. However, since we are interested in the large-scalestructure of the configurations, we shall only discuss thetotal pair correlation described by Eq. (4). Note that the r g (r ,t ) r g (r ,t ) t=10 t=10 t=10 t=10 t=10 -1 t=10 t=10 ρ =0.6T=0.1 FIG. 6: Time evolution of the pair correlation function g ( r, t )for a quench at ρ = 0 . T = 0 .
1, and various logarithmi-cally separated times. The inset shows a zoom of the smallamplitude oscillation around g ( r, t ) ≈ pair correlation function is frequently measured in col-loidal experiments using confocal microscopy techniquesand hence is a quantity that is experimentally accessible.In Fig. 6 we show examples of the time evolution of g ( r, t ) for a quench to ρ = 0 . T = 0 .
1. By definition, g ( r, t ) is proportional to the probability to find a parti-cle at distance r from a particle located at the origin.Therefore g ( r, t ) describes for short distances the localamorphous structure of the dense phase. Accordingly,it shows a pronounced first peak corresponding to inter-particle distances (occurring at r ≈ . g ( r, t )rapidly converges to unity if r increases, for the hetero-geneous phase separating systems it shows oscillationsaround unity even at large distances, as shown in Fig. 6.If t is large, the first oscillation below 1 represents theaverage distance between a particle taken at random in adense domain to a neighboring gas region. This physicalinterpretation suggests that a possible quantitative defi-nition of the average domain size, L ( t ), can be obtainedfrom g ( L ( t ) , t ) = 1. The data shown in Fig. 6, however,indicate that the oscillations of g ( r, t ) around unity atlarge r have a rather small amplitude, in particular if t islarge (see Inset in Fig. 6). Our simulations have shownthat the above definition is physically sensible, in that itcoincides well with the typical domain size seen in thesnapshots [56]. However, we also noticed that this mea-surement is prone to very large statistical fluctuations,since the amplitude of the oscillations in g ( r, t ) is verysmall. As a result, the use of this method to follow the q -1 S ( q ,t ) t=10 t=10 t=10 t=13.2t=4.37t=1.45t=39.8 ρ =0.6T=0.1 FIG. 7: Time evolution of the static structure factor, Eq. (5),for a quench at ρ = 0 . T = 0 .
1, and various logarithmi-cally separated times. The lines are fits using Eq. (6). time evolution of the phase separating systems does re-quire a very large numerical effort since highly accuratepair correlation functions must be measured.Additional visual inspections indicate that the largedistance behavior of g ( r, t ) is in fact strongly influencedby a small number of very large domains found in the sys-tem, whose statistical properties strongly fluctuate fromone run to another. Therefore, although this two-pointfunction is a simple structural correlation which corre-sponds also to the quantity usually analyzed in theo-retical calculations, our work suggests that, at least forparticle-based numerical simulations, it does not repre-sent the most practical choice to determine the averagedomain size. B. Structure factor
Since the static structure factor is the Fourier trans-form of the pair correlation function [11], it carries apriori the same physical content. It is, however, moreeasily accessible to experiments using for instance lightscattering techniques. It is defined as S ( q , t ) = 1 N N (cid:88) j =1 N (cid:88) k =1 (cid:104) exp[ i q · ( r j ( t ) − r k ( t ))] (cid:105) . (5)As for the pair correlation function, we perform a spher-ical average to obtain S ( q, t ) = S ( | q | , t ). In Fig. 7, weshow the results for the structure factor for the sameset of parameters as for the pair correlation function pre-sented in Fig. 6. An advantage of S ( q, t ) over the pair cor-relation function is that the local structure of the densephase at short-distance and the inhomogeneous bicontin-uous structure present at large scales show up in S ( q, t )at very different wavevectors and thus can easily be stud-ied independently: While the local structure appears asa sharp peak near q ≈ π/σ ≈
6, the large domains at larger scale produce a signal at much lower wavevectors,and it is this low- q range which is shown in Fig. 7.As reported in previous work on spinodal decomposi-tion in fluids [6, 48], the static structure factor develops atlow wavevector a peak whose maximum intensity, S (cid:63) ( t ),grows and whose peak position, q (cid:63) ( t ), moves to lower q as time increases. The peak position directly revealsthe typical domain size in the phase separating system, L ( t ) ≈ π/q (cid:63) ( t ), and the q -dependence near the peak canbe adjusted using the following empirical formula [59]: S ( q, t ) = S (cid:63) q/q (cid:63) ) q/q (cid:63) ) , , (6)where the numerical factors are chosen such that S ( q = q (cid:63) , t ) = S (cid:63) . This formula interpolates in a simple mannerbetween the expected quadratic behavior at low q , S ( q (cid:28) q (cid:63) , t ) ∝ q , and Porod’s law describing the structure ofthe interfaces at larger q , namely S ( q (cid:29) q (cid:63) , t ) ∝ q − , ina three dimensional space. We have used Eq. (6) to fitthe data shown in Fig. 7, which gives us direct access tothe growing length scale L ( t ) characterizing the spinodaldecomposition.We have found that this approach is much more re-liable to obtain a quantitative estimate of the averagedomain size than using the function g ( r, t ), since thelarge-scale signal is better resolved in Fourier than inreal space. Also, we have noticed that fluctuations areless pronounced in S ( q, t ) than in g ( r, t ), which impliesa reduced numerical effort. The drawback of this sim-ple measurement of the domain size is however readilyobserved in Fig. 7. Since with increasing time the peakposition shifts to lower q , for large t the peak appears atthe border of the accessible wavevector range, which isbounded at low q by the system size, i.e. q ≥ πL . This isparadoxical at first sight, because the snapshots shownin Fig. 5 seem to indicate that even at large times theaverage domain size remains quite a bit smaller than thebox size.The reason for this is that the behavior of the structurefactor is, just as for the pair correlation function, stronglydominated by the largest domains in the system. There-fore, despite the several advantages mentioned above forthe structure factor, it suffers from the practical draw-back that accurate measurements of this two-point corre-lation function require system sizes that are considerablylarger than the typical domain size. To circumvent thisdifficulty, we have turned to a slightly more complicatedobservable, as we now describe. C. The chord length distribution
The definition of the coarse-grained density field, givenby Eq. (3), allows to locate the position of the interfacesseparating the two phases in the course of the phase sep-aration process. This information can then be used to0
FIG. 8: Chords are defined by the intersection of straightlines with the interfaces in the phase separating system. Inthis bidimensional example the red segments belong to thedense phase (yellow areas) and their length give the chordlength of the dense phase. measure the distribution of the domain size in the bicon-tinuous structures.To this end, we measure the so-called chord length dis-tribution (see Fig. 8) [60, 61]. We define chords by twoconsecutive intersections of a straight line with the in-terfaces present in the system. In practice, we measurechords along the three axis of the lattice used to coarse-grain the density, and measure the length (cid:96) of the seg-ments belonging either to the gas or to the dense phase.For this one has of course to take into account the peri-odic boundary conditions.By repeating this measurement over the entire latticeused to determine the coarse-grained density, we obtainthe distribution of chord lengths P ( (cid:96) ), either for chordsin the gas phase, or for chords in the dense phase. Rep-resentative results for the time evolution of these dis-tributions after a quench to ρ = 0 . T = 0 . (cid:96) = L/ L ( t ) = (cid:90) ∞ d(cid:96) P ( (cid:96), t ) (cid:96). (7)In addition, although equivalent at long times, we findthat the distribution of chord length in the gas phaseyields more accurate results at short times than the oneof the dense phase, see Fig. 9. Therefore, in the restof this paper we shall use Eq. (7) for the chord lengthdistribution in the gas phase as our quantitative deter-mination of the average domain size characterizing our l -6 -5 -4 -3 -2 -1 P ( l ,t ) ρ =0.6T=0.1 t = . t = . t = . t = . t=120.2 t=10 t=10 gas (a) l -6 -5 -4 -3 -2 -1 P ( l ,t ) t=10 t=0.83 t=7.59 t=13.2 t=10 t=39.8 t=10 ρ =0.6T=0.1liquid (b) FIG. 9: Time evolution of the chord length distributionsmeasured in the gas phase (a) and in the dense phase (b)after a quench to ρ = 0 . T = 0 . phase separating structures. V. TEMPORAL EVOLUTION
In this subsection, we analyze the temporal evolutionof the phase separation process, and study how the ki-netics depends on the state point chosen for the quench,focusing in particular on the influence of temperature.
A. Finite size effects
One of the central question we wish to answer iswhether or not the phase separation kinetics is arrestedat sufficiently low temperatures. In the previous sectionswe have shown that decreasing T does indeed lead to astrong slowing down of the relaxation dynamics. Beforeone starts to characterize this slowing down in a morequantitative manner it is, however, important to recallthat this relaxation dynamics does depend to some ex-tent on the size of the system and hence one has to checkthe influence of these finite size effects. In particular, wefind that coarsening stops earlier if the systems size issmall. The reason for this is that the interfaces are frus-trated by the periodic boundary conditions which con-1 t L ( t ) N=49KL
S(q) (t) N=300KN=1M ρ =0.6T=0.1 N=8K
FIG. 10: Influence of system size on the growth of the averagedomain size for a quench to ρ = 0 . T = 0 .
1, and systemsize between N = 8 · and N = 10 . The data for L S ( q ) ( t ) =2 π/q (cid:63) are deduced from the analysis of the structure factorshown in Fig. 7. Note that this data has been multiplied by0.38 in order to match the length scale L ( t ) at intermediatetimes. strain and slow down their motion. This remark is alsoexperimentally relevant for studies of phase separation inconfined geometries [62]. Therefore, before concluding onthe possibility of kinetically arrested phase separations, itis important to make sure that our results do not dependcrucially on the chosen system size.To this end, we have performed a systematic study ofthe influence of a finite system size on the growth of theaverage domain size. Some of these results are presentedin Fig. 10 for a quench to ρ = 0 . T = 0 .
1. Theseresults confirm that the average domain length reachedat a given time after the quench increases with increasingthe system size. However, we find that this effect doesnot influence the results in a strong manner. For the par-ticular case shown in Fig. 10, the domain size increasesby about 10 % when N increases by more than 2 ordersof magnitude. Furthermore we find no N dependencewithin the error bars for the final sizes shown in Fig. 10.Therefore we have decided to perform most of our studiesusing N = 3 · , as a compromise between a very largesystem, and a broad time window in which the dynamicscan be probed.Another finding documented in Fig. 10 is that thegrowth of the length scale extracted from the chordlength distribution, Eq. (7), matches the one obtainedfrom the length scale extracted from the dominantwavevector q (cid:63) in the static structure factor, Eq. (6). Inthe graph we have multiplied the latter by a constantfactor 0.38 and find that at short and intermediate timesthe two curves do indeed track each other. Thus one canconclude that the chord length distribution represents anefficient and accurate way of extracting the average do-main size in phase separating systems. t L ( t ) T=0.015T=0.05T=0.1 T=0.2T=0.3T=0.4 t T=0.6T=0.5 t t ρ =0.6 (a) t L ( t ) ρ =0.3 ρ =0.2 ρ =0.4 ρ =0.6 ρ =0.8 (b) FIG. 11: a) Influence of temperature on the domain growthfor ρ = 0 .
6. Various power-laws are included as well. b)Influence of the density on the domain growth for T = 0 . B. Growth of domain size
We now study how the density and temperature forthe quench influence the kinetics of the phase separa-tion process. Our numerical results are summarized inFig. 11.We first discuss the influence of temperature for a givendensity, ρ = 0 .
6, see Fig. 11a. For a relatively shallowquench, T = 0 .
5, the data shows an upward curvature inthis double logarithmic representation before it saturatesat long times at a system size dependent value, indicat-ing complete phase separation. For relatively early times,the domain growth is approximated well by L ( t ) ∼ t / ,crossing over to a faster growth at larger times beforethe saturation. Such an apparent square root time de-pendence has been found in earlier work on the liquid-gasphase separation [6, 48, 49]. Its physical interpretation isthat it represents an effective power-law growth [7] inter-polating between two regimes, L ( t ) ∼ t / at short timesfollowed by L ( t ) ∼ t at longer times, which are theoret-ically expected power-laws controlling spinodal decom-position at early and late times [5]. These regimes haveindeed been separately observed in specifically dedicatedsimulations [9]. Whereas the first regime corresponds to asurface-tension driven domain growth limited by particlediffusion, the second one is observed when hydrodynam-2ics becomes relevant. We note that our data in Fig. 11are consistent with this scenario, but do not exhibit con-vincing power-law regimes over broad and distinct timewindows.The situation for deeper quenches at lower tempera-ture is more unusual. For T ≤ . t − dependencecrosses over to one that is significantly slower, one can-not argue that the exponent 0.5 is related to a cross-overbehavior to the hydrodynamic regime. In fact, due tothe very large viscosity of the liquid at low T it must beexpected that hydrodynamics ceases to play a role [13].What is a somewhat surprising is that at long times thegrowth is even slower than the usual t / law. This in-dicates that at low temperatures surface tension is nolonger the main mechanism that drives the coarseningprocess when the domain size becomes large and that in-stead a different coarsening regime sets in. This resultis in agreement with the snapshots presented in Fig. 3that show that at low temperatures the interface can berather rough (see also Fig. 13).When temperature becomes very small, T ≤ .
2, thedata in Fig. 11 indicates that the domain growth at longtimes is not well described by a power-law dependence,as the curves appear to be bent in this log-log repre-sentation. This indicates that at these low temperaturesthe growth is logarithmic, a functional form that is foundquite often in glassy systems. For very low temperatures, T = 0 .
015 (which is about 1 % of the critical tempera-ture), the domain size ceases to grow at long times andbecomes nearly constant within our statistical accuracy.The observation of very slow domain growth at lowtemperature is quite generic, as demonstrated in Fig. 11b,where the density is varied for a constant low tempera-ture T = 0 .
1. The data for ρ = 0 . ρ = 0 . -1 t S u r f a c e ( t ) T=0.015T=0.03 T=0.05T=0.1T=0.2T=0.3T=0.4 ρ =0.6 (a) t -7.0-6.9-6.8-6.7-6.6-6.5 E p ( t ) T=0.015T=0.03T=0.05 T=0.1 T = . T = . T = . ρ =0.6 (b) FIG. 12: Time evolution of the surface of the interface (a)and potential energy (b) during the phase separation at ρ =0 . that the reported observation of gels formed by “kinet-ically arrested” spinodal decomposition in colloidal sys-tems [26, 28] is likely only an approximation (albeit aphysically relevant one) resulting from the combinationof deep quenches and short observation timescales. C. Energy density and area of interfaces
Before discussing the details of the microscopic dynam-ics at low temperature, we present two additional observ-ables that are useful to characterize the dynamics of thesystem: The energy density and area of interfaces.The time evolution of these two quantities is shownin Fig. 12 for the density ρ = 0 . t − dependence of these two quan-tities is of course very different: The energy density de-3creases rapidly since it is dominated by the bulk behav-ior of the fluid which must thermalize at low temper-ature after the sudden quench from the high tempera-ture. By contrast, because the system must first createa large quantity of interfaces right after the quench froma homogeneous initial state, the amount of interfaces isnon-monotonous in that it first increases rapidly at earlytimes, before the coarsening process starts and the sur-face decreases again, see Fig. 12a.For times t >
10 the evolution of these quantities ismore similar, and follows from similar physical consider-ations. At large times, the area of the interface decreasesslowly as a result of the coarsening process which elimi-nates small domains and generates larger ones. The en-ergy density has a more complicated behavior because itreceives contributions from both the energetically costlyinterfaces as well as from the bulk of the dense domains.Since both contributions decrease slowly with time, theenergy density also decreases with increasing t . If oneassumes that interfaces dominate the time dependenceat long times, then the energy density should display atime dependence that is very similar to the one of thefraction of interfaces, as confirmed by the data in Fig. 12at high T . However, since we have seen that at low tem-peratures the size of the domains is not governed by thesurface tension (cf. discussion on Fig. 11) one can expectthat at low T the time dependence of the energy and ofthe surface are not the same.The time dependence of the surface as well as of theenergy depends strongly on temperature in that the re-laxation becomes slower, in agreement with our find-ings regarding the size of the domains. Figure 12ashows that the amount of surface at a given (large)time increases monotonically with decreasing tempera-ture, showing that at low T the system has more smallerdomains and their surface is rougher. At the lowest tem-peratures the time dependence of the surface is compat-ible with a logarithmic decay, in agreement with our re-sults on the growth of the domain sizes.The increasing fraction of small domains at low tem-peratures has an implication on the value of the energy ata given (large) time: Since the energy of the dense phasedecreases with T , one finds that for intermediate tem-peratures the energy at large t decreases if T is lowered.However, at even lower T , the system starts to have somany small domains that are rough and that cost energy,that the overall energy starts to increase again, leadingto a non-monotonic T − dependence of that quantity (seeFig. 12b). VI. INTERMITTENT DYNAMICS AT LOWTEMPERATURES
In this section we provide evidence that the qualitativechange in the growth law at low temperatures is also ac-companied by a qualitative evolution of the microscopicmechanisms driving the coarsening dynamics.
A. How domains coarsen
If temperature is not very low, the microscopic dy-namics of coarsening is well understood. At early timesof the spinodal decomposition, a bicontinuous structureemerges rapidly, which is characterized by a well-definedlength scale that reflects the intrinsic instability of the ho-mogeneous system after the quench into the coexistenceregion. Moreover this bicontinuous structure is charac-terized by curved interfaces that store a large amountof potential energy. In this case surface tension is themain driving force for the phase separation process thatfollows the spinodal instability, and domains coarsen inorder to reduce the curvature of the interfaces and theirtotal area. At the microscopic scale, this process can pro-ceed because particles can easily move within the densephase in response to this driving force.At low temperature, we observe that surface tensionbecomes unable to advance the coarsening process, be-cause the dense phase is now an (aging) glassy mate-rial that has a very high viscosity and is in fact visco-elastic [22, 63]. As a consequence, surface tension is nolonger able to relax in a significant manner the curvedinterfaces formed during the phase separation process.In a recent experimental work on attractive col-loids [33], it has been found that particles located nearthe interface of the bicontinuous structure have a mobil-ity which is larger than the one of particles in the bulkof the domains. Therefore it has been concluded thatsurface-enhanced mobility provides an important contri-bution to the dynamics for deep quenches [33]. We haveinvestigated whether this effect is also relevant for oursimulations. First of all we point out that, due to en-ergetic considerations, it is more favorable for the sys-tem to put the A particles at the interface (since theA-B interaction is stronger than the A-A interaction).Thus such a micro-segregation would a priori give riseto an enhanced mobility of the particles at the interface.However, despite this we have not found that for deepquenches particles at the interfaces are significantly moremobile than the ones inside the dense phase. This resultis in fact consistent with earlier studies of the presentbinary Lennard-Jones in inhomogeneous geometries [64].While particles near surfaces are in general indeed moremobile than particles in the bulk at equivalent thermo-dynamic conditions, one should also notice that this dy-namical difference is usually only relevant in the narrowrange of temperatures which corresponds to the inter-val between bulk and surface glass temperatures, wherebulk diffusion is already arrested while surface diffusionis not. The quench depth should therefore be specificallyadjusted to have conditions for which surface-enhanceddiffusion is as effective as it seems to be the case in theexperiments of Ref. [33].Although we do find in our low-temperature simula-tions that slow coarsening persists, this domain growthis driven neither by surface tension nor by interface-enhanced particle diffusion. Instead, the complex bi-4
FIG. 13: Time series showing the breaking of two thin do-mains in a quench at ρ = 0 . T = 0 .
1. Only a small fraction(16 %) of the entire system is shown for clarity, and about3000 particles with the largest mobility are highlighted, theyare located near the breaking points. The first breaking oc-curs between t = 50 and 1100 and the second one between t = 1300 and 2100. continuous structure formed at long times continues toevolve by intermittent breaking of thin necks, which inturn allows the structure to relax further [38, 65]. Weillustrate this process in the time series of Fig. 13, wherethe most mobile particles over the considered time win-dow are highlighted. Visual inspection shows that forthis low temperature, T = 0 .
1, particles located at theinterface of a dense domain are nearly arrested. A sec-ond observation is that the interfaces are relatively rough,which confirms that surface tension is no longer efficientat low temperatures, and that it is unable to relax inter-faces that are very curved. A third observation is thatthe time evolution of the overall structure occurs whenthin domains suddenly break, which occurs twice in thisspecific time series, first near t ≈ t ≈ B. Intermittent dynamics in space and time
In this subsection we provide further evidence for aqualitative change in the microscopic dynamics betweenshallow and deep quenches. In the previous section wehave argued that at low T coarsening proceeds by in-termittent domain breaking. However, even for shallowquenches domains grow and thin domains can break. Themain difference between the two situations is that at lowtemperatures, mobility is highly localized in space andtime, and domain breaking appears very rarely, whereasat high T no such localization is observed.This difference is illustrated in Figs. 14 and 15, wherethe 1 % most mobile particles over different time inter-vals are marked by spheres, and their displacements arerepresented by a vector. The rest of the particles are5 FIG. 14: Most mobile particles (1 % of all particles, shownas spheres) and their displacements (shown as cylinder) for agiven time interval, t ∈ [38 . , .
5] after a quench to ρ = 0 . T = 0 . represented by small dots. The time intervals are chosenso that the images are taken for comparable evolutionsof the average domain size. As a consequence, the timesare much shorter at the high temperature (Fig. 14) thanat the lower one (Fig. 15), and time windows becomebroader with increasing time in Fig. 15.At high temperature ( T = 0 .
5, Fig. 14) the most mo-bile particles appear anywhere in the system (both atinterfaces and in the bulk domains where mobility ishigh) and their displacements are essentially uncorre-lated. This picture is representative of shallow quenchesand indicates that dynamic heterogeneity at the particlescale is not very relevant for ordinary phase separations.This behavior is in strong contrast with the one foundat low temperature, T = 0 . FIG. 15: Most mobile particles (1 % of all particles, shownas spheres) and their displacements (shown as cylinder) forthree given time intervals, t ∈ [209 , t ∈ [1905 , t ∈ [5754 , ρ = 0 . T = 0 .
1. The small dots are the remaining 99% of theparticles.
VII. RELATION BETWEEN BINODAL ANDGLASS LINES
In this last section, we provide more details and dis-cussion about the determination and behavior of the co-existence line at low temperatures in the phase diagramof Fig. 2.Why is this an issue? In the phase diagram of Fig. 2,the binodal line crosses the glass transition line near thepoint ρ ≈ .
15 and T ≈ .
35. This implies that fortemperatures lower than T = 0 .
35, we cannot determinethe coexistence line using standard equilibrium simula-tions, and we have to rely on nonequilibrium protocolsto extend the binodal down to T = 0. Therefore the low-temperature extension of the coexistence line changes na-ture at low temperature since it is not uniquely definedanymore but is instead dependent on the protocol.Because of this difficulty, two experimental groups haveinvestigated this issue in more detail [26, 28]. In bothcases the measurements proceed as follows: The systemis first quenched into the coexistence region where phaseseparation takes place, and becomes nearly arrested onexperimental timescales. Then the volume occupied bythe dense phase is determined experimentally, from whichthe density is determined. Different methods have beenused to measure the evolution of this density as a functionof quench depth, and two qualitatively distinct resultswere found, as mentioned above in Sec. II C.For the case of our numerical simulations we have fol-lowed a similar approach to determine the coexistenceline. After quenching the system to a given state pointwe determined the density of each phase using the coarse-grained density field described in Sec. II B. Following thisprocedure, we obtained for each configuration both thevolume occupied by the dense phase as well as the num-ber of particles it contains, from which the density iseasily deduced. Note that for a given temperature, thedensity of the dense phase will depend in principle onboth the time spent since the quench, and on the to-tal density of the system at which the quench has beenperformed.In Fig. 16a we demonstrate that after the quench themeasured density converges very rapidly to its asymp-totic value. This is an important result since it showsthat the heterogeneous nature of the configurations doesnot preclude a quantitative determination of the densityof the glass phase. In other words, the time dependenceof the glass density is not an issue. Therefore, for a givendensity, here ρ = 0 .
6, we can vary the quench depthand obtain the temperature evolution of the glass den-sity, which we can include in the ρ − T phase diagram, -1 t ρ g l a ss ( t ) T=0.1 T=0.2 T=0.3T=0.4 T=0.5 ρ =0.6 (a) ρ T gas+liquid liquid glassgas+glass b i noda l g l a s s t r a n s i t i o n li n e ρ =0.8 ρ =0.6 ρ =0.4 ρ =0.2 Vollmayr et al. bulk density ( ρ =0.6) (b) FIG. 16: a) Density of the dense phase during phase separa-tion at ρ = 0 . see Fig. 16b. In this representation, we focus directlyon the relevant temperature region below the intersec-tion with the glass transition line. We observe that thecoexistence line determined using quenches at ρ = 0 . ρ = 0 . ρ = 0 .
8. We find that all densities produce very sim-ilar coexistence lines, with differences in density of about2 % between the two extremes, the larger quench densityproducing a smaller glass density.Finally we compare the results for the coexistence linewith a very different numerical approach. In Ref. [71], the7present Lennard-Jones binary mixture was used to studythe influence of cooling rates on the structure of the glass.These quenches were done at constant pressure, whichwas zero. Thus, the produced homogeneous glass config-urations were adjusting their densities to maintain a zeropressure, and these densities were recorded numerically.We have included the temperature evolution of these den-sities in Fig. 16b as well and we see that they comparevery well with our determination of the coexistence re-gion. This is expected because the dense phase in ourphase separating systems coexists with a gas phase withvanishing pressure. This comparison seems to confirmour finding that the coexistence line does not exhibit areentrant behavior as a result of the crossing of the glasstransition line.
VIII. SUMMARY AND CONCLUSION
We have used large-scale molecular dynamics simula-tions to study the influence of a temperature quench onthe liquid-gas phase separation kinetics in a Lennard-Jones fluid, and therefore the competition between thephase separation kinetics and the glass transition occur-ring at low temperature in bulk liquids. This representstherefore an example of a viscoelastic phase separation.Our main finding is the observation that the phaseseparation kinetics changes qualitatively with decreas-ing temperature: The microscopic dynamics evolves froma diffusive motion driven by surface tension for shallowquenches, to a qualitatively different coarsening regime in which the dynamics becomes strongly intermittent,spatially heterogeneous and thermally activated at lowtemperature, leading to logarithmically slow growth ofthe typical domain size.The microscopic description of the coarsening pro-cess occurring in our simulations at low temperatures,which results from the intermittent release of mechan-ical constraints, is strongly reminiscent of the physicalscenario put forward to explain experimental and simu-lation results obtained in a broad variety of soft glassymaterials for which unusual aging dynamics has been re-ported [24, 68, 69, 72]. In future work, it would be in-teresting to compare time correlation functions measuredin numerical simulations such as ours to the outcome oflight scattering experiments performed in soft glassy ma-terials in their aging regime. Such studies would allow toobtain a better understanding how these soft glass ma-terials are related to the gel-like structures investigatedhere.
Acknowledgments
We thank D. Reichman for initially stimulating thiswork, and B. Coasne, T. Gibaud, and P. Royall for use-ful discussions. The research leading to these results hasreceived funding from the European Research Councilunder the European Union’s Seventh Framework Pro-gramme (FP7/2007-2013) / ERC Grant agreement No306845. W. Kob acknowledges support from the InstitutUniversitaire de France. [1] R. C. Desai and R. Kapral,
Dynamics of self-organizedand self-assembled structures (Cambridge UniversityPress, Cambridge, 2009).[2] M. C. Cross and P. C. Hohenberg, Rev. Mod. Phys. ,851 (1993).[3] A. J. Bray, Adv. Phys. , 357 (1994).[4] P. C. Hohenberg and B. I. Halperin, Rev. Mod. Phys ,435 (1977).[5] H. Furukawa, Adv. Phys. , 703 (1985).[6] S. W. Koch, R. C. Desai, and F. F. Abraham, Phys. Rev.A , 2152 (1983).[7] S. Bastea and J. L. Lebowitz, Phys. Rev. Lett. , 3776(1995).[8] M. Laradji, S. Toxvaerd, and O. G. Mouritsen, Phys.Rev. Lett. , 2253 (1996).[9] V. M. Kendon, M. E. Cates, I. Pagonabarraga, J.-C. De-splat, and P. Bladon, J. Fluid Mech. , 147 (2001).[10] S. Ahmad, S. K. Das and S. Puri, Phys. Rev. E ,031140 (2012).[11] J. P. Hansen and I. R. McDonald, Theory of Simple Liq-uids (Elsevier, Amsterdam, 1986).[12] T. A. Witten,
Structured fluids (Oxford University Press,Oxford, 2004).[13] H. Tanaka, J. Phys.: Condens. Matter , R207 (2000).[14] H. Tanaka and Y. Nishikawa, Phys. Rev. Lett. , 078103 (2005).[15] V. Testard, L. Berthier, and W. Kob, Phys. Rev. Lett. , 125702 (2011).[16] D. Sappelt and J. J¨ackle, Europhys. Lett. , 13 (1997).[17] S. K. Danchinov, Y. D. Shibanov, and Y. K. Godovsky,Colloid Polym. Sci. , 234 (1999).[18] F. Varrato, L. Di Michele, M. Belushkin, N. Dorsaz, S.H. Nathan, E. Eiser, and G. Foffi, Proc. Nat. Ac. Science , 19155 (2012).[19] K. Binder and W. Kob, Glassy materials and disorderedsolids (World Scientific, Singapore, 2011).[20] L. Berthier and G. Biroli, Rev. Mod. Phys. , 587(2011).[21] J.-P. Bouchaud, L.F. Cugliandolo, J. Kurchan, and M.M´ezard, Physica A , 243 (1996).[22] W. Kob and J.-L. Barrat, Phys. Rev. Lett. , 4581(1997).[23] L. Buisson, L. Bellon, and S. Ciliberto, J. Phys.: Con-dens. Mat. , S1163 (2003).[24] L. Cipelletti and L. Ramos, J. Phys.: Condens. Matter , R253 (2005).[25] S. Manley, H. M. Wyss, K. Miyazaki, J. C. Conrad, V.Trappe, L. J. Kaufman, D. R. Reichman, and D. A.Weitz, Phys. Rev. Lett. , 238302 (2005).[26] P. J. Lu, E. Zaccarelli, F. Ciulla, A. B. Schofield, F. Sciortino, and D. A. Weitz, Nature , 499 (2008).[27] E. Zaccarelli, P. J. Lu, F. Ciulla, D. A. Weitz, and F.Sciortino, J.Phys.: Condens. Matt. , 494242 (2008).[28] F. Cardinaux, T. Gibaud, A. Stradner, and P. Schurten-berger, Phys. Rev. Lett. , 118301 (2007).[29] T. Gibaud and P. Schurtenberger, J. Phys.: Condens.Matt. , 32220 (2009).[30] T. Gibaud, F. Cardinaux, J. Bergenholtz, A. Stradner,and P. Schurtenberger, Soft Matter , 857 (2011).[31] A. I. Campbell, V. J. Anderson, J. S. van Duijneveldt,and P. Bartlett, Phys. Rev. Lett. , 208301 (2005).[32] L. J. Teece, M. A. Faers, and P. Bartlett, Soft Matter ,1341 (2011).[33] I. Zhang, C. P. Royall, M. A. Faers, and P. Bartlett, Softmatter , 2076 (2013)[34] G. Foffi, C. de Michele, F. Sciortino, and P. Tartaglia, J.Chem. Phys. , 224903 (2003).[35] P. Charbonneau and D. Reichman, Phys. Rev. E ,050401(R) (2007).[36] P. Charbonneau and D. Reichman, Phys. Rev. E ,011507 (2007).[37] E. Zaccarelli, J. Phys.: Condens. Matter , 323101(2007).[38] H. Tanaka and T. Araki, EPL , 58003 (2007).[39] F. Sciortino and E. Zaccarelli, Curr. Op. Sol. State Mat.Scien. , 246 (2011).[40] C. P. Royall and S. R. Williams, J. Phys. Chem. B ,7288 (2011).[41] E. R. Dufresne, H. Noh, V. Saranathan, S. G. J. Mochrie,H. Cao, and R. O. Prum, Soft Matter , 1792 (2009).[42] J. D. Forster, H. Noh, S. F. Liew, V. Saranathan, C. F.Schreck, L. Yang, J-G Park, R. O. Prum, C. S. O’Hern, S.G. J. Mochrie, H.Cao, and E. R. Dufresne, Adv. Mater. , 2939 (2010).[43] K. van Gruijthuijsen, V. Herle, R. Tuinier, P. Schurten-berger, and A. Stradner, Soft Matter , 1547 (2012).[44] W. Ma, C. Yang, X. Gong, K. Lee, and A. J. Heeger,Adv. Funct. Mater. , 1617 (2005).[45] J. Peet, A. J. Heeger, and G. C. Bazan, Acc. Chem. Res. , 1700 (2009).[46] J. S. Moon, J. K. Lee, S. Cho, J. Byun, and A. J. Heeger,Nano. Lett. , 230 (2009).[47] A. Onuki and S. Puri, Phys. Rev. E , 1331 (1999).[48] R. Yamamoto and K. Nakanishi, Phys. Rev. B , 14958(1994); Phys. Rev. B , 2715 (1995).[49] E. Velasco and S. Toxvaerd, Phys. Rev. Lett. , 388(1993). [50] B. D. Butler, H. J. M. Hanley, D. Hansen, and D. J.Evans, Phys. Rev. Lett. , 4468 (1995).[51] J. Sabin, A. E. Bailey, G. Espinosa, and B. J. Frisken,Phys. Rev. Lett. , 195701 (2012).[52] W. Kob and H. C. Andersen, Phys. Rev. Lett. , 1376(1994); Phys. Rev. E , 4626 (1995); Phys. Rev. E ,4134 (1995).[53] S. Plimpton, J. Comp. Phys. , 1 (1995).[54] S. Majumder and S. K. Das, EPL , 46002 (2011).[55] S. Ahmad, S. K. Das and S. Puri, Phys. Rev. E ,040107 (2010).[56] V. Testard, Etude par simulations num´eriques del’influence de la transition vitreuse sur la s´eparation dephase liquide-gaz , Th`ese de l’Universit´e Montpellier 2(2011).[57] L. Berthier and G. Tarjus, Phys. Rev. E , 031502(2010).[58] S. Sastry, Phys. Rev. Lett. , 590 (2000).[59] H. Furukawa, Physica A , 497 (1984).[60] P. Levitz, Adv. Colloid Interface Sci. , 71 (1998).[61] L. D. Gelb and K. E. Gubbins, Langmuir , 2097 (1998).[62] P. S. Sarangapani, Y. Yu, J. Zhao, and Y. Zhu, Phys.Rev. E , 061406 (2008).[63] D. El Masri, L. Berthier, and L. Cipelletti, Phys. Rev. E , 031503 (2010).[64] R. Malshe, M. D. Ediger, L. Yu, and J. J. de Pablo, J.Chem. Phys. , 194704 (2011).[65] T. Koyama, T. Araki, and H. Tanaka, Phys. Rev. Lett. , 065701 (2000).[66] M. Tsamados, A. Tanguy, C. Goldenberg, and J.-L. Bar-rat, Phys. Rev. E , 026112 (2009).[67] Slow Relaxations and Nonequilibrium Dynamics in Con-densed Matter , edited by J.-L. Barrat, J. Dalibard, M.Feigelman, and J. Kurchan (Springer, Berlin, 2003).[68] L. Cipelletti, S. Manley, R. C. Ball, and D. A. Weitz,Phys. Rev. Lett. , 2675 (2000).[69] A. Duri, D. A. Sessoms, V. Trappe, and L. Cipelletti,Phys. Rev. Lett. , 085702 (2009).[70] Dynamical heterogeneities in glasses, colloids and granu-lar materials , Eds. L. Berthier, G. Biroli, J.-P. Bouchaud,L. Cipelletti, W. van Saarloos (Oxford University Press,Oxford, 2011).[71] K. Vollmayr, W. Kob, and K. Binder, J. Chem. Phys. , 4714 (1996).[72] M.-A. Suarez, N. Kern, E. Pitard, and W. Kob, J. Chem.Phys.130