# Third and fifth harmonic responses in viscous liquids

Samuel Albert, M. Michl, P. Lunkenheimer, A. Loidl, P. M Déjardin, F.

Ladieu

► To cite this version:

Samuel Albert, M. Michl, P. Lunkenheimer, A. Loidl, P. M Déjardin, et al. Third and fifth harmonic responses in viscous liquids. R. Richert. Nonlinear Dielectric Spectroscopy, Springer, 2018, Advances in Dielectrics, 978-3-319-77574-6. 10.1007/978-3-319-77574-6_7. hal-01826296

HAL Id: hal-01826296

https://hal.science/hal-01826296

Submitted on 5 Jul 2018

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Third and fifth harmonic responses in viscous liquids

S. Albert, M. Michl, P. Lunkenheimer, A. Loidl, P. M. Déjardin, and F. Ladieu

## 1 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 Fig. 1 displaying the third harmonic

S. Albert

SPEC, CEA, CNRS, Université Paris-Saclay, CEA Saclay Bat 772, 91191 Gif-sur-Yvette Cedex, France. e-mail: samuel.albert@ens-lyon.fr

M. Michl

Experimental Physics V, Center for Electronic Correlations and Magnetism, University of Augsburg, 86159 Augsburg, Germany. e-mail: marion.michl@physik.uni-augsburg.de

P. Lunkenheimer

Experimental Physics V, Center for Electronic Correlations and Magnetism, University of Augsburg, 86159 Augsburg, Germany. e-mail: peter.lunkenheimer@physik.uni-augsburg.de

A. Loidl

Experimental Physics V, Center for Electronic Correlations and Magnetism, University of Augsburg, 86159 Augsburg, Germany. e-mail: alois.loidl@physik.uni-augsburg.de

P. M. Déjardin

LAMPS Université de Perpignan Via Domitia - 52 avenue Paul Alduy - 66860 Perpignan Cedex, France. e-mail: dejardip@univ-perp.fr

F. Ladieu

SPEC, CEA, CNRS, Université Paris-Saclay, CEA Saclay Bat 772, 91191 Gif-sur-Yvette Cedex, France. e-mail: francois.ladieu@cea.fr
<figure><img src="image_2.png" /><figcaption>**Fig. 1** From Refs. [1, 2, 3, 6]. Third-harmonic susceptibilities of very different types of glasses approaching their glass transition. (A): In the Ag:Mn spin glass [1], the static value of $\chi_3^{(3)}$ diverges when approaching the critical temperature $T_c \approx 2.94$ K [1]. (B): Similar arguments are used to rationalize the third-harmonic dielectric susceptibility of an orientational glass [2]. (C): In glycerol [3, 4], the modulus of the -dimensionless- cubic susceptibility $\chi_3^{(3)}$ has a peak as function of frequency, which increases anomalously upon cooling. (D): Strain-stress experiment in the colloidal system studied in Refs. [5, 6]. When increasing the volumic density $\phi$ the increasing peak of $Q_0 = |\chi_3^{(3)}|/\chi_1$ reveals that any shear strain connects the system to a non equilibrium steady state -see [5, 6]-. In all these four examples $\chi_3^{(3)}$ unveils informations about the nature of the glassy state that cannot be obtained by studying the linear susceptibility $\chi_{lin}$.</figcaption></figure>

