Fluctuating Phases and Fluctuating Relaxation Times in Glass Forming Liquids

The presence of fluctuating local relaxation times, $\tau_r(t)$ has been used for some time as a conceptual tool to describe dynamical heterogeneities in glass-forming systems. However, until now no general method is known to extract the full space and time dependent $\tau_r(t)$ from experimental or numerical data. Here we report on a new method for determining the local phase field, $\phi_r(t) = \int^{t} dt'/\tau_r(t')$ from snapshots $\{r(t_i)\}_{i=1...M}$ of the positions of the particles in a system, and we apply it to extract $\phi_r(t)$ and $\tau_r(t)$ from numerical simulations. By studying how the phase field depends on the number of snapshots, we find that it is a well defined quantity. By studying fluctuations of the phase field, we find that they describe heterogeneities well at long distance scales.


I. INTRODUCTION
As glass-forming liquids enter the supercooled state, relaxation processes gradually but dramatically slow down, and at the same time they become non-exponential [1][2][3]. One possible way of understanding this behavior would be to think of the relaxation function F (t) as the sum of a large number of exponentially decaying contributions, each with different relaxation times; with each contribution corresponding to a particular region in the system [2,3]. This is one aspect of the picture of dynamical heterogeneity, in which mesoscopic regions relax differently from each other and from the bulk [1][2][3][4][5][6][7][8][9][10][11][12]. In this framework, it is assumed that the differences between regions are not frozen in but dynamical, i.e. as the system evolves, the local relaxation times τ r (t) evolve, so that "fast" regions become "slow", and viceversa. This same picture has also been invoked to explain other phenomena occuring in the same regime, such as the breakdown of the Stokes-Einstein relation between the diffusion coefficient and the shear viscosity. Direct evidence for dynamical heterogeneity has been found in particle tracking experiments in glassy colloidal systems [5][6][7] and granular systems [9], and in numerical simulations [10,11].
The origin of dynamical heterogeneity is still uncertain, although some explanations have been proposed [12]. One of the proposed explanations is that dynamical heterogeneity emerges as a consequence of dynamical constraints, which produce non trivial structure in the space of (space-time) particle trajectories [13,14]. Another proposal is that as a liquid is supercooled it eventually undergoes a Random First Order Transition in which the liquid freezes into a mosaic of aperiodic crystals [15,16], which give rise to the heterogeneity that is observed. A third approach postulates that dynamical heterogeneity emerges from the presence of soft (Goldstone) modes associated with a broken continuous symmetry under time reparametrizations t → φ(t) [17][18][19][20][21][22][23][24][25][26][27].
A quantitative description of dynamical heterogeneity in terms of the presence of locally fluctuating relaxation times τ r (t) has in principle some strong advantages. One of them is its simplicity and intuitive appeal. Another one is that the basic quantities that appear in this description are intrinsically instantaneous, as opposed to other common descriptions for which the basic quantities describe changes in the system over a finite time interval, which complicates their interpretation. Despite those advantages, however, such a description has proved elusive. For example, even the question of experimentally determining the relationship between the lifetime τ ex of regions of heterogeneous dynamics and the bulk α-relaxation time τ α has proved controversial, with some results indicating that τ ex /τ α ∼ 1, and other results indicating that τ ex /τ α ≫ 1 [2].
In the present work, our goal is twofold. On one hand, we establish a connection between the presence of time reparametrization fluctuations and the presence of fluctuating local relaxation times τ r (t). On the other hand, we introduce a new method to extract the actual values of those relaxation times τ r (t) from experimental or numerical data, and test this method on data from numerical simulations.
Our paper is organized as follows. In Sec. II we discuss the connection between the presence of time reparametrization fluctuations and a description of the system in terms of fluctuating local relaxation times. In Sec. III we explain our proposed method to extract local relaxation times from numerical or experimental data. In Sec. IV we discuss the details of the numerical simulations that we have used to test the method, and we present the results of those tests. Finally, in Sec. V, we summarize our results and present some remarks and ideas for further work.

