Third and fifth harmonic responses in viscous liquids

We review the works devoted to third and fifth harmonic susceptibilities in glasses, namely and . We explain why these nonlinear responses are especially well adapted to test whether or not some amorphous correlations develop upon cooling. We show that the experimental frequency and temperature dependences of and of have anomalous features, since their behavior is qualitatively different to that of an ideal dipolar gas, which is the high temperature limit of a fluid. Most of the works have interpreted this anomalous behavior as reflecting the growth, upon cooling, of amorphously ordered domains, as predicted by the general framework of Bouchaud and Biroli (BB). We explain why most—if not all—of the challenging interpretations can be recast in a way which is consistent with that of BB. Finally, the comparison of the anomalous features of and of shows that the amorphously ordered domains are compact, i.e. the fractal dimension df  is close to the dimension d of space. This suggests that the glass transition of molecular liquids corresponds to a new universality class of critical phenomena.

Abstract. We review the works devoted to third and fifth harmonic susceptibilities in glasses, namely χ (3) 3 and χ (5) 5 . We explain why these nonlinear responses are especially well adapted to test whether or not some amorphous correlations develop upon cooling. We show that the experimental frequency and temperature dependences of χ (3) 3 and of χ (5) 5 have anomalous features, since their behavior is qualitatively dierent to that of an ideal dipolar gas, which is the high temperature limit of a fluid. Most of the works have interpreted this anomalous behavior as reflecting the growth, upon cooling, of amorphously ordered domains, as predicted by the general framework of Bouchaud and Biroli (BB). We explain why most-if not all-of the challenging interpretations can be recast in a way which is consistent with that of BB. Finally, the comparison of the anomalous features of χ (5) 5 and of χ (3) 3 shows that the amorphously ordered domains are compact, i.e. the fractal dimension d f is close to the dimension d of

Why measuring harmonic susceptibilities? Some facts and an oversimplified argument
Most of our everyday materials are glasses, from window glasses to plastic bottles, and from colloids to pastes and granular materials. Yet the formation of the glassy state is still a conundrum and the most basic questions about the nature of the glassy state remain unsolved, e.g. it is still hotly debated whether glasses are genuine solids or merely hyperviscous liquids. Over the past three decades, the notion evolved that higher-order harmonic susceptibilities are especially well suited to unveil the very peculiar correlations governing the glass formation, yielding information that cannot be accessed by monitoring the linear response. This is illustrated in figure 1 displaying the third harmonic cubic susceptibility χ (3) 3 defined in section 2.1 for four very dierent kinds of glasses [1][2][3][4][5][6]. In the case of spin glasses [1,7] see figure 1(A), it was discovered in the eighties that χ (3) 3 diverges at the spin glass transition temperature T SG , revealing the long range nature of the spin glass amorphous order emerging around T SG . Here the expression 'amorphous order' corresponds to a minimum of the free energy realized by a configuration which is not spatially periodic. Similar non-linear susceptibility experiments have been performed by Hemberger et al [2] on an orientational glass former. In orientational glasses, electric dipolar or quadrupolar degrees of freedom undergo a cooperative freezing process without long-range orientational order [8]. As illustrated in figure 1(B), the divergence of |χ (3) 3 | is not accompanied by any divergence of the linear susceptibility |χ 1 |. We shall show in equations (1) and (2) that this is intimately related to the very notion of amorphous ordering. For structural glasses, e.g. glycerol, it was discovered [3,4] less than 10 years ago that |χ (3) 3 (ω, T )| has a hump close to the α relaxation frequency f α , the height of this hump increasing anomalously upon cooling-by anomalous increase we just mean that this increase has to be attributed to correlations between molecules since this increase does not arise without correlations, e.g. in a dipolar gas. A hump of |χ (3) 3 | has also been recently discovered in a colloidal glass [5,6], in the vicinity of the β relaxation frequency f β , revealing that any shear strain connects the system to a non equilibrium steady state, see [5,6]. Of course, as detailed balance does not hold in colloids, the comparison of colloidal glasses with spin glasses, orientational glasses, and structural glasses cannot be quantitative. However, the four very dierent kinds of glasses of figure 1 have the common qualitative property that nonlinear cubic responses unveil new information about the glassy state formation.
Let us now give an oversimplified argument explaining why nonlinear responses should unveil the correlations developing in glasses. We shall adopt the dielectric  [1], the value of χ (3) 3 is plotted in logarithmic scale as a function of frequency, for various values of the reduced temperature (T − T c )/T c . Note that the static value of χ (3) 3 increases by two orders of magnitude when T varies by less than 8%: this can be fitted by a critical law where the static value of χ (3) 3 diverges when approaching the critical temperature T c 2.94 K, see [1] for more details about the values of the critical exponents. Reprinted figure with permission from [1], copyright (1986) by the American Physical Society. (B): similar arguments are used to rationalize the third-harmonic dielectric susceptibility of an orientational glass [2]. Reprinted figure with permission from [2], copyright (1996) by the American Physical Society. (C): in glycerol [3,4], the modulus of the dimensionless cubic susceptibility X (3) 3 has a peak as function of frequency, which increases anomalously upon cooling: the expression 'anomalous increase' means that this increase is absent in a system without any correlations e.g. in a dipolar gas. Therefore this increase can be attributed to correlations between molecules, even though a genuine divergence is not observed-see text for an explanation of the dierence of the nonlinear responses between spin glasses and molecular glasses. (D): strain-stress experiment in the colloidal system studied in [5,6]. When increasing the volumic density φ the increasing peak of Q 0 = |χ (3) 3 /χ 1 | reveals that any shear strain connects the system to a non equilibrium steady state, see [5,6]. In all these four examples χ (3) J. Stat. Mech. (2019) 124003 language adapted to this review devoted to supercooled liquids where detailed balance holds, and consider a static electric field E st applied onto molecules carrying a dipole moment µ dip . At high temperature T the system behaves as an ideal dipolar gas and its polarization P is given by: where a d is the molecular d-dimensional volume, L d is the suitable Langevin function expressing the thermal equilibrium of a single dipole in dimension d, and where the numerical prefactors of the linear, third, and fifth order responses correspond to the case d = 3. Assume now that upon cooling some correlations develop over a characteristic lengthscale , i.e. molecules are correlated within groups containing N corr = ( /a) d f molecules, with d f the fractal dimension characterizing the correlated regions. The characteristic lengthscale can be directly deduced from the spatial decrease of the suitable correlation function-see section 3. Therefore the corresponding domains must be regarded as being statistically independent from each other, and one can use equation (1), provided that we change the elementary volume a d by that of a domain, namely a d ( /a) d , as well as the molecular dipole µ dip by that of a domain, namely µ dip ( /a) (d f /2) . Here, the exponent d f /2 expresses the amorphous ordering within the correlated regions, i.e. the fact that the orientation of the correlated molecules looks random in space, allowing the use of the law of large numbers-should we want to describe a para-ferro transition we would use the exponent d f (and not d f /2) to express the dipole of a correlated domain. We obtain in the case of amorphous ordering: which shows that the larger the order k of the response, the stronger the increase of the response when increases. As d f d, equation (2) shows that the linear response never diverges with : it is always, for any , of the order of µ 2 dip /(a d k B T ). This can be seen directly in equation (2) in the case d f = d ; while for d f < d one must add to equation (2) the polarization arising from the uncorrelated molecules not belonging to any correlated region. This insensitivity of the linear response to directly comes from the amorphous nature of orientations that we have assumed when rescaling the net dipole of a domain by using the power d f /2. By contrast, in a standard para-ferro transition one would use instead a power d f to rescale the moment of a domain, and we would find that the linear response diverges with as soon as d f > d/2 which is the standard result close to a second order phase transition. For amorphous ordering, the cubic response is thus the lowest order response diverging with , as soon as d f > d/2. This is why cubic responses as well as higher order responses are ideally suited to test whether or not amorphous order develops in supercooled liquids upon cooling. https://doi.org/10.1088/1742-5468/ab371e For spin-glasses, the above purely thermodynamic argument is enough to relate the divergence of the static value of χ (3) 3 see figure 1(A) to the divergence of the amorphous order correlation length . For structural glasses this argument must be complemented by some dynamical argument, since we have seen on figure 1(C) that the anomalous behavior of χ (3) 3 takes place around the relaxation frequency f α . This has been done, on very general grounds, by the predictions of Bouchaud and Biroli (BB), who anticipated [9] the main features reported in figure 1(C). Because BB's predictions will be used to interpret the experimental behavior of nonlinear responses in molecular glass formers, they will be explained, and tested by a comparison to experiments, only in section 3, after having reviewed in section 2 the experimental features of nonlinear responses in molecular glass formers. Section 3 will be followed by section 4 where more specific approaches of the cubic response of glass forming liquids will be explained and compared to experiments. Beyond their apparent diversity, we shall show that these more specific approaches can be unified by the fact that in all of them, N corr is a key parameter even though it is sometimes implicit. Finally, the appendix contains some additional material for the readers aiming at deepening their understanding of this field of high harmonic responses.
The rationale for not devoting section 2 to theoretical considerations but instead to the review of the main experimental features of third and fifth harmonic susceptibilities is twofold. Firstly, because of the generality of equation (2), we are led to anticipate that χ 3 and χ 5 have common anomalous features that can be interpreted as reflecting the evolution of and thus of N corr upon cooling. Secondly, the nonlinear response of a system without any correlations i.e. a dipolar gas is exactly known: this will serve as a benchmark, and any qualitative dierence between this benchmark and the experimental behavior will be attributed to correlations between molecules. To make the comparison easier, we shall define dimensionless nonlinear responses in such a way that they are independent on the temperature T in a dipolar gas, which will allow a thorough discussion of the T dependence observed at finite frequencies in molecular glass formers. Overall, the data and considerations reported in section 2 have been chosen to elicit directly the eect of correlations in molecular glasses on a general basis, and therefore the content of section 2 is meant to gather facts allowing to test not only the known theoretical approaches but also the future ones.