cubic susceptibility $\chi_3^{(3)}$ -defined in Section 2.1- for four very different kinds of glasses [1, 2, 3, 4, 5, 6]. In the case of spin glasses [1, 7] -see Fig. 1(A)-, it was discovered in the eighties that $\chi_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 Fig. 1(B), the divergence of $|\chi_3^{(3)}|$ is not accompanied by any divergence of the linear susceptibility $|\chi_1|$.
We shall show in Eqs. (1)-(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 $|\chi_3^{(3)}(\omega, T)|$ has a hump close to the $\alpha$ relaxation frequency $f_\alpha$, the height of this hump increasing anomalously upon cooling. A hump of $|\chi_3^{(3)}|$ has also been recently discovered in a colloidal glass [5, 6], in the vicinity of the $\beta$ relaxation frequency $f_\beta$, 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 different kinds of glasses of Fig. 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 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 $\mu_{dip}$. At high temperature $T$ the system behaves as an ideal gas and its polarization $P$ is given by:

$$P = \frac{\mu_{dip}}{a^d} \mathcal{L}_d \left( \frac{\mu_{dip} E_{st}}{k_B T} \right) \\ \approx \frac{1}{3} \frac{\mu_{dip}}{a^d} \left( \frac{\mu_{dip} E_{st}}{k_B T} \right) - \frac{1}{45} \frac{\mu_{dip}}{a^d} \left( \frac{\mu_{dip} E_{st}}{k_B T} \right)^3 + \frac{2}{945} \frac{\mu_{dip}}{a^d} \left( \frac{\mu_{dip} E_{st}}{k_B T} \right)^5 + \dots \quad (1)$$

where $a^d$ is the molecular $d$-dimensional volume, $\mathcal{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 $\ell$, i.e. molecules are correlated within groups containing $N_{corr} = (\ell/a)^{df}$ molecules, with $df$ the fractal dimension characterizing the correlated regions. Because these domains are independent from each other, one can use Eq. (1), provided that we change the elementary volume $a^d$ by that of a domain - namely $a^d(\ell/a)^d$, as well as the molecular dipole $\mu_{dip}$ by that of a domain -namely $\mu_{dip}(\ell/a)^{(d_f/2)}$. Here, the exponent $df/2$ expresses the amorphous ordering within the correlated regions, i.e. the fact that the orientation of the correlated molecules looks random in space. We obtain:

$$\frac{P}{\mu_{dip}/a^d} \approx \frac{1}{3} \left(\frac{\ell}{a}\right)^{d_f-d} \left(\frac{\mu_{dip}E_{st}}{k_B T}\right) - \frac{1}{45} \left(\frac{\ell}{a}\right)^{2d_f-d} \left(\frac{\mu_{dip}E_{st}}{k_B T}\right)^3 + \\ + \frac{2}{945} \left(\frac{\ell}{a}\right)^{3d_f-d} \left(\frac{\mu_{dip}E_{st}}{k_B T}\right)^5 + \dots \quad (2)$$

which shows that the larger the order $k$ of the response, the stronger the increase of the response when $\ell$ increases. As $d_f \le d$, Eq. (2) shows that the linear response
never diverges with $\ell$: it is always, for any $\ell$, of the order of $\mu_{dip}^2/(a^d k_B T)$. This can be seen directly in Eq. (2) in the case $d_f = d$; while for $d_f < d$ one must add to Eq. (2) the polarization arising from the uncorrelated molecules not belonging to any correlated region. This insensitivity of the linear response to $\ell$ 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 $\ell$ 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 $\ell$, 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.

For spin-glasses, the above purely thermodynamic argument is enough to relate the divergence of the static value of $\chi_3^{(3)}$ -see Fig 1-(A)- to the divergence of the amorphous order correlation length $\ell$. For structural glasses this argument must be complemented by some dynamical argument, since we have seen on Fig. 1-(C) that the anomalous behavior of $\chi_3^{(3)}$ takes place around the relaxation frequency $f_\alpha$. This has been done, on very general grounds, by the predictions of Bouchaud and Biroli (BB), who anticipated [9] the main features reported in Fig. 1-(C). BB's predictions will be explained in Section 3. Before, we shall review in Section 2 the main experimental features of third and fifth harmonic susceptibilities. Because of the generality of Eq. (2) and of BB's framework, we anticipate that $\chi_3$ and $\chi_5$ have common anomalous features that can be interpreted as reflecting the evolution of $\ell$ -and thus of $N_{corr}$- upon cooling. The end of the chapter, Section 4, will be devoted to more specific approaches of the cubic response of glass forming liquids. Beyond their apparent diversity, we shall show that they can be unified by the fact that in all of them, $N_{corr}$ is a key parameter -even though it is sometimes implicit-. The Appendix contains some additional material for the readers aiming at deepening their understanding of this field of high harmonic responses.

## 2 Experimental behavior of third and fifth harmonic susceptibilities

### 2.1 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 :

$$P(t) = \sum_{m=0}^{\infty} P_{2m+1}(t) \qquad (3)$$
where because of the $E \to -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:

$$ \frac{P_{2m+1}(t)}{\varepsilon_0} = \int_{-\infty}^{\infty} \dots \int_{-\infty}^{\infty} \chi_{2m+1}(t-t'_1, \dots, t-t'_{2m+1}) E(t'_1) \dots E(t'_{2m+1}) dt'_1 \dots dt'_{2m+1} \quad (4) $$

Because of causality $\chi_{2m+1} \equiv 0$ whenever one of its arguments is negative. For a field $E(t) = E \cos(\omega t)$ of frequency $\omega$ and of amplitude $E$, it is convenient to replace $\chi_{2m+1}$ by its $(2m+1)$-fold Fourier transform and to integrate first over $t'_1, \dots, t'_{2m+1}$. Defining the one-fold Fourier transform $\phi(\omega)$ of any function $\phi(t)$ by $\phi(\omega) = \int \phi(t)e^{-i\omega t} dt$ (with $i^2 = -1$) and using $\int e^{-i(\omega_1-\omega)t} dt = 2\pi\delta(\omega_1-\omega)$, where $\delta$ is the Dirac delta function, one obtains the expression of $P_{2m+1}(t)$. This expression can be simplified by using two properties: (a) the fact that the various frequencies $\omega_\lambda$ play the same role, which implies $\chi_{2m+1}(-\omega, \omega, \dots, \omega) = \chi_{2m+1}(\omega, -\omega, \dots, \omega)$; (b) the fact that $\chi_{2m+1}$ is real in the time domain implying that $\chi_{2m+1}(-\omega, \dots, -\omega)$ is the complex conjugate of $\chi_{2m+1}(\omega, \dots, \omega)$. By using these two properties, we obtain the expression of all the $P_{2m+1}(t)$, and in the case of the third order polarisation this yields:

$$ \frac{P_3(t)}{\varepsilon_0} = \frac{1}{4} E^3 |\chi_3^{(3)}(\omega)| \cos(3\omega t - \delta_3^{(3)}(\omega)) + \frac{3}{4} E^3 |\chi_3^{(1)}(\omega)| \cos(\omega t - \delta_3^{(1)}(\omega)) \quad (5) $$

where we have set $\chi_3(\omega, \omega, \omega) = |\chi_3^{(3)}(\omega)|e^{-i\delta_3^{(3)}(\omega)}$, and $\chi_3(\omega, \omega, -\omega) = |\chi_3^{(1)}(\omega)|e^{-i\delta_3^{(1)}(\omega)}$. Similarly, for the fifth-order polarisation, we obtain:

$$ \frac{P_5(t)}{\varepsilon_0} = \frac{1}{16} E^5 |\chi_5^{(5)}(\omega)| \cos(5\omega t - \delta_5^{(5)}(\omega)) + \frac{5}{16} E^5 |\chi_5^{(3)}(\omega)| \cos(3\omega t - \delta_5^{(3)}(\omega)) + \\ + \frac{10}{16} E^5 |\chi_5^{(1)}(\omega)| \cos(\omega t - \delta_5^{(1)}(\omega)) \quad (6) $$

where, we have set $\chi_5(\omega, \omega, \omega, \omega, \omega) = |\chi_5^{(5)}(\omega)|e^{-i\delta_5^{(5)}(\omega)}$, and similarly $\chi_5(\omega, \omega, \omega, \omega, -\omega) = |\chi_5^{(3)}(\omega)|e^{-i\delta_5^{(3)}(\omega)}$ as well as $\chi_5(\omega, \omega, \omega, -\omega, -\omega) = |\chi_5^{(1)}(\omega)|e^{-i\delta_5^{(1)}(\omega)}$.

For completeness we recall that the expression of the linear polarisation $P_1(t)$ is $P_1(t)/\varepsilon_0 = E|\chi_1(\omega)|\cos(\omega t - \delta_1(\omega))$ where we have set $\chi_1(\omega) = |\chi_1(\omega)|e^{-i\delta_1(\omega)}$. 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 $\omega$. The only exception to this simplification is in Fig. 11 (see below) where for convenience the linear susceptibility is denoted $\chi_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 $\chi_3^{(3)}$ and/or the fifth-harmonic fifth-order susceptibility $\chi_5^{(5)}$ to which this chapter is devoted. In
Section 2.3 and Figs. 8-9, we shall briefly compare $\chi_3^{(3)}$ with other cubic susceptibilities, namely $\chi_3^{(1)}$ already defined in Eq. (5) as well as $\chi_{2;1}^{(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_3^{(1)}$ of the cubic polarization, which now contains two terms:

$$ \frac{P_3^{(1)}(t)}{\varepsilon_0} = \frac{3}{4} |\chi_3^{(1)}(\omega)| E^3 \cos(\omega t - \delta_3^{(1)}(\omega)) + 3 |\chi_{2,1}^{(1)}(\omega)| E_{st}^2 E \cos(\omega t - \delta_{2,1}^{(1)}(\omega)) \quad (7) $$

where we have defined $|\chi_{2,1}^{(1)}(\omega)| \exp(-i\delta_{2,1}^{(1)}(\omega)) = \chi_3(0, 0, \omega)$.

For any cubic susceptibility – generically noted $\chi_3$ – or for any fifth-order susceptibility – generically noted $\chi_5$ – the corresponding dimensionless susceptibility $X_3$ or $X_5$ is defined as :

$$ X_3 = \frac{k_B T}{\varepsilon_0 \Delta \chi_1^2 a^3} \chi_3, \quad X_5 = \frac{(k_B T)^2}{\varepsilon_0^2 \Delta \chi_1^3 a^6} \chi_5 \quad (8) $$

where $\Delta\chi_1$ is the “dielectric strength”, i.e. $\Delta\chi_1 = \chi_{lin}(0) - \chi_{lin}(\infty)$ where $\chi_{lin}(0)$ is the linear susceptibility at zero frequency and $\chi_{lin}(\infty)$ 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 great advantage to be both dimensionless and independent of the field amplitude.

## 2.2 Frequency and temperature dependence of third harmonic susceptibility

In this section we review the characteristic features of $\chi_3^{(3)}$ both as a function of frequency and temperature. We separate the effects at equilibrium above $T_g$ and those recorded below $T_g$ in the out-of-equilibrium regime.
<figure><img src="image_3.png" /><figcaption>**Fig. 2** Third-order harmonic component of the dielectric susceptibility of propylene carbonate [10]. Spectra of $|\chi_3^{(3)}|E^2$ are shown for various temperatures measured at a field of 225 kV/cm. The yellow-shaded plane indicates the plateau arising in the trivial regime.</figcaption></figure>

### 2.2.1 Above $T_g$

In the $\alpha$ regime:

Fig. 2 shows the modulus $|\chi_3^{(3)}|$ for propylene carbonate [10]. 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_\alpha \le 0.05$, a plateau is observed as indicated by the shaded area in Fig. 2, i.e. $|\chi_3^{(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 gas, $|\chi_3^{(3)}|$ has a plateau below the relaxation frequency and monotonously falls to zero as one increases the frequency. Because the observed plateau in Fig. 2 is reminiscent to the ideal 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 $\chi_3$ are “simple” –see Section 7–, but that the glassy correlations do not change qualitatively the shape of $\chi_3^{(3)}$ in this range. Physically, an ideal gas of dipoles corresponds to the high-$T$ limit of a fluid. This is why it is a useful benchmark which allows to distinguish the “trivial” features and those involving glassy correlations.

2. When rising the frequency above $0.05f_\alpha$ one observes for $|\chi_3^{(3)}|$ a hump for a frequency $f_{peak}/f_\alpha \simeq c$ where the constant $c$ does not depend on $T$ and weakly depends on the liquid (e.g., $c \simeq 0.22$ for glycerol and $c \simeq 0.3$ for propylene carbonate).
bonate). This hump is followed by a power law decrease $|\chi_3^{(3)}| \sim f^{-\beta_3}$ where $\beta_3 < 1$ is close [3] to the exponent governing the decrease of $|\chi_1|$ above $f_\alpha$. 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 $|\chi_1|$ of the supercooled liquids. This is why this hump has been termed the “glassy contribution” to $\chi_3$. On a more quantitative basis, the proportionality of $f_{\text{peak}}$ and of $f_\alpha$ has been observed for $f_\alpha$ ranging from 0.01 Hz to 10 kHz -above 10 kHz the measurement of $\chi_3^{(3)}$ is obscured by heating issues, see [11] and Section 7-.

<figure><img src="image_4.png" /><figcaption>**Fig. 3** (a) Modulus, real, and imaginary part of the third-order dielectric susceptibility $\chi_3^{(3)}$ (times $E^2$) of 1-propanol at 120 K as measured with a field of 468 kV/cm [16]. The solid lines were calculated according to Refs. [17]. (b) Same for glycerol at 204 K and 354 kV/cm [16].</figcaption></figure>

The consistency of the above considerations can be checked by comparing the third-order susceptibility of canonical glass formers to that of monohydroxy alcohols. The linear dielectric response of the latter is often dominated by a Debye relax-
ation 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 [12]. This process represents an idealised Debye-relaxation case as it lacks
the heterogeneity-related broadening found for other glass formers. Moreover, cor-
relations or cooperativity should not play a significant role for this process, because
cluster-cluster interactions can be expected to be rare compared to the intermolec-
ular interactions governing the α relaxation in most canonical glass formers [13].
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
dipole gas" case mentioned above. The monohydroxy alcohol 1-propanol is espe-
cially well suited to check this notion because here transitions between different
chain topologies, as found in several other alcohols affecting the nonlinear response
[14, 15], do not seem to play a role [15]. Figure 3(a) shows the frequency-dependent
modulus, real, and imaginary part of χ₃⁽³⁾E² for 1-propanol at 120 K [13, 16]. In-
deed, no hump is observed in |χ₃⁽³⁾|(ν) as predicted for a non-cooperative Debye
relaxation. The solid lines were calculated according to Refs. [17], accounting for
the expected trivial polarization-saturation effect. 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 monohy-
droxy 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 $\chi_3^{(3)}$ spectra related to the conventional $\alpha$ relaxation of canonical
glass formers exhibit strong deviations from the trivial response, just as expected in
the presence of molecular correlations. As an example, Fig. 3(b) shows the modu-
lus, real, and imaginary part of $\chi_3^{(3)}E^2$ of glycerol at 204 K. Again the lines were
calculated assuming the trivial nonlinear saturation effect only [17]. Obviously, this
approach is insufficient 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 response. This mirrors the fact that, on long time scales, the
liquid flow smoothes out any glassy correlations.

When varying the temperature, two very different behaviors of $\chi_3^{(3)}$ are observed:

1. In the plateau region the weak temperature dependence of $\chi_3^{(3)}$ is easily captured by converting $\chi_3^{(3)}$ into its dimensionless form $X_3^{(3)}$ by using Eq. (8): one observes [3, 4] that in the plateau region $X_3^{(3)}$ does not depend at all on the temperature. Qualitatively this is important since in an ideal gas of dipoles $X_3^{(3)}$ does also not depend on temperature, once plotted as a function of $f/f_\alpha$. This reinforces the “trivial” nature of the plateau region, i.e. the fact that it is not qualitatively affected by glassy correlations.

2. In the hump region, $|X_3^{(3)}(f/f_\alpha)|$ increases upon cooling, again emphasizing the “anomalous” –or “non trivial”– behavior of the glassy contribution to $\chi_3^{(3)}$. This increase of the hump of $|X_3^{(3)}|$ has been related to that of the apparent
<figure><img src="image_5.png" /><figcaption>**Fig. 4** From Ref. [18]: For several glass formers, $N_{corr}(T)$ as extracted from the hump of $|X_3^{(3)}|$ (left axis) closely follows $E_{act}(T)$, deduced from the temperature dependence of the $\alpha$-relaxation time [10] (right axis). The abbreviations stand for propylene carbonate (PCA), 3-fluoroaniline (FAN), 2-ethyl-1-hexanol (2E1H), cyclo-octanol (c-oct), and a mixture of 60% succinonitrile and 40% glutaronitrile (SNGN).</figcaption></figure>

activation energy $E_{act}(T) \equiv \partial \ln \tau_{\alpha}/\partial (1/T)$ -see refs. [10, 18]- as well as to $T\chi_T \equiv |\partial \ln \tau_{\alpha}/\partial \ln T|$ [3, 4, 19, 20]. Note that because the experimental temperature interval is not so large, the temperature behavior of $E_{act}$ and of $T\chi_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} \sim N_{corr}(T)$ dates back to the work of Adam and Gibbs [21]; while another series of papers [22, 23] proposed a decade ago that $N_{corr} \propto T\chi_T$. Fig. 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 different fragilities, including a plastic crystal, where only the orientational degrees of freedom experience the glass transition [24].

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 $\alpha$ peak [25]. Figure 5(a) shows loss spectra of glycerol, measured at low and high fields up to 671 kV/cm [26, 27], where the excess wing is indicated by the dashed lines. (It should be noted that the difference of these loss curves for high and low fields is directly related to the cubic susceptibility $\chi_3^{(1)}$, defined in Eq. (5) [16].) As already reported in the seminal paper by Richert and Weinstein [28], in Fig. 5(a) at the right flank of the $\alpha$-relaxation peak a strong field-induced increase of the dielectric loss is found while no significant field dependence is detected at its low-
frequency flank. In Ref. [28] it was pointed out that these findings are well consistent with the heterogeneity-based box model (see Section 4.3). However, as revealed by Fig. 5(a), remarkably in the region of the excess wing no significant nonlinear effect is detected. Time-resolved measurements, later on reported by Samanta and Richert [29], revealed nonlinearity effects in the excess-wing region when applying the high field for extended times of up to several 10000 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 differs from that of the $\alpha$ relaxation.

To check whether weaker nonlinearity in the excess-wing region is also revealed in higher-harmonic susceptibility measurements, Fig. 5(b) directly compares the modulus of the linear dielectric susceptibility of glycerol at 191 K to the third-order susceptibility $|\chi_3^{(3)}|$ (multiplied by $E^2$) [30]. (We show $|\chi_1|$ corrected for $\chi_{1,\infty} = \epsilon_\infty - 1$ caused 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 found in $|\chi_3^{(3)}(\nu)|$. 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 $\alpha$ relaxation.

The excess wing is often regarded as the manifestation of a secondary relaxation process, partly superimposed by the dominating $\alpha$-relaxation [31, 32]. 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 [33, 34]. Moreover, in a recent work [35] it was pointed out that the small or even absent nonlinear effects in the excess-wing region can also be consistently explained within the framework of the coupling model [34], where the excess wing is identified with the so-called "nearly constant loss" caused by caged molecular motions.

### 2.2.2 Below $T_g$

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 transition temperature $T_g$ has been crossed. The mechanism of aging is still a matter of controversy [36, 37, 38, 39, 40], owing to the enormous theoretical and numerical difficulties inherent to out-of-equilibrium processes. Experimentally, a few clear cut results have been obtained in spin glasses [41] where it was shown, by using nonlinear techniques, that the increase of the relaxation time $\tau_\alpha$ 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 [42] to the interpretation given previously in the experiments of Ref. [41].
<figure><img src="image_6.png" /><figcaption>**Fig. 5** (a) Dielectric loss of glycerol measured at fields of 14 kV/cm (open symbols) and 671 kV/cm (closed symbols) shown for four temperatures [27]. The solid lines were measured with 0.2 kV/cm [26]. The dashed lines indicate the excess wing. (b) Open triangles: Absolute values of $\chi_1$ (corrected for $\chi_{1,\infty} = \epsilon_\infty - 1$) at 14 kV/cm for glycerol at 191 K. Closed triangles: $\chi_3^{(3)}E^2$ at 565 kV/cm [30]. 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 $\chi_3^{(3)}(\nu)$.</figcaption></figure>

In structural glasses, the aging properties of the linear response have been reported more than one decade ago [43, 44]. More recently, the aging properties of $\chi_3^{(3)}$ were reported in glycerol [45] and its main outputs are summarized in Figs. 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 effect is the increase of
<figure><img src="image_7.png" /><figcaption>**Fig. 6** From [45]. During the aging of glycerol -at $T_g$ – 8K- the increase of $\tau_\alpha$ with the aging time $t_a$ is measured by rescaling the aging data -symbols- of $\chi_1''$ -right axis- onto the equilibrium data -solid black line-. The corresponding scaling fails for $X_3^{(3)}(f, t_a)$ -left axis- revealing the increase of $N_{corr}$ during aging. See [45] for details about the quantity $z(t_a)/z(T)$ which is involved in the left axis but varies by less than 2% during aging-.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Fig. 7** From [45]. The values of $\delta = N_{corr}(t_a)/N_{corr}(eq)$ extracted from Fig. 6 show the increase of $N_{corr}$ during aging. Inset: different theories are tested gathering equilibrium and aging experiments.</figcaption></figure>

the relaxation time $\tau_\alpha$ with $t_a$. In Ref. [45] $\tau_\alpha$ increases by a factor $\simeq 6$ between the arrival at $T_w$ -i.e. $t_a = 0$- and the finally equilibrated situation reached for $t_a \gg \tau_{\alpha,eq}$ where $\tau_\alpha$ is equal to its equilibrium value $\tau_{\alpha,eq}$ -and no longer evolves with $t_a$. This variation of $\tau_\alpha$ 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 $\chi''(f, t_a)$. This is summarized in Fig. 6 for 5 different frequencies: when plotted as a function of $f/f_\alpha(t_a) \equiv 2\pi f \tau_\alpha(t_a)$, the aging values of $\chi''(f, t_a)$ -symbols- are nicely rescaled onto the equilibrium values $\chi''(f, eq)$ -continuous line- measured when $t_a \gg \tau_{\alpha,eq}$.
The most important experimental result is that this scaling fails for $|X_3^{(3)}(f, t_a)|$ as shown by the left axis of Fig. 6: For short aging times, the difference 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 Fig. 7 where the variation of $\delta = 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 $|\chi_3^{(3)}|$ measured at equilibrium above $T_g$ and by translating the $\tau_\alpha(t_a)$ in terms of a fictive temperature $T_{fict}(t_a)$ which decreases during aging, finally reaching $T_w$ when $t_a \gg \tau_{\alpha,eq}$. This yields the continuous line in Fig. 7, which fairly well captures the data drawn from the aging of $\chi_3^{(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 Fig. 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 $\chi_3^{(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 Fig. 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 different predictions by $\ln(\tau_\alpha/\tau_0) = YN_{corr}^{\psi/3}/(k_B T)$ with $Y \sim T$; $\psi = 3/2$ for Random First Order Transition theory (RFOT) [46] while $Y \sim 1$; $\psi = 1$ for the numerical approach of Ref. [47], Fig. 7 is designed to test these two predictions -see Ref. [45] for details-: it shows that both of them are consistent with experiments -contrary to another prediction relying onto a critical relation $\tau_\alpha \propto N_{corr}^z$, yielding an unrealistic large value of $z \sim 20$ to account for the experiments-.

## 2.3 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 $\chi_3^{(3)}$ as well as the first-harmonic cubic susceptibilities $\chi_3^{(1)}$ and $\chi_{2;1}^{(1)}$ introduced in Eq. (7). We remind that $\chi_{2;1}^{(1)}$ corresponds to the case where a static field $E_{st}$ is superimposed to the ac field $E \cos(\omega t)$.

Figs. 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 $|\chi_3^{(1)}|$ and of $|\chi_{2;1}^{(1)}|$ are the same as those of $|\chi_3^{(3)}|$: when increasing the frequency, one first observes a low frequency plateau, followed by a hump in
<figure><img src="image_9.png" /><figcaption>**Fig. 8** From [48]: For glycerol and $f_\alpha \approx 2$ Hz, modulus -top panel- and phase -bottom panel- of the three cubic susceptibilities defined in Eqs. (5) and (7). The salient features of the three cubic susceptibilities are similar, which strongly suggests a common physical origin -see text-. Dotted lines are $Arg[X_{2;1}^{(1)}] + \pi$ or $+2\pi$ and support Eqs. (9) and (10).</figcaption></figure>

the vicinity of $f_\alpha$ and then by a power law decrease $\sim f^{-\beta_3}$. The most important differences 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 $|\chi_3^{(1)}|$ and for $|\chi_{2;1}^{(1)}|$ what we have already seen for $|\chi_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 $E_{act}(T) = \partial \ln \tau_\alpha / \partial (1/T)$ which in this $T$ range is very close to $T\chi_T = |\partial \ln \tau_\alpha / \partial \ln T|$ [3, 4, 10, 19, 20, 48].

2. The phases of the three cubic susceptibilities basically do not depend explicitly on temperature, but only on $u = f/f_\alpha$, 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
<figure><img src="image_10.png" /><figcaption>Fig. 9 From [48]. Same representation as in Fig. 8 but for propylene carbonate.</figcaption></figure>

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):

$$ \operatorname{Arg}[X_3^{(1)}] \approx \operatorname{Arg}[X_{2,1}^{(1)}] + \pi \quad \text{for } f/f_\alpha \ge 0.5; \qquad (9) $$

$$ \operatorname{Arg}[X_3^{(1)}] \approx \operatorname{Arg}[X_3^{(3)}] \quad \text{for } f/f_\alpha \ge 5 \qquad (10) $$

which are quite non trivial relations.

3. In the phase of $\chi_3^{(1)}$ of propylene carbonate (Fig. 9), a jump of $\pi$ is observed which is accompanied by the indication of a spikelike minimum in the modulus -see [48] for more details-. A similar jump may also be present in glycerol (Fig. 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 $\chi'_1$, while around the hump $\chi'_1$ is enhanced. At the frequency of the jump, both effects
compensate and this coincides with a very low value of the imaginary part of
$X_3^{(1)}$.

## 2.4 Frequency and temperature dependence of fifth harmonic susceptibility

In this section, we first explain why measuring $\chi_5^{(5)}$ is interesting for a better understanding of the glass transition. We then see the characteristic features of $\chi_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 \text{ Hz} \le f_\alpha \le 10 \text{ 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 $\tau_\alpha$. 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 $\tau_\alpha$: extrapolating the above result, e.g., $|X_3| \propto |\partial \ln \tau_\alpha / \partial \ln T|$ and assuming a VFT law for $\tau_\alpha$, one concludes that the experimental characteristic times corresponding to an increase of $|X_3|$ by one order of magnitude is $0.1 \text{ ms} \le \tau_\alpha \le 10^{18} \text{ 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 [49], 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 Eq. (2) that, apart from $\chi_1$ which is blind to glassy correlations, all other responses $\chi_{k \ge 3}$ grow as power laws with the amorphous ordering length $\ell$: $\chi_3 \propto (1/a)^{2d_f-d}$ and $\chi_5 \propto (1/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 $\chi_5 \propto \chi_3^{\mu(d_f)}$ where $\mu(d_f) = (3d_f - d)/(2d_f - d)$.

Measuring $\chi_5$ is of course extremely difficult, because, for the experimentally available electric fields, one has the hierarchy $|\chi_1|E \gg |\chi_3|E^3 \gg |\chi_5|E^5$. However this was done in Ref. [50] and we shall now briefly review the corresponding results.
<figure><img src="image_11.png" /><figcaption>**Fig. 10** From Ref. [50]. Measured values of |χ₅⁽⁵⁾| for glycerol - upper panel- and propylene carbonate - lower panel- (the spheres and cubes in the upper panel indicate results from two different experimental setups). The hump lies at the same frequency as for |χ₃⁽³⁾| and has significantly stronger variations in frequency and in temperature, see Figs 11 and 12. The arrows indicate the peak positions fα in the dielectric loss. The yellow-shaded planes indicate the plateau arising in the trivial regime.</figcaption></figure>

### 2.4.2 Characteristic features of the fifth order susceptibility

The modulus |χ₅⁽⁵⁾| of glycerol and propylene carbonate [50] can be seen in Fig. 10 as a function of frequency and temperature. Similarly to what has been seen in section 2.2 on |χ₃⁽³⁾|, the frequency dependence can be separated in two domains (see also Fig. 11):

1. For very low reduced frequencies ($f/f_α$ ≤ 0.05), there is a plateau (indicated by the yellow-shaded planes in Fig. 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
<figure><img src="image_12.png" /><figcaption>Fig. 11 From Ref. [50]. For glycerol, comparison of the fifth, third and linear susceptibilities -the latter is noted $|\chi_1^{(1)}|$. The hump for $|\chi_5^{(5)}|$ is much stronger than that of $|\chi_3^{(3)}|$. The dashed lines are the trivial contribution -see [50] for details.</figcaption></figure>

expected from a high temperature liquid of dipoles, depicted by the "trivial" $X_k^{(k)}$ curves represented as dotted lines in Fig. 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 $|\chi_3^{(3)}|$ in both glycerol and propylene carbonate. Again one finds that, for the five temperatures where the peak is studied, $f_{peak}/f_\alpha = 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 Fig. 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^{(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_3^{(3)}|$:

1. In the plateau region at $f/f_\alpha \le 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^{(5)}$ from $\chi_5^{(5)}$ -see Eq. (8)- is extremely efficient to remove all trivial temperature dependences. As the trivial behavior depends on frequency -see the dashed lines of Fig. 11-, the "singular" parts of $X_3$ and of $X_5$ are obtained as follows:

$$X_{3,sing.}^{(3)} \equiv X_{3}^{(3)} - X_{3,trivial}^{(3)}, \qquad X_{5,sing.}^{(5)} \equiv X_{5}^{(5)} - X_{5,trivial}^{(5)} \qquad (11)$$
<figure><img src="image_13.png" /><figcaption>**Fig. 12** From Ref. [50] Temperature evolution of the singular parts of fifth and third order re-
sponses. All quantities are normalized at a given temperature -namely 207 K for glycerol, upper
panel; and 164 K for propylene carbonate, bottom panel-. This allows to determine the exponent
$\mu$ relating $|X_5|$ and $|X_3|^\mu$ and to conclude that the amorphously ordering domains are compact -see
text-. The hatched areas represent the uncertainty on $\mu$.</figcaption></figure>

which correspond in Fig. 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<sup>(5)</sup><sub>5,sing.</sub> (f<sub>peak</sub>)| is compared to that of |X<sup>(3)</sup><sub>3,sing.</sub> (f<sub>peak</sub>)|^μ where μ is an exponent that is determined experimentally by looking for the best overlap of the two series of data in Fig. 12 -see [50] for details-. This leads us to values of μ = 2.2 ± 0.5 in glycerol and μ = 1.7 ± 0.4 in propylene carbonate. Therefore, within experimental uncertainties, results for |X<sup>(3)</sup><sub>3</sub>| and |X<sup>(5)</sup><sub>5</sub>| would seem to advocate a value of μ ≈ 2. With μ = (3d<sub>f</sub> - d)/(2d<sub>f</sub> - d) as seen in Eq. (2) -see also Eq. (13) below-, this corresponds to a fractal dimensions of d<sub>f</sub> ≈ 3.
## 3 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.

### 3.1 Bouchaud-Biroli's predictions

#### 3.1.1 General considerations about $\chi_{2k+1}$

To illustrate the general relations existing between the susceptibility $\chi_{2k+1}$ and the correlation function of order $2k+2$ -with $k \ge 0$- in a system at thermal equilibrium, let us consider a sample, submitted to a constant and uniform magnetic field $h$, containing $N$ spins with an Hamiltonian $\mathcal{H}$ that depends on the spin configuration “c”. The elementary relations of statistical physics yield the magnetisation $M \equiv \sum_i <S_i>/ (Na^3)$ where $a^3$ is the elementary volume and where the thermal average $<S_i>$ is obtained with the help of the partition function $Z = \sum_c \exp(-\beta \mathcal{H} + \beta h \sum_k S_k)$ by writing $<S_i> = \sum_c S_i \exp(-\beta \mathcal{H} + \beta h \sum_k S_k)/Z$ with $\beta = 1/(k_B T)$. The linear response $\chi_1 \equiv (\partial M/\partial h)_{h=0}$ is readily obtained:

$$Na^3 \chi_1 = \frac{1}{\beta Z} \left(\frac{\partial^2 Z}{\partial h^2}\right)_{h=0} - \frac{1}{\beta} \left(\frac{\partial Z}{Z \partial h}\right)_{h=0}^2 = \beta \left(\sum_{i1,i2} \langle S_{i1}S_{i2} \rangle - (\sum_{i1} \langle S_{i1} \rangle)^2\right) \quad (12)$$

which shows that the linear response is related to the connected two-point correlation function. Repeating the argument for higher-order responses -e.g. $\chi_3 \propto (\partial^3 M/\partial h^3)_{h=0}$-, one obtains that $\chi_{2k+1}$ is connected to the $(2k+2)$ points correlation function -e.g., $\chi_3$ is connected to a sum combining $\langle S_{i1}S_{i2}S_{i3}S_{i4} \rangle$, $\langle S_{i1}S_{i2}S_{i3} \rangle \langle S_{i4} \rangle$, $\langle S_{i1}S_{i2} \rangle \langle S_{i3}S_{i4} \rangle$ etc...-.

#### 3.1.2 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 $\chi_1$ which does not diverge at $T_{SG}$: this
is because the various terms of $\sum_{i1;i2} <S_{i1}S_{i2}>$ cancel since half of them are positive and the other half are negative. By contrast the cubic susceptibility $\chi_3$ contains a term $\sum_{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 $\chi_3$ when decreasing $T$ towards $T_{SG}$, as already illustrated in Fig. 1-(A). By adding a standard scaling assumption close to $T_{SG}$ one can account for the behavior of $\chi_3$ at finite frequencies, i.e. one easily explains that $\chi_3$ is frequency independent for $\omega\tau_\alpha \le 1$, and smoothly tends to zero at higher frequencies. Finally, similar scaling arguments about correlation functions easily explain the fact that the stronger $k \ge 1$ the more violent the divergence of $\chi_{2k+1}$ in spin glasses, as observed experimentally by Levy et al [51].

### 3.1.3 The glass forming liquids case

The case of glass forming liquids is of course different from that of spin glasses for some obvious reasons (e.g. molecules have both translational and rotational degrees of freedom). As it has been well established that rotational and translational degrees of freedom are well coupled in most of liquids, it is tempting to attempt a mapping between spin glasses and glass forming liquids by replacing the spins $S_i$ by the local fluctuations of density $\delta\rho_i$ or by the dielectric polarisation $p_i$. As far as nonlinear responses are concerned, this mapping requires a grain of salt because (a) there is no frozen-in disorder in glass forming liquids, and (b) there is a nonzero value of the molecular configurational entropy $S_c$ around $T_g$.

The main physical idea of BB's work [9] is that these difficulties have an effect which is important at low frequencies and negligible at high enough frequencies:

1. Provided $f \ge f_\alpha$, 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 Ref. [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 $\chi_3(\omega)$. 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, 50]:

$$X_{2k+1}^{\text{glass}}(f, T) = [N_{\text{corr}}(T)]^{\alpha_k} \times \mathcal{H}_k \left(\frac{f}{f_\alpha}\right) \quad \text{with} \quad \alpha_k = (k+1) - d/d_f \quad (13)$$

where the scaling functions $\mathcal{H}_k$ do not explicitly depend on temperature, but depend on the kind of susceptibility that is considered, i.e. $X_3^{(1)}$, $X_3^{(3)}$ or $X_{2,1}^{(1)}$ in
the third order case $k=1$. We emphasize that in Ref. [9] the amorphously ordered
domains were assumed to be compact, i.e. $d_f = d$, yielding $\alpha_1 = 1$ i.e. $X_3 \propto N_{corr}$.
The possibility of having a fractal dimension $d_f$ lower than the spatial dimension
$d$ was considered in Ref. [50] 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 \propto N_{corr}^2$.

2. In the low frequency regime $f \ll f_\alpha$, relaxation has happened everywhere in the system, destroying amorphous order [52] and the associated anomalous response to the external field and $\mathcal{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 effectively independent molecules.

Due to the definition adopted in Eq. (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_{2k+1}^{glass}$ should exceed the trivial contribution, leading to hump-shaped non-linear susceptibilities, peaking at $f_{peak} \sim f_\alpha$, where the scaling function $\mathcal{H}_k$ reaches its maximum.

## 3.2 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 $|\chi_3^{(3)}|, |\chi_3^{(1)}|, |\chi_{2;1}^{(1)}|$ and of $|\chi_5^{(5)}|$ have a humped shape in frequency, contrary to $|\chi_1|$.

2. Due to the fact that $\mathcal{H}_k$ does not depend explicitly on $T$, the value of $f_{peak}/f_\alpha$ should not depend on temperature, consistent with the experimental behavior.

3. Because of the dominant role played by the glassy response for $f \ge 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 func-
tions, it is natural that the behavior of their modulii and phases are quantitatively
related especially at high frequency where the "trivial" contribution can be ne-
glected, consistent with Eqs. (9)-(10) –see below for a more quantitative argu-
ment in the context of the so-called "Toy model"– [53].

Having shown that BB’s prediction is consistent with experiments, the temper-
ature variation of $N_{corr}$ can be drawn from the increase of the hump of $X_3$ upon
cooling. It has been found [3, 4, 10, 19, 20] 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\chi_T$, which was
proposed in Refs. [22, 23] as a simplified estimator of $N_{corr}$ in supercooled liq-
uids. The convergence of these different estimates, that rely on general, model-free
theoretical arguments, is a strong hint that the underlying physical phenomenon is
indeed the growth of collective effects 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 $\ell$ of amorphously ordered domains is (much) larger than the molecular size $a$. This naturally explains the similarities of the cubic responses in microscopically very different liquids such as glycerol and propylene carbonate, as well as many other liquids [10, 20]. Indeed the microscopic differences are likely to be wiped out for large $\ell \propto N_{corr}^{1/df}$, much like in usual phase transitions.

### 3.3 Static vs dynamic length scale? $\chi_3$ and $\chi_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 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 Ref. [50] 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 [50]: 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) [56] 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 [50]. This is not the case for theories (such as RFOT [46] or Frustration theories [57]) where a non-trivial thermodynamic critical point drives the glass transition: in this case the incipient amorphous order allows to account [50] for the observed features of $X_3$ and $X_5$. This is why it was argued in [48, 50] 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 [55] –this statement will be reinforced in section 4.2.

## 4 More specific models for harmonic susceptibilities

We now review the various other approaches that have been elaborated for the non-linear 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.

## 4.1 Toy and pragmatical models

The “Toy model” has been proposed in Refs. [19, 58] as a simple incarnation of the BB mechanism, while the “Pragmatical model” is more recent [59, 60]. Both models start with the same assumptions: (i) each amorphously ordered domain is compact and contains $N_{corr}$ molecules, which yields a dipole moment $\propto \sqrt{N_{corr}}$ and leads to an anomalous contribution to the cubic response $X_3^{glass} \propto N_{corr}$; (ii) there is a crossover at low frequencies towards a trivial cubic susceptibility contribution $X_3^{triv}$ 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 $\delta$, favoring one well over the other. The most important difference between the Toy and the Pragmatical model comes from the description of the low-frequency crossover, see Refs. [58] and [60] for more details.

On top of $N_{corr}$ and $\delta$, the Toy model uses a third adjustable parameter, namely the frequency $f^*$ below which the trivial contribution becomes dominant. In Ref. [58], both the modulus and the phase of $X_3^{(3)}(\omega, T)$ and of $X_3^{(1)}(\omega, T)$ in glycerol were well fitted by using $f^* \simeq f_\alpha/7$, $\delta = 0.6$ and, for $T = 204$ K, $N_{corr} = 5$ for $X_3^{(3)}$ and $N_{corr} = 15$ for $X_3^{(1)}$. Fig. 13 gives an example of the Toy model prediction for $X_3^{(3)}$ in glycerol. Besides, in Ref. [19], the behavior of $X_{2,1}^{(1)}(\omega, T)$ in glycerol was fitted with the same values of $\delta$ and of $f^*$ but with $N_{corr} = 10$ (at a slightly different temperature $T = 202$ K). Of course, the fact that a different 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 [19, 58] by using the very same values of $f^*/f_\alpha$ and of $\delta$. This means that Eqs. (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 Eqs. (9) and (10).

Let us mention briefly the Asymmetric Double Well Potential (ADWP) model [61], which is also about species living in a double well of asymmetry energy $\Delta$, 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 [61, 62] only when the reduced asymmetry $\delta = \tanh(\Delta/(2k_B T))$ is close to a very specific value, namely $\delta_c = \sqrt{1/3}$, where $X_3$ vanishes at zero frequency due to the compensation of its several terms. However, at the fifth order [62] this compensation happens for two values of $\delta$ very different from $\delta_c$: as a result the model cannot predict
<figure><img src="image_14.png" /><figcaption>**Fig. 13** From Ref. [58]. Fit of the values of $X_3^{(3)}$ measured in glycerol -symbols- at 204 K by using the Toy model with $N_{corr} = 5$, $\delta = 0.6$ and $f^* \approx f_\alpha/7$. The prediction of the Toy model is given by the two thick solid lines (main panel for the modulus of $X_3^{(3)}$ and inset for its phase).</figcaption></figure>

a hump happening both for the third and for the fifth order in the same parametric regime, contrarily to the experimental results of Ref. [50]. This very recent calculation of fifth order susceptibility [62] 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 [58] an anomalous contribution $X_{2k+1}^{glass} \sim [N_{corr}]^k$, provided that $N_{corr}$ is large enough, the magnitude of this contribution is much larger than that of the small trivial contribution $X_{2k+1}^{triv.} \sim 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 \ll f_\alpha$. As for the right side of the peak, it comes from the fact that $|X_{2k+1}| \to 0$ when $f \gg f_\alpha$ for the simple reason that the supercooled liquid does not respond to the field at very large frequencies.

## 4.2 Entropic effects

A contribution to nonlinear responses was recently calculated by Johari in Refs [63, 64] 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 $\chi_{2;1}^{(1)}$ experiments in Refs [65, 66, 67, 68] -see however Ref. [69] for a case where the agreement is not as good-. It was then extended to pure ac experiments -and thus to $\chi_3^{(3)}$- in Refs. [70, 71]. The relation between Johari's idea and $N_{corr}$ was made in Ref. [48].
### 4.2.1 When a static field $E_{st}$ is applied

Let us start with the case of $\chi_{2;1}^{(1)}$ experiments, i.e. with the case where a static field $E_{st}$ is superimposed onto an ac field $E \cos(\omega t)$. In this case, there is a well defined variation of entropy $[\delta S]_{E_{st}}$ induced by $E_{st}$, which, for small $E_{st}$ and a fixed $T$, is given by:

$$[\delta S]_{E_{st}} \approx \frac{1}{2} \epsilon_0 \frac{\partial \Delta \chi_1}{\partial T} E_{st}^2 a^3, \quad (14)$$

where $a^3$ is the molecular volume. Eq. (14) holds generically for any material. However, in the specific case of supercooled liquids close enough to their glass transition temperature $T_g$, a special relation exists between the molecular relaxation time $\tau_\alpha$ and the configurational contribution to the entropy $S_c$. This relation, first anticipated by Adam and Gibbs [21], can be written as:

$$\ln \frac{\tau_{\alpha}(T)}{\tau_{0}} = \frac{\Delta_{0}}{T S_{c}(T)} \quad (15)$$

where $\tau_0$ is a microscopic time, and $\Delta_0$ is an effective energy barrier for a molecule. The temperature dependence of $TS_c(T)$ quite well captures the temperature variation of $\ln(\tau_\alpha)$, at least for a large class of supercooled liquids [72].

Following Johari [63, 64] let us now assume that $[\delta S]_{E_{st}}$ is dominated by the dependence of $S_c$ on field, -see the Appendix of Ref. [48] for a further discussion of this important physical assumption-. Combining Eqs. (14) and (15), one finds that a static field $E_{st}$ produces a shift of $\ln(\tau_\alpha/\tau_0)$ given by:

$$[\delta \ln \tau_{\alpha}]_{E_{st}} = - \frac{\Delta_0}{T S_c^2} [\delta S]_{E_{st}} \quad (16)$$

As shown in Ref. [48] this entropic effect gives a contribution to $X_{2,1}^{(1)}$, which we call $J_{2,1}^{(1)}$ after Johari. Introducing $x = \omega\tau_\alpha$, the most general and model-free expression of $J_{2,1}^{(1)}$ reads:

$$J_{2,1}^{(1)} = - \frac{k_B \Delta_0}{6 S_c^2} \left[ \frac{\partial \ln(\Delta \chi_1)}{\partial T} \right] \left[ \frac{\partial \chi_{lin}}{\partial \ln x} \right] \propto \frac{1}{S_c^2} \quad (17)$$

where $\chi_{lin}$ is the complex linear susceptibility.

Eq. (17) deserves three comments:

1. $|J_{2,1}^{(1)}|$ has a humped shaped in frequency with a maximum in the region of $\omega\tau_\alpha \approx 1$, because of the frequency dependence of the factor $\propto \partial\chi_{lin}/\partial \ln x$ in Eq. (17).

2. The temperature variation of $J_{2,1}^{(1)}$ is overwhelmingly dominated by that of $S_c^{-1}$ because $S_c \propto (T - T_K)$ -with $T_K$ the Kauzmann temperature-.

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 state and $N_{corr}$ would diverge-. In other words, there is a relation between $S_c^{-1}$ and $N_{corr}$, which yields [48]:

$$J_{2,1}^{(1)} \propto N_{corr}^q, \qquad (18)$$

where it was in shown in Ref. [48] 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 [48].

b. the exponent $q$ lies in the $[1/3; 3/2]$ interval [48] when one uses the RFOT and plays with its two critical exponents $\Psi$ and $\theta$. Notably, taking the “recommended RFOT values” $-\Psi = \theta = 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 [48] for a refined discussion-.

### 4.2.2 When a pure ac field $E \cos(\omega t)$ is applied

Motivated by several works [65, 66, 67, 68] showing that Johari’s reduction of entropy fairly well captures the measured $\chi_{2;1}^{(1)}$ in various liquids, an extension of this idea was proposed in Refs. [70, 71] for pure ac experiments, i.e. for $\chi_3^{(3)}$ and $\chi_3^{(1)}$. This has given rise to the phenomenological model elaborated in Refs. [70, 71] where the entropy reduction depends on time, which is nevertheless acceptable in the region $\omega\tau_\alpha \le 1$ where the model is used. Fig. 14 shows the calculated values for $|\chi_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 Ref. [71], the hump displayed in Fig. 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 effects, 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 effects can be calculated by using the Adam-Gibbs relation. Once combined, these two assumptions give a contribution to $\chi_{2;1}^{(1)}$ reasonably well in agreement with the measured values in several liquids [65, 66, 67, 68]. An extension to $\chi_3^{(3)}$ is even possible, at least in the region $\omega\tau_\alpha \le 1$ and fairly well accounts for the measured hump of $|\chi_3^{(3)}|$ in glycerol [70, 71] -a figure similar to Fig. 11 for $|\chi_5^{(5)}(\omega)/\chi_5^{(5)}(0)|$ is even obtained in Ref. [71]-. As shown in Eq. (18), this entropy contribution 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, Eq. (18) supports the interpretation that the various cubic susceptibilities $\chi_3$ are related to static amorphous correlations, as discussed in Section 3.3.
<figure><img src="image_15.png" /><figcaption>Fig. 14 From Ref. [71]. The model elaborated in Refs. [70, 71] includes three contributions - entropy reduction, Box model, and trivial-. It predicts for $|\chi_3^{(3)}|$ the solid lines which account very well for the measured values in glycerol in frequency and in temperature. The peak of $|\chi_3^{(3)}|$ arises because of the entropy reduction effect (noticed “sing. $T_{fic.}$”) which completely dominates the two other contributions in the peak region, as shown by the inset.</figcaption></figure>

## 4.3 Box model

### 4.3.1 Are nonlinear effects related to energy absorption ?

The “Box model” is historically the first model of nonlinear response in supercooled liquids, designed to account for the Nonresonant Hole Burning (NHB) experiments [73]. When these pioneering experiments were carried out, a central question was whether the dynamics in supercooled liquids is homogeneous or heterogeneous. In the seminal ref. [73] it was reported that when applying a strong ac field $E$ of angular frequency $\omega$, the changes in the dielectric spectrum are localised close to $\omega$ and that they last a time of the order of $1/\omega$. These two findings yield a strong qualitative support to the heterogeneous 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 $\tau$ and obeying the Debye dynamics. The distribution of the various $\tau$'s is chosen to recover the measured non Debye spectrum by adding the various linear Debye susceptibilities $\chi_{1,dh} = \Delta\chi_1/(1 - i\omega\tau)$ of the various domains. For the nonlinear response, the Box model assumes that it is given by the Debye linear equation in which $\tau(T)$ is replaced by $\tau(T_f)$ where the fictive temperature $T_f = T + \delta T_f$ is governed by the constitutive equation -see e.g. [28, 77]-:

$$c_{dh} \frac{\partial(\delta T_f)}{\partial t} + \kappa \delta T_f = \frac{1}{2} \varepsilon_0 \chi_{1,dh}'' \omega E^2 \quad (19)$$
with $c_{dh}$ the volumic specific heat of the DH under consideration, $\kappa$ the thermal conductance (divided by the DH volume $v$) between the DH and the phonon bath, $\tau_{therm} = c/\kappa$ the corresponding thermal relaxation time. In Eq. (19), only the constant part of the dissipated power has been written, omitting its component at $2\omega$ which is important only for $\chi_3^{(3)}$ -see e.g. [77]-. From Eq. (19) one easily finds the stationary value $\delta T_f^*$ of $\delta T_f$ which reads:

$$ \delta T_f^* = \frac{\tau_{therm} \epsilon_0 \Delta \chi_1 E^2}{\tau 2c_{dh}} \frac{\omega^2 \tau^2}{1 + \omega^2 \tau^2} \quad (20) $$

As very clearly stated in the seminal Ref. [73] because the DH size is smaller than 5nm, the typical value of $\tau_{therm}$ is at most in the nanoseconds range: this yields, close to $T_g$, a vanishingly small value of $\tau_{therm}/\tau$, which, because of Eq. (20), gives fully negligible values for $\delta T_f^*$. The choice of the Box model is to increase $\tau_{therm}$ by orders of magnitude by setting $\tau_{therm} = \tau$, 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 efficiency*: it allows to account reasonably well for the NHB experiments [73] and thus to bring a strong support to the heterogeneous character of the dynamics in supercooled liquids. Since the seminal ref. [73], some other works have shown [28, 74, 75, 76] that the Box model efficiently accounts for the measured $\chi_3^{(1)}(f > f_\alpha)$ in many glass forming liquids. It was shown also [77] that the Box model is not able to fit quantitatively the measured $X_3^{(3)}$ (even though some qualitative features are accounted for), and that the Box model only provides a vanishing contribution to $X_{2,1}^{(1)}$ – see [19].

The key choice $\tau_{therm} = \tau$ made by the Box model has two important consequences for cubic susceptibilities: it implies a) that $\chi_3^{(1)}$ mainly comes from the energy absorption (since the source term in Eq. (19) is the dissipated power) and b) that $\chi_3^{(1)}$ does not explicitly depend on the volume $v = N_{corr}a^3$ of the DH’s (see [28, 77]). However, alternative models of nonlinear responses are now available [58, 60] where, instead of choosing $\tau_{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 $\tau_c$ in Appendix 3). By using this molecular physics approach one obtains that $\chi_3^{(1)}$ is governed by $N_{corr}$ and not by energy absorption. For $\chi_3^{(1)}$, writing loosely $P_3^{(1)} \approx \partial P_1/(\partial \ln \tau)\delta \ln \tau$, one sees that the pivotal quantity is the field induced shift of the relaxation time $\delta \ln \tau$. Comparing the Box model (BM) and, e.g., the Toy model (TM), one gets respectively:

$$ \delta \ln \tau_{BM} \approx -\frac{1}{2} \chi_T \frac{\epsilon_0 \Delta \chi_1 E^2}{c_{dh}} ; \quad \delta \ln \tau_{TM} \approx -\frac{3}{2} \frac{N_{corr}}{T} \frac{\epsilon_0 \Delta \chi_1 E^2}{k_B/a^3} \quad (21) $$
where we remind our definition $\chi_T = |\partial \ln \tau_\alpha / \partial T|$ and where the limit $\omega \tau \gg 1$, relevant for the $\chi_3^{(1)}(f > f_\alpha)$ was taken in the Box model, while the simplest case (symmetric double well with a net dipole parallel to the field) was considered for the Toy model. Eq. (21) deserves two comments:

1. one sees that the two values of $\delta \ln \tau$ are similar provided $N_{corr}$ and $T\chi_T$ are proportional – which is a reasonable assumption as explained above and in Refs. [4, 22, 23]. Taking reasonable values of this proportionality factor, it was shown in Ref. [48] that $\chi_3^{(1)}(f > f_\alpha)$ is *the same* in the two models. This sheds a new light on the efficiency of the Box model and on consequence b).

2. Let us shortly discuss consequence a). In the Toy model, $\delta \ln \tau$ 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 $\propto \sqrt{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 $\delta \ln \tau$ because dissipation depends only on $\chi_1''$, which in the Toy model does not depend on $N_{corr}$. In the Toy model, as in the Pragmatical model [60] and the Diezemann model [61], the heating is neglected because at the scale of a given DH it is vanishingly small as shown above when discussing $\tau_{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 [11]. This purely exogeneous effect can be precisely calculated by solving the heat propagation equation, see e.g. ref.[11] and Appendix 2, 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 Refs. [70, 71], 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_{2,1}^{(1)}|$ and of $|X_3^{(3)}|$ would be mainly due to the entropy effect, contrarily to the hump of $|X_3^{(1)}|$ which would be due to the Box model contribution. As noted in Ref. [48], this means that very different 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 Eqs. (9) and (10): why should $X_3^{(1)}$ and $X_3^{(3)}$ have the same phase at high frequencies if their physical origin is different?

This is why it was emphasized in Ref. [48] that there is no reason for such a similarity if the growth of $X_3^{(1)}$ and $X_3^{(3)}$ are due to independent mechanisms. Because entropic effects have been related to the increase of $N_{corr}$ -see Eqs. (17) and (18)-, everything becomes instead very natural if the Box model is recasted in a frame-
work where $X_3^{(1)}$ is related to the glassy correlation volume. As evoked above, a first step in this direction was done in Ref. [48] where it was shown that the Box model prediction for $X_3^{(1)}$ at high frequencies is identical to the above Toy model prediction, provided $N_{corr}$ and $T\chi_T$ are proportional. In all, it is argued in Ref. [48] that the only reasonable way to account for the similarity of all three cubic susceptibilities, demonstrated experimentally in Figs. 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.

## 5 Conclusions

We have reviewed in this chapter 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 experimental findings are finally consistent with the existence of an underlying 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 [79, 5] 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 differences 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 Kinetically Constrained Model family will be extremely welcome to go beyond the general arguments given in Refs. [50]. Moreover, it would be very interesting to access $\chi_3$ (and $\chi_5$) in molecular liquids at higher temperatures, closer to the Mode Coupling Transition temperature $T_{MCT}$, and/or for frequencies close to the fast $\beta$ process where more complex, fractal structures with $d_f < d$ may be anticipated [80, 81]. This will require a joined effort 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 non-linear 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.

**ACKNOWLEDGEMENTS** We thank C. Alba-Simionesco, Th. Bauer, U. Buchenau, A. Coniglio, G. Johari, K. Ngai, R. Richert, G. Tarjus, and M. Tarzia for interesting discussions. The work in Saclay has been supported by the Labex RTRA grant Aricover and by the Institut des Systèmes Complexes ISC-PIF. The work in Augsburg was supported by the Deutsche Forschungsgemeinschaft via Research Unit FOR1394.

# 6 Appendix 1: making sure that exogeneous effects are negligible

We briefly explain how the nonlinear effects reported here have been shown to be -mainly- free of exogeneous effects:

1. The global homogeneous heating of the samples by the dielectric energy dissipated by the application of the strong ac field $E$ was shown to be fully negligible for $X_3^{(3)}$ as long as the inverse of the relaxation time $f_\alpha$ is $\le 1$ kHz, see Ref. [11]. Note that these homogeneous heating effects contribute much more to $X_3^{(1)}$: to minimize them, one can either keep $f_\alpha$ below 10 Hz [4], and/or severely limit the number $n$ of periods during which the electric field is applied -see, e.g., [28, 78]).

2. The contribution of electrostriction was demonstrated to be safely negligible in Refs. [74, 4], both by using theoretical estimates and by showing that changing the geometry of spacers does not affect $X_3^{(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., Figs. 8 and 9-. On the one hand it was shown that the ion heating contribution is fully negligible in $X_{2,1}^{(1)}$ (see Ref. [19]), on the other hand it is well known that ions affect the linear response $\chi_1$ at very low frequencies (say $f/f_\alpha \le 0.05$): this yields an upturn on the out-of-phase linear response $\chi_1''$, which diverges as $1/\omega$ instead of vanishing as $\omega$ in an ideally pure liquid containing only molecular dipoles. This may be the reason why most of the $\chi_3$ measurements are reported above $0.01 f_\alpha$: at lower frequencies the nonlinear responses is likely to be dominated by the ionic contribution.
## 7 Appendix 2: Trivial third and fifth harmonic susceptibilities

As explained in the main text, in the long time limit -i.e. for $f/f_\alpha \ll 1$, the liquid flow destroys the glassy correlations, making each molecule effectively 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 $\tau_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 $\Delta\chi_1/(1 - i\omega\tau_D)$, hence the subscript “Debye” in the Eq. (22) below. By using Refs. [17], and following the definitions given in the main text, as well as Eqs. (5)-(8) above, one gets for the dimensionless nonlinear responses of such an ideal gas, setting for brevity $x = \omega\tau_D$:

$$ 
\begin{aligned}
X_{3,Debye}^{(3)} &= \left(\frac{-3}{5}\right) \frac{3 - 17x^2 + ix(14 - 6x^2)}{(1 + 9x^2)(9 + 4x^2)(1 + 9x^2)} \\
X_{5,Debye}^{(5)} &= \frac{432(72 - 2377x^2 - 1979x^4 + 2990x^6)}{1680(1 + x^2)(4 + x^2)(9 + 4x^2)(1 + 9x^2)(9 + 16x^2)(1 + 25x^2)} \\
&\quad +i \frac{432x(246 - 737x^2 - 1623x^4 + 200x^6)}{560(1 + x^2)(4 + x^2)(9 + 4x^2)(1 + 9x^2)(9 + 16x^2)(1 + 25x^2)}
\end{aligned}
\quad (22) 
$$

In Ref. [50] the trivial response combined the above $X_{k,Debye}^{(k)}$ with a distribution $\mathcal{G}(\tau)$ of relaxation times $\tau$ chosen to account for the linear susceptibility of the supercooled liquid of interest. In Refs [58, 19] a slightly different modelization was used since $\mathcal{G}(\tau)$ was replaced by the Dirac delta function $\delta(\tau - \tau_\alpha)$, i.e. $\tau_D$ was simply replaced by $\tau_\alpha$ for the cubic trivial susceptibilities.

## 8 Appendix 3: Derivation of the Toy model from Langevin Fokker-Planck considerations

In this section we shall rederive the phenomenological Toy model of Ladieu et al. [58] 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 Ref. [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 $\tau_\alpha$, such interactions lose their coherence and the typical behavior of the liquid is that of an ideal gas. Before and around $\tau_\alpha$, these interactions manifest themselves in a frequency range $\omega \approx 1/\tau_\alpha$. Thus, stricto sensu, our modelling of this specific pro-
cess 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 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 $\mu_d = \mu \sqrt{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 $\mathbf{q}$ at temperature T, viz.

$$ \mathbf{q}(t) = \{q_1(t), \dots, q_n(t)\} $$

Inside each DH, each elementary constituent is assumed to interact via a multidimensional interaction potential $V_{int}(\mathbf{q})$ that possesses a double-well structure with minima at $\mathbf{q}_A$ and $\mathbf{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.

$$ \dot{q}_i = -\frac{1}{\zeta} \frac{\partial V_T}{\partial q_i}(\mathbf{q}, t) + \Xi_i(t) \quad (23) $$

where $\zeta$ is a generalized friction coefficient, $V_T = V_{int} + V_{ext}$, $V_{ext}$ is the potential energy of externally applied forces and the generalized forces $\Xi_i(t)$ have Gaussian white noise properties, namely

$$ \overline{\Xi_i(t)} = 0, \overline{\Xi_i(t) \Xi_j(t')} = \frac{2kT}{\zeta} \delta_{ij} \delta(t-t') \quad (24) $$

Thus, the dynamics of a DH is represented by the stochastic differential equations (23) and (24), which are in effect 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(\mathbf{q}, t)$ to find the system in state $\mathbf{q}$ at time $t$ which corresponds to Eqs. (23) and (24), namely

$$ \begin{aligned} \frac{\partial W}{\partial t}(\mathbf{q}, t) &= \frac{1}{2\tau_c} \nabla \cdot [\nabla W(\mathbf{q}, t) + \beta W(\mathbf{q}, t) \nabla V_T(\mathbf{q}, t)] \\ &= L_{FP}(\mathbf{q}, t) W(\mathbf{q}, t) \end{aligned} \quad (25) $$

where $2\tau_c = \zeta/(kT)$ is the characteristic time of fluctuations, $\nabla$ is the del operator in $\mathbf{q}$ space, and $L_{FP}(\mathbf{q}, t)$ is the Fokker-Planck operator. We notice that Eq. (25) may also be written
$$ \frac{\partial W}{\partial t}(\mathbf{q}, t) = \frac{1}{2\tau_c} \nabla \cdot \left\{ e^{-\beta V_T(\mathbf{q}, t)} \nabla \left[ W(\mathbf{q}, t) e^{\beta V_T(\mathbf{q}, t)} \right] \right\} \quad (26) $$

Now we use the transformation [84]

$$ \phi(\mathbf{q}, t) = W(\mathbf{q}, t) e^{\beta V_T(\mathbf{q}, t)} \quad (27) $$

so that Eq. (26) becomes

$$ \begin{aligned} \frac{\partial \phi}{\partial t}(\mathbf{q}, t) - \beta \frac{\partial V_T}{\partial t}(\mathbf{q}, t) \phi(\mathbf{q}, t) &= \frac{1}{2\tau_c} e^{\beta V_T(\mathbf{q}, t)} \nabla \cdot \left\{ e^{-\beta V_T(\mathbf{q}, t)} \nabla \phi(\mathbf{q}, t) \right\} \\ &= L_{FP}^{\dagger}(\mathbf{q}, t) \phi(\mathbf{q}, t) \end{aligned} \quad (28) $$

where $L_{FP}^{\dagger}(\mathbf{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 Eq. (28) means that $W$ is near its equilibrium value, so restricting further calculations to low frequencies, $\omega\tau_c \ll 1$ (quasi-stationary condition). Hence, Eq. (28) now reads

$$ \frac{\partial \phi}{\partial t}(\mathbf{q}, t) \approx L_{FP}^{\dagger}(\mathbf{q}, t) \phi(\mathbf{q}, t) \quad (29) $$

Now, the interpretation of the Fokker-Planck equation (25) (or equally well the Langevin equations (23)) with time-dependent potential in terms of usual population equations with time-dependent rate coefficients has a meaning, since now Eq. (27) 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) = \rho_0 \mu_d \int \cos \vartheta(\mathbf{q}) W(\mathbf{q}, t) d\mathbf{q} \quad (30) $$

where $\rho_0$ is the number of DH per unit volume, and $\vartheta(\mathbf{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 Eq. (30)

$$ P(t) = \rho_0 \mu_d \left[ \int_{\text{well A}} \cos \vartheta(\mathbf{q}) W(\mathbf{q}, t) d\mathbf{q} + \int_{\text{well B}} \cos \vartheta(\mathbf{q}) W(\mathbf{q}, t) d\mathbf{q} \right] \quad (31) $$

Now, it is known from the Kramers theory of chemical reaction rates [84] that at sufficiently large energy barriers, most of the contributions of the integrands come from the minima of the wells, therefore we have

$$ P(t) \approx \rho_0 \mu_d \left[ \cos \vartheta(\mathbf{q}_A) \int_{\text{well A}} W(\mathbf{q}, t) d\mathbf{q} + \cos \vartheta(\mathbf{q}_B) \int_{\text{well B}} W(\mathbf{q}, t) d\mathbf{q} \right] \quad (32) $$
Now, the integrals represent the relative populations $x_i(t) = n_i(t)/N$, $i = A, B$ in each well (we assume that $W(\mathbf{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 \qquad (33)$$

Thus, Eq. (32) reads

$$P(t) \approx \rho_0 \mu_d [\cos \vartheta(\mathbf{q}_A) x_A(t) + \cos \vartheta(\mathbf{q}_B) x_B(t)] \qquad (34)$$

We assume now for simplicity that $\vartheta(\mathbf{q}_B) = \pi - \vartheta(\mathbf{q}_A)$, so that

$$P(t) \approx \rho_0 \mu_d \cos \vartheta(\mathbf{q}_A) [x_A(t) - x_B(t)] \qquad (35)$$

Finally, since $\rho_0 = N/V$ where $V$ is the volume of the polar substance made of DH only, we obtain

$$P(t) \approx \frac{\mu_d \cos \vartheta(\mathbf{q}_A)}{N v_{DH}} [n_A(t) - n_B(t)] \qquad (36)$$

where $v_{DH}$ is the volume of a DH. This is the definition of the polarization in the Toy model.

In order to determine the polarization (36), we need to calculate the dynamics of $n_i(t)$. From the conservation law -Eq. (33)-, we have

$$\dot{x}_A(t) = -\dot{x}_B(t) = \int_{\text{well A}} \frac{\partial W}{\partial t}(\mathbf{q},t) d\mathbf{q} \qquad (37)$$

By using the Fokker-Planck equation (26) and limiting well A to a closed generalized bounding surface constituting the saddle region $\partial A$, we have by Gauss's theorem

$$\dot{x}_A(t) = -\dot{x}_B(t) = \frac{1}{2\tau_c} \int_{\partial A} e^{-\beta V_T(\mathbf{q},t)} \nabla\phi(\mathbf{q},t) \cdot \mathbf{v}_{\mathbf{q}} dS_{\mathbf{q}} \qquad (38)$$

where $v_{\mathbf{q}}$ is the outward normal to the bounding surface and $dS_{\mathbf{q}}$ is a generalized surface element of the bounding surface, and where we have used Eq. (27). Now, we follow closely Coffey et al. [85] and introduce the crossover function $\Delta(\mathbf{q},t)$ via the equation

$$\phi(\mathbf{q},t) = \phi_A(t) + [\phi_B(t) - \phi_A(t)] \Delta(\mathbf{q},t) \qquad (39)$$

where $\Delta(\mathbf{q},t) = 0$ if $\mathbf{q} \in \text{well A}$ while $\Delta(\mathbf{q},t) = 1$ if $\mathbf{q} \in \text{well B}$ and exhibits strong gradients in the saddle region $\partial A$ allowing the crossing from A to B (and vice-versa) by thermally activated escape. By combining Eqs. (38) and (39), we have immediately
$$ \dot{x}_A(t) = -\dot{x}_B(t) = \frac{\phi_B(t) - \phi_A(t)}{2\tau_c} \oint_{\partial A} e^{-\beta V_T(\mathbf{q},t)} \nabla[\Delta(\mathbf{q},t)] \cdot \mathbf{v}_{\mathbf{q}} dS_{\mathbf{q}} \quad (40) $$

Now,

$$ x_i(t) = \phi_i(t) x_i^s(t), \quad x_i^s(t) = \int_{\text{well i}} W_s(\mathbf{q}, t) d\mathbf{q} \quad (41) $$

where $W_s(\mathbf{q},t)$ is a normalized solution of the Fokker-Planck equation

$$ L_{FP}(\mathbf{q}, t) W_s(\mathbf{q}, t) = 2 \tau_c \frac{\partial W_s}{\partial t}(\mathbf{q}, t) \approx 0 \quad (42) $$

because the frequencies we are concerned with are very small with respect to the inverse thermal fluctuation time $\tau_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

$$ x_A^s(t) + x_B^s(t) = 1 \quad (43) $$

Using Eqs. (41) and (43), we may easily show that [85]

$$ \phi_B(t) - \phi_A(t) = \left( \frac{1}{x_A^s(t)} + \frac{1}{x_B^s(t)} \right) [x_B(t) x_A^s(t) - x_A(t) x_B^s(t)] \quad (44) $$

By combining Eqs. (38) and (44), we readily obtain

$$ \dot{x}_A(t) = -\dot{x}_B(t) = \Gamma(t) (x_B(t) x_A^s(t) - x_A(t) x_B^s(t)) \quad (45) $$

where the overall time-dependent escape rate $\Gamma(t)$ is given by [85]

$$ \Gamma(t) = \frac{1}{2\tau_c} \left( \frac{1}{x_A^s(t)} + \frac{1}{x_B^s(t)} \right) \oint_{\partial A} e^{-\beta V_T(\mathbf{q},t)} \nabla\phi(\mathbf{q},t) \cdot \mathbf{v}_{\mathbf{q}} dS_{\mathbf{q}} \quad (46) $$

Finally, by setting

$$ \Pi_{AB}(t) = \Gamma(t) x_B^s(t), \quad \Pi_{BA}(t) = \Gamma(t) x_A^s(t) \quad (47) $$

we arrive at the population equations

$$ \dot{n}_A(t) = -\dot{n}_B(t) = -\Pi_{AB}(t) n_A(t) + \Pi_{BA}(t) n_B(t) \quad (48) $$

The obtaining of a more explicit formula for the various rates involved in Eq. (47) is not possible, due to the impossibility to calculate the surface integral in Eq. (46) explicitly, in turn due to the fact that $V_T$ is not known explicitly. Then, the rates in Eqs. (47) and (48) are estimated using Arrhenius's formula. All subsequent derivations regarding the Toy model of Ladieu et al. [58] follow immediately and will not be repeated here due to lack of room and straightforward but laborious algebra.
References

1. L.-P. Lévy, and A. T. Oglieski *Nonlinear dynamic susceptibilities at the spin-glass transition of Ag:Mn*, Phys. Rev. Lett. **57**, 3288 (1986).

2. J. Hemberger, H. Ries, A. Loidl, and R. Böhmer, *Static Freezing Transition at a Finite Temperature in a Quasi-One-Dimensional Deuteron Glass*, Phys. Rev. Lett. **76**, 2330 (1996).

3. C. Crauste-Thibierge, C. Brun, F. Ladieu, D. L'Hôte, G. Biroli, and J.-P. Bouchaud, *Evidence of Growing Spatial Correlations at the Glass Transition from Nonlinear Response Experiments*, Phys. Rev. Lett. **104**, 165703 (2010).

4. C. Brun, F. Ladieu, D. L'Hôte, M. Tarzia, G. Biroli, and J.-P. Bouchaud, *Nonlinear dielectric susceptibilities: Accurate determination of the growing correlation volume in a supercooled liquid*, Phys. Rev. B **84**, 104204 (2011).

5. J. M. Brader, M. Siebenbürger, M. Ballauff, K. Reinheimer, M. Wilhelm, S. J. Frey, F. Weysser, and M. Fuchs, *Nonlinear response of dense colloidal suspensions under oscillatory shear: Mode-coupling theory and Fourier transform rheology experiments*, Phys. Rev. E **82**, 061401 (2010).

6. R. Seyboldt, D. Merger, F. Coupette, M. Siebenbürger, M. Ballauff, M. Wilhelm, and M. Fuchs, *Divergence of the third harmonic stress response to oscillatory strain approaching the glass transition*, Soft Matter, **12**, 8825 (2016).

7. K. Binder and A. P. Young, *Spin glasses: experimental facts, theoretical concepts, an dopen questions*, Rev. Mod. Phys. **58**, 801 (1986).

8. U. T. Höchli, K. Knorr, and A. Loidl, *Orientational glasses*, Adv. Phys. **39**, 405 (1990).

9. J.-P. Bouchaud, and G. Biroli, *Nonlinear susceptibility in glassy systems: A probe for cooperative dynamical length scales*, Phys. Rev. B **72**, 064204 (2005).

10. Th. Bauer, P. Lunkenheimer, and A. Loidl, *Cooperativity and the Freezing of Molecular Motion at the Glass Transition*, Phys. Rev. Lett. **111**, 225702 (2013).

11. C. Brun, C. Crauste-Thibierge, F. Ladieu, and D. L'Hôte, *Study of the heating effect contribution to the nonlinear dielectric response of a supercooled liquid*, J. Chem. Phys. **133**, 234901 (2010).

12. C. Gainaru, R. Meier, S. Schildmann, C. Lederle, W. Hiller, E. A. Rössler, and R. Böhmer, *Nuclear-magnetic-resonance measurements reveal the origin of the Debye process in mono-hydroxy alcohols*, Phys. Rev. Lett. **105**, 258303 (2010).

13. Th. Bauer, M. Michl, P. Lunkenheimer, and A. Loidl, *Nonlinear dielectric response of Debye, α, and β relaxation in 1-propanol*, J. Non-Cryst. Solids **407**, 66 (2015).

14. L. P. Singh and R. Richert, *Watching hydrogen-bonded structures in an alcohol convert from rings to chains*, Phys. Rev. Lett. **109**, 167802 (2012).

15. L. P. Singh, C. Alba-Simionesco, and R. Richert, *Dynamics of glass-forming liquids. XVII. Dielectric relaxation and intermolecular association in a series of isomeric octyl alcohols*, J. Chem. Phys. **139**, 144503 (2013).

16. P. Lunkenheimer, M. Michl, Th. Bauer, and A. Loidl, *Investigation of nonlinear effects in glassy matter using dielectric methods*, Eur. Phys. J. Special Topics **226**, 3157 (2017).

17. W.T. Coffey, and B.V. Paranjape, Proc. Roy. Ir. Acad. **78 A**, 17 (1978); J. L. Déjardin and Yu. P. Kalmykov, *Nonlinear dielectric relaxation of polar molecules in a strong ac electric field: Steady state response*, Phys. Rev. E **61**, 1211 (2000).

18. M. Michl, Th. Bauer, and A. Loidl, *Nonlinear Dielectric Spectroscopy in a Fragile Plastic Crystal*, J. Chem. Phys. **144**, 114506 (2016).

19. D. L'Hôte, R. Tourbot, F. Ladieu, and P. Gadige, *Control parameter for the glass transition of glycerol evidenced by the static-field-induced nonlinear response*, Phys. Rev. B **90**, 104202 (2014).
20. R. Casalini, D. Fragiadakis, and C. M. Roland, *Dynamic correlation length scales under isochronal conditions*, J. Chem. Phys. **142**, 064504 (2015).

21. G. Adam and J. H. Gibbs, *On the Temperature Dependence of Cooperative Relaxation Properties in Glass-Forming Liquids*, J. Chem. Phys. **43**, 139 (1965).

22. L. Berthier, G. Biroli, J.-P. Bouchaud, L. Cipeletti, D. El Masri, D. L'Hôte, F. Ladieu, and M. Pierno, *Direct Experimental Evidence of a Growing Length Scale Accompanying the Glass Transition*, Science **310**, 1797 (2005).

23. C. Dalle-Ferrier, C. Thibierge, C. Alba-Simionesco, L. Berthier, G. Biroli, J.-P. Bouchaud, F. Ladieu, D. L'Hôte, and G. Tarjus, *Spatial correlations in the dynamics of glassforming liquids: Experimental determination of their temperature dependence*, Phys. Rev. E **76**, 041510 (2007); L. Berthier, G. Biroli, J.-P. Bouchaud, W. Kob, K. Miyazaki, and D. R. Reichman, *Spontaneous and induced dynamic fluctuations in glass formers*. I. General results and dependence on ensemble and dynamics, J. Chem. Phys. **126**, 184503 (2007).

24. R. Brand, P. Lunkenheimer, and A. Loidl, *Relaxation dynamics in plastic crystals*, J. Chem. Phys. **116**, 10386 (2002).

25. P. Lunkenheimer, U. Schneider, R. Brand, and A. Loidl, *Glassy dynamics*, Contemp. Phys. **41**, 15 (2000).

26. U. Schneider, P. Lunkenheimer, R. Brand, and A. Loidl, *Dielectric and far-infrared spectroscopy of glycerol*, J. Non-Cryst. Solids **235-237**, 173 (1998).

27. Th. Bauer, P. Lunkenheimer, S. Kastner, and A. Loidl, *Nonlinear dielectric response at the excess wing of glass-forming liquids*, Phys. Rev. Lett. **110**, 107603 (2013).

28. R. Richert and S. Weinstein, *Nonlinear Dielectric Response and Thermodynamic Heterogeneity in Liquids*, Phys. Rev. Lett. **97**, 095703 (2006).

29. S. Samanta and R. Richert, *Limitations of heterogeneous models of liquid dynamics: Very slow rate exchange in the excess wing*, J. Chem. Phys. **140**, 054503 (2014).

30. Th. Bauer, *Nichtlineare dielektrische Spektroskopie zum Nachweis von Kooperativitt und Heterogenitt in glasbildenden Flissigkeiten* (Verlag Dr. Hut, Munich, 2015).

31. U. Schneider, R. Brand, P. Lunkenheimer, and A. Loidl, *Excess wing in the dielectric loss of glass formers: A Johari-Goldstein β relaxation?*, Phys. Rev. Lett. **84**, 5560 (2000).

32. A. Döß, M. Paluch, H. Sillescu, and G. Hinze, *From strong to fragile glass formers: Secondary relaxation in polyalcohols*, Phys. Rev. Lett. **88**, 095701 (2002).

33. M. Beiner, H. Huth, and K. Schrter, *Crossover region of dynamic glass transition: general trends and individual aspects*, J. Non-Cryst. Solids **279**, 126 (2001).

34. K. L. Ngai, *An extended coupling model description of the evolution of dynamics with time in supercooled liquids and ionic conductors*, J. Phys.: Condens. Matter **15**, S1107 (2003).

35. K. L. Ngai, *Interpreting the nonlinear dielectric response of glass-formers in terms of the coupling model*, J. Chem. Phys. **142**, 114502 (2015).

36. L. Berthier and G. Biroli, *Theoretical perspective on the glass transition and amorphous materials*, Rev. Mod. Phys. **83**, 587 (2011).

37. A. Q. Tool, *Relation between inelastic deformability and thermal expansion of glass in its annealing range*, J. Am. Ceram. Soc. **29**, 240 (1946).

38. O. S. Narayanaswamy, *Model of structural relaxation in glass*, J. Am. Ceram. Soc. **54**, 491 (1971).

39. C. T. Moynihan, A. J. Easteal, M. A. Bolt, and J. Tucker, *Dependence of fictive temperature of glass on cooling rate*, J. Am. Ceram. Soc. **59**, 12 (1976).

40. V. Lubchenko and P. G. Wolynes, *Theory of aging in structural glasses*, J. Chem. Phys. **121**, 2852 (2004).

41. F. Bert, V. Dupuis, E. Vincent, J. Hammann, and J.-P. Bouchaud, *Spin anisotropy and slow dynamics in spin glasses*, Phys. Rev. Lett. **92**, 167203 (2004).

42. F. Baity-Jesi et al., *Matching Microscopic and Macroscopic Responses in Glasses*, Phys. Rev. Lett. **118**, 157202 (2017).

43. R. L. Leheny and S. R. Nagel, *Frequency-domain study of physical aging in a simple liquid*, Phys. Rev. B **57**, 5154 (1998).

44. P. Lunkenheimer, R. Wehn, U. Schneider, and A. Loidl, *Glassy aging dynamics*, Phys. Rev. Lett. **95**, 055702 (2005).
45. C. Brun, F. Ladieu, D. L'Hôte, G. Biroli, and J.-P. Bouchaud, Evidence of growing spatial correlations during the aging of glassy glycerol, Phys. Rev. Lett. **109**, 175702 (2012).

46. M. Dzero, J. Schmalian, and P. G. Wolynes, in *Structural glasses and supercooled Liquids: theory, experiment, and applications*, P. G. Wolynes, V. Lubchenko, Eds. (Wiley, 2012), pp. 193-222; G. Biroli, and J.-P. Bouchaud, in *Structural glasses and supercooled liquids: theory, experiment, and applications*, P. G. Wolynes, V. Lubchenko, Eds. (Wiley, 2012), pp. 31-114; G. Biroli, and J.-P. Bouchaud, *The Random First-Order Transition Theory of Glasses: a critical assessment*, arXiv:0912.2542v1.

47. C. Cammarota, A. Cavagna, G. Gradenigo, T. S. Grigera and P. Verrocchio, *Numerical determination of the exponents controlling the relationship between time, length, and temperature in glass-forming liquids*, J. Chem. Phys. **131**, 194901 (2009).

48. P. Gadige, S. Albert, M. Michl, Th. Bauer, P. Lunkenheimer, A. Loidl, R. Tourbot, C. Wiertel-Gasquet, G. Biroli, J.-P. Bouchaud, and F. Ladieu, *Unifying different interpretations of the nonlinear response in glass-forming liquids*, Phys. Rev. E **96**, 032611 (2017).

49. Y. Kimura, S. Hara, and R. Hayakawa, *Nonlinear dielectric relaxation spectroscopy of ferroelectric liquid crystals*, Phys. Rev. E **62**, R5907-R5910 (2000).

50. S. Albert, Th. Bauer, M. Michl, G. Biroli, J.-P. Bouchaud, A. Loidl, P. Lunkenheimer, R. Tourbot, C. Wiertel-Gasquet, and F. Ladieu, *Fifth-order susceptibility unveils growth of thermodynamic amorphous order in glass-formers*, Science **352**, 1308 (2016).

51. L. P. Levy, *Critical-dynamics of metallic spin glasses*, Phys. Rev. B **38**, 4963-4973 (1988).

52. This comes from the non zero value of the configurational entropy in supercooled liquids. This configurational entropy does not exist in the case of spin glasses, hence there is a finite probability for an amorphously ordered domain to come back to a given configuration even after having relaxed. As a result, instead of a hump for $X_3$ around $f_\alpha$, there is a plateau for $X_3(f \le f_\alpha) \simeq X_3(0)$, the amplitude of which diverges when approaching the spin glass transition -see L. Lévy, Phys. Rev. B **38**, 4963 (1988).

53. It would be worth to develop further the MCT analysis of [54] to obtain detailed analytical predictions on the phases.

54. M. Tarzia, G. Biroli, A. Lefèvre, J.-P. Bouchaud, *Anomalous nonlinear response of glassy liquids: General arguments and a mode-coupling approach*, J. Chem. Phys. **132**, 054501 (2010).

55. The growth of $\chi_3$ and $\chi_5$ can be induced by purely dynamical correlation, as it is the case in MCT theory and in the MCT regime. However, in the activated dynamics regime, the only theory that can explain a growth of $\chi_3$ and $\chi_5$ are the ones in which dynamical correlations are driven by static ("point-to-set") correlations.

56. D. Chandler, and J. P. Garrahan, *Dynamics on the Way to Forming Glass: Bubbles in Space-Time*, Annu. Rev. Phys. Chem. **61**, 191-217 (2010).

57. G. Tarjus, S. A. Kivelson, Z. Nussinov, and P. Viot, *The frustration-based approach of supercooled liquids and the glass transition: a review and critical assessment*, J. Phys: Cond. Matt. **17**, R1143-R1182 (2005).

58. F. Ladieu, C. Brun, D. L'Hôte, *Nonlinear dielectric susceptibilities in supercooled liquids: A toy model*, Phys. Rev. B **85**, 184207 (2012).

59. U. Buchenau, *Retardation and flow at the glass transition*, Phys. Rev. E, **93**, 032608 (2016).

60. U. Buchenau, *Pragmatical access to the viscous flow of undercooled liquids.*, Phys. Rev. E **95**, 062603 (2017).

61. G. Diezemann, *Nonlinear response theory for Markov processes: Simple models for glassy relaxation*, Phys. Rev. E **85**, 051502 (2012).

62. G. Diezemann, *Nonlinear response theory for Markov processes II: Fifth-order response functions*, Phys. Rev. E **96**, 022150 (2017).

63. G. P. Johari, *Effects of electric field on the entropy, viscosity, relaxation time, and glass-formation*, J. Chem. Phys. **138**, 154503 (2013).

64. G. P. Johari, *Effects of electric field on thermodynamics and ordering of a dipolar liquid*, J. Chem. Phys. **145**, 164502 (2016).

65. S. Samanta, and R. Richert, *Dynamics of glass-forming liquids. XVIII. Does entropy control structural relaxation times?*, J. Chem. Phys. **142**, 044504 (2015).
66. A. R. Young-Gonzales, S. Samanta, and R. Richert, *Dynamics of glass-forming liquids*. *XIX*. *Rise and decay of field induced anisotropy in the non-linear regime*, J. Chem. Phys. **143**, 104504 (2015).

67. B. Riechers, K. Samwer, and R. Richert, *Structural recovery in plastic crystals by time-resolved non-linear dielectric spectroscopy*, J. Chem. Phys. **142**, 154504 (2015).

68. S. Samanta and R. Richert, *Electrorheological Source of Nonlinear Dielectric Effects in Molecular Glass-Forming Liquids*, J. Phys. Chem. B **120**, 7737 (2016).

69. S. Samanta, Yamamuro, and R. Richert, *Connecting thermodynamics and dynamics in a supercooled liquid: Cresolphthalein-dimethylether*, Thermochim Acta **636**, 57 (2016).

70. R. Richert, *Non-linear dielectric signatures of entropy changes in liquids subject to time dependent electric fields*, J. Chem. Phys. **144**, 114501 (2016).

71. P. Kim, A. R. Young-Gonzales, and R. Richert, *Dynamics of glass-forming liquids*. *XX*. *Third harmonic experiments of non-linear dielectric effects versus a phenomenological model*, J. Chem. Phys. **145**, 064510 (2016).

72. R. Richert and C. A. Angell, *Dynamics of glass-forming liquids*. *V*. *On the link between molecular dynamics and configurational entropy*, J. Chem. Phys. **108**, 9016 (1998).

73. B. Schiener, R. Böhmer, A. Loidl, and R. V. Chamberlin, *Nonresonant Spectral Hole Burning in the Slow Dielectric Response of Supercooled Liquids*, Science **274**, 752 (1996); B. Schiener, R. V. Chamberlin, G. Diezemann, and R. Böhmer, *Nonresonant dielectric hole burning spectroscopy of supercooled liquids*, J. Chem. Phys. **107**, 7746 (1997).

74. S. Weinstein and R. Richert, *Nonlinear features in the dielectric behavior of propylene glycol*, Phys. Rev. B **75**, 064302 (2007).

75. L.-M. Wang and R. Richert, *Measuring the Configurational Heat Capacity of Liquids*, Phys. Rev. Lett. **99**, 185701 (2007).

76. A. Khalife, U. Pathak, and R. Richert, *Heating liquid dielectrics by time dependent fields*, Eur. Phys. J. B, **83**, 429, (2011).

77. C. Brun, C. Crauste-Thibierge, F. Ladieu, and D. L'Hôte, *Third harmonics nonlinear susceptibility in supercooled liquids: A comparison to the box model*, J. Chem. Phys. **134**, 194507 (2011).

78. see Supplementary Information of [Th. Bauer, P. Lunkenheimer, S. Kastner, and A. Loidl, Phys. Rev. Lett. **110**, 107603 (2013).].

79. E. Lippiello, F. Corberi, A. Sarracino, and M. Zannetti, *Nonlinear response and fluctuation-dissipation relations*, Phys. Rev. E **78**, 041120 (2008); E. Lippiello, F. Corberi, A. Sarracino, and M. Zannetti, *Nonlinear susceptibilities and the measurement of a cooperative length*, Phys. Rev. E **77**, 212201 (2008).

80. G. Biroli, J.-P. Bouchaud, K. Miyazaki, and D. R. Reichman, *Inhomogeneous Mode-Coupling Theory and Growing Dynamic Length in Supercooled Liquids*, Phys. Rev. Lett. **97**, 195701 (2006).

81. J. D. Stevenson, J. Schmalian, and P. G. Wolynes, *The shapes of cooperatively rearranging regions in glass-forming liquids*, Nature Physics **2**, 268-274 (2006).

82. P.M. Dejardin and F. Ladieu, unpublished work.

83. H. Risken, *The Fokker-Planck Equation*, 2nd Edition (Springer, Berlin, 1989).

84. H. A. Kramers, *Brownian Motion in a field of force and the diffusion theory of chemical reactions*, Physica **7**, 284 (1940).

85. W.T. Coffey, D.A. Garanin, and D.J. McCarthy, *Crossover Formulas in the Kramers Theory of Thermally Activated Escape Rates: Application to spin systems*, Adv. Chem. Phys. **117**, 483 (2001).