II. CONNECTION BETWEEN FLUCTUATING LOCAL RELAXATION TIMES AND TIME REPARAMETRIZATION FLUCTUATIONS
At temperatures close to the calorimetric glass transition temperature T g , the dynamics of glass forming liquids is dramatically different from the high temperature dynamics [2,28]. At higher temperatures the relaxation is exponential, a direct consequence of Brownian-like motion. By contrast, at lower temperatures the relaxation is exponential only at the early times (with fast dynamics) and is followed by a plateau at intermediate times. This plateau is attributed to the caging effect that results from having each particle temporarily trapped by its neighbors. At longer timescales, as the cages start to break and the particles move distances of the order of the inter-particle distances, the relaxation function starts to decay. Unlike what happens at high temperatures, the final decay, i.e. the long time slow dynamics, is usually not exponential and in many cases it is described well by the Kohlrausch-Williams-Watts (KWW) "stretched exponential" function where we are using some correlation function C(t− t ′ ) between the states of the system at times t and t ′ as an example of a relaxation function. Here A and β are constants and β is sometimes called the "stretching exponent" or the "Kohlrausch exponent". Two extreme scenarios were proposed to explain the non-exponential behavior [2,28]. In one of them the liquid is heterogeneous, the relaxation is exponential in each small region, there is a spatial distribution of relaxation times P (τ ), and the relaxation function is given by where C r (t, t ′ ) is the local two-time correlation in a small region around point r and τ r is the corresponding relaxation time. In the other extreme scenario, the dynamics of the liquid is homogeneous, the distribution of the relaxation times is a delta function, and the relaxation in all regions is described by a unique function g(x), which is non-exponential: Eqs. (2) and (3) can be generalized to allow for the simultaneous presence of both heterogeneous relaxation times and non-exponential local relaxation functions [29]: where P g (τ ) is the probability density of relaxation times compatible with the local relaxation function g(x). If we allow the time difference t− t ′ to be long enough, we cannot assume that each local region is characterized by a single relaxation time τ r , because τ r could in principle have fluctuated over this time interval. To take this effect into account, let us divide the time interval [t ′ , t) into n − 1 sub-intervals {[t i−1 , t i )} i=2,...,n such that for each sub-interval the fluctuations in τ r (t) are negligible. Hence, for long time intervals, the local correlation is given by where t n = t and t 1 = t ′ . In the limit of infinitesimal time intervals we get t i − t i−1 → dt and the sum becomes an integral, .
In this general expression for C r (t, t ′ ), the relaxation is allowed to be locally non-exponential, and the relaxation time is allowed to fluctuate both in space and time. Another way of looking at local correlations is to consider first the exponential relaxation case where If we allow τ to fluctuate in space and time but postulate that the same differential equation still holds, then we have where C r (t, t ′ ) and τ r (t) are defined as before and we have used the initial condition C 0 = lim t→t ′+ C(t, t ′ ). Eq. (8) can be extended to describe the more general case of non-exponential relaxation. For example, we can modify the right hand-side of the differential equation in Eq.
Here the function G[C r (t, t ′ )] represents the effect of non-exponential relaxation. The solution of the new differential equation, which represents the long-time limit of the local correlation function, is given by where φ r (t) ≡ t dt ′ τ r (t ′ ) , and C(x) is a monotonous decreasing function that satisfies C ′ (x) = −G[C(x)] and 0 ≤ C(x) ≤ 1.
A completely different way of deriving Eq. (9) is based on the presence of time reparametrization symmetry, and we describe it next. Time reparametrization symmetry refers to an invariance under transformations of the time variable t → φ(t), where φ(t) is a continuous and monotonically increasing function of t. This symmetry has been shown to be present in the long-time dynamics of mean field spin glass models [30,31] and of short range spin glass models [17,19,21,25]. Numerical studies in spin glasses [19] and structural glasses [20,26] have given strong evidence supporting the presence of this symmetry. The symmetry has also been used, from a different point of view, to determine properties of symmetric observables whose full dynamical properties are hard to compute [32,33]. The symmetry is spontaneously broken by correlations and relaxation functions. We can illustrate the breaking of the symmetry by considering the correlation function , which is only possible if the correlation is independent of time. But we know that in the long time limit the correlation function of a glass forming liquid decays from its plateau value, therefore the symmetry is broken by C(t − t ′ ). The broken reparametrization symmetry is expected to give rise to Goldstone modes [17][18][19][20][21][22][23][24][25][26][27] associated to smooth variations in space of the reparametrized variable, that is t → φ r (t). The spatio-temporal fluctuations of φ r (t) imply that different regions of a sample relax differently from each other, that the relaxation of the different regions is advanced (or retarded) with respect to each other and the bulk, and that "advanced" and "retarded" regions can switch roles in the course of the relaxation. Hence, the spatio-temporal fluctuations induced by the broken reparametrization symmetry provide a possible explanation for dynamical heterogeneities in glass forming and glassy systems. In addition, the broken reparametrization symmetry leads to recovering Eq. (9) if we ignore longitudinal fluctuations [17][18][19]. Longitudinal fluctuations should be suppressed by coarse graining in larger regions because at long distances they are less correlated than transverse fluctuations [34]. Finally, it should be noticed that there is a gauge symmetry under time-independent shifts of the reparametrization variable φ r (t), that is in Eq. (9), the local correlation is unchanged by the transformation φ r (t) → φ r (t) + ρ r . Consequently, when we analyze data we cannot determine an absolute φ r (t). To work with gaugeinvariant quantities, we either have to go back to studying two-time quantities like φ r (t) − φ r (t ′ ) or we must work with time derivatives of φ r (t). By analogy to the spin glass case, we expect that time reparametrization symmetry, if present, can only be exact in the long time limit, and that at any finite timescale there are symmetry breaking terms, which are expected to become weaker at longer times. In particular, if there is a finite relaxation time in the system, it will provide a cutoff for long timescales, thus precluding the symmetry to show its full effect. By contrast, for lower temperatures, the dynamics at longer timescales can be explored, and the effects of the symmetry could in principle manifest themselves more clearly. Based on this discussion, we want to test the hypothesis that at low temperatures and for large coarse-graining regions, the long time correlations are well described by Eq. (9). This is equivalent to the statement that transverse fluctuations (Goldstone modes) are much stronger than longitudinal fluctuations. This interplay between temperature and region sizes implies that the crossover from dominant longitudinal fluctuation to dominant transverse fluctuations should occur for smaller coarse-graining regions for lower temperatures and for larger coarse-graining regions for higher temperatures. In particular, at temperatures moderately below the mode coupling critical temperature T c , we expect that for long timescales the coarse graining size for which the transverse fluctuations start to be dominant will be smaller than at temperatures moderately above T c .

III. METHOD TO EXTRACT PHASES AND INSTANTANEOUS RELAXATION RATES FROM NUMERICAL OR EXPERIMENTAL DATA
Our method for the extraction of φ r (t) from local two-time correlations is inspired by the work in Refs. [26,27]. The authors of Refs. [26,27] use a functional form g global (x) = q EA exp −|x| β , where q EA and β are fitting parameters, to fit global correlations in molecular dynamics simulations of particle systems and polymer systems. They find functions φ(t) for the different systems such that the data are described by As discussed in Sec. II, local fluctuations can be described in terms of Goldstone modes such that In principle, the functions g global (x) and g(x), respectively describing global two-time correlations and local two-time correlations, could have completely different forms. For example, as discussed in Sec. II, some of the initial motivation for the picture of dynamical heterogeneity came from the idea that non-exponential relaxation in a macroscopic sample is due to the combined effect of local exponential relaxations with different relaxation times. In this picture, g global (x) ∝ exp(−|x| β ) and g(x) ∝ exp(−|x|). In Refs. [26,27], as a simplifying assumption, the two functions were taken to be equal. In the present work, we focus exclusively on the local function g(x). However, we will consider different coarse graining sizes when determining C r (t, t ′ ). Each of these coarse graining sizes will give rise to a different form for g(x). In practice, to simplify the determination of g(x), we will restrict it to be a member of a family of functional forms g(x; α), parametrized by a vector α of p components, and for each coarse graining size both the parameter vector α and the fluctuating phases φ r (t) will be determined by fitting C r (t, t ′ ) according to Eq. (12). For technical reasons, we impose the condition that g be an even function, g(−x) = g(x).
Let's consider a data set containing "snapshots" of the relevant degrees of freedom of the system, taken at times {t i } i=1,...,M . In our case the positions of the particles are recorded at each time step. For a given coarse graining region centered at a point r in the sample, our data points will be the M (M − 1)/2 two-time local correlations C r (t i , t j ) calculated from the recorded data. Considering for the moment a fixed parameter vector α, we have M fitting parameters {φ r (t i )} i=1,...,M , thus giving us a fitting problem with ∼ M/2 data points per fitting parameter. However, since the function g(x) is nonlinear, this becomes a nonlinear fitting problem with a large number M of fitting parameters. In order to make this problem manageable, we want to convert it into a linear fitting problem. To do this, we expand the rhs of Eq. (12) to linear order in the local fluctuations where φ(t) is a global phase, still to be determined, and the local fluctuations of the phase are given by At this point, we have simplified the problem by converting it into a linear fitting problem, at the price of introducing the extra phase variables {φ(t i )} i=1,··· ,M . Before we describe the method for determining the phases, we define three quantities to measure fluctuations in the two-time local correlations. The total fluctuations are defined by and contain complete information about the space dependence of the local two time correlations. We decompose the total fluctuations as the sum of a transverse component δC T r (t, t ′ ) and a longitudinal component The transverse component is associated to time reparametrization fluctuations, and it is defined by the linear term in Eq. (13), namely To satisfy Eq. (16), the longitudinal component is defined by In what follows we explain how to perform a least squares fit of the local correlation C r (t, t ′ ) by the rhs of Eq. (13).
In the case of numerical data, we normally have not only data for different coarse grained regions inside the sample, but also data for different independent runs. If we have N R independent simulation runs for a system of volume V , we can imagine juxtaposing the configurations for all the runs, thus creating a larger volume V = N R V . Assuming that each coarse graining region B r centered around point r has a volume V cg , the total number of non-overlapping regions per run is ω = [V /V cg ], and the total number of non-overlapping regions including all runs is The fit we want to perform corresponds to minimizing the quantity with respect to all phase fluctuations and with respect to α. Here is the residual corresponding to the region centered at r, is the number of degrees of freedom in the fit of the data corresponding to one particular coarse graining region, and r denotes a sum over non-overlapping coarse graining regions. From now on we consider only times belonging to the set {t i } i=1,··· ,M , and we simplify our notation by indicating times as subindices, i.e. δφ i r ≡ δφ r (t i ), δC L ji r ≡ δC L r (t j , t i ), and similarly for all other one-and two-time variables.
To minimize E, we observe that for a fixed parameter vector α, the determination of the local fluctuations {δφ i r } for a specific region centered at r can be performed by separately minimizing the corresponding ǫ ({δφ i r } i=1,··· ,M ; α, {C ji r } 1≤i<j≤M ) for that particular region. For this reason, instead of minimizing with respect to all variables in one step, the problem becomes significantly simplified if we perform an iterated minimization, In other words, we first keep α fixed and separately minimize with respect to the phase fluctuations in each coarse graining region, and then minimize the result with respect to α. The minimization with respect to the phases, at fixed α, is performed in two steps. In the first step we determine global phases. In the second step, we use the global phases from the first step to determine local phases. Next, we describe these steps in detail, but without derivations. The details of the derivations can be found in Appendix A. A.
Step One: determining global phases at fixed α The global phases {φ i } i=1,...,M define the point around which the Taylor expansion in Eq. (13) is performed. The only requirement on them is to be close enough to the values of φ i r such that the expansion is accurate to first order in the fluctuations. In particular, we only need to determine them to within an error of the order of the fluctuations in the φ i r . There is more than one possible way to satisfy these conditions, and for simplicity we do it by choosing them to be the phases that best represent the global correlations for the given value of α. We define the quantityǭ where C ji ≡ C(t j , t i ), and g (1) (φ j − φ i ; α) is the first order Taylor expansion of g(φ j − φ i ; α) with respect to the phase difference, taken about To minimizeǭ we impose the condition that all derivatives ofǭ with respect to φ k must be zero. Additionally, we fix the gauge with the condition Thus we obtain the matrix equation The solution to Eq. (29) is the set of global phases {φ i } i=1,...,M which minimizesǭ for the given α. We use the global phases as input in the second step, which we describe next.

B.
Step Two: determining the local phases at fixed α In the second step, we minimize ǫ({δφ i r }; α) with respect to the local phase fluctuations {δφ 1 r , · · · , δφ M r }, for each coarse graining region, while keeping α fixed at the values used in the first step. By Taylor expanding g(φ i r − φ j r ; α) about the global phase differences, φ i − φ j , imposing the condition that all derivatives of ǫ({δφ i r }; α) with respect to the phase fluctuations {δφ 1 r , · · · , δφ M r } should be zero, and fixing the gauge with the condition we obtain the matrix equation where The solution to this matrix equation is the vector δ φ r = (δφ 1 r , · · · , δφ M r ) containing the local fluctuating phase differences that minimize ǫ for the given value of the parameter vector α and for the region B r . By averaging the minimum values of ǫ({δφ i r }; α) for all regions, we obtain ǫ φ ( α).

C. Minimization with respect to α
The optimum value of α is obtained by numerically minimizing ǫ φ ( α) with respect to it. Once the optimum α has been obtained, the solutions of Eq. (33) for that specific value of α and for each region B r give the optimal choice for the fluctuating phases δφ i r .

IV. TESTING THE METHOD WITH DATA FROM NUMERICAL SIMULATIONS OF GLASS-FORMING MODEL SYSTEMS
We tested the method on data from classical Molecular Dynamics simulations of two glass forming systems, each one consisting of an 80 : 20 binary mixture of 1000 particles interacting via purely repulsive Weeks-Chandler-Andersen (WCA) potentials [24][36]. The systems were quenched from an initial temperature T i = 5.0 much higher than the Mode Coupling critical temperature T c = 0.263 ± 0.01 [24]. In [24], T c was determined by measuring the α-relaxation time τ α (T ) and fitting it with a power law as a function of temperature, τ α (T ) = τ 0 Tc T −Tc γ , with γ = 2.1 ± 0.7. One of the systems was quenched to a temperature T = 0.23 ∼ 0.9T c and it did not reach equilibrium during the simulation. We call this system the aging system. The other system was quenched to a temperature T = 0.29 ∼ 1.1T c , and it reached equilibrium within the duration of the simulation. In our studies we use the data corresponding to times after the system reached equilibrium and we call this system the equilibrium system. After the quench, the systems were allowed to evolve for times much longer than the typical vibrational periods of the particles. Snapshots of the systems were taken periodically during the evolution and the positions of the particles were stored.
To probe a system, we divide it into coarse graining boxes B r , centered at the space points r. The local correlation for each box is defined by [20] Here N B r is the number of particles in the coarse graining box B r at time t ′ , and | q| is chosen to be the value given by the location of the main peak of the static structure factor S(q), i.e. q = 7.2. C r (t, t ′ ) measures the extent to which the configuration in the coarse graining region has changed between the times t ′ and t. If the particles initially in the region have not moved much then C r (t, t ′ ) is larger than 1/2, and if the particles have moved a significant distance then C r (t, t ′ ) ≪ 1. In the first case, we say the the region is "slow", and in the second, that it is "fast". To implement the method, we choose the fit function g (x; α) to be a stretched exponential In this case the vector of parameters is α ≡ (q EA , β), where q EA is the plateau value of the correlation function and β is the stretching exponent. We begin the analysis by determining the optimal values of α at different coarse graining sizes; the results are summarized in Fig. (1) for both systems. Each one of the curves is drawn using two  colors: one part is drawn in black and the other in a color other than black. The part that is drawn in black is the part where longitudinal fluctuations, i.e. the fit residuals, account for more than 40% of the total variance of the fluctuations. We observe that the optimal β decreases with increasing coarse graining volume and the optimal q EA increases with coarse graining volume. Both quantities approach constant values at high coarse graining volumes. The trend on q EA appears to be an artifact of the fact that fluctuations are stronger for smaller coarse graining volumes. Since fluctuations that increase C r (t, t ′ ) above a value given by the Debye-Waller factor are unlikely, there is a bias for fluctuations to reduce, rather than increase, the value of C r (t, t ′ ) for times when C r (t, t ′ ) is high. Since the largest values that C r (t, t ′ ) takes are associated to the optimal q EA , this bias will push q EA down, particularly when fluctuations are stronger due to the small size of the coarse graining regions. The trend on β is reminiscent of one of the initial motivations for the proposal that dynamics must be heterogeneous, i.e. the idea [2,3] that in systems for which the bulk relaxation is non-exponential, local regions have exponential relaxation -β ≈ 1 -with heterogeneous timescales, and this translates into a lower exponent β for the system as a whole. However, as we will discuss below, for smaller coarse graining volumes the majority of the fluctuations are captured by the fit residuals and not by the fitting function itself, so it is unclear whether or not this trend in the optimal value of the fit parameter β has any physical significance.

A. Fluctuations
As we mentioned in Sec. II we expect longitudinal fluctuations to be less correlated at long distances, and consequently to be more strongly suppressed by coarse graining, than transverse fluctuations [34]. It is not immediately obvious how to quantify the strength of the fluctuations, since they are functions of the position and of two times. However, in this context it becomes natural to think of them as vectors in an Euclidean vector space, with the inner product defined by and the Euclidean norm defined by Here M is the number of snapshots, { r} are the centers of each one of the ω non-overlapping coarse-graining boxes B r in the volume V of the system, and · · · denotes an average over thermal fluctuations, which in the case of quantities obtained from numerical simulations corresponds to an average over independent simulation runs. The only slightly unusual aspect in these definitions is the prefactor 2/ωM (M − 1), which was included so that the square of the norm is an "intensive" quantity, which gives the average over and therefore they satisfy the Pythagorean condition In Fig. (2) we show how the strengths of the fluctuations vary with the average number of particles N B r per coarse graining box. In the figure, the two top panels correspond to the aging system and the two bottom ones correspond to the equilibrium system. For each system, the left panel shows the magnitudes ||δC||, ||δC T ||, and ||δC L || of the total fluctuations, their transverse component and their longitudinal component; and the right panel shows the ratios ||δC T ||/||δC||, and ||δC L ||/||δC||. Due to the Pythagorean condition, the sum of the squares of these two ratios is unity in all cases.
We expect that as the coarse graining is increased, all fluctuations will be suppressed by the effect of averaging. Indeed, we find that the magnitudes of all three kinds of fluctuations decrease as we increase the coarse graining volume, both for the aging system and for the equilibrium system. The magnitude of the longitudinal fluctuations decreases at a faster rate than the total fluctuations while the magnitude of the transverse fluctuations decreases at a slower rate. Indeed at small coarse graining volumes the ratio ||δC T ||/||δC|| is smaller than the ratio ||δC L ||/||δC||, but the first one increases and the second one decreases as the coarse graining volume increases. For the aging system the two ratios cross at N B r ≈ 30 and for the equilibrium system they cross at N B r ≈ 100. This is consistent with our expectation that longitudinal fluctuations should be more strongly suppressed by coarse graining than transverse fluctuations. We also observe that transverse fluctuations are dominant for a wider range of coarse graining sizes in the case of the aging system than in the case of the equilibrium system. This is indeed what should be expected if the transverse fluctuations are the Goldstone fluctuations associated with time reparametrization symmetry. Since time reparametrization symmetry is a long time asymptotic effect, and in the case of the equilibrium system, which is at higher temperature than the aging system, the relaxation time provides a cutoff for long timescales, we expect that   For both systems, the first ratio is smaller than the second one at small coarse grainings, and larger than the second one at large coarse grainings. The two ratios cross at NB r ≈ 30 for the aging system and at NB r ≈ 100 for the equilibrium system. Error bars are shown, with some of the bars not clearly visible because they are the same size as the symbols or smaller.
the transverse fluctuations associated with this symmetry will manifest themselves less strongly in the equilibrium system than in the aging one. From now on, we consider in detail one coarse graining volume for each system. We choose in each case the smallest coarse-graining volume for which transverse fluctuations capture 60% of the total variance of the fluctuations, i.e. ||δC T || 2 ≥ 0.6 ||δC|| 2 . This corresponds to N B r = 125 for the aging system and N B r = 216 for the equilibrium system. Table I contains a summary of the fitting parameters for those coarse graining volumes.
Since the number of snapshots M used in the analysis cannot be increased indefinitely, it is important to test that the results are robust with respect to changes in M , and that meaningful results can be obtained even for relatively small values of M . There are two aspects to this question: one is to show that the phases φ r (t) are well defined, by showing that their determination is robust with respect to changes in the number of snapshots M , and the other is to show that the magnitudes of the transverse, longitudinal and total fluctuations are not singular as a function of M . With regards to the first aspect, we find that as the number of snapshots is increased, φ r (t) quickly converges to a relatively smooth function. Fig. (3) shows that φ r (t) changes very little if the number of snapshots M used to compute it is 7 or larger in the case of the aging system or if it is 24 or larger in the case of the equilibrium system. It is noticeable that the minimum number M st of snapshots needed for a stable determination of φ r (t) is comparable to the total variation ∆φ of the global phase over the interval considered: in the aging case M st ≈ 7 and ∆φ ≈ 4, while in the equilibrium case M st ≈ 24 and ∆φ ≈ 18. It is tempting to speculate that at least one or two configurations per relaxation time are needed in order to capture a minimum of information about the α relaxation in the system, and this would lead to the prediction that M st /∆φ ∼ 1-2. However at this point the data we have available are insufficient to draw any definite conclusions.  We now consider the issue of the smoothness of the dependence of the magnitudes of the different kinds of fluctuations on the number M of snapshots. In Fig. (4) we plot the magnitude ||δC|| of the total fluctuations, the magnitude ||δC T || of the transverse fluctuations, and the magnitude ||δC L || of the longitudinal fluctuations as functions of 1/M . We find that all three fluctuation magnitudes are gently increasing functions of M , and that each of them appears to approach a constant as M → ∞. For the specific coarse graining sizes shown in the figure, the longitudinal fluctuations are always weaker than the transverse ones, but the opposite case is found for small enough coarse graining sizes.