Definitions
When submitted to an electric field E(t) depending on time t, the most general expression of the polarisation P (t) of a dielectric medium is given by a series expansion: where because of the E → −E symmetry, the sum contains only odd terms, and the (2m + 1)-order polarisation P 2m+1 (t) is proportional to E 2m+1 . The most general expression of P 2m+1 (t) is given by, see e.g. [1] or equation (20) of [10]: Because of causality χ 2m+1 ≡ 0 whenever one of its arguments is negative. Even though the right hand side of equation (4) looks quite involved, it can be considerably simplified by performing exactly the same operations as for the standard linear response: with the obvious extension for the (2m + 1)-fold Fourier transform of any function φ(t 1 , t 2 , ...t 2m+1 ). Note that the Fourier modes exist both for positive and negative angular frequencies. In the physically important case where φ(t) has only real values, one shows that φ(ω) is the complex conjugate of φ(−ω).
In the case of the third order polarisation, this yields, see [11] as well as appendix D: where we have set χ 3 (ω, ω, ω) = |χ , and χ 3 (ω, ω, −ω) = |χ . Similarly, for the fifth-order polarisation, we obtain: where, we have set , and similarly . For completeness we recall that the expression of the linear polarisation P 1 (t) is In the linear case, we often drop the exponent indicating the harmonic, since the linear response P 1 (t) is by design at the fundamental angular frequency ω. The only exception to this simplification is in figure 11 (see below) where for convenience the linear susceptibility is denoted χ (1) 1 .
Up to now we have only considered nonlinear responses induced by a pure ac field E, allowing to define the third harmonic cubic susceptibility χ and/or the fifthharmonic fifth-order susceptibility χ (5) 5 to which this paper is devoted. In section 2.3 and figures 8 and 9, we shall briefly compare χ (3) 3 with other cubic susceptibilities, namely χ (1) 3 already defined in equation (5) as well as χ (1) 2;1 that we introduce now. This supplementary cubic susceptibility is one of the new terms arising when a static field E st is superimposed on top of E. Because of E st , new cubic responses arise, both for even and odd harmonics. For brevity, we shall write only the expression of the first harmonic part P (1) 3 of the cubic polarization, which now contains two terms: P (1) 3 (ω)) + 3|χ 2,1 (ω)) (7) where we have defined |χ (1) 2,1 (ω)| exp (−iδ (1) 2,1 (ω)) = χ 3 (0, 0, ω). For any cubic susceptibility-generically noted χ 3 -or for any fifth-order susceptibility-generically noted χ 5 -the corresponding dimensionless susceptibility X 3 or X 5 is defined as: is the linear susceptibility at zero frequency and χ lin (∞) is the linear susceptibility at a high frequency where the orientational mechanism has ceased to operate. Note that X 3 as well as X 5 have the three important advantages: (i) they are dimensionless, (ii) they are independent of the field amplitude, (iii) they are independent of temperature in an ideal gas of dipoles [19], which justifies the use of X 3 and X 5 to elicit the possible temperature dependence of the eect of correlations in the nonlinear responses of molecular glass formers. Before moving to the experimental behavior of third and fifth order responses, we justify the series expansion made in equations (1)-(3) by two kinds of experimental arguments. Firstly, in practice it turns out that, with the typical fields used experimentally, the linear response is much larger than the cubic responses (typically by two J. Stat. Mech. (2019) 124003 orders of magnitude), which are themselves typically two orders of magnitude larger than the fifth order responses-this is why experiments are so dicult. This strong hierarchy between the successive orders of the response in the field gives, in our view, a justification for the series expansion made in equations (1)-(3). In such a series expansion, the coecients of the various terms E 2m+1 are derivatives of order (2m + 1) of the P (E(t)) curve: these derivatives are thus calculated at zero field. This strongly suggests that the various nonlinear suceptibilities reported hereafter are related to the physics taking place at zero field. This statement about the validity of the series expansion is reinforced by a second (and more specific) fact: most (if not all) of the nonlinear susceptibilities reported hereafter have been extracted after having thoroughly checked, see e.g. [3,11,54], that the associated nonlinear signal follows the expected ∝ E 2m+1 dependence for its modulus with a field independent value of the phase. Again this allows to conclude that the measured nonlinear susceptibilities are related to physics at zero field. We remind the reader that, from the strictly theoretical point of view, the theorems expressing the susceptibilities in terms of the appropriate correlation function at zero field only exist, for an arbitrary frequency, in the case of the linear response. For nonlinear susceptibilities, such theorems only exist for the static case ω = 0, see section 3.

Frequency and temperature dependence of third harmonic susceptibility
In this section we review the characteristic features of χ 3 both as a function of frequency and temperature. We separate the eects at equilibrium above T g and those recorded below T g in the out-of-equilibrium regime.

