{'fraction_non_alphanumeric': 0.033603333920291134, 'fraction_numerical': 0.006661971004284792, 'mean_word_length': 4.427683975788468, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Critical behaviour of large scale dynamical heterogeneities in glasses: a complete theory

In this talk I will present a complete theory for the behaviour of large-scale dynamical heterogeneities in glasses. Following the work arXiv:1001.1746 I will show that we can write a (physically motivated) simple stochastic differential equation that is potentially able to explain the behaviour of large scale dynamical heterogeneities in glasses. It turns out that this behaviour is in the same universality class of the dynamics near the endpoint of a metastable phase in a disordered system, as far as reparametrization invariant quantities are concerned. Therefore Large scale dynamical heterogeneities in glasses have many points in contact with the Barkhausen noise. Numerical verifications of this theory have not yet done, but they are quite possible.


Introduction
In this talk I will give a unified picture of the dynamic heterogeneities that arise in the slow dynamics of glassy systems. I will base my presentation on the recent comprehensive study of the dynamics near the mode coupling transition [1], most of it still being unpublished [2], that is at the basis of most of the analytic conclusions. A very important role in the arguments will be played by the results that have been obtained in [3] on the numerical simulations of glasses forming mixtures and by the extreme large scale simulations of spin glasses that have been done using the Janus machine [4,5,6].
It is refreshing to see how many different pieces of the puzzle fit together and how this new picture helps us to answer to old questions and to formulate new questions. As we shall see many progresses can be done if we concentrate our attention on reparametrization invariant quantities, i.e. on the relations among quantities that are obtained by eliminating the time parametrically. We present here a complete theoretical approach that is able to predict many of the observed quantities.
2 Let us eliminate the time from the plots!

The importance of being reparametrization invariant
The very strong slowing down in glassy material is a very impressive feature. Often most of the attention is devoted to the study of the time dependence of the observables. However this may be not wise.
It is important to understand first what happens and to postpone the study of when it happens. It would be very difficult to understand when, if we do not know what. In plain words let us suppose that we have a movie of the evolution of a glass. We would like to have an understanding of what we see on the various pictures without having to know the exact time at which each picture is taken. We stress the importance of having a wellcharacterized synchronic description before attempting a diachronic one 1 . From this point of view we can get a greater insight if we consider two observables, e.g. A(t) and C(t) and we plots A(t) versus C(t), parametrically in time. Sometimes we say that the corresponding function A(C) is reparametrization invariant, because it would not change if we redefine the time in a diffe! rent way.
When the dynamics become very slow, reparametrization invariant quantities are much simpler to be understood theoretically and they have a much higher degree of universality. Indeed when the dynamics become very slow the actual time dependence becomes irrelevant in deducing the relation among different quantities: we can do adiabatic-like approximations where the systems goes through quasi-equilibrium states. The actual time dependence is much harder to find out: for example in the case of a system where the dynamics is dominated by barriers crossing, the actual time dependence is related to the distribution of the barriers that is a rather complex quantity to compute.
Of course in a dully world, where there is only one quantity that we can measure at a given time, reparametrization invariant quantities would be trivial. Fortunately there are many observations that can be done at the same time, and we have many reparametrization invariant quantities to study. In spite of the fact that reparametrization invariant quantities have already been noticed in the past, their importance is not widely appreciated. We will present now some examples that illustrate their importance and simplicity.
Of course we want (and eventually we will be able) to understand the time dependence of various quantities, but this very difficult task may be successful only after that we establish a firm command of the behaviour of reparametrization invariant quantities.

Fluctuation dissipation relations in off-equilibrium systems.
Fluctuations dissipations relations in off-equilibrium systems are one of the most powerful tools to investigate the properties of the energy landscape of glassy systems. They are well known and I will just recall the definitions and the main results [7,8].
During aging we define a two times correlation function C(t w , t w + t) and a relaxation function R(t w , t w + t). For a given t w one defines a function R tw (C) by eliminating parametrically the time: The fluctuation dissipation relations tell us that the limit is well defined. Moreover the function R(C) (that is experimentally observable) encodes interesting properties of the (free-)energy landscape. In particular it is know that: • The form the function R(C) can be deduced by general arguments: for example in structural glasses we expect that R(C) has slope 1 for C > C * and R(C) has slope m for C < C * , where C * is the position or the plateau in the two times correlation function C(t w , t w + t).
• The function R(C) does not change if the same system is put in an off-equilibrium situations using a method that is different from a temperature jump (as done in aging): e.g. it is subject to a small time dependent force, for example the very gentle stirring of a fluid.
• Even in infinite range models, where all static quantities can be computed analytically, also the function R(C) can be computed analytically in spite of the fact that (at the present moment) we are unable to compute C(t w , t w + t) for large t w .
It is crucial to note that it is possible to relate the function R(C) to a purely equilibrium quantity (i.e. the probability distribution of the overlaps P (q)). Indeed this is the main point of this talk: it ought to be possible to compute reparametrization invariant quantities without having to solve the dynamical equations, but computing the properties of a simple model (in the best of the possible worlds an equilibrium model).
In a similar way we can introduce local fluctuation-dissipation relations that have interesting properties that cannot be described here for lack of space [9].

Dynamical heterogeneities in mean field theory
We consider a very simple case: a mean field model [10,11,12,13,14] (i.e. a random p-spin Ising model on the Bethe lattice) that has a standard dynamical transition at a temperature T d . The correlation time of a system of N spins diverges when N → ∞ and T → T d . For a given system we can define the correlation functions as where the overline denotes the average over t ′ and q(t ′ , t ′ + t) is the (normalized) overlap between a configuration at time t and another configuration at time t ′ . The quantity C 2 N,T (t)−(C N,T (t)) 2 tells us how strongly the correlations fluctuates in time: this is the simplest example of dynamical heterogeneity.
The behaviour of the function C N,T (t) is rather complex: • For T > T d the function C T (t) ≡ lim N →∞ C N,T (t) can be studied using mode-coupling type equations. Exactly at T d for large t we have that C T d (t) ≈ C * + Bt −b (for a review see for example [15]) .
• Exactly at T d the function C T d ,N (t) has a plateau at C * for large N and it decrease to zero when the time goes to infinity slower and slower when N goes to infinity.
• In the mean field case (where the mode coupling equations are correct) one can distinguish a short time β regime where the C ≈ C * + Bt −b and a later large time regime α regime where C ≈ C * − At a + O(t 2a ). The α regime disappear (in mean field theory) below the dynamical temperature.
We stress that at the present moment we have a very little command on the form of the function C N,T (t), we only know the value C * . However in the critical region it is possible to make detailed predictions on the behaviour of the function C 2 N,T (C) defined as These predictions will be described in details in [2]. Let me present the simplest and more surprising one (coming from numerical analysis or from the analytic study of the differential equation (9)). If we define we find If we define the dynamical susceptibility as χ(C) = C 2 d (C) − C 2 we find: The corresponding Binder cumulant can be computed and it is given by At the point where χ(C) is maximum we have C = C * /2 and B = 1. The previous equations have a very simple meaning: in the long time limit for large systems at T d the quantity C(t ′ , t ′ + t) may have (with probability one) two values (0 or C * ). The result is very simple in spite of the rather complex dependence of the correlation function on time and N .
It is interesting that we can obtain the results for the function C 2 d (C) by considering the following differential equation, where we have introduced a fictitious time s: with the boundary condition C η (0) = 1, where W (C) has the shape given in fig. (1), i.e. minimum at the origin and an horizontal flex at C * and η is a random number whose variance goes to zero as a power of N . If we are interested to compute only the function C 2 d (C), the precise N dependence of η is not relevant. We define where · means average over the noise η. If we eliminate the fictitious time s in a parametric way, we find the result in equation (6). The learned reader may be quite surprised: equation (9) is quite different from a mode coupling equation and it does not have all the complications of the mode coupling equation (memory kernel, non-locality in time and all other paraphernalia). This is the main point we are making in this talk: as far as reparametrization invariant quantities are concerned, the two equations are equivalent, in spite of the enormous difference in complexity. Equation (9) is enough to get a physical understanding of the behaviour of the fluctuations.

Two approaches to finite volume analysis
We would like to spend a few words on what happens in a finite volume. In the infinite volume limit the correlation C(t ′ , t ′ + t) does not depends on t ′ , i.e. it is an intensive quantity that has fluctuations proportional to 1/V . When we go to finite volume systems, the fluctuations may be very strong.
We would to point out that in a finite volume there are two different procedures to make the averages: • The procedure that has been described previously: we compute the averages over t ′ of the observables A(t ′ , t ′ + t) and of the correlation C(t ′ , ′ t + t) and we plot parametrically in time A(t) versus C(t).
• For each t ′ we compute At the end we do the average over t ′ of A(t ′ , C) in order to get A(C).
There are advantages and disadvantages in both procedures. It is important to stress that the results may be quite different depending on the chosen procedure. It is likely that in some cases it would be convenient to analyze the data using both approaches.
3 Dynamical heterogeneities in short range models: a naive approach

Overlaps and correlation function
Let us consider a system where the approach at equilibrium is very slow and the correlation functions takes a very long time to decay. The correlation functions have a plateau that may become very long. Quite often the Figure 2: The correlation C(t w , t w + t) in spin glasses below T c during aging: the different curves are taken at different times t w (data taken from [4]). The crosses indicate the extrapolation at infinite waiting time t w , The horizontal line is the asymptotic position of the plateau (estimated using other information).
parameter that controls the length of the plateau is the temperature (near a dynamical critical point) or the waiting time.
Given two configurations that we denote (σ and τ ), we define q σ,τ (x) as the similarity (overlap) of the two configurations in the region of space around the point x. Usually the definition is such that In the spin case q σ,τ (x) is just the product of the spins of two configurations at point x. For fluids different definitions are possible: we can take q σ,τ (x) = 1, if a region around x of size a has the same particle content in the configurations σ and τ ; otherwise, if the particle content is different, q σ,τ (x) = 0. In the simplest case, where there is one kind of particles we can just take q(x) = σ(x)τ (x), where σ(x) and τ (x) are the smeared densities around the point x.
In the aging regime the correlation function C(t w , t w + t) is given by where the overline denotes the space average over the point y and we have used the shorthand notation: We are interested in computing the average of the connected correlation functions These correlations can be used to trace the properties of heterogeneities in the dynamics [16,17,18,19]. Physically the strong slowing down observed in glasses is connected to the existence of cooperatively rearranging regions. These regions manifest themselves as large heterogeneities in the dynamics that imply the existence of long range correlations in the function G(t w , t w + t; x).
Our aim is to get precise predictions on the form of the correlations, on the behaviour of the correlations length, on critical exponents, upper critical dimensions and all that.
According to the previous discussions we eliminate parametrically the time t (the new clock becomes C): we consider the correlation function Our aim is to characterize in some way the expectation values of local quantities in the ensemble where t w goes to infinite at fixed C. It may be convenient to indicate the average as · D C .

A bold conjecture
Let us do something completely different: we take two replicas σ and τ . The first one is a standard equilibrium configuration while the second is a configuration constrained to have an overlap q with the first one. We restrict the statistical sum to all the configurations τ that have an overlap q with σ. In other words we define the following expectation values: where the average over σ is implicit and the subscript R stands for replica. At this stage it is convenient to introduce the so called replica potential W (q) [14,21,22,23,24]. We consider an equilibrium configuration σ. We consider a restricted partition function Z σ (q) where we sum over all the configurations that are at an overlap q with the configuration σ; The restricted partition function Z σ (q) is the denominator in equation (16). Finally where the overline here denotes the average over the different equilibrium configurations σ. Therefore W (q) is β times the increase in free energy when we restrict the sum to those configurations that have an overlap q with a generic configuration σ. In other word exp(−N W (q)) is the typical probability of finding a configuration with overlap q.
The naive conjecture is that in the region of C that corresponds to large times (i.e. those where dC dt is very small), the two previously defined ensemble do coincide for q = C: where in the first case A is the product of factors of the form q(t w , t w + t; x), while in the second case it is the product of factors of the form q σ,τ (x). This is a far reaching conjecture whose validity should strongly simplify the theoretical analysis and the interpretation of the data. It tames all the difficulties of the dynamics and it reduces the analysis to an equilibrium computation. Figure 4: Left panel: the correlation function at distance r as function of q 2 or C 2 , full line equilibrium data, points dynamical data [6]; right panel the correlation length for heterogeneities ζ, from the dynamical data as function C 2 at different t w [4].

A Favorable case: spin glasses
We consider the case of a three dimensional spin glass at zero magnetic field. In this case the replica potential W (q) is symmetric (W (q) = W (−q)) and it satisfies the following relations: It is interesting to look to reparametrization invariant quantities. In this case it should be possible to do an explicit check of the correctness of equation (19). Unfortunately the direct comparison is not so simple because both the t w dependence in the dynamics and the dependence on the size of the box L in the equilibrium data are strong and the extrapolations to infinite t w and infinite L are not trivial. We can use the time-space equivalence principle [6] that state that the finite volume corrections and the finite time corrections are similar for a box of size L * (t): But how to determine L * (t)? We could do a fit in order to match the correlations functions, however we can proceed in a niftier way.
In the spin glass case quenched disorder is present: we can study the correlations induced by disorder. This can be done by simulating the dynamics starting from the two different configurations. We can consider the overlaps among two different configurations both at time t w , i.e. q(σ 1 (t w ), σ 2 (t w ), x) and define their correlation function C 4 (t w , x). According to general principles we should have that The space-time equivalence principle implies that where · R 0 L denotes the expectation values computed in a box of size L. One finds in the low temperature phase (i.e. at a temperature T ≈ 0.6T c ): • Also in the infinite time limit the function C 4 (t w , x) goes to zero: lim tw→∞ C 4 (t w , x) ∝ x −α with α ≈ 0. 4 for large x.
• At finite time the function C 4 (t w , x) goes to zero at large distances exponentially (or faster) and one can define a correlations length ξ(t). 2 The correlation length as a function of time is shown in fig. (3).
• It is possible to satisfy equation (22) with an L * (t w ) ≈ 3.7ξ(t w ).
Using these results one can just determine the correlation length ξ(t w ) from the C 4 (t w , x) correlation and check directly equation (20). The agreement is excellent, see fig. (4), left panel. If we plot the correlation length for the dynamic heterogeneities we find the results in the right panel of fig. (4). One can see that in the β-region (i.e. high C 2 ) the correlation length becomes independent from t w while in the α-region (i.e. low C 2 ) the correlations length increases with t w . The interested reader should look to the original papers in order to know what we have learned on the behaviour of dynamical heterogeneities in this case. Here we want to stress that the naive approach to the computation of the heterogeneities works very well; this is likely due to the vanishing of the replica potential W (q); as we shall see later, in structural glasses the situation is more complex.

Near the dynamical transition 4.1 The behaviour of the replica potential
We have introduced the replica potential W (q) in the previous section. In mean field models near the dynamical transition a secondary minimum appears [14]. The value of W (q * ) at the secondary minimum is the configurational entropy Σ C . Indeed the probability of a given valley is exp(−N W (q * )) and their number is exp(N W (q * )).
In the dynamics we expect naively that there is slowing down in the correlations near the dynamic transition, as discussed in details in section (2.3). Moreover we should have that for T < T d , the equilibrium correlation function does not decays to zero, i.e. lim t→∞ C(t) = C * = 0. However in short range model the convexity of the free energy implies that the correct formula of the potential is given by the convex envelope of W (q), see the insert in fig. (5). It is easy to check that the Maxwell construction implies a constant W ′ (q) in a wide region, that correspond to phase separation in configuration space.
Numerical simulations near the dynamical point for a binary mixture of three dimensional soft sphere clearly show the presence of a Maxwell phenomenon: this is an important test of the basic assumptions of the theory.
One finds that there are two different regimes: • The β-regime, C(t) > C * . We are approaching the saddle point. Heterogeneities are due to fluctuations in the shape of the barrier. This is a quasi-equilibrium situation. The Dynamics gives the same results of the constrained replica approach and the equivalence should be perfect quite near to the plateau.
• The α regime, C(t) < C * . Heterogeneities are due to variations in the escape time due to fluctuations in the shape of the barrier. Dynamics does not give the same results of the constrained replica approach: there is no Maxwell construction in the dynamics, therefore it seems that we do not have any way to use the dynamic-replica connection in this region.
It is interesting to note that the mode coupling equations for the frozen correlations in the region T < T d have been derived in a purely equilibrium setting in a recent remarkable paper [20], adding more evidence for the dynamic-replica connection in the β-regime.

A world of caution: two different correlation functions for dynamical heterogeneities
It is usually not recognized that we can define two different correlation functions for studying the behaviour of dynamical heterogeneities. They display a quite different behaviour. In order to be clear we denote by σ the configuration at time t ′ (the initial configuration) and by τ the configuration at time t ′ + t. We also indicate by an overline the average of the τ configurations and by · σ the average over the τ configurations at fixed σ (we assume that the dynamics is stochastic, as it happens for the Langevin equation or for Monte Carlo). The standard correlation function is However we can define a different correlation function The correlation function G 1 (x) measures the connected correlations of τ at fixed σ, i.e. the heterogeneities induced by the noise, while the correlation function G 0 (x) takes care also of the correlations induced by the fluctuations in the initial configuration. Neglecting loops near the critical point the theoretical prediction (in momentum space) are [21,25,26,2]: where ξ is the correlation length. Unfortunately there are very few measurements done for the second correlation function G 1 , but, especially in the theoretical analysis, one has to take care of its existence. It is also important to note the presence of a double pole in G 0 and not in G 1 .
As we shall see later the situation is very similar to that of a ferromagnetic material in presence of a random field. In this case one can define two correlation functions that are characterized by different exponents [27,28].

A proof of the replica-dynamics connection in perturbation theory
Here we present the hard analytic results that form the backbone of this talk. We can now stop to wave hands and to make analogical argument and we can start to look to the bare formulae. Here I will only present the results; the actual derivations can be found in [1,2]. These analytic results concern mainly the β regime and not the late decay in the α regime.
We will consider for simplicity a short range p-spin model (or any other model with a simple mean field theory). It is crucial to note that we can construct two quite different perturbative expansions: • The perturbative expansion in the dynamics: the starting point (i.e. the zero-order approximation) is the solution of the mode coupling equations in the infinite range model. This perturbation can be constructed for example by using the Martin-Siggia-Rose formalism [29,2].
• We can study the model with one replica constrained to stay at a fixed overlap from a reference configuration. In this case we can just construct the standard loop expansion for the two coupled replicas in the framework of the replica formalism. Essentially we work with an action of an order parameter matrix q ab with a, b = 1, . . . n in the limit of n = 1 replicas. The general diagrammatic rules of [30] can be considerably simplified due to the n = 1 limit and display a great similarity to the Random Field Ising Model [27].
At all order in perturbation theory in the loop expansion one can prove that the leading terms near C * do coincide. More precisely in the limit T → T d and C near to C * (and bigger than C * ), the two theories only differ by subleading terms. Therefore in perturbation theory they have the same critical exponents.
There are many more detailed results; the main ones are: • The correlations (neglecting loops) are given by the previous expressions (25).
• The upper critical dimension is 8: the previous expressions (25) for the correlations should be valid for dimensions greater than 8.
• The theory can be related to a stochastic differential equation that has a clear physical meaning (we will explore this point later on). This equivalence can be shown diagrammatically at all orders and also in a non-perturbative way similar to [31].
• Perturbative dimensional reduction is valid, i.e. the theory is equivalent to a non-random theory in dimensions D − 2.
• The sum of the perturbative expansion is ill-defined (sic erat in fatis). Perturbation theory indeed misses the existence of the Maxwell construction.
I will concentrate my attention on the equilibrium stochastic differential equation. Let us present one of the most striking results: if we consider the leading terms in the perturbative expansion in the low temperature phase, just below T d , we can consider the following stochastic differential equation: where ω(x) is a Gaussian short range noise, (i.e. ω(x)ω(y) = gδ(x − y)). The quantity g plays the role of the coupling constant 3 . Let us consider the overlap profile at large times. In the region T < T d we define In perturbation theory q ∞ (t ′ , x) is a random (correlated) variable, that depends on the initial conditions because in the region below T d the dynamics is no more ergodic and the memory of the initial condition is persistent. We can prove that in perturbation theory the statistical properties of the overlap profile q ∞ (t ′ , x) are just the same of the solution of the previous stochastic differential equation, if g is chosen in a suitable way. Therefore one finds for the correlations of the heterogeneities at large time the following relations: The physics is quite clear. The choice of the initial conditions in the dynamics induce a fluctuating (pointdependent) linear term in the potential (that induces local fluctuations of the critical temperature) and the effects of these fluctuations are the dominant one.
The previous results can be generalized to the late time correlations. Indeed we can introduce a fictitious time s and write the following evolution equation: with the boundary conditions at time s = 0 given by q(x, 0) = 1. We can argue that in the whole β region (always near to the critical point), if we limit ourselves to reparametrization invariant quantities, the previous equation gives the correct results. It is interesting to note that the previous equation just describes the naive decay of the overlap profile in a system with a potential We remark that the fictitious time s in the previous equation should not be identified with the real time: the dependence of q(s), if we solve equation (29), is quite different from the behaviour of q(t) in mode-coupling theory and there should be no confusion among the two time-like variables s and t.

Beyond Perturbation theory
It can be shown that eq. (26) starts to have many solutions when we approach the dynamical temperature. In presence of multiple solutions (as it happens for the Random Field Ising Model) everything becomes more complex (multiple solutions cannot be seen in perturbation theory so that this problem does not affect the perturbative analysis that we have presented above). Therefore perturbation theory is ill-defined. This is a problem also for the equilibrium differential equation (26) because in presence of multiple solutions we have to make a choice and select one of the solutions.
The dynamic formulation, eq. (29) is free from these ambiguity and we bet that it gives the correct answer (the question cannot be decided in perturbation theory). Therefore in the region where the equation has many solutions, the solution relevant for the dynamics is uniquely identified by We are near to the end of our journey. We still have to compute the critical exponent and the critical behaviour of the correlation of q * h (x). This can be easily done in 8 − ǫ dimensions, but in dimension 3 we fear that only a numerical simulation of equation (29) could give the correct results. The simulation is quite easy and there should be no problem in extracting the correct results.
What happens in the more difficult α region? We conjecture that, if we consider only reparametrization invariant quantities, the previous stochastic evolution equation gives the correct results for reparametrization invariant quantities, both in the β and in the α region. Indeed the differential equation (29) make sense also in the α region: we know already that it gives the correct results in the mean field case, see eq. (6).
If this conjecture is correct we have a very simple and unified picture for describing the heterogeneities both in the β and in the α regions. The behaviour in the α region should be very similar to the one described in the recent paper [33], i.e. the nucleation of regions of small overlap that invade the space. A hand waving argument suggest that that in the fully fledged α region, the overlap susceptibility is proportional to ξ 3 (as observed firstly in [18] at low temperatures and recently [35] also near the dynamical transition), however numerical simulations are needed to check the validity of the assumptions.
Our approach has many points in common with the one of [33]: these authors say that the dynamics of glasses near the dynamical phase transition has many points in common with the process of decay of the unstable phase near the end of the metastable phase. Our point is more radical: when we consider the case where some space disorder is present, the two physical situations belongs to the same universality class, if we restrict our study to reparametrization invariant quantities.
Up to now we have not included activated processes. The previous conjecture may be naturally generalized by introducing a more general evolution equation that can be used to compute the properties of reparametrization invariant quantities when activated processes are present below T c : where η(x, s) is just a thermal noise. Apparently the physics is right. This equation seems to predict that the crossing of the barrier is going through a sequence of noise induced avalanches [32]. The previous equation is quite appealing because of its simplicity. It would be very interesting to carefully compare the predictions of these equations with the experimental data, also in order to understand if something is missing. We finally consider one of the few presentations of numerical data in a reparametrization invariant way: i.e. the evaluation of the Binder parameter of the overlap in finite volume [35]( see fig. (6)) at the point of maximum of the dynamical susceptibility. The function B(V /ξ 3 ) is computable by solving numerically equation (29). However we can already see without solving numerically the equation that there is a qualitative agreement with the theoretical predictions. Indeed in the small volume limit we should get the same result of the mean field theory (see section (2.3): a bimodal (symmetric) distribution of the overlap and therefore B = 1, i.e.B = −2/3 in the normalization of figure (6). One the other end for large volume we expect a Binder parameter decaying as the inverse of the volume.
It is remarkable that equation (32) should also describe the large scale behaviour of the approach to equilibrium in a wide class of systems where avalanches and crackling noise is present [36,37,38], for example in ferromagnetic systems where quenched randomness is present and Barkhausen noise is present. The subject has been widely studied and the reader should look to the original papers [36,37,38]. We only remark that region that is relevant for the glassy transition is the small randomness region, where one has eventually to take care of the activated processes. Indeed, if activated processes are neglected and we consider equation (29) a first order transition is present [36,37,38].

Conclusions
We have written various equations, that can be very easily solved numerically, that should be enough to compute the properties of the dynamical heterogeneities, if we limit ourselves to the sector of reparametrization invariant quantities. The results in the β region are very solid, while they become less solid when we enter in the α region and we consider activated process.
We have proposed a single stochastic differential equation that allow us to obtain a unified picture of the heterogeneities related to the glass transition. This equation also describes the dynamics near the endpoint of a metastable phase in a disordered system. Therefore, as far as reparametrization invariant quantities are concerned, heterogeneities in glassy systems have many points in contact with the Barkhausen noise [36,37,38] and they should belong to the same universality class.
It would be very interesting to compare the predictions of the theory (including the various critical exponents) with the experimental and numerical data, with the aim of understanding the points of contacts between theory and experiments and to understand the possible discrepancies. In order to reach this goal it would be very important if the (experimental or numerical) data were presented also in reparametrization invariant form, where the time is eliminated in a parametric way. This would help not only the comparison with the present theory, but also among different systems.