B. Local Phases and Local Relaxation Times
Up to this point we have presented evidence in favor of the statement that the method we are introducing leads to a robust determination of the transverse and longitudinal components δC T and δC L of the fluctuations and of the fluctuating phases φ r (t). We now use the method to show some examples of the phases φ r (t) and their time derivativesφ r (t) ≡ dφ r (t)/dt. The time derivatives are particularly interesting, not only because they are gauge invariant, but also because, according to Eq. (10), they can in principle be interpreted as local relaxation rates or inverse local relaxation times τ −1 r (t). In Fig. (5) we show the time evolution of the local phases φ r (t) at a fixed point r in the sample for different thermal histories, for both the aging system and the equilibrium system. For comparison, we also show the evolution of the global phase φ(t), and we indeed find that for different thermal histories the local phase φ r (t) gently fluctuates around the global value.
In Figs. (6) and (7) we show the time evolution of the time derivativeφ r (t) of the local phase at the same point in space and for the same thermal histories as in the previous figure. In Fig. (6),φ r (t) is rescaled so that the instantaneous global relaxation rate is unity, and results are shown both for the aging and the equilibrium cases. In the case of the equilibrium system, the rescaling is done approximately, by multiplying the time derivative of the phase by the α-relaxation rate τ α ≡ 1/dφ/dt, where · · · denotes a time average, and the fluctuating values ofφ r (t) are plotted as functions of the rescaled time t/τ α , which is approximately equal to the global phase plus a constant. In the case of the aging system, the rescaling is done by dividing by dφ/dt, and the scaled values ofφ r (t) are plotted as functions of the global phase φ(t). In both cases the time derivative of the local phase fluctuates strongly, often becoming more than twice as large as the global relaxation rate dφ/dt, or even becoming negative. It should be pointed out that at the times whenφ r (t) becomes negative, its interpretation as a local relaxation rate 1/τ r (t) becomes problematic. However, from a more general point of view, it is not surprising that transient fluctuations in a small region of the sample may give rise to changes that make the configuration of the region temporarily become closer to what it was in the past, thus making the two time correlation C r (t, t ′ ) increase over time instead of decreasing, with the effect that, for some time at least,φ r (t) becomes negative. In Fig. (7) we plot the time derivative of the local phase in the aging system, without rescaling it, as a function of time in Lennard-Jones units. For times longer than 2.5 Lennard-Jones time units, the global relaxation rate dφ/dt decreases as a function of time, which is a direct consequence of the aging in the system, namely the fact that the dynamics becomes slower as the relaxation progresses. In this figure we can notice again the strong fluctuations of the time derivative of the local phase with respect to the global value, but we also notice that both the typical value and the fluctuations ofφ r (t) decrease together with the decrease of dφ/dt. Also shown in the same figure is a fit of the global relaxation rate by the form dφ/dt = t t0 α , where t 0 ≈ 3.4 × 10 −4 and α ≈ −0.91, a result that is in agreement with the results in Refs. [26,27]. it is plotted as a function of the global phase φ(t). In the right panel, corresponding to the equilibrium system with NB r = 216 and M = 94, the rescaling ofφ r is performed approxmately by multiplying it with the α-relaxation rate τα ≡ 1/dφ/dt, and the rescaled time derivative is plotted as a function of the rescaled time t/τα, which is approximately equal to the global phase φ(t) plus a constant. Error bars are smaller than the symbols.