Above T g .
In the α regime. Figure 2 shows the modulus |χ (3) 3 | for propylene carbonate [12]. It is an archetypical example of what has been measured in glass forming liquids close to T g . For a given temperature one distinguishes two domains: 1. For very low frequencies, f /f α 0.05, a plateau is observed as indicated by the shaded area in figure 2, i.e. |χ 3 | does not depend on frequency. This is reminiscent of the behavior of an ideal gas of dipoles where each dipole experiences a Brownian motion without any correlation with other dipoles. In such an ideal dipolar gas, |χ 3 | has a plateau below the relaxation frequency and monotonously falls to zero as one increases the frequency. Because the observed plateau in figure 2 is reminiscent to the ideal dipolar gas case, it has sometimes [3,4] been called the 'trivial' regime. What is meant here is not that the analytical expressions of the various χ 3 are 'simple'-see section 7-,but that the glassy correlations do not change qualitatively the shape of χ 2. When rising the frequency above 0.05f α one observes for |χ (3) 3 | a hump for a frequency f peak /f α c where the constant c does not depend on T and weakly depends on the liquid (e.g. c 0.22 for glycerol and c 0.3 for propylene carbonate). This hump is followed by a power law decrease |χ [3] to the exponent governing the decrease of |χ 1 | above f α . Qualitatively, this hump is important since it exists neither in the cubic susceptibility of an ideal gas of dipoles nor in the modulus of the linear response |χ 1 | of the supercooled liquids. This is why this hump has been termed the 'glassy contribution' to χ 3 . On a more quantitative basis, the proportionality of f peak and of f α has been observed for f α ranging from 0.01 Hz to 10 kHz (above 10 kHz the measurement of χ (3) The consistency of the above considerations can be checked by comparing the thirdorder susceptibility of canonical glass formers to that of monohydroxy alcohols. The linear dielectric response of the latter is often dominated by a Debye relaxation process, which is commonly ascribed to the fact that part of the molecules are forming chain-like hydrogen-bonded molecule clusters with relatively high dipolar moments [14]. This process represents an idealised Debye-relaxation case as it lacks the heterogeneity-related broadening found for other glass formers. Moreover, correlations or cooperativity should not play a significant role for this process, because cluster-cluster interactions can be expected to be rare compared to the intermolecular interactions governing the α relaxation in most canonical glass formers [15]. Thus, this relaxation process arising from rather isolated dipolar clusters distributed in a liquid matrix can be expected to represent a good approximation of the 'ideal dipolar gas' case mentioned above. The monohydroxy alcohol 1-propanol is especially well suited to check this notion because here transitions between dierent chain topologies, as found in several other alcohols aecting the nonlinear response [16,17], do not seem to play a role [17]. Figure 3(a) shows the frequency-dependent modulus, real, and imaginary 3 E 2 for 1-propanol at 120 K [15,18]. Indeed, no hump is observed in |χ (3) 3 |(ν) as predicted for a non-cooperative Debye relaxation. The solid lines were calculated according to [19], accounting for the expected trivial polarization-saturation eect. Indeed, the spectra of all three quantities are reasonably described in this way. In the calculation, for the molecular volume an additional factor of 2.9 had to be applied to match the experimental data, which is well consistent with the notion that the Debye relaxation in the monohydroxy alcohols arises from the dynamics of clusters formed by several molecules.
In marked contrast to this dipole-gas-like behavior of the Debye relaxation of 1-propanol, the χ 3 spectra related to the conventional α relaxation of canonical glass formers exhibit strong deviations from the trivial response, just as expected in the presence of molecular correlations. As an example, figure 3(b) shows the modulus, real, and imaginary part of χ 3 E 2 of glycerol at 204 K. Again the lines were calculated assuming the trivial nonlinear saturation eect only [19]. Obviously, this approach is insucient to provide a reasonable description of the experimental data. Only the detection of plateaus in the spectra arising at low frequencies agrees with the calculated trivial 3 (times E 2 ) of 1-propanol at 120 K as measured with a field of 468 kV cm −1 [18]. The solid lines were calculated according to [19]. (b) Same for glycerol at 204 K and 354 kV cm −1 [18]. Reproduced from [18] (2017) (© EDP Sciences and Springer-Verlag GmbH Germany 2017). With permission of Springer.

J. Stat. Mech. (2019) 124003
response. This mirrors the fact that, on long time scales, the liquid flow smoothes out any glassy correlations.
When varying the temperature, two very dierent behaviors of χ 3 are observed: 1. In the plateau region the weak temperature dependence of χ 3 is easily captured by converting χ by using equation (8): one observes [3,4] that in the plateau region X does not depend at all on the temperature. Qualitatively this is important since in an ideal gas of dipoles X does also not depend on temperature, once plotted as a function of f /f α . This reinforces the 'trivial' nature of the plateau region, i.e. the fact that it is not qualitatively aected by glassy correlations.

In the hump region, |X
3 ( f /f α )| increases upon cooling, again emphasizing the 'anomalous'-or 'non trivial'-behavior of the glassy contribution to χ (3) 3 . This increase of the hump of |X (3) 3 | has been related to that of the apparent activation energy E act (T ) ≡ ∂ ln τ α /∂(1/T ), see [12,20], as well as to T χ T ≡ |∂ ln τ α /∂ ln T | [3,4,21,22]. Note that because the experimental temperature interval is not so large, the temperature behavior of E act and of T χ T is extremely similar. Both quantities are physically appealing since they are related to the number N corr (T ) of correlated molecules: the line of thought where E act ∼ N corr (T ) dates back to the work of Adam and Gibbs [23]; while another series of papers [24,25] proposed a decade ago that N corr ∝ T χ T . Figure 4 illustrates how good is the correlation between the increase of the hump of |X (3) 3 |, left axis, and E act (T ). This correlation holds for 5 glass formers, of extremely dierent fragilities, including a plastic crystal, where only the orientational degrees of freedom experience the glass transition [28]. It should be noted that, within the framework described above, the case of 1-propanol seems to be somewhat inconsistent: as revealed by figure 3(a), the modulus of its χ 3 does not exhibit a hump and can be rather well fitted without invoking cooperativity. Therefore, for its Debye relaxation one would not expect any deviations from Arrhenius behavior caused by cooperativity. However, while the Debye-relaxation time of 1-propanol exhibits a weaker temperature dependence than the α relaxation [26], it still deviates from Arrhenius behavior. For the monohydroxy alcohol 2-ethyl-1-hexanol (2E1H) the situation is dierent and a hump is detected in |χ (3) 3 | [12]. As revealed by figure 4 (squares and corresponding line), N corr (T ) of this alcohol is consistent with the weak non-Arrhenius behavior of its Debye-relaxation [27], mirrored by a weakly varying E act (T ). It seems that the molecule clusters, commonly believed to cause the Debye relaxation in monohydroxy alcohols [14], at least to some extent interact with each other in 2E1H. For 1-propanol, one may speculate that its 'trivial' nonlinear response arising from polarization saturation may be significantly stronger compared to the cooperativity contribution in |χ

J. Stat. Mech. (2019) 124003
In the excess wing regime. In the dielectric-loss spectra of various glass formers, at high frequencies the excess wing shows up, corresponding to a second, shallower power law at the right flank of the α peak [29]. Figure 5(a) shows loss spectra of glycerol, measured at low and high fields up to 671 kV cm −1 [30,31], where the excess wing is indicated by the dashed lines. (It should be noted that the dierence of these loss curves for high and low fields is directly related to the cubic susceptibility χ (1) 3 , defined in equation (5), [18]). As already reported in the seminal paper by Richert and Weinstein [32], in figure 5(a) at the right flank of the α-relaxation peak a strong field-induced increase of the dielectric loss is found. In [32], no significant field dependence is detected at its low-frequency flank. Similar results were also obtained for other glass formers [31,75]. In [32], it was pointed out that these findings are well consistent with the heterogeneity-based box model (see section 4.3). However, as revealed by figure 5(a), remarkably in the region of the excess wing no significant nonlinear eect is detected. Time-resolved measurements, later on reported by Samanta and Richert [33], revealed nonlinearity eects in the excess-wing region when applying the high field for extended times of up to several 10 000 cycles. Anyhow, the nonlinearity in this region seems to be clearly weaker than for the main relaxation and the nonlinear behavior of the excess wing diers from that of the α relaxation.
To check whether weaker nonlinearity in the excess-wing region is also revealed in higher-harmonic susceptibility measurements, figure 5(b) directly compares the modulus of the linear dielectric susceptibility of glycerol at 191 K to the third-order susceptibility |χ by the ionic and electronic polarizability, whose contribution in the modulus strongly superimposes the excess wing.) While the linear response exhibits a clear signature of the excess wing above about 100 Hz (dashed line), no trace of this spectral feature is  3 (ν)|. Thus, we conclude that possible nonlinearity contributions arising from the excess wing, if present at all, must be significantly weaker than the known power-law decay of the third-order susceptibility at high frequencies, ascribed to the nonlinearity of the α relaxation.
The excess wing is often regarded as the manifestation of a secondary relaxation process, partly superimposed by the dominating α-relaxation [35,36]. Thus the weaker nonlinearity of the excess wing seems to support long-standing assumptions of the absence of cooperativity in the molecular motions that lead to secondary relaxation processes [37,38]. Moreover, in a recent work [39] it was pointed out that the small or even absent nonlinear eects in the excess-wing region can also be consistently explained within the framework of the coupling model [38], where the excess wing is identified with the so-called 'nearly constant loss' caused by caged molecular motions. 3 E 2 at 565 kV cm −1 [34]. The solid lines indicate similar power laws above the peak frequency for both quantities. The dashed line indicates the excess wing in the linear susceptibility at high frequencies, which has no corresponding feature in Below T g , the physical properties are aging, i.e. they depend on the time t a elapsed since the material has fallen out of equilibrium, i.e. since the glass trans ition temperature T g has been crossed. The mechanism of aging is still a matter of controversy [40][41][42][43][44], owing to the enormous theoretical and numerical diculties inherent to out-of-equilibrium processes. Experimentally, a few clear cut results have been obtained in spin glasses [45] where it was shown, by using nonlinear techniques, that the increase of the relaxation time τ α with the aging time t a can be rather convincingly attributed to the growth of the number N corr of correlated spins with t a . Very recently extremely sophisticated numerical simulations have been carried out by the so called Janus international collaboration, yielding, among many other results, a strong microscopic support [46] to the interpretation given previously in the experiments of [45].
In structural glasses, the aging properties of the linear response have been reported more than one decade ago [47,48]. More recently, the aging properties of χ 3 were reported in glycerol [49] and its main outputs are summarized in figures 6 and 7. A glycerol sample previously well equilibrated at T g + 8 K was quenched to a working temperature T w = T g − 8 K and its third harmonic cubic susceptibility was continuously monitored as a function of t a . The dominant eect is the increase of the relaxation time τ α with t a . In [49] τ α increases by a factor 6 between the arrival at T w , i.e. t a = 0, and the finally equilibrated situation reached for t a τ α,eq where τ α is equal to its equilibrium value τ α,eq and no longer evolves with t a . This variation of τ α with the aging time t a can be very accurately deduced from the shift that it produces on the imaginary part of the linear response χ ( f , t a ). This is summarized in figure 6 for 5 dierent frequencies: when plotted as a function of f /f α (t a ) ≡ 2πf τ α (t a ), the aging values of χ ( f , t a ) symbols are nicely rescaled onto the equilibrium values χ ( f , eq), continuous line, measured when t a τ α,eq . The most important experimental result is that this scaling fails for |X  . During the aging of glycerol at T g − 8 K the increase of τ α with the aging time t a is measured by rescaling the aging data symbols of χ 1 -right axis-onto the equilibrium data-solid black line. The corresponding scaling fails for X (3) 3 (f , t a ) -left axis-: this can be interpreted as revealing the increase of N corr during aging (see [49] for more details, and for the definition of the quantity z(t a )/z(T) which is involved in the left axis but varies by less than 2% during aging). aging times, the dierence between aging data (symbols) and equilibrium values (continuous line) is largest. This has been interpreted as an increase of N corr with the aging time t a . This increase of N corr (t a ) towards its equilibrated value N corr (eq) is illustrated in figure 7 where the variation of δ = N corr (t a )/N corr (eq) is plotted as a function of t a . It turns out to be independent of the measuring frequency, which is a very important self consistency check.
The increase of N corr during aging can be rather well captured by extrapolating the N corr (T ) variation obtained from the growth of the hump of |χ (3) 3 | measured at equilibrium above T g and by translating the τ α (t a ) in terms of a fictive temperature T fict (t a ) which decreases during aging, finally reaching T w when t a τ α,eq . This yields the continuous line in figure 7, which fairly well captures the data drawn from the aging of χ 3 . Because this extrapolation roughly agrees with the aging data, one can estimate that the quench from T g + 8 K to T w = T g − 8 K corresponds to a doubling of N corr,eq . The approximately 10% increase reported in figure 7 is thus the long time tail of this increase, while the first 90% increase cannot be measured because it takes place during the quench.
Beyond the qualitative result that N corr increases during aging, these χ 3 (t a ) data can be used to test quantitatively some theories about the emergence of the glassy state. By gathering, in the inset of figure 7, the equilibrium data symbols lying in the [1;1.3] interval of the horizontal axis and the aging data translated in terms of T fict (t a ) symbols lying in the [2;2.3] interval, one extends considerably the experimental temperature interval, which puts strong constraints onto theories. Summarizing two dierent predictions by ln( figure 7 is designed to test these two predictions (see [49] for details): it shows that both of them are consistent with experiments contrary to another prediction relying onto a critical relation τ α ∝ N z corr , yielding an unrealistic large value of z ∼ 20 to account for the experiments.

Strong similarities between third and first cubic susceptibilities
We now come back to equilibrium measurements i.e. above T g and compare the behavior of the third-harmonic cubic susceptibility χ (3) 3 as well as the first-harmonic cubic susceptibilities χ (1) 3 and χ (1) 2;1 introduced in equation (7). We remind that χ (1) 2;1 corresponds to the case where a static field E st is superimposed to the ac field E cos(ωt). Figures 8 and 9 show the modulus and the phases of the three cubic susceptibilities for glycerol and for propylene carbonate.
1. For the modulus: at a fixed temperature, the main features of the frequency dependence of |χ (1) 3 | and of |χ (1) 2;1 | are the same as those of |χ 3 |: when increasing the frequency, one first observes a low frequency plateau, followed by a hump in the vicinity of f α and then by a power law decrease ∼f −β 3 . The most important dierences between the three cubic susceptibilities are the precise location of the hump and the absolute value of the height of the hump. As for the temperature dependence one recovers for |χ (1) 3 | and for |χ (1) 2;1 | what we have already seen for |χ (3) 3 |: once put into their dimensionless forms X 3 the three cubic susceptibilities do not depend on T in the plateau region, at variance with the region of the hump where they increase upon cooling typically as 2. The phases of the three cubic susceptibilities basically do not depend explicitly on temperature, but only on u = f /f α , through a master curve that depends only on the precise cubic susceptibility under consideration. These master curves have the same qualitative shape as a function of u in both glycerol and propylene carbonate. We note that the phases of the three cubic susceptibilities are related to each other. In the plateau region all the phases are equal, which is expected because at low frequency the systems responds adiabatically to the field. At higher frequencies, we note that for both glycerol and propylene carbonate (expressing the phases in radians): which are quite non trivial relations.
3. In the phase of χ 3 of propylene carbonate (figure 9), a jump of π is observed which is accompanied by the indication of a spikelike minimum in the modulus, see [52] for more details. A similar jump may also be present in glycerol ( figure  8). This jump in the phase happens at the crossover between the T-independent 'plateau' and the strongly T-dependent hump. More precisely in the 'plateau' region one observes a reduction of the real part of the dielectric constant χ 1 , while around the hump χ 1 is enhanced. At the frequency of the jump, both eects compensate and this coincides with a very low value of the imaginary part of X In this section, we first explain why measuring χ (5) 5 is interesting for a better understanding of the glass transition. We then see the characteristic features of χ (5) 5 as a function of frequency and temperature.

2.4.1.
Interest in the fifth-order susceptibility. In the previous sections, we have seen that the increase of the hump of |X 3 | upon cooling has been interpreted as reflecting that of the correlation volume N corr a 3 . However in practice, this increase of N corr remains modest (typically it is an increase by a factor 1.5) in the range 0.01 Hz f α 10 kHz where the experiments are typically performed. Physically this may be interpreted by the fact that an increase of N corr changes the activation energy, yielding an exponentially large increase of the relaxation time τ α . Now if one demands, as in standard critical phenomena, to see at least a factor of 10 of increase of |X 3 | to be able to conclude on criticality, one is lead to astronomical values of τ α : extrapolating the above result, e.g. |X 3 | ∝ |∂ ln τ α /∂ ln T | and assuming a VFT law for τ α , one concludes that the experimental characteristic times corresponding to an increase of |X 3 | by one order of magnitude is 0.1 ms τ α 10 18 s. This means experiments lasting longer than the age of the universe.
This issue of astronomical time scales can be circumvented by using a less commonly exploited but very general property of phase transitions: close to a critical point all the responses diverge together [53], since the common cause of all these divergences is the growth of the same correlation length. Showing that all the responses of order k behave as a power law of the first diverging susceptibility is another way of establishing criticality. For glasses, we have seen in equation (2) that, apart from χ 1 which is blind to glassy correlations, all other responses χ k 3 grow as power laws with the amorphous ordering length : χ 3 ∝ (l/a) 2d f −d and χ 5 ∝ (l/a) 3d f −d . Therefore, assuming that the main cause for the singular responses appearing in the system is the development of correlations, there should be a scaling relation between the third and fifth order responses, namely one should observe Measuring χ 5 is of course extremely dicult, because, for the experimentally available electric fields, one has the hierarchy |χ 1 |E |χ 3 |E 3 |χ 5 |E 5 . However this was done in [54] and we shall now briefly review the corresponding results. 1. For very low reduced frequencies (f /f α 0.05), there is a plateau (indicated by the yellow-shaded planes in figure 10) where the reduced response X 5 depends neither on frequency nor on temperature. In this plateau, the behavior of the supercooled liquid cannot be qualitatively distinguished from the behavior expected from a high temperature liquid of dipoles, depicted by the 'trivial' X (k) k curves represented as dotted lines in figure 11. 2. At higher frequencies, we can observe a hump of |X (5) 5 | that remarkably occurs at the same peak frequency f peak as in |χ (3) 3 | in both glycerol and propylene car- bonate. Again one finds that, for the five temperatures where the peak is studied, f peak /f α = c, where the constant c does not depend on T and weakly changes with the liquid. This peak is much sharper for |X (5) 5 | than for |X (3) 3 |: this is clearly evidenced by figure 11 where the linear, cubic and fifth-order susceptibilities are compared, after normalisation to their low-frequency value. This shows that the anomalous features in the frequency dependence are stronger in |X (5) 5 | than in |X 3 |: this may be regarded as a sign of criticality since close to a critical point, the larger the order k of the response, the stronger the anomalous features of X k .
A second, and more quantitative indication of incipient criticality is obtained by studying the temperature dependence of |X (5) 5 | and by comparing it with that of |X 1. In the plateau region at f /f α 0.05, the value of |X (5) 5 | does not depend on the temperature. This shows that the factor involved in the calculation of the dimensionless X 5 , see equation (8), is extremely ecient to remove all trivial temperature dependences. As the trivial behavior depends on frequency see the dashed lines of figure 11, the 'singular' parts of X 3 and of X 5 are obtained as follows: 3,trivial , X 5,sing. ≡ X which correspond in figure 11 to a complex subtraction between the measured data symbols and the trivial behavior dashed lines. 2. Around the hump, the temperature behavior of |X (5) 5,sing. ( f peak )| is compared to that of |X 3,sing. ( f peak )| µ where µ is an exponent that is determined experimentally by looking for the best overlap of the two series of data in figure 12, see [54] for details. This leads us to values of µ = 2.2 ± 0.5 in glycerol and µ = 1.7 ± 0.4 in Figure 11. For glycerol, comparison of the fifth, third and linear susceptibilities, the latter is noted |χ (1) 1 |. The hump for |χ (5) 5 | is much stronger than that of |χ (

Testing Bouchaud-Biroli's predictions as well as the general theories of the glass transition
Having shown the experimental data for the nonlinear responses, we now move to the interpretation part and start with Bouchaud-Biroli's approach (BB), which is the most general one. The more specific and/or phenomenological approaches of nonlinear responses will be detailed in section 4.
The linear response χ 1 ≡ (∂M/∂h) h=0 is readily obtained: which shows that the linear response is related to the connected two-point correlation function. Repeating the argument for higher-order responses e.g. χ 3 ∝ (∂ 3 M/∂h 3 ) h=0 , one obtains that χ 2k+1 is connected to the (2k + 2) points correlation function, e.g.

The spin glass case.
Spin glasses are characterized by the fact that there is frozen disorder, i.e. the set of the interaction constants {J i;j } between two given spins S i and S j is fixed once and for all, and has a random sign-half of the pairs of spins are coupled ferromagnetically, the other half antiferromagnetically. Despite the fact that the system is neither a ferromagnet, nor an antiferromagnet, upon cooling it freezes, below a critical temperature T SG , into a solid long range ordered state called a spin glass state. This amorphous ordering is not detected by χ 1 which does not diverge at T SG : this is because the various terms of i1;i2 S i1 S i2 cancel since half of them are positive and the other half are negative. By contrast the cubic susceptibility χ 3 contains a term i1;i2 S i1 S i2 2 which does diverge since all its components are strictly positive: this comes from the fact that the influence S i1 S i2 of the polarization of spin S i1 on spin S i2 may be either positive or negative, but it has the same sign as the reverse influence S i2 S i1 of spin S i2 on spin S i1 . This is why the amorphous ordering is directly elicited by the divergence of the static value of χ 3 when decreasing T towards T SG , as already illustrated in figure 1(A). By adding a standard scaling assumption close to T SG one can account for the behavior of χ 3 at finite frequencies, i.e. one easily explains that χ 3 is frequency independent for ωτ α 1, and smoothly tends to zero at higher frequencies. Finally, similar scaling arguments about correlation functions easily explain the fact that the stronger k 1 the more violent the divergence of χ 2k+1 in spin glasses, as observed experimentally by Levy et al [55]. The main physical idea of BB's work [9] is that these diculties have an eect which is important at low frequencies and negligible at high enough frequencies: 1. Provided f f α , i.e. for processes faster than the relaxation time, one cannot distinguish between a truly frozen glass and a still flowing liquid. If some amorphous order is present in the glass forming system, then non-trivial spatial correlations should be present and lead to anomalously high values of non-linear susceptibilities: this holds for very general reasons (e.g. the Langevin equation for continuous spins which is used in [9] needs not to specify the detailed Hamiltonian of the system) and comes from an analysis of the most diverging term in the four terms contributing to χ 3 (ω). If the amorphous correlations extend far enough to be in the scaling regime, one can neglect the subleading terms and one predicts that the nonlinear susceptibilities are dominated by the glassy correlations and given by [9,54]: where the scaling functions H k do not explicitly depend on temperature, but depend on the kind of susceptibility that is considered, i.e. X 3 , X 3 or X 2,1 in the third order case k = 1. We emphasize that in [9] the amorphously ordered domains were assumed to be compact, i.e. d f = d, yielding α 1 = 1 i.e. X 3 ∝ N corr . The possibility of having a fractal dimension d f lower than the spatial dimension d was considered in [54] where the fifth order response was studied. As already shown in section 2.4.2, the experimental results were consistent with d f = d, i.e. X 5 ∝ N 2 corr . 2. In the low frequency regime f f α , relaxation has happened everywhere in the system, destroying amorphous order [55] 5 and the associated anomalous response to the external field and H k (0) = 0. In other words, in this very low frequency regime, every molecule behaves independently of others and X 2k+1 is dominated by the 'trivial' response of eectively independent molecules.
Due to the definition adopted in equation (8), the trivial contribution to X 2k+1 should not depend on temperature (or very weakly). Hence, provided N corr increases upon cooling, there will be a regime where the glassy contribution X glass 2k+1 should exceed the trivial contribution, leading to hump-shaped non-linear susceptibilities, peaking at f peak ∼ f α , where the scaling function H k reaches its maximum. https://doi.org/10.1088/1742-5468/ab371e J. Stat. Mech. (2019) 124003

Experimentally testing BB's predictions
We now briefly recall why all the experimental features reported in section 2 are well accounted for by BB's prediction: 1. The modulus of both the third order susceptibilities |χ 2;1 | and of |χ (5) 5 | have a humped shape in frequency, contrary to |χ 1 |.
2. Due to the fact that H k does not depend explicitly on T, the value of f peak /f α should not depend on temperature, consistent with the experimental behavior.
3. Because of the dominant role played by the glassy response for f f peak , the T-dependence of |X 2k+1 | will be much stronger above f peak than in the trivial low-frequency region.
4. Finally, because non-linear susceptibilities are expressed in terms of scaling functions, it is natural that the behavior of their modulii and phases are quantitatively related especially at high frequency where the 'trivial' contribution can be neglected, consistent with equations (9) and (10)-see below for a more quantitative argument in the context of the so-called 'Toy model 6 '.
Having shown that BB's prediction is consistent with experiments, the temperature variation of N corr can be drawn from the increase of the hump of X 3 upon cooling. It has been found [3,4,12,21,22] that the temperature dependence of N corr inferred from the height of the humps of the three X 3 's are compatible with one another, and closely related to the temperature dependence of T χ T , which was proposed in [24,25] as a simplified estimator of N corr in supercooled liquids. The convergence of these dierent estimates, that rely on general, model-free theoretical arguments, is a strong hint that the underlying physical phenomenon is indeed the growth of collective eects in glassy systems-a conclusion that will be reinforced by analyzing other approaches in section 4. Let us again emphasize that the BB prediction relies on a scaling argument, where the correlation length of amorphously ordered domains is (much) larger than the molecular size a. This naturally explains the similarities of the cubic responses in microscopically very dierent liquids such as glycerol and propylene carbonate, as well as many other liquids [12,22]. Indeed the microscopic dierences are likely to be wiped out for large ∝ N 1/d f corr , much like in usual phase transitions.

Static versus dynamic length scale? χ 3 and χ 5 as tests of the theories of the glass transition
We now shortly discuss whether N corr , as extracted from the hump of |X 3 |, must be regarded as a purely dynamical correlation volume, or as a static correlation volume. This ambiguity arises because theorems relating (in a strict sense) nonlinear responses to high-order correlation functions only exist in the static case, and that supplementary arguments are needed to interpret the humped shape of X 3 (and of X 5 ) observed

J. Stat. Mech. (2019) 124003
experimentally. In the original BB's work [9] it was clearly stated that N corr was a dynamical correlation volume since it was related to a four point time dependent correlation function. This question was revisited in [54] where it was argued that the experimental results could be accounted for only when assuming that N corr is driven by static correlations. This statement comes from an inspection of the various theories of the glass transition [54]: as we now briefly explain, only the theories where the underlying static correlation volume is driving the dynamical correlation volume are consistent with the observed features of nonlinear responses.
As a first example, the case of the family of kinetically constrained models (KCMs) [57] is especially interesting since dynamical correlations, revealed by, e.g. four-point correlation functions, exist even in the absence of a static correlation length. However in the KCM family, one does not expect any humped shape for nonlinear responses [54]. This is not the case for theories (such as RFOT [50] or Frustration theories [58]) where a non-trivial thermodynamic critical point drives the glass transition: in this case the incipient amorphous order allows to account [54] for the observed features of X 3 and X 5 . This is why it was argued in [52,54] that, in order for X 3 and X 5 to grow, some incipient amorphous order is needed, and that dynamical correlations in strongly supercooled liquids are driven by static ('point-to-set') correlations 7 -this statement will be reinforced in section 4.2.

More specific models for harmonic susceptibilities
We now review the various other approaches that have been elaborated for the nonlinear responses of glass forming liquids. We shall see that most of them, if not all, are consistent with BB's approach since they involve N corr as a key implicit or explicit parameter.

Toy and pragmatical models
The 'Toy model' has been proposed in [21,59] as a simple incarnation of the BB mechanism, while the 'Pragmatical model' is more recent [60,61]. Both models start with the same assumptions. (i) Each amorphously ordered domain is compact and contains N corr molecules: because amorphous order amounts to drawing the orientation of the N corr dipoles at random in space, one gets, as explained in section 1, a net dipole moment ∝ √ N corr . This leads to an anomalous contribution to the cubic response X glass 3 ∝ N corr . (ii) There is a crossover at low frequencies towards a trivial cubic susceptibility contribution X triv 3 which does not depend on N corr . More precisely, in the 'Toy model' each amorphously ordered domain is supposed to live in a simplified energy landscape, namely an asymmetric double-well potential with a dimensionless asymmetry δ, favoring one well over the other. The most important dierence between the Toy and the Pragmatical model comes from the description of the low-frequency crossover, see [59] and [61] for more details. https://doi.org/10.1088/1742-5468/ab371e

J. Stat. Mech. (2019) 124003
On top of N corr and δ, the Toy model uses a third adjustable parameter, namely the frequency f * below which the trivial contribution becomes dominant. In [59], both the modulus and the phase of X 3 (ω, T ) in glycerol were well fitted by using f * f α /7, δ = 0.6 and, for T = 204 K, N corr = 5 for X (3) 3 and N corr = 15 for X (1) 3 . Figure 13 gives an example of the Toy model prediction for X (3) 3 in glycerol. Besides, in [21], the behavior of X (1) 2,1 (ω, T ) in glycerol was fitted with the same values of δ and of f * but with N corr = 10 (at a slightly dierent temperature T = 202 K). Of course, the fact that a dierent value of N corr must be used for the three cubic susceptibilities reveals that the Toy model is oversimplified, as expected. However, keeping in mind that the precise value of N corr does not change the behavior of the phases, we note that the fit of the three experimental phases is achieved [21,59] by using the very same values of f * /f α and of δ. This means that equations (9) and (10) are well accounted for by the Toy model by choosing two free parameters. This is a quantitative illustration of how the BB general framework does indeed lead to strong relations between the various non-linear susceptibilities, such as those contained in equations (9) and (10).
Let us mention briefly the asymmetric double well potential (ADWP) model [62], which is also about species living in a double well of asymmetry energy ∆, excepted that two key assumptions of the Toy and Pragmatical models are not made: the value of N corr is not introduced, and the crossover to trivial cubic response is not enforced at low frequencies. As a result, the hump for |X (3) 3 | is predicted [62,63] only when the reduced asymmetry δ = tanh(∆/(2k B T )) is close to a very specific value, namely δ c = 1/3, where X 3 vanishes at zero frequency due to the compensation of its several terms. However, at the fifth order [63] this compensation happens for two values of δ very dierent from δ c : as a result the model cannot predict a hump happening both for the third and for the fifth order in the same parametric regime, contrarily to the  [54]. This very recent calculation of fifth order susceptibility [63] reinforces the point of view of the Toy and Pragmatical models, which do predict a hump occurring at the same frequency and temperature due to their two key assumptions (N corr and crossover to trivial nonlinear responses at low frequencies). This can be understood qualitatively: because the Toy model predicts [59] an anomalous contrib ution X glass 2k+1 ∼ [N corr ] k , provided that N corr is large enough, the magnitude of this contrib ution is much larger than that of the small trivial contribution X triv. 2k+1 ∼ 1, and the left side of the peak of |X 2k+1 | arises just because the Toy model enforces a crossover from the large anomalous response to the small trivial response at low frequencies f f α . As for the right side of the peak, it comes from the fact that |X 2k+1 | must decrease when raising the frequency far above f α for a very general reason: the higher the frequency, the more dicult it is for the supercooled liquid to follow the alternation of the field. During its natural time scale τ α , the field reverses back and forth many times, which eectively diminishes its eect on the system. This diminution reinforces as the frequency increases, which naturally explains why the magnitude of all the responses at all orders decrease with increasing f in the range f f α . In practice, other degrees of fredom, not involved in the glass transition, may also add up to the measured di electric response: this why, e.g. the real part of the linear susceptibility χ 1,∞ at very high frequency does not vanish but instead tends to n 2 opt − 1 where n opt is the optical index.

Entropic eects
A contribution to nonlinear responses was recently calculated by Johari in [64,65] in the case where a static field E st drives the supercooled liquid in the nonlinear regime. Johari's idea was positively tested in the corresponding χ (1) 2;1 experiments in [66][67][68][69], see however [70] for a case where the agreement is not as good. It was then extended to pure ac experiments, and thus to χ (3) 3 , in [71,72]. The relation between Johari's idea and N corr was made in [52].

When a static field E st is applied.
Let us start with the case of χ (1) 2;1 experiments, i.e. with the case where a static field E st is superimposed onto an ac field E cos(ωt). In this case, there is a well defined variation of entropy [δS] Est induced by E st , which, for small E st and a fixed T, is given by: where a 3 is the molecular volume. Equation (14) holds generically for any material. However, in the specific case of supercooled liquids close enough to their glass trans ition temperature T g , a special relation exists between the molecular relaxation time τ α and the configurational contribution to the entropy S c . This relation, first anticipated by Adam and Gibbs [23], can be written as: where τ 0 is a microscopic time, and ∆ 0 is an eective energy barrier for a molecule. The temperature dependence of TS c (T) quite well captures the temperature variation of ln(τ α ), at least for a large class of supercooled liquids [73]. Following Johari [64,65] let us now assume that [δS] Est is dominated by the dependence of S c on field,-see the appendix of [52] for a further discussion of this important physical assumption. Combining equations (14) and (15), one finds that a static field E st produces a shift of ln(τ α /τ 0 ) given by: As shown in [52] this entropic eect gives a contribution to X 2,1 , which we call J 2,1 after Johari. Introducing x = ωτ α , the most general expression of J (1) 2:1 reads: where χ lin is the complex linear susceptibility. Equation (17) deserves three comments: 2,1 | has a humped shaped in frequency with a maximum in the region of ωτ α 1, because of the frequency dependence of the factor ∝ ∂χ lin /∂ ln x in equation (17). 3. The smaller S c , the larger must be the size of the amorphously ordered domains -in the hypothetical limit where S c would vanish, the whole sample would be trapped in a single amorphously ordered sate and N corr would diverge. In other words, there is a relation between S −1 c and N corr , which yields [52]:

The temperature variation of J
where it was in shown in [52] that: a. the exponent q lies in the [2/3; 2] interval when one combines the Adam-Gibbs original argument with general constraints about boundary conditions [52].
b. the exponent q lies in the [1/3; 3/2] interval [52] when one uses the RFOT and plays with its two critical exponents Ψ and θ. Notably, taking the 'recommended RFOT values' Ψ = θ = 3/2 for d = 3 gives q = 1, which precisely corresponds to BB's prediction. In this case, entropic effects are a physically motivated picture of BB's mechanism-see [52] for a refined discussion.

4.2.2.
When a pure ac field E cos(ωt) is applied. Motivated by several works [66][67][68][69] showing that Johari's reduction of entropy fairly well captures the measured χ 3 and χ 3 . This has given rise to the phenomenological model elaborated in [71,72] where the entropy reduction depends on time, which is nevertheless acceptable in the region ωτ α 1 where the model is used. Figure 14 shows the calculated values for |χ (3) 3 | at three temperatures for glycerol. The calculation fairly well reproduces the hump of the modulus observed experimentally (the phase has not been calculated). As very clearly explained in [72], the hump displayed in figure 14 comes directly from the entropic contribution and not from the two other contributions included in the model (namely the 'trivial' or 'saturation' contribution, and the Box model contribution see section 4.3 below).
Summarizing this section about entropy eects, we remind the two main assumptions made by Johari: (i) the field-induced entropy variation mainly goes into the configurational part of the entropy; (ii) its eects can be calculated by using the Adam-Gibbs relation. Once combined, these two assumptions give a contribution to χ (1) 2;1 reasonably well in agreement with the measured values in several liquids [66][67][68][69]. An extension to χ 3 is even possible, at least in the region ωτ α 1 and fairly well accounts for the measured hump of |χ (3) 3 | in glycerol [71,72], a figure similar to figure 11 for |χ (5) 5 (ω)/χ (5) 5 (0)| is even obtained in [72]. As shown in equation (18), this entropy contrib ution to cubic responses is related to N corr , which is consistent with the general prediction of BB. Additionally, because S c is a static quantity, equation (18) supports the interpretation that the various cubic susceptibilities χ 3 are related to static amorphous correlations, as discussed in section 3.3. Figure 14. The model elaborated in [71,72] includes three contributions entropyreduction, Box model, and trivial. It predicts for |χ (3) 3 | the solid lines which account very well for the measured values in glycerol in frequency and in temperature. The peak of |χ (3) 3 | arises because of the entropy reduction eect (noticed 'sing. T fic. ') which completely dominates the two other contributions in the peak region, as shown by the inset.  [74]. When these pionneering experiments were carried out, a central question was whether the dynamics in supercooled liquids is homogeneous or heterogeneous. In the seminal [74] it was reported that when applying a strong ac field E of angular frequency ω, the changes in the dielectric spectrum are localised close to ω and that they last a time of the order of 1/ω. These two findings yield a strong qualitative support to the heterogenous character of the dynamics, and the Box model was designed to provide a quantitative description of these results. Accordingly, the Box model assumes that the dielectric response comes from 'domains' that will be later called dynamical heterogeneities (DH), each domain being characterized by its dielectric relaxation time τ and obeying the Debye dynamics. The distribution of the various τ 's is chosen to recover the measured non Debye spectrum by adding the various linear Debye susceptibilities χ 1,dh = ∆χ 1 /(1 − iωτ ) of the various domains. For the nonlinear response, the Box model assumes that it is given by the Debye linear equation in which τ (T ) is replaced by τ (T f ) where the fictive temperature T f = T + δT f is governed by the constitutive equation, see e.g. [32,78]: with c dh the volumic specific heat of the DH under consideration, κ the thermal conductance (divided by the DH volume v) between the DH and the phonon bath, τ therm = c/κ the corresponding thermal relaxation time. In equation (19), only the constant part of the dissipated power has been written, omitting its component at 2ω which is important only for χ 3 , see e.g. [78]. From equation (19) one easily finds the stationary value δT f of δT f which reads: As very clearly stated in the seminal [74] because the DH size is smaller than 5 nm, the typical value of τ therm is at most in the nanoseconds range: this yields, close to T g , a vanishingly small value of τ therm /τ, which, because of equation (20), gives fully negligible values for δT f . The choice of the Box model is to increase τ therm by orders of magnitude by setting τ therm = τ, expanding onto the intuition that this is a way to model the 'energy storage' in the domains. The main justification of this choice is its eciency: it allows to account reasonably well for the NHB experiments [74] and thus to bring a strong support to the heterogeneous character of the dynamics in supercooled liquids. Since the seminal [74], some other works have shown [32,[75][76][77]] that the Box model eciently accounts for the measured χ 3 ( f > f α ) in many glass forming liquids. It was shown also [78] that the Box model is not able to fit quantitatively the measured X (even though some qualitative features are accounted for), and that the Box model only provides a vanishing contribution to X (1) 2,1 -see [21].
The key choice τ therm = τ made by the Box model has two important consequences for cubic susceptibilities: it implies (a) that χ (since the source term in equation (19) is the dissipated power) and (b) that χ (1) 3 does not explicitly depend on the volume v = N corr a 3 of the DH's (see [32,78]). However, alternative models of nonlinear responses are now available [59,61] where, instead of choosing τ therm , one directly resolves the microscopic population equations, which is a molecular physics approach, and not a macroscopic law transferred to microscopics. The population equations approach is equivalent to solving the relevant multidimensional Fokker-Planck equation describing the collective tumbling dynamics of the system at times longer than the time between two molecular collisions (called τ c in appendix C). By using this molecular physics approach one obtains that χ (1) 3 is governed by N corr and not by energy absorption. For χ (1) 3 we summarize the argument given in the end of the appendix of [52]: writing loosely P (1) 3 ≈ ∂P 1 /(∂ ln τ )δ ln τ , one sees that the pivotal quantity is the field induced shift of the relaxation time δ ln τ (more precisely δ ln τ is the only contribution contained in the Box model, while the Toy model generally contains also another contribution). Comparing the Box model (BM) and, e.g. the Toy model (TM), one gets respectively: where we remind our definition χ T = |∂ ln τ α /∂T | and where the limit ωτ 1, relevant for the χ (1) 1. one sees that the two values of δ ln τ are similar provided N corr and T χ T are proportional-which is a reasonable assumption as explained above and in [4,24,25]. Taking reasonable values of this proportionality factor, it was shown in [52] that χ (1) 3 ( f > f α ) is the same in the two models. This sheds a new light on the eciency of the Box model and on consequence (b). 2. Let us shortly discuss consequence (a). In the Toy model, δ ln τ directly expresses the field induced modification of the energy of each of the two wells modeling a given DH. It comes from the work produced by E onto the DH and this is why it involves N corr : the larger this number, the larger the work produced by the field because the net dipole of a DH is ∝ √ N corr and thus increases with N corr . It is easy to show that the dissipation, i.e. the 'energy absorption', is not involved in δ ln τ because dissipation depends only on χ 1 , which in the Toy model does not depend on N corr . In the Toy model, as in the Pragmatical model [61] and the Diezemann model [62], the heating is neglected because at the scale of a given DH it is vanishingly small as shown above when discussing τ therm . Of course, at the scale of the whole sample, some global heating arises for thick samples and/or high frequencies because the dissipated power has to travel to the electrodes which are the actual heat sinks in dielectric experiments [13]. This purely exogeneous eect can be precisely calculated by solving the heat propagation equation, see e.g. [13] and appendix B, and must not be confused with what was discussed in this section.

4.3.2.
Gathering the three measured cubic susceptibilities. As explained above, in [71,72], the three experimental cubic susceptibilities have been argued to result from a superposition of an entropic contribution and of an energy absorption contribution coming from the Box model (plus a trivial contribution playing a minor role around the peaks of the cubic susceptibilities). More precisely, the hump of |X (1) 2,1 | and of |X (3) 3 | would be mainly due to the entropy eect, contrarily to the hump of |X (1) 3 | which would be due to the Box model contribution. As noted in [52], this means that very dierent physical mechanisms would conspire to give contributions of the same order of magnitude, with phases that have no reason to match as they do empirically, see equations (9) and (10): why should X (1) 3 and X 3 have the same phase at high frequencies if their physical origin is dierent? This is why it was emphasized in [52] that there is no reason for such a similarity if the growth of X (1) 3 and X (3) 3 are due to independent mechanisms. Because entropic eects have been related to the increase of N corr , see equations (17) and (18), everything becomes instead very natural if the Box model is recasted in a framework where X (1) 3 is related to the glassy correlation volume. As evoked above, a first step in this direction was done in [52] where it was shown that the Box model prediction for X (1) 3 at high frequencies is identical to the above Toy model prediction, provided N corr and T χ T are proportional. In all, it is argued in [52] that the only reasonable way to account for the similarity of all three cubic susceptibilities, demonstrated experimentally in figures 8 and 9, is to invoke a common physical mechanism. As all the other existing approaches, previously reviewed, relate cubic responses to the growth of the glassy correlation volume, reformulating the Box model along the same line seems to be a necessity.

Conclusions
We have reviewed in this paper the salient features reported for the third and fifth harmonic susceptibilities close to the glass transition. This is a three decades long story, which has started in the mid-eighties as a decisive tool to evidence the solid, long range ordered, nature of the spin glass phase. The question of whether this notion of 'amorphous order' was just a curiosity restricted to the somehow exotic case of spin glasses remained mostly theoretical until the seminal work of Bouchaud and Biroli in 2005. This work took a lot from the spin glass physics, and by taking into account the necessary modifications relevant for glass forming liquids, it has anticipated all the salient features discovered in the last decade for the three cubic susceptibilities X 3 . This is why, in most of the works, the increase of the hump of X 3 upon cooling has been interpreted as reflecting that of the glassy correlation volume. Challenging alternative and more specific interpretations have been proposed, but we have seen that most-if not all of them-can be recasted into the framework of BB. The avenue opened by BB's prediction was also used to circumvent the issue of exponentially long time scales which are the reason why the nature of the glass transition is still debated: this is how the idea of comparing the anomalous features of X 3 and of X 5 has arisen. The exper imental findings are finally consistent with the existence of an underlying J. Stat. Mech. (2019) 124003 thermodynamic critical point, which drives the formation of amorphously ordered compact domains, the size of which increases upon cooling. Last we note that this field of nonlinear responses in supercooled liquids has been inspiring both theoretically [5,79] and experimentally, e.g. for colloidal glasses: the very recent experiments [6] have shed a new light on the colloidal glass transition and shown interesting dierences with glass forming liquids.
All these progresses open several routes of research. On the purely theoretical side, any prediction of nonlinear responses in one of the models belonging to the kin etically constrained model family will be extremely welcome to go beyond the general arguments given in [54]. Moreover, it would be very interesting to access χ 3 (and χ 5 ) in molecular liquids at higher temperatures, closer to the mode coupling trans ition temper ature T MCT , and/or for frequencies close to the fast β process where more complex, fractal structures with d f < d may be anticipated [80,81]. This will require a joined eort of experimentalists to avoid heating issues and of theorists to elicit the nature of nonlinear responses close to T MCT . Additionally, one could revisit the vast field of polymers by monitoring their nonlinear responses, which should shed new light onto the temperature evolution of the correlations in these systems. Therefore there is likely much room to deepen our understanding of the glass transition by carrying out new experiments about nonlinear susceptibilities.

J. Stat. Mech. (2019) 124003
2. The contribution of electrostriction was demonstrated to be safely negligible in [4,75], both by using theoretical estimates and by showing that changing the geometry of spacers does not aect X 3 . 3. As for the small ionic impurities present in most of liquids, we briefly explain that they have a negligible role, except at zero frequency where the ion contribution might explain why the three X 3 's are not strictly equal, contrarily to what is expected on general grounds see, e.g. figures 8 and 9. On the one hand it was shown that the ion heating contribution is fully negligible in X 2,1 (see [21]), on the other hand it is well known that ions aect the linear response χ 1 at very low frequencies (say f /f α 0.05): this yields an upturn on the out-of-phase linear response χ 1 , which diverges as 1/ω instead of vanishing as ω in an ideally pure liquid containing only molecular dipoles. This may be the reason why most of the χ 3 measurements are reported above 0.01f α : at lower frequencies the nonlinear responses is likely to be dominated by the ionic contribution.

Appendix B. Trivial third and fifth harmonic susceptibilities
As explained in the main text, in the long time limit i.e. for f /f α 1, the liquid flow destroys the glassy correlations, making each molecule eectively independent of others. This is why we briefly recall what the nonlinear responses of an ideal gas of dipoles are, where each dipole is independent of others, and undergoes a Brownian rotational motion of characteristic time τ D due to the underlying thermal reservoir at temperature T. The linear susceptibility of such an ideal gas of dipoles is given by the Debye susceptibility ∆χ 1 /(1 − iωτ D ), hence the subscript 'Debye' in the equation (B.1) below. By using [19], and following the definitions given in the main text, as well as equations (5)-(8) above, one gets for the dimensionless nonlinear responses of such an ideal dipolar gas, setting for brevity x = ωτ D : In [21,59] a simplified calculation was used for the cubic trivial susceptibilities for which G (τ ) was replaced by the Dirac delta function δ(τ − τ α ), i.e. τ D was simply replaced by τ α . A more consistent calculation, yielding slightly dierent results, was made in [54] where the trivial response at any order was obtained by combining the above X (k) k,Debye with a distribution G (τ ) of relaxation times τ chosen to account for the linear susceptibility of the supercooled liquid of interest. This was done because in supercooled liquids, the frequency dependence of χ 1 is qualitatively similar to that of a dipolar gas: therefore, one can always account for χ 1 (ω) as resulting from the mere addition J. Stat. Mech. (2019) 124003 of trivial Debye linear responses. This is not the case for nonlinear responses X 2m+1>1 because of the qualitative dierence arising in the frequency dependence of |X 2m+1>1 | in supercooled liquids and in the Debye case. This does not mean that the linear response is trivial in all respects: for example its characteristic frequency f α = 1/(2πτ α ) has a super-Arrhenius behavior which has no counterpart in the Debye case.

Appendix C. Derivation of the Toy model from Langevin Fokker-Planck considerations
In this section we shall rederive the phenomenological Toy model of Ladieu et al [59] starting from the Langevin-Fokker-Planck equation, which is the starting point of Bouchaud and Biroli when they illustrate their general theoretical ideas in the last part of [9]. We shall idealize the supercooled state of a liquid as follows. At high temperatures, the liquid is made of molecules the interactions between which are completely negligible. On cooling, the molecules arrange themselves in groups, called 'dynamical heterogeneities' (DH), between which there are no interactions. Inside a typical group, specific intermolecular interactions manifest themselves dynamically, by which we mean that in a time larger than a characteristic time τ α , such interactions lose their coherence and the typical behavior of the liquid is that of an ideal dipolar gas. Before and around τ α , these interactions manifest themselves in a frequency range ω ≈ 1/τ α . Thus, stricto sensu, our modelling of this specific process pertains to the behavior of the various dielectric responses of a DH, linear and nonlinear, near this frequency range. This indeed implies that information regarding the 'ideal dipolar gas' phase must be added to fit experimental data. It may be shown on fairly general grounds that either for linear and non-linear responses, such extra information simply superposes onto the specific behavior that has been alluded to above [82]. Now, we consider that (a) a given DH has a given size at temperature T , (b) that a DH is made of certain mobile elements that do interact between themselves, (c) that there are no interactions between DHs, (d) that the dipole moment of a DH is µ d = µ √ N corr , and (e) that all constituents of a DH are subjected to Brownian motion.
In order to translate the above assumptions in mathematical language, we assign to each constituent of a DH a generalized coordinate q i (t), so that each DH is described by a set of generalized coordinates q at temperature T , viz.
Inside each DH, each elementary constituent is assumed to interact via a multidimensional interaction potential V int (q) that possesses a double-well structure with minima at q A and q B , and are sensitive both to external stresses and thermal agitation. The equations of motion may be described by overdamped Langevin equations with additive noise, viz.
where ζ is a generalized friction coecient, V T = V int + V ext , V ext is the potential energy of externally applied forces and the generalized forces Ξ i (t) have Gaussian white noise properties, namely Thus, the dynamics of a DH is represented by the stochastic dierential equations (C.1) and (C.2), which are in eect the starting point of the Bouchaud-Biroli theory, as stated above. A totally equivalent representation of these stochastic dynamics is obtained by writing down the Fokker-Planck equation [83] for the probability density W (q, t) to find the system in state q at time t which corresponds to equations (C.1) and (C.2), namely where 2τ c = ζ/ (kT ) is the characteristic time of fluctuations, ∇ is the del operator in q space, and L FP (q, t) is the Fokker-Planck operator. We notice that equation (C.3) may also be written ∂W ∂t (q, t) = 1 2τ c ∇ · e −βV T (q,t) ∇ W (q, t) e βV T (q,t) .

(C.4)
Now we use the transformation [84] φ (q, t) = W (q, t) e βV T (q,t) (C.5) so that equation (C.4) becomes ∂φ ∂t (q, t) − β ∂V T ∂t (q, t) φ (q, t) = 1 2τ c e βV T (q,t) ∇ · e −βV T (q,t) ∇φ (q, t) where L † FP (q, t) is the adjoint Fokker-Planck operator [83]. Next, we make the first approximation in our derivation, namely, we assume that the time variation of V T is small with respect to that of W. If the time dependence of V T is contained in, say, the application of a time-varying uniform AC field only, this implies immediately that neglecting the second term in the left hand side of equation (C.6) means that W is near its equilibrium value, so restricting further calculations to low frequencies, ωτ c << 1 (quasi-stationary condition). Hence, equation (C.6) now reads Now, the interpretation of the Fokker-Planck equation (C.3) (or equally well the Langevin equation (C.1)) with time-dependent potential in terms of usual population equations with time-dependent rate coecients has a meaning, since now equation (C.5) means detailed balancing. The polarization of an assembly of noninteracting DH in the direction of the applied field may then be defined as P (t) = ρ 0 µ d cos ϑ (q) W (q, t) dq (C.8)

J. Stat. Mech. (2019) 124003
where ρ 0 is the number of DH per unit volume, and ϑ (q) is the angle a DH dipole makes with the externally applied electric field. Because of the double-well structure of the interaction potential, we may equally well write equation (C.8) cos ϑ (q) W (q, t) dq . (C.9) Now, it is known from the Kramers theory of chemical reaction rates [84] that at suciently large energy barriers, most of the contributions of the integrands come from the minima of the wells, therefore we have Now, the integrals represent the relative populations x i (t) = n i (t) /N , i = A, B in each well (we assume that W (q, t) is normalized to unity), where n i (t) is the number of DH states in well i, and N the total number of DH. At any time t, we have the conservation law x A (t) + x B (t) = 1. (C.11) Thus, equation (C.10) reads We assume now for simplicity that ϑ (q B ) = π − ϑ (q A ), so that Finally, since ρ 0 = N/V where V is the volume of the polar substance made of DH only, we obtain where υ DH is the volume of a DH. This is the definition of the polarization in the Toy model. In order to determine the polarization (C.14), we need to calculate the dynamics of n i (t). From the conservation law equation (C.11), we havė By using the Fokker-Planck equation (C.4) and limiting well A to a closed generalized bounding surface constituting the saddle region ∂A, we have by Gauss's theoreṁ x A (t) = −ẋ B (t) = 1 2τ c ∂A e −βV T (q,t) ∇φ (q, t) · ν q dS q (C. 16) where ν q is the outward normal to the bounding surface and dS q is a generalized surface element of the bounding surface, and where we have used equation (C.5). Now, we follow closely Coey et al [85] and introduce the crossover function ∆ (q, t) via the equation where ∆ (q, t) = 0 if q ∈ well A while ∆ (q, t) = 1 if q ∈ well B and exhibits strong gradients in the saddle region ∂A allowing the crossing from A to B (and vice-versa) by thermally activated escape. By combining equations (C. 16) and (C.17), we have immediatelyẋ Now, where W s (q, t) is a normalized solution of the Fokker-Planck equation because the frequencies we are concerned with are very small with respect to the inverse thermal fluctuation time τ c and because the time-dependent part of the potential V T is much smaller than other terms in it at any time. We have Using equations (C. 19) and (C.21), we may easily show that [85] φ B (t) − φ A (t) = 1 x s A (t) The obtaining of a more explicit formula for the various rates involved in equation (C.25) is not possible, due to the impossibility to calculate the surface integral in equation (C.24) explicitly, in turn due to the fact that V T is not known explicitly. Then, the rates in equations (C. 25) and (C.26) are estimated using Arrhenius's formula. All subsequent derivations regarding the Toy model of Ladieu et al [59] follow immediately and will not be repeated here due to lack of room and straightforward but laborious algebra.

J. Stat. Mech. (2019) 124003
Appendix D. Derivation of equation (5) As an example of the use of the nonlinear formalism, we derive equation (5).