V. SUMMARY AND CONCLUSIONS
In this work we have shown that a quantitative description of glassy relaxation in terms of local fluctuating phases and relaxation times is possible. On one hand, starting from the case of exponential relaxation we presented a line of phenomenological arguments that shows how a local phase function φ r (t) can emerge in the description of the data, and how its time derivative can, under certain circumstances, be interpreted as a local fluctuating relaxation rate 1/τ r (t). On the other hand, we showed that in a theoretical framework that postulates the presence of a broken symmetry under time reparametrizations, it is expected that as the dynamics becomes more glassy, fluctuations should be dominated by Goldstone modesi.e. transverse fluctuations -which are described naturally in terms of the same local fluctuating phases φ r (t). Besides establishing a connection between the two points of view, we have also derived a practical method for extracting the phases φ r (t) from numerical or experimental data.
We have applied this method to two numerical simulation datasets, one corresponding to a system undergoing aging at a temperature slightly below the mode coupling critical temperature T c and another one corresponding to a system in equilibrium at a temperature slightly above T c . We have found that in both cases the results are robust with respect to changes in the number M of position snapshots used in the analysis, as long as M is not too small.
The time reparametrization symmetry framework predicts that for larger coarse graining regions and lower temperatures the transverse fluctuations should be dominant, and indeed, this is what we found in our results. For both systems, the ratio between the magnitude of the transverse fluctuations and the magnitude of the total fluctuations grows monotonously with the coarse graining size. For the smallest coarse graining sizes, longitudinal fluctuations dominate, but the transverse fluctuations dominate for large coarse graining sizes. Also, as expected, the range of coarse graining sizes for which transverse fluctuations dominate is wider in the system that is in contact with a heat reservoir at a lower temperature than in the system that is in equilibrium at a higher temperature.
We have shown more detailed results for the fluctuating phase φ r (t) and its time derivativeφ r (t), for a coarse graining size such that the transverse fluctuations capture slightly more than 60% of the total fluctuations. For this coarse graining size, we found that both in the equilibrium system and in the aging system, the time derivative of the local phase fluctuates rather strongly over timescales of the order of the instantaneous global relaxation times or shorter. These fluctuations are strong enough that the time derivative often takes negative values. The presence of those negative values, although not surprising, makes an interpretation ofφ r (t) as a relaxation rate somewhat problematic. It is tempting to speculate that this interpretation may be more cleanly applicable for large enough coarse graining sizes, for which the total fluctuations are expected to be weak enough thatφ r (t) should always remain positive. Since our method involves fitting the local two-time correlations C r (t, t ′ ) by the expression g[φ r (t) − φ r (t ′ )], one of our results is the functional form for g(x). In our fits we modeled this function as a stretched exponential g(x) = q EA exp(−|x| β ). Phenomenological arguments have been made that would indicate that the relaxation should be exponential for small enough coarse graining regions. If this was the case, we should find the fitting parameter β approaching unity for smaller coarse graining sizes. Our results in this respect are inconclusive: although the fitting parameter β does approach unity as the coarse graining size is reduced, and even becomes larger than unity in the case of the aging system, this happens for a range of coarse graining sizes in which the fluctuations are mostly longitudinal, or in other words most of the fluctuations are "explained" by the fitting residuals and not by the fit itself.
One of the ways in which the study of glasses is more difficult than the study of other condensed matter systems is that in the case of glasses there are no simple observables that allow probing glassiness directly. For example, if we want to probe magnetic ordering, a one-time observable, the magnetization S( r, t) , is enough to detect the presence of ordering, and magnetic fluctuations can be probed by the two-point function S( r, t) · S( r ′ , t ′ ) . For the case of crystalline ordering, the ordering can be probed by the one-time Fourier transformed particle density ρ( k, t) ≡ N −1 j exp[i k · r j (t)] , and correlations can be probed by the dynamic structure factor S( k, In the case of structural glasses, by contrast, the presence of a glass state is normally detected by probing the freezing of the density fluctuations, which already requires a two-time observable such as the intermediate scattering function F s ( q, t, t ′ ) or the mean square displacement ∆(t, t ′ ); and the corresponding fluctuations are quantified by a 4-point correlation function such as G 4 ( r, t − t ′ ) or its Fourier transform S 4 ( q, t, t ′ ), which correspond to taking products of two-time observables probing different spatial points in the system. The description of glassy dynamics that we have presented here, although it is extracted from the two-time quantities that probe freezing in the system, has the advantage that it is expressed in terms of one time quantities, the local phases φ r (t), or, more properly, their gauge invariant time derivativesφ r (t). This should make it possible, in principle, to probe correlations in the dynamics by computing two-point functions, such as φ r (t)φ r ′ (t ′ ) , as opposed to four-point functions. One way of thinking about this is that this advantage in simplicity is gained by coarse graining both in space and in time. The coarse graining in space is explicit, and has been discussed extensively in this work; the coarse graining in time happens in the sense that the fitting procedure at the heart of the method focuses on extracting information about the slow part of the dynamics and ignores the fast part. On one hand, it could be argued that this coarse graining comes with two kinds of costs: one is the loss of some information about processes that are "fast" in space and time, and the other is the fact that some extra steps have been introduced in the analysis of the data. On the other hand, it is possible that this loss of information and these extra steps in the analysis may make it easier to observe more cleanly the essential features of the dynamics.
At this point we have introduced a new general analysis tool to study glassy dynamics. Since the starting point of this analysis is the determination of local two-time correlations, it can be applied to any experiment or numerical simulation data set which provides enough information to compute those correlations. This includes not only results from numerical simulations, but also results from visualization experiments [5][6][7]9] in which the particle positions are determined as a function of time. More generally, there are other kinds of experiments, such as photon correlation imaging [35], which provide data for local two-time correlations, and in this case too this type of analysis is applicable.
The investigation presented here leaves many questions open to be addressed. Those questions include the determination of the basic statistical properties of the fluctuating quantityφ r (t), such as its probability distribution and its correlation functions, the study of how those properties vary with changes in the control parameter, namely either the temperature or the volume fraction, and also the study of how they change across different systems exhibiting glass-like dynamics. In particular, the correlation function φ r (t)φ r ′ (t ′ ) encodes important information about dynamical heterogeneity: in the equal-time t = t ′ case it gives very direct information about the spatial extent of the heterogeneous regions, and in the equal-point r = r ′ case it gives information about the persistence in time of those regions. Another set of questions that would be worth investigating is the relationship between the present description of the dynamics, and descriptions based on other probes such as the dynamical susceptibility χ 4 or the various point-to-set correlation functions. Another set of questions that could be addressed is a more detailed study of the form of the function g(x), and in particular the question of the possible presence of simple exponential relaxation at relatively small coarse graining scales.
[36] Both kinds of particles have the same mass M. The number density is ρ = 1.204. The repulsive potential is given by V ab (r) = 4ǫ ab σ ab r 12 − σ ab r 6 + 1 θ(r cutoff ab − r). Here a, b ∈ {A, B} denote the particle type for the two particles that are interacting, σ ab and ǫ ab respectively denote the distance and energy scales for the interaction between a particle of species a and a particle of species b, θ(x) is the Heaviside theta function, and r cutoff ab = 2 1/6 σ ab is the cutoff distance for the potential. The Boltzmann constant is taken to be unity and the (Lennard-Jones) units of length, energy and time are given by σAA, ǫAA and τLJ ≡ (σ 2 AA M/48ǫAA) 1/2 .