Stanimir S Iliev

Nina S Pesheva

Vadim S Nikolayev
email: vadim.nikolayev@cea.fr

Contact angle hysteresis and pinning at periodic defects in statics

Keywords: number(s): 68.03.Cd, 05.90.+m, 68.08.Bc

This article deals with the theoretical prediction of the wetting hysteresis on nonideal solid surfaces in terms of the surface heterogeneity parameters. The spatially periodical chemical heterogeneity is considered. We propose precise definitions for both the advancing and the receding contact angles for the Wilhelmy plate geometry. It is well known that in such a system, a multitude of metastable states of the liquid meniscus occurs for each different relative position of the defect pattern on the plate with respect to the liquid level. As usual, the static advancing and receding angles are assumed to be a consequence of the preceding contact line motion in the respective direction. It is shown how to select the appropriate states among all metastable states. Their selection is discussed. The proposed definitions are applicable to both the static and the dynamic contact angles on heterogeneous surfaces. The static advancing and receding angles are calculated for two examples of periodic heterogeneity patterns with sharp borders: the horizontal alternating stripes of a different wettability (studied analytically) and the doubly periodic pattern of circular defects on a homogeneous base (studied numerically). The wetting hysteresis is determined as a function of the defect density and the spatial period. A comparison with the existing results is carried out.

I. INTRODUCTION

There are yet many open questions in the wetting statics and dynamics [1] in spite of the strong interest in these topics by both the scientific community and industry. One of the open issues concerns the wetting on heterogeneous and/or rough solid surfaces. The surface heterogeneities that we call hereafter "defects" are well known to be responsible for the stick-slip dynamics of the contact line (CL) and the nonuniqueness of the static contact angle (CA), which is a phenomenon responsible for the appearance of the wetting hysteresis. It is ubiquitous for all real surfaces [2] and is thus an important problem.

The main objective of the theoretical description of wetting is a prediction of behavior on macroscopic length scale (e.g., the average meniscus curvature or drop size) in terms of the heterogeneity parameters, defined on a much smaller length scale.

The Cassie-Baxter relation for the averaged cosine of the static CA was one of the first attempts of such an averaging. However, the introduction of this unique averaged CA is not sufficient for characterizing heterogeneous surfaces due to the existence of wetting hysteresis. One needs to define other parameters to characterize the wetting on a heterogeneous and/or rough surfaces. The static advancing θ a and receding θ r contact angles (ACA and RCA) are used for this purpose. However, the definitions in the literature for the wetting hysteresis are not consistent (see Ref. [3], and references therein) and lead to ambiguous results.

The appearance of different values of the CA is related to the existence of both a multitude of metastable states of the liquid meniscus for a given position of the defect pattern on the plate with respect to the liquid and the different possible positions of the defect pattern on the plate with respect to the liquid, over which the averaging should be performed. A wetting hysteresis theory should give answers to the following questions: (i) how to determine the CA distribution along the CL for a specific system realization, (ii) how to select the advancing and receding CAs for the specific system realization from the set of all possible values of averaged CAs, (iii) how to choose from different suggested ways of averaging the CA (see Sec. III A), and (iv) how to average the CAs over the surface, due to the existence of many different possible relative positions of the defect pattern on the plate with respect to the liquid meniscus. The answer to these questions is not unique and there is still an ongoing discussion. Further work is needed to elucidate these problems.

An effort towards a more consistent theoretical description of hysteresis, i.e., definitions and clear procedure for their application is thus necessary. The objective of this study is to suggest precise expressions for these angles that can be used in both the static and the dynamic (for small capillary and Reynolds numbers) cases. We limit ourselves here to the study of chemical heterogeneity in statics; however, the approach, we use here is more general and is applicable for studying the system in contact with arbitrary heterogeneous and/or rough surfaces.

A large number of studies of CA hysteresis on heterogeneous surfaces exist. The hysteresis was calculated for different regular and random patterns of the solid surface by the continuum media energy minimization (see, e.g., [4][5][6][7][8][9][10]), lattice-Boltzmann method [11], phase-field models [12], and molecular dynamics approaches [13], to mention just a few.

We focus here on the Wilhelmy balance geometry (Fig. 1) used conventionally for the CA determination. The static ACA measurement in such a system may be performed by fixing the vertical position of the defect pattern on the plate with respect to the liquid level far from the plate, after having moved the plate downwards into the liquid pool. Similarly, the RCA may be obtained after the upward motion. The different possible relative positions of the defect pattern on the plate with respect to the liquid level far from the plate (see Fig. 1) are characterized by the values of a parameter H (detailed further in the text). For each H , there is a variety of metastable meniscus shapes (that we call states) that may be characterized by the CL height h (cf. Fig. 1), averaged along the CL (the angle brackets are used hereafter to denote such a "horizontal" averaging). It is well known [4,[14][15][16] that among them the most important for the hysteresis determination are the two states with the highest and the lowest possible values of h , beyond which no metastable meniscus state exists. We denote them with the "+" and "-" superscripts, respectively, so that the metastable states exist for h - h h + . The ACA is related to the "-" state, while RCA is defined by the "+" state. These states are the most important for the determination of the CA hysteresis.

The derived RCA and ACA expressions are then applied to two examples of heterogeneous pattern. For the first of them, the horizontal alternating stripes of a different wettability, the hysteresis is obtained analytically. For the second, the doubly periodic pattern of circular defects on a homogeneous base, two methods are used and their results compared: the direct numerical simulation method and the semianalytical small slope approximation. We compare our results with the existing hysteresis determination methods (in particular, with [17]) and consider the wetting hysteresis dependence on the concentration of periodically distributed defects.

II. PROBLEM FORMULATION

A flat heterogeneous solid plate is immersed vertically into a liquid pool. The plate vertical position may be varied, so that the defect pattern position H with respect to the liquid level in the pool is controlled. One of the plate faces (we do not consider the other) is described with the Cartesian coordinates (y,z) where the y axis is horizontal and the z axis is directed upward, as shown in Fig. 1, with zero coinciding with the liquid level far from the moving plate. The liquid free surface is denoted by lg and its surface tension by γ lg . The CL, which the liquid meniscus forms with the solid surface, is denoted by L and its height by h. We consider here the case when the free liquid surface lg is a single-valued function of x 0 and y: z = f (x,y); h = f (0,y). A CL point is then described by the vector R = [y,h]. The distance between the plate and the container wall, opposite to the plate, is denoted by d x , where d x l c [l c = √ γ lg /(ρ lρ g )g is the capillary length; the density of the liquid and of the ambient gas are denoted by ρ l and ρ g , respectively; g is the gravity acceleration]. The width of the plate is denoted by d y .

In the past, many studies (see, e.g., [9]) used the approach of Joanny and de Gennes [18] developed for smooth variation of solid surface energy. We focus here on the case on the "mesatype" defects when the plate surface consists of two kinds of domains (i) (i = 1,2) of different wettability with sharp borders between them. It is assumed that both domains are homogeneous, with no intrinsic wetting hysteresis. Inside the domain (i) , the solid-liquid interface energy γ sl (y,z) = γ (i) sl and the solid-gas surface energy γ sg (y,z) = γ (i) sg . These surface energies correspond to the Young CAs,

cos θ (i) eq = γ (i) sg -γ (i) sl γ lg . ( 1
)
Without any loss of generality, it is assumed in the following that θ (1) eq < θ (2) eq .

(

The domain pattern is assumed to be periodical, with spatial periods λ y and λ z in y and z directions, respectively (Fig. 1).

It is assumed that d y λ y . Here we concentrate on the case of mesoscopic scale heterogeneities. Because of their large length scale, thermal effects [1] are not expected to influence the hysteresis so that the zero-temperature case may be considered. The free energy U of the system is thus [START_REF] Landau | Fluid Mechanics[END_REF]

U = lg γ lg d + sl γ sl d + sg γ sg d + l ρ l gzd + g ρ g gzd , ( 3
)
where the liquid and gas volumes are denoted l and g , respectively. For a given choice of the surface pattern, the free energy (3) is uniquely determined by the position of the liquid meniscus and thus U is a functional of f (x,y). The meniscus metastable states (including those related to the advancing and receding contact angles) correspond to the local minima of the free energy (3). The minimization should thus be performed to answer the question (i) posed in Sec. I. For mesa type defects, the minimization leads [6] to the following conditions for the equilibrium local CA θ ( R), which is the meniscus slope at x = 0. For the points R ∈ L that belong to the interiors of domains, the following condition is valid:

θ ( R) =
θ (1) eq , if R ∈ (1) , θ (2) eq , if R ∈ (2) .

(

) 4
For the points R b ∈ L at the borders between the domains, the following inequality should be satisfied θ (1) eq θ ( R b ) θ (2) eq .

(

In other words, when the CL attains a point R b at a border between the domains, the CL remains pinned as long as the local CA satisfies the double inequality (5). In this case the equilibrium may occur only when the liquid situates at the side of the defect border, corresponding to the domain with the better wettability (domain 1).

III. DEFINITIONS OF ADVANCING AND RECEDING CONTACT ANGLES

Consider now the questions (iii) and (iv) of Sec. I, i.e., the CA averaging over the solid surface. It is more convenient to discuss separately the averaging over the horizontal (y) and the vertical (z) directions; cf. Fig. 1.

A. Horizontal averaging

After an equilibrium meniscus state is found for a fixed plate position, one can define the averaged CA. There are many possible ways to define averaged angles. One can, e.g., introduce the angle averaged along the CL L (or along the width of the plate), R ∈ L,

1 |L| L θ ( R)dl, ( 6
)
where |L| is the CL length. Other definitions of average angle have been also considered by different researchers. For example, a projection angle is introduced in [6] and an apparent angle is defined in [START_REF] Brandon | [END_REF] (see also [1]). Such definitions are basically similar to that of David and Neumann [9], who define the average angle θ by the expression

θ = arcsin 1 - h 2 2l 2 c , (7)
deduced from the equilibrium averaged over y CL height h via the well known relation between the CA θ and the CL height h on an ideally homogeneous wall [START_REF] Landau | Fluid Mechanics[END_REF]:

h = l c 2(1 -sin θ ). (8)
In the present study, we adopt another definition of averaging the CA, related to the values of the RCA and ACA, measured experimentally by the Wilhelmy balance method [21,22]. The average angle is deduced from the vertical component of the force, exerted by the liquid on the Wilhelmy plate, more specifically, from its "capillary" part, related to the presence of the CL. Note that we consider only one side of the plate (Fig. 1), while the force is exerted on both sides. The force per unit length of the CL L is γ lg cos θ ( R) m( R) [23], where the unit vector m( R), defined at R ∈ L, is normal to the CL L, is tangent to the liquid-solid interface, and points outward from the liquid. The z component of the total force per unit width of the plate is thus [10,22]

F = γ lg |L| L cos θ ( R) m( R) • zdl, ( 9
)
where z is the unit vector directed along the z axis. A change of variables allows us to write this as

F = γ lg d y d y 0 cos θ (y,h(y))dy. (10)
It is evident now that the ratio F /γ lg defines the cosine of the CA averaged over y:

cos θ = 1 d y d y 0 cos θ (y,h(y))dy.
The average angle is obtained as the inverse cosine of this quantity. For the case of periodical defects, θ and h are periodical functions of y and it is sufficient to average over one period λ y ,

cos θ = 1 λ y λ y 0 cos θ (y,h(y))dy. (11)
The CA hysteresis definition [1] as a difference of the cosines of the advancing and receding CAs gives one more reason to average the cosine of the local CAs and not the CAs themselves.

B. Vertical averaging

For a pattern whose wettability varies in the vertical direction, the "±" states that are crucial for the wetting hysteresis determination depend also on the relative position H (in the vertical direction) of the heterogeneity pattern with respect to the liquid level. Therefore, to obtain the wetting hysteresis in this case one needs to consider also the averaging along the vertical direction.

For spatially periodical patterns, the CL heights h + and h - are periodical functions not only of y but also of H ; h ± (y,H + λ z ) = h ± (y,H ) for any H , where λ z is the vertical period of the defect pattern. The meniscus states and the corresponding CLs thus need to be determined only for H ∈ [0,λ z ).

For the CLs with h ± = h ± (y,H ), one can calculate the corresponding CA values and their horizontally averaged cosines,

cos θ + = 1 λ y λ y 0 cos θ (y,h + (y,H ))dy, ( 12a
) cos θ -= 1 λ y λ y 0 cos θ (y,h -(y,H ))dy. ( 12b
)
It is evident that cos θ -and cos θ + are the upper and the lower bounds of cos θ , respectively. They can be used to obtain the advancing and receding CAs [4]. One yet needs to find out how. The most straightforward choices [4] that provide the upper and lower bounds for θ (y,z) are

cos θ r, * = min 0 H <λ z cos θ + , cos θ a, * = max 0 H <λ z cos θ -. ( 13
)
Apart from (13), in Schwartz and Garoff [4] another definition is discussed of wetting hysteresis, the straight contact line (SCL) approximation defined as cos θ a,SCL = max

0 z<λ z 1 λ y λ y 0 cos θ (y,z)dy, cos θ r,SCL = min 0 z<λ z 1 λ y λ y 0 cos θ (y,z)dy, (14)
where θ (y,z) is given by Eq. ( 4).

We suggest here a different definition (more specifically, a way of averaging) that has the advantage of being applicable in both the static and dynamic cases (in the limiting case of small capillary and Reynolds numbers), since it is based on our experience in the dynamic case. More precisely, since the static advancing and receding contact angles are defined as resulting from the preceding CL motion in the respective direction (caused by the plate motion), let us consider first the dynamic case by introducing the plate velocity u. The plate velocity is assumed to be positive if the plate motion is directed upward. Then one can write H = ut and express the CL height as h = h(y,t). Experimentally, the dynamic apparent CA is obtained by averaging over some time P ,

cos θ app = 1 λ y P P 0 dt λ y 0 cos θ (y,h(y,t))dy. (15)
When the defect pattern is periodic (in the vertical direction), so is the CL dynamics after an initial transient time.

A parametric dependence of h on u is assumed here. Such a definition does not contain the time explicitly and may thus be applied also to the static case where u → 0. This expression is close to the surface averaging of the CA cosine that was used by many authors; see e.g., [24]. By applying this expression to the states h ± discussed above, one obtains finally the expressions for the advancing and receding CAs:

cos θ r = 1 λ y λ z λ z 0 dH λ y 0 cos θ (y,h + (y,H ))dy, (17a) cos θ a = 1 λ y λ z λ z 0 dH λ y 0 cos θ (y,h -(y,H ))dy. (17b)
From Eqs. ( 13) and (17), it is quite clear that θ r, * > θ r and θ a < θ a, * . We consider the use of expressions (17) more appropriate than that of (13). Indeed, since the static configuration appears necessarily after the CL motion, the advancing and receding contact angles must satisfy the definitions θ a = θ eff (u → -0) and θ r = θ eff (u → +0) that result in (17). Note that such a "dynamic" definition corresponds to the advancing and receding CAs, calculated in numerical approaches that use dynamic (see, e.g., [25]) or pseudodynamic [9] algorithms. Equation (17) is not the only possible way for defining the ACA and RCA. For example, for horizontal averaging instead of Eq. ( 10) one can use Eqs. ( 6) and (7). If one uses different

h (2) (1) 0 H h cos (1) cos (a) (b)
(1)

h (1) 0 H h (2) h h cos (1) cos (2)
cos FIG. 2. (Color online) The highest possible metastable value h + of the CL height (black bold solid line) and the respective CA (red bold dashed line) are shown as functions of the plate position H for the horizontal stripe pattern for one period λ. This corresponds to the receding case. Two possible situations occur which differ by the relation between h and λ (2) : (a) h λ (2) and (b) h < λ (2) . A part of the period, where the CL belongs to the domain (1) , is shaded.

possible averaging also in the vertical direction, then one obtains nine different ways for obtaining the ACAs and RCAs and these all lead to different values of the CA hysteresis.

The importance of the solutions h ± (y,H ) (for a given plate displacement H ) is not limited to the expressions (17) only.

From their definition, it is apparent that the states h ± (y,H ) may be considered as corresponding to situations, where H varies with time (H = ut) and u approaches zero [u → +0 for h + (y,H ) and u → -0 for h -(y,H )]. They present a clear demonstration of the CL stick-slip motion (see, e.g., Figs. 2 and3

below).

One should note also that the generalization of the expressions (17) to nonperiodic patterns can be obtained by letting both λ y and λ z go to infinity (in practice, to the respective dimensions of the plate).

IV. HYSTERESIS CALCULATION METHODS

The free energy functional Eq. ( 3), as pointed out in the previous sections, depends on the plate vertical position H . One needs to obtain the equilibrium meniscus states f ± (x,y) that correspond to the CL heights h ± for each

H ∈ [0,λ z ) [this concerns question (ii) in the Introduction]. (2) h (1) 0 H h (1) h h (2) (1) cos (2) cos (a)
(2)

h 0 H h (1) h (2) h cos (2) cos (1) cos (b)
FIG. 3. (Color online) The lowest possible metastable value h - of the CL height (black bold solid line) and the respective CA (red bold dashed line) are shown as functions of the plate position H for the horizontal stripe pattern for one period λ. This corresponds to the advancing case. Two situations arise: (a) h λ (1) (black slanted dotted line (the dotted line partially coincides with the solid line), a nonperiodic part of the curve h -(H ) that appears because of the choice of the plate position corresponding to H = 0) and (b) h < λ (1) . A part of the period, where the CL belongs to the domain (1) , is shaded.

From these solutions, the distributions of the local CAs θ (y,h ± (y,H )) are found. They are used in Eqs. (17) to find the CA hysteresis.

At fixed H , we find the meniscus states f ± (x,y) numerically, with an iterative minimization procedure of the energy functional (3). To find the h + state that corresponds to the local minimum with the largest h value, one needs to start iterations from a state with the CL position h + 0 (y) with h + 0 h + ; otherwise, the minimization would bring the system into a wrong local minimum. The opposite inequality should be satisfied for the initial meniscus state h - 0 (y) used for the "-" state search. The states with the constant CL heights h (1) and h (2) , where

h (i) = l c 2 1 -sin θ (i) eq , ( 18
)
may be used for h + 0 (y) and h - 0 (y), respectively. Indeed, the above criteria on h ± 0 are satisfied because from the conditions (4) and ( 5) it is evident that for any CL position h(y), h (2) h(y) h (1) . The initial meniscus shapes corresponding to h (i) are obtained with a well known analytical solution of a capillary rise along a partially wetted vertical wall [START_REF] Landau | Fluid Mechanics[END_REF].

Every iteration step (let us say, number k) of the minimization procedure consists of two parts.

(i) First, a minimum of the free energy (3) is found for a fixed CL obtained on the iteration k -1 with CL height h k-1 = f k-1 (0,y). This results in a meniscus shape f k (x,y) satisfying the Laplace equation. The resulting local CA distribution θ k ( R) along the CL L ( R = [y,h k (y)]) may be different from its equilibrium counterpart θ ( R) given by the conditions (4) and ( 5).

(ii) In the second part of the iteration step we obtain the new position of the CL height h k by shifting the points of the CL vertically in the direction, decreasing the system free energy. The direction of the CL change, which decreases the system energy, is determined by considering the variation δ R of the point R belonging to the CL which leads to the variation of the system energy

δU = γ lg [cos θ ( R) -cos θ k ( R)](δ R ˙ z). ( 19
)
The iterative process stops when there is no movement of the CL that lowers the energy of the system.

In this work we obtain results using two different methods for the first part of the iteration step. The first method (method A) provides the numerical minimization of the energy functional. This is realized as follows. A part of the surface z = f (x,y) defined on a set of width λ y and length d x = 10l c is approximated by a set of triangles with N i × N j node points r ij = [x i ,y j ,f (x i ,y j )]; i = 1, . . . ,N i ; j = 1, . . . ,N j . The CL is represented by the points r 1j = [0,y j ,h(y j )],j = 1, . . . ,N j . The periodic boundary conditions f (x i ,y 1 ) = f (x i ,y N j ),i = 1, . . . ,N i are imposed in the y direction to simulate periodical solutions. The meniscus shape is determined by the analysis of the variation δU of the free energy (3) calculated for small virtual nodal displacements δ r ij . The number of points N j , which approximate the CL, controls the precision with which the CL and CA are determined. High precision of the determination of the CA value is achieved by comparing with solutions obtained with a larger N j value. The set N j = 60,80,100,120 and the value N i = 1000 are used. Convergence is assumed to be achieved when the cosine of ACA (or RCA) differs from the respective value obtained for the smaller N j belonging to the above set by less than 0.01.

The second method (method B) uses an asymptotic solution [4] of the Laplace equation, obtained for small meniscus slopes and arbitrary h k (y) function. In other words, an analytical approximation is used instead of the rigorous numerical minimization. The corresponding local CA value may be expressed as [22] cos

θ k ( R) = - 1 2λ y ∞ n=-∞ l -2 c + π 2 n 2 λ 2 y × λ y -λ y dy h k (y ) cos πn(y -y ) λ y . ( 20
)
Instead of using Eq. ( 19), one finds δh k from the pseudodynamic equation

δh k (y) = [cos θ ( R) -cos θ k ( R)]δt, ( 21
)
where δt is a small pseudo time step.

Note that such a theory reduces to that of Joanny and de Gennes [18] in the limit λ y /l c → 0, i.e., for the flat liquid wedge.

V. RESULTS AND DISCUSSION

A. Periodic stripe pattern

The horizontal stripes of different surface energy present the simplest heterogeneity pattern where wetting hysteresis occurs [9,16]. The alternating stripes of widths λ (1) ,λ (2) , characterized by the CA values θ (1) eq < θ (2) eq , respectively, are considered. The CL heights h ± (H ) (the y dependence is absent in this case) are periodic functions of period λ z = λ (1) + λ (2) ≡ λ. H = 0 is assumed to correspond to such a plate position, where the CL of height h = h (1) coincides with the stripe border, for which the domain (2) is above the domain (1) . To find h + (H ), consider the increase of H that corresponds to the upward displacement of the plate, i.e., the CL receding. It is evident that the meniscus will remain at equilibrium (i.e., θ = θ (1) eq and h = h (1) ) as long as H < λ (1) , shown by the solid line in Fig. 2. At H = λ (1) , the CL slips to a new position with a smaller height and θ becomes larger. If the stripe (2) is wide enough, i.e., λ (2) > h ≡ h (1) h (2) , the CL attains the equilibrium height h (2) during the slip [Fig. 2(b)]. If, on the contrary, h λ (2) , the CL stops its slip at the domain border [Fig. 2(a)] at the height h = h (1) λ (2) and becomes pinned to that border because the condition (5) of pinning holds. A similar pinning occurs in the case λ (2) > h at H = λh (when the CL attains the stripe border). The CL remains pinned while the inequality (5) holds, i.e., while h h (1) . Note that the CL does not situate at all inside the stripe (2) in the case of Fig. 2(a).

The H dependence of the CL height is equivalent to the time dependence introduced as H = ut in the limit u → 0; cf. Sec. III. For H belonging to regions where h(H ) = h (i) (a horizontal line), the meniscus is at equilibrium and CL moves freely inside the stripe (i) ; θ = θ (i) eq . The CL pinning at a stripe border corresponds to slanted (at 45 • ) h(H ) line portions. The slip corresponds to its vertical portions.

The determination of h -(H ) is similar. Again, we start from the same initial position with h = h (1) (Fig. 3). The direction of plate displacement is now downward (advancing CL), and H < 0. Since h (2) < h (1) , the CL remains pinned at the border until its height reaches h (2) , which occurs at H =h. From that moment on, θ = θ (2) eq and h = h (2) as long as H > -λ (2)h. At H = -λ (2) h, CL slips to a new position. If the stripe (1) is wide enough (λ (1) h), the equilibrium height h = h (1) and the CA θ = θ (1) eq are both attained [Fig. 3(b)]. In the opposite case [Fig. 3(a)], the CL becomes pinned at the domain border when its height attains h = λ (1) + h (2) < h (1) and the new period begins. In this case the CL never again attains its initial height. A part of the h -(H ) line between λ (1) + h (2) and h (1) does not belong to the periodical curve and for this reason is shown as dotted in Fig. 3(a). Note, that in this case the CL never lies inside the domain (1) .

The dependence of the cosine of the CA on H can be obtained from the dependence on H of the CL height h(H ) by the help of (8). The respective cos(θ )(H ) dependencies are sketched in Figs. 2 and3 with red dashed lines.

By comparing Figs. 2 and 3 one can check that the condition h -(H ) h + (H ) is verified (provided both functions are reduced to the same period).

Note, that the simulation from Fig. 2 of David and Neumann [9], where λ (1) = λ (2) < h, corresponds to the cases of Figs. 2(a) and 3(a).

Consider now the wetting hysteresis. With Eqs. (17), one can obtain by integration the following analytical expressions for RCA and ACA: cos θ r = 1 λ λ (1) cos θ (1) eq + max(0,λ (2) h)

× cos θ (2) eq - l c 6 2 3/2 1 + sin θ (1) eq 3/2 -4 - 1 l 2 c [max(0, h -λ (2) ) + h (2) ] 2 3/2 , ( 22a
)
cos θ a = 1 λ λ (2) cos θ (2) eq + max(0,λ (1) h)

× cos θ (1) eq + l c 6 2 3/2 1 + sin θ (2) eq 3/2 -4 - 1 l 2 c [min(0,λ (1) -h) + h (1) ] 2 3/2 . ( 22b
)
Since one has |h (i) | l c √ 2, the above expressions are always real valued. For comparison, the expressions ( 13) and ( 14) result in θ r,SCL = θ r, * = θ (1) eq and θ a,SCL = θ a, * = θ (2) eq independently of the pattern parameters. These values are independent of the density of the domain boundaries (where, in fact, the pinning occurs). The advancing and receding CAs calculated with Eqs. ( 22) as functions of the surface coverage p (2) = λ (2) /λ by (2) stripes are shown in Fig. 4 for the following choice of parameters: θ (2) eq = 80 • , θ (1) eq = 30 • , and λ/ l c = 0.3,0.1,0.01. The values of these parameters are the same as in Fig. 8 below. One can see that in the limit λ → 0, the ACA tends to θ (1) eq while RCA tends to θ (2) eq . The spatially averaged (Cassie) CA,

cos θ C = p (i) cos θ (i) eq + (1 -p (i) ) cos θ (3-i) eq , ( 23
)
is also shown on this figure (i = 2). It always situates between the receding and the advancing CAs. Normally, the hysteresis should vanish when the defect coverage tends to zero, i.e., for a homogeneous plate. This corresponds to two cases p (2) → 0 ( (1) -covered plate) and p (2) → 1 ( (2) -covered plate). The hysteresis does not, however, vanish in both these limits in Fig. 4. This happens because the stripe pattern is special. The above limits correspond to the infinitely thin defect stripes. The pinning in this case occurs because the whole CL remains pinned at the defect boundary. This becomes evident if one considers, e.g., Fig. 3 for the FIG. 4. (Color online) Cosines of the advancing and receding CAs for θ (2) eq = 80 • and θ (1) eq = 30 • as functions of the surface coverage p (2) by (2) stripes for λ/ l c = 0.3,0.1,0.01. advancing case. In the limit λ 2 → 0 (i.e., p → 0), the CL stick-slip behavior persists. It is shown below that such a pinning does not occur for y-dependent defect patterns.

B. Doubly periodic pattern of circular defects on a homogeneous base

As an example of heterogeneity pattern, where the CA θ depends not only on the vertical, but also on the horizontal coordinate, a doubly periodic defect pattern (see Fig. 1), studied previously in [22], is considered here. It consists of circles of radius a ordered on a square lattice with a distance λ y = λ z ≡ λ between the centers of two neighboring circles. One of the sides of the square lattice is parallel to the liquid level. For such a pattern, we have chosen the plate position H = 0 to be the one, where the centers of one of the rows of circular defects coincides with the liquid level far from the plate. From now on all the lengths are considered to be dimensionless, in the units l c . In the absence of analytical solution for such a problem, a numerical analysis is performed with the two methods of Sec. IV.

Contact line pinning and hysteresis

The metastable CL states are studied in this section in the small slope approximation (method B), and for its validation a comparison is made with the method A results. In order to do that, the system parameters are chosen to be the same as in [22]: The circular defects (domains (2) , θ (2) eq = 110 • ) are less wettable than the rest of the plate (domain (1) with θ (1) eq = 70 • ). The period of heterogeneity pattern is λ = 0.3 and the radius of circular defects is a = 0.1.

First, at a fixed plate position H , the meniscus states with the highest and the lowest possible averaged CL heights h ± (y,H ) are analyzed. During the minimization process, the CL passes through several rows of defects before the equilibrium state is attained.

An example for h + (y,H ), obtained by method A, is shown in Figs. 5 [START_REF] Landau | Fluid Mechanics[END_REF].

(a) and 5(c). (for 20 values of H ,

H i = iλ/20, i = 0,1, . . . ,
As mentioned in Sec. IV, the straight CL of height h = h (1) needs to be chosen for the initial state for the iterative procedure for every H to find h + (y,H ). Since the CL height at the circular spots should be lower than outside them, the CL starts to slide down during the iterations until the equilibrium state is reached. The CL reaches the lower border of the defect and sticks to it. This situation is common for several heights H i , i = 0, . . . ,9 (where H i = iλ/20), since a range of CAs is allowed at the border of defect according to the condition (5); see Fig. 6. After slip of the CL to the next range of defects, the same behavior is observed for i = 17,18,19. When the position of heterogeneity pattern with respect to the liquid level is such that the initial CL belongs entirely to the domain (1) , the initial CL is in equilibrium and θ (y,H i ) = θ (1) eq . This is the case for i = 10, . . . , 16.

The situation is different for the states with the lowest average height h -(y,H i ) shown in Figs. 5(b) and 5(d) (here H i = -iλ/20). During the minimization process, the parts of the initial CL (chosen to be of height h = h (2) ) lying outside the circular defects move toward the height h (1) and reach θ = θ (1) eq . At the same time, the parts of the CL inside defects maintain the CA θ = θ (2) eq . Unlike the h + case, there are no parts of CL pinned at the border of the defects; a part of the CL lies inside defects.

The spatial variation (along the y axis) of θ reflects the above described CL behavior. In the h -(advancing) case, a part of the CL always lies inside defects so that θ sharply varies at their borders from θ (1) eq to θ (2) eq . In the h + (receding) case, a part of CL is pinned on the defect borders and θ smoothly according to (5). Its variation is shown in Fig. 6 for one period λ of the the CL shapes corresponding to Figs. 5(b) and5(d).

With the one can obtain the averaged cosines of the CAs. These quantities are shown in Fig. 7 as functions of the plate vertical position H with respect to the liquid level. The results of the small-slope approximation are shown also in the figure with lines: the solid lines for the receding and the dashed lines for the advancing. The symbols show the numerical results by method A: circles for the advancing CAs and squares for the receding CAs.

The numerical results, obtained by method A, are compared to those of the method B. One can see that the agreement is very good. Note that, similarly to Figs. 2 and3, Fig. 7 shows the variation of cos θ for advancing and receding motion (where H should be associated with ut) in the limit u → 0. One can compare the behavior for the stripe and the circle patterns. The receding curve in Fig. 7 is directly comparable with Fig. 2(a), however, with a shift caused by the different choices of the H = 0 position. Since the circle diameter is similar to λ (2) of the stripe pattern and the CL initial height in the circle pattern is about h (1) ≈ λ + 0.5a, the curves in Fig. 2(a) should be shifted to the left by about λ (1) + 0.25λ (2) to be compared. Such a shift makes both cases look completely similar. When H increases, first pinning (corresponding to the slanted parts) occurs both in the receding curve of Fig. 7 and in Fig. 2(a). The pinning corresponds to the curves 0-9 in Figs. 5(a) and 5(c). Then the CL becomes straight in Figs. 5(a) and (5 c), which corresponds to the horizontal parts of both the receding curve in Fig. 7 and Fig. 2(a). Slip (corresponding to the vertical parts) occurs in Fig. 7 and in shifted Fig. 2(a) almost at the same position.

For the advancing case, the difference between the stripe and circle patterns is larger. Similarly to the receding case, the curves of Fig. 3(a) should be shifted by about 0.75λ (2) to the right. The advancing curve in Fig. 7 is quite different from that in Fig. 3(a). In the circle case, cos θ -> cos θ (2) during the whole period. In other words, the CL never becomes straight [cf. Figs. 5(b) and 5(d)] and no horizontal portion in the cos θ -(H ) curve exists. This occurs because there is no possibility for a straight CL to situate entirely inside the (2) domain.

From Figs. 5(c) and 5(d), one can deduce that the pinning occurrence at the domain border is independent of the direction of the CL motion and depends only on the side of the border at which the liquid situates. Indeed, the pinning for both advancing and receding motion occurs at a border where the liquid situates at the more wettable side.

One can check in Fig. 7 that the inequalities cos θ (2) cos θ - cos θ + cos θ (1) are satisfied, as they should be according to the definition of the "±" states.

One needs to note that for all meniscus shapes, the average cosine of the CA [Eq. (11)] coincides with high precision with the cosine of the average CA, obtained by the averaged height of the CL [Eq. (7)], the difference between the two values being less than 0.01.

The wetting hysteresis interval [θ r ,θ a ] = [73.58 • ,92.7 • ] given by Eqs. (17) can now be compared to that given by Eq. ( 13) [θ r, * ,θ a, * ] = [70 • ,96.5 • ]. The SCL approximation (14) results in [θ r,SCL ,θ a,SCL ] = [70 • ,90 • ]. One can see that the hysteresis values given by Schwartz and Garoff [4] are considerably off those obtained with Eqs. (17) compatible with the dynamics history of the CL motion. All the above numbers have been obtained with method A of Sec. IV. On the other hand, the accuracy of the small-slope approximation (method B) results in [θ r ,θ a ] = [73.6 • ,91.2 • ]. These data provide a good accuracy: The hysteresis is within 2 • from the method A results.

Wetting hysteresis parametric analysis

In this section we discuss the dependence of the wetting hysteresis magnitude on the surface coverage of the plate by the defects for the circle pattern and compare it to the stripe pattern; cf. Fig. 4. We take advantage of the possibility to calculate with the method A the menisci of arbitrary (high) slopes and use the values θ (2) eq = 80 • and θ (1) eq = 30 • so that the results are comparable with those of the previous works [9,10]. Similarly to the previous section, the case of less wettable defects on more wettable base is considered first. The surface coverage is defined in this case as p (2) = πa 2 /λ 2 . The calculations have been performed for p (2) < 0.75: For a > λ/2 (p > π/4) the circles overlap and the pattern topology changes.

The cosines of the advancing and receding CAs calculated with Eq. ( 17) for both patterns are shown in Fig. 8 for different but fixed λ values and varying a. The tendency is similar FIG. 8. (Color online) Cosines of the advancing and receding CAs for the doubly periodic pattern of circular defects as functions of the defect surface coverage for fixed values of the spatial period λ = 0.3,0.1,0.01. The equilibrium CA is θ (2) eq = 80 • on the circles and θ (1) eq = 30 • on the base. The prediction of Raj et al. [17] is also shown with open circles. at fixed a and varying λ. For a more complete picture, the Cassie"s angle (23) (i = 2) is also shown on this figure with a dotted line.

The RCA value remains close to θ (1) eq and the cosine of the ACA decreases with p (2) . As a result, the hysteresis amplitude increases with the surface coverage, which agrees qualitatively with the experimental and other numerical [10,[25][26][27] data. The difference of the cosines of the receding and advancing CAs obtained by the help of Eqs. ( 9) and ( 10) is not a linear or quadratic function of the defect concentration. This results is in contrast to the case of randomly heterogeneous surfaces, where we obtain quadratic dependence on the concentration of defects [10].

Note that the defect coverage p (2) cannot be considered as a unique parameter, characterizing hysteresis; the hysteresis depends also on the heterogeneity period λ. It is, however, likely that the latter dependence will disappear in the limit λ → 0, i.e., for periods of the heterogeneity patterns much smaller than the macroscopic liquid surface curvature. This tendency is suggested by the expression (22) for the cosine of the local CA that saturates in the limit λ/ l c → 0.

At this stage, a comparison can be performed with the semiempirical model of Raj et al. [17]. This model considers the CL behavior for the same circle pattern. Raj et al. proposed that the receding CL stays always between the less wettable defects, which corresponds to θ r,Raj = θ (1) eq . As for the advancing CL, it was argued to belong partly (on average, the CL length 2a per period λ) to the defects. This results in [note that there are typos in Eq. ( 8) of Raj et al. [17] that is supposed to correspond to their Fig. 5(b); numbers 1 and 2 should be interchanged] cos θ a,Raj = 2a λ cos θ (2) eq + 1 -2a λ cos θ (1) eq .

(

) 24
This result is shown with open circles in Fig. 8. One can see that Eq. ( 24) agrees with the rigorous numerical data in the limit λ → 0 [the model ( 24) is independent of this parameter]. The same tendency was obtained [10] for the case of square shaped defects. The above data are obtained for less wettable defects on a more wettable base. One can consider the opposite case of the domain (1) formed by the circular defects with θ (1) eq = 30 • on the base (2) with θ (2) eq = 80 • . The defect surface coverage in this case is p (1) = πa 2 /λ 2 . The ACA and RCA results are inverted with respect to the former case. The ACA remains close to θ (2) eq , while the cosine of RCA grows with p (1) . The hysteresis vanishes in the limit p (1) → 0. The Cassie's angle, given by Eq. ( 23) with i = 1, situates again between the ACA and RCA. An example of the calculated RCA variation with p (1) for λ = 0.03 is shown in Fig. 9. Raj et al. [17] provide the following expression for this case: cos θ r,Raj = 2a λ cos θ (1) eq + 1 -2a λ cos θ (2) eq . (25) This dependence is also shown in Fig. 9. The agreement with the numerical data, obtained for a small λ, is good.

FIG. 9. Cosine of as function of the surface coverage, calculated for λ = 0.03 for the doubly pattern of circular with more wettable circles with θ (1) eq = 30 • and less wettable base with θ (2) eq = 80 • . The prediction of Raj et al. [17] is also shown with open circles.

VI. CONCLUSION

In this work we suggest a general approach for determination of the advancing and receding contact angles and discuss how to find them from known properties of the substrate heterogeneity. Our expressions for the static angles account for the preceding CL motion history and result in the values measurable with the Wilhelmy balance. We have shown that when determining the ACA and RCA it matters how one averages the CA, first at fixed position of the plate pattern with respect to the liquid level and second over the many different possible relative positions of the defect pattern on the plate with respect to the liquid meniscus.

Two examples of periodic heterogeneity patterns are used to illustrate the wetting hysteresis determination. One of them is the simplest pattern in which the hysteresis occurs, the alternating horizontal stripes of different wettability. The analytical expressions for the advancing and receding contact angles are obtained for it. The analysis of the other, doubly periodic pattern of circular defects, has been carried out by methods A and B. Their comparison shows the validity of the small-slope approximation that needs much less computation time. A comparison with other existing expressions for the wetting hysteresis has also been carried out. It is confirmed that the determination [10] of the CA hysteresis using the averaged CL height or through the averaged cosine of the CA lead to the same results.

The pinning on both the stripe and the doubly periodic patterns has been discussed in some detail. Their similarities and differences have been revealed. It has been shown that the pinning occurrence at a border between the domains of different wettability is independent of the direction of the CL motion. The pinning occurs at the border where the liquid situates at the more wettable side.

Next, the hysteresis was analyzed as a function of the defect density. This function also differs from the dependence, found for the randomly distributed defects, and is close to the dependence, found for the square checkerboard defect pattern. It was shown that for the case when the pattern period is much smaller than the average liquid meniscus radius of curvature, the wetting hysteresis is well described by the semiempirical expressions of Raj et al. [17].

Cassie's angle is one of the parameters, characterizing the heterogeneous surfaces, however, it is not sufficient for characterizing the heterogeneous surfaces. Our results displayed in Figs. 4, 8, and 9 show that the ACAs and the RCAs [calculated by the help of Eq. ( 17)] remain different from Cassie's angle (its value is always between the ACA and the RCA). This holds true also in the limit of minuscule defects. We note here again that these results are valid for heterogeneous surfaces with abrupt changes of the surface tension.

Although the definitions for RCA and ACA were studied on the example of two heterogeneity patterns, there is no limitation for studying the rough surfaces in Cassie's regime. The application to this case is straightforward. However, the application of this approach to studying CA hysteresis on rough surfaces in Wentzel's regime, as well as the application of the same approach for studying the wetting hysteresis in dynamics, needs further development and will be a subject of future research.

FIG. 1 .

 1 FIG. 1. (Color online) Schematic drawing of the Wilhelmy plate geometry. The vertical plate has surface defects of different wettability, shown with darker color.

  The temporal period is P = λ z /|u|. By performing a change of variables one obtains cos θ eff = 1 y,h(y,H ))dy.

FIG. 5 .

 5 FIG. 5. (Color online) Contact line shapes h ± (y,H i ) calculated for 20 different positions i = 0,1, . . . ,19 of the heterogeneity pattern; the index i is the curve parameter. Panels (a) and (c) correspond to h + (y,H i ) with the highest h (receding case, plate moves upwards, H i = iλ/20) and (b) and (d) correspond to h -(y,H i ) with the lowest h (advancing case, plate moves downwards, H i = -iλ/20). The CLs in panels (a) and (b) are in the reference system of the liquid pool. The CLs in panels (c) and (d) are in the reference system of the plate; however, in panel (c) every CL height h i (and the respective defects) is shifted by -H i and in panel (d) by H i . The defect borders (dashed lines) are shown too.

FIG. 6 .

 6 FIG. 6. (Color online) Variation of the local CA along the CL for the receding CL shapes shown in Figs. 5(a) and 5(c).

FIG. 7 .

 7 FIG. 7. (Color online) Horizontally averaged cosines of CAs corresponding to the receding and advancing cases. The lines come from the small-slope approximation (method B); the symbols show the method A results. Each symbol corresponds to a in Figs. 5. The numbering is left to right for the receding and right to left for the advancing curve.


D L'hôte

R Tourbot

F Ladieu

P Gadige

Control parameter for the glass transition of glycerol evidenced by the static-field-induced nonlinear response

Keywords: number(s): 64.70.kj, 77.22.-d, 89.75.Da

By studying a nonlinear susceptibility on supercooled glycerol, we show that applying a static field E st increases the glass transition temperature T g by an amount quadratic in E st . This has important consequences: (i) it reinforces the relation between the two paths put forward in the last years to unveil the dynamical correlation volume close to T g ; (ii) it clarifies the interpretation of nonlinear measurements; (iii) it yields a new control parameter of the glass transition, which paves the way for experiments deepening our understanding of glasses.

I. INTRODUCTION

Many materials of industrial importance are glasses, from window glasses and plastic bottles to emulsions, foams, dense assemblies of grains, and so on. Understanding the formation of these amorphous solids, the so-called glass and jamming transitions, is still a challenge, despite a large research effort over the last decades [START_REF]Dynamical Heterogeneities in Glasses, Colloids and Granular Media[END_REF]. All these systems exhibit a remarkable "viscous slowing down" phenomenon, i.e., a dramatic increase of the relaxation time τ α when varying the control parameter by a modest amount. Although there is some evidence that jamming and glass transitions correspond to different physical phenomena [START_REF] Ikeda | [END_REF], it was shown [3][4][5] in both cases that, close enough to the transition, relaxation occurs by groups of N corr dynamically correlated particles, with some groups being faster than others [hence the name dynamical heterogeneities (DH) for these groups]. However, the underlying mechanism of such a heterogeneous dynamics remains controversial [3][4][5]. For supercooled liquids, according to some approaches [6,7] the drainage of free volume when density increases is the main reason for the strong increase of τ α when lowering the temperature T towards the glass transition temperature T g . At the other extreme, another scenario relies on thermally activated processes on a constant-density energy landscape [8][9][10]. This has motivated significant experimental efforts [11][12][13][14][15] where τ α has been monitored both as a function of temperature T and pressure , to assess which is the dominant control parameter of the glass transition.

Here we report the existence of a new control parameter of the glass transition. Namely, we show that applying a static electrical field E st onto a prototypical glass forming liquid-glycerol-is equivalent to increasing T g by an amount δT g (E st ) ∝ E 2 st . Even though such an effect is anticipated in the framework of mode coupling theory (MCT) [16], this is the first experimental evidence of the influence of a static electric field on T g . This is not surprising owing to the smallness of the effect that we report here (δT g 0.3 mK for E st = 1 MV/m). Even if this is much smaller than the shifts of T g reported when varying , we shall show that the T g (E st ) results reported here have important physical consequences-in the same way as the T g ( ) dependence * deceased † francois.ladieu@cea.fr gives a deeper insight about the glass transition [11][12][13][14][15]. First, our results establish a quantitative link between the two independent methods put forward theoretically to relate N corr to macroscopic observables [17][18][19]. This link is of great importance since it deepens the overall consistency of the recent progress made about the characterization of the DH's. Second, our results put a severe constraint on the phenomenological models of the nonlinear response, which clarifies the interpretation of nonlinear measurements. Third, the fact that E st is a control parameter opens the way for new experiments, deepening our understanding of glasses.

The paper is organized as follows. We first define and explain in Sec. II how we measure the nonlinear susceptibility χ (1) 2;1 associated to a static field. Then, in Sec. III, we show the behavior of χ (1) 2;1 and why it unveils the above-mentioned shift of T g with E st . Section IV is devoted to the discussion of the main physical consequences of our experiments. Finally, in Sec. V, we give some perspectives and summarize our work.

II. EXPERIMENTS

A. Definitions

In a closed cell, ultrapure glycerol was put between stainless steel electrodes and submitted to an electric field E. The resulting polarization P was overwhelmingly dominated by its linear component P Lin . We recall that, for an ac field E ac cos(ωt) oscillating as a function of the time t with the frequency f = ω/(2π ), the complex linear susceptibility χ Lin = χ Liniχ Lin (with i 2 = -1) is defined by P Lin = 0 E ac [χ Lin cos(ωt) + χ Lin sin(ωt)] where 0 is the vacuum dielectric constant. At a given T , χ Lin is peaked at a frequency f α , which yields the relaxation time τ α ≡ 1/(2πf α ). As detailed in Refs. [20][21][22] (see also Fig. 1), in our "two samples bridge" P Lin cancels very efficiently. This allows to study the nonlinear polarization δP ≡ P -P Lin at any f , at variance with other setups [23][24][START_REF] Herweg | For the first works about nonlinear dielectric effects[END_REF]. For an isotropic liquid submitted to an external electric field E(t), the most general expression of the nonlinear polarization δP (t) is given by a series expansion in E(t) [note that no even terms are allowed because of the symmetry with respect to field reversal E(t) → -E(t)]:

δP (t) 0 = ∞ -∞ χ 3 (t -t 1 ,t -t 2 ,t -t 3 ) × E(t 1 )E(t 2 )E(t 3 )dt 1 dt 2 dt 3 + • • • . ( 1
)
FIG. 1. (Color online) Inset: Representation of our "two samples bridge": the source voltage V s (t) is the sum of a constant voltage V st and of an ac voltage V ac cos(ωt). The left (right) branch contains the thick (thin) sample of impedance Z thick (Z thin ). With a lock-in amplifier, we measure the 1ω component δV 1ω of the differential voltage V m between the two middle points of each arm of the bridge, after having tuned the values of the impedances z thin and z thick so as to suppress accurately the linear voltage coming from the linear response P Lin . Main panel: Example of the rms values of δV 1ω (V ac ,V st ) measured in the stationary regime (see text) by varying V st and keeping constant V ac . To get accurate values of the quantity of interest V 1ω (V ac ,V st ) = δV 1ω (V ac ,V st ) -δV 1ω (V ac ,0), we measure the stationary regime value of δV 1ω (V ac ,0) before switching on and after switching off V st . We average these two values of δV 1ω (V ac ,0) and subtract this average to δV 1ω (V ac ,V st ) to get the sought V 1ω (V ac ,V st ) which is directly proportional to χ (1) 2;1 . One observes on the graph that V 1ω (V ac ,V st ) has a modulus proportional to V 2 st (the thick solid line is such a quadratic dependence) and a phase independent on V st , in good agreement with Eq. (2). Note that these two key features are not observed for the values of δV 1ω (V ac ,V st )-a single value of δV 1ω (V ac ,0) is plotted on the graph at an arbitrary value of V st = 1.1V-. The data plotted here were obtained in the following conditions: the thickness of the thin (thick) sample is 16.8 μm (46.2 μm); T = 197 K; f = 0.8 Hz 3f α , where f α is the relaxation frequency; z thin 3.56 M ; z thick 9.9 M ; V ac = 7.84 V rms , S = 2.27 cm 2 where S is the surface of each sample. To make the connection between measured voltage and polarization, one has to use the standard complex impedance relating voltage to the intensity I , as well as the relation I = S∂P /∂t, see, e.g., Ref. [20].

For the fields used here, δP is overwhelmingly dominated by the cubic response (∝E 3 ), this is why the higher-order terms proportionnal to E 5 ,E 7 , . . . , are not written explicitly in Eq. [START_REF]Dynamical Heterogeneities in Glasses, Colloids and Granular Media[END_REF]. All the nonlinear cubic responses come from Eq. ( 1) and can be characterized by the triple Fourier transform kernel χ 3 (ω 1 ,ω 2 ,ω 3 ). For a pure ac field E ac , δP contains two terms studied in Refs. [20][21][22]. The first one oscillates at 1ω and is given by δP 1ω (E ac ) = 3 4 0 E 3 ac |χ (1) 3 | cos(ωtδ (1) 3 ) where we have noticed the corresponding cubic susceptibility [22] by χ (1) 3 instead of the more heavy notation χ 3 (-ω,ω,ω). The second one ∝χ 3 (ω,ω,ω)E 3 ac oscillates at 3ω, and similarly [20,21], we call the corresponding cubic susceptibility χ (3) 3 . In the present work, a static field E st was applied on top of an ac field E ac . This yields a series of new cubic terms in δP , e.g., terms proportional to E 2 ac E st or to E ac E 2 st . Here we focus on the new cubic term P 1ω , which appears due to the static field in the 1ω component of δP . Since P 1ω ≡ δP 1ω (E ac ,E st ) -δP 1ω (E ac ,0), where we recall that δP 1ω (E ac ,0) is the term proportional to χ (1) 3 mentionned above, one gets

P 1ω = 3 0 E ac E 2 st χ (1) 2;1 cos ωt -δ (1) 2;1 , ( 2
)
where χ (1) 2;1 stands for χ 3 (0,0,ω), and -δ (1) 2;1 is its phase.

B. Measurement of χ (1)

2;1

In our experiment E ac is kept constant, and thanks to our bridge setup described in Fig. 1, we monitor directly the first harmonics component δP 1ω of δP for various values of E st . We then form the difference P 1ω = δP 1ω (E ac ,E st ) -δP 1ω (E ac ,0). As a result P 1ω is directly proportional to the sought χ (1) 2;1 . More precisely, for each given values of T ,f , we first measure the value of δP 1ω (E ac ,0), then the static field is switched on to a given E st and we wait for a time τ st , long enough to let all the transients vanish, before measuring δP 1ω (E ac ,E st ). Then we switch off E st and measure again the stationary regime value of δP 1ω (E ac ,0), average it with the value measured before switching on E st , subtract this average to δP 1ω (E ac ,E st ), and obtain the sought P 1ω . We then do again this sequence with the opposite sign of E st : P 1ω is measured again, with the same delay of τ st . This inversion allows both to check that our measurements do not depend on the sign of E st , as expected from Eq. ( 2), and to avoid any spurious electrode charging effects coming from any slow drift of the small amount of free ions in the samples. For a significant number of points, by increasing the value of τ st by a factor of 4, we checked that τ st was long enough so as to have no influence on our measurements. This is very important since it shows that our results are really in the stationary regime, i.e., they are not affected by the transient effects coming from the switching on of the static field. In practice, f τ st and f α τ st were typically much larger than 10, the minimal values. Finally, we checked that P 1ω was proportional to E ac . Our "two samples bridge" setup has already been depicted in Refs. [20][21][22], but for the sake of clarity, and owing to the fact that χ (1) 2;1 has never been reported before, all the key information is given in the caption of Fig. 1, hence its unusual length. Typically, in the thinnest sample E ac = 0.65 MV/m, and 0 E st 4.1 MV/m. We emphasize that, in Fig. 1, the analysis of P 1ω as a function of E st yields the quadratic dependence expected from Eq. ( 2), giving χ (1) 2;1 with a ±4% error bar.

III. RESULTS

In the previous works [26][27][28] devoted to the effects of pressure , the shift δT g ( ) of T g is obtained by superimposing the polarization obtained at two different pressures. We follow this rule here: Stating that E st shifts T g by δT g (E st ) is equivalent to fulfill the equation P 1ω (E ac ,E st ,T ) = P 1ω [E ac ,0,T -δT g (E st )]. Note that, as pressure experiments are done with a pure ac field E ac , δT g ( ) is the solution of P (E ac , ,T ) = P [E ac ,0,T -δT g ( )]. Here the corresponding equation must be written on the odd harmonics only since E st generates even harmonics (0ω, 2ω) which do not exist when E st = 0. This is why we solve the equation on P 1ω , FIG. 2. (Color online) Comparison of χ (1) 2;1 (symbols) and of -κ(∂χ Lin /∂T ) (lines, shifted horizontally by a factor λ eff 0.80). Left axis: Moduli (full symbols or plain lines). From top to bottom T = 197, 202, 211, and 218 K, which amounts, at E st = 0, to f α = 0.27, 2.1, 60, and 520 Hz. Right axis: Phases (open symbols or dashed lines). Here λ eff = 0.80 ± 0.04, κ = (1.18 ± 0.03) × 10 -16 Km 2 /V 2 . Symbols are larger than the error bar: ±4% for |χ (1) 2;1 | and ±4 degrees for its phase.

the similar one on the 3ω component of P requires to measure five order responses, i.e., a more sensitive setup. For symmetry reasons we set δT g = 3κE 2 st with κ a constant. As we shall find δT g (E st )

T , and as P is overwhelmingly dominated by P Lin , Eq. ( 3) below is obtained by the following: (i) lowest-order expansion of P Lin (E ac ,0,T -3κE 2 st ); (ii) identification of the expansion with Eq. ( 2). This is how we obtain χ (1) 2;1 (f,T ) = -κ

∂χ Lin ∂T f/λ eff ,T ,E st =0 , ( 3
)
where the introduction of the effective parameter λ eff -equal to 1 in the simplest case-is justified below. Figure 2 shows that, for our data, Eq. ( 3) holds, as soon as f/f α 0.2, with λ eff 0.80 and κ 1.18 × 10 -16 Km 2 /V 2 . Figure 2 corresponds to a T range where τ α varies by more than three decades: in this range, κ and λ eff may slightly increase FIG. 3. (Color online) Inset: Qualitative picture of the trivial, singular, and global χ (1) 2;1 . Main panel: Evolution of T g with E st by the quantity 3κE 2 st . The size of the points is larger than the ±2.5% uncertainty on κ.

FIG. 4. (Color online)

Inset: Measured dimensionless susceptibility X (1) 2;1 (same T s as in Fig. 2). Main panel: T dependencies of various quantities Z normalized to 1 at 202 K. when lowering T but this possible effect is extremely tiny since it is close to the experimental uncertainty (namely ±2.5% on κ and ±5% on λ eff ). We emphasize that our data obey Eq. ( 3) both for the moduli and for the phases, which is highly nontrivial. By contrast, keeping λ eff reasonably close to 1, one finds [21,22] that Eq. ( 3) holds neither for χ (1) 3 nor for χ (3) 3 . Let us now briefly explain why for our data λ eff is not strictly equal to 1. We first recall, see Fig. 2, that Eq. ( 3) fails for our data at low values of f/f α , i.e., for time scales larger than τ α , over which the liquid flows: this destroys the glassy correlations, and makes any molecule effectively independent of others [START_REF]Dynamical Heterogeneities in Glasses, Colloids and Granular Media[END_REF]4]. When f f α , we thus expect a "trivial" physics, resembling that of an ideal gas: the corresponding |χ (1) 2;1 | triv is depicted as a dashed line [29] in the inset of Fig. 3. It is qualitatively different from the "singular" |χ (1) 2;1 | sing [thin solid line in the inset of Fig. 3] that one obtains when the glassy correlations dominate. Increasing f from f f α , we thus expect a crossover between |χ (1) 2;1 | triv and |χ (1) 2;1 | sing , yielding the "global" |χ (1) 2;1 |: the crossover distorts |χ (1) 2;1 | even close to f α , which may slightly shift its peak-see the thick solid line in the inset of Fig. 3. As this crossover has little impact on ∂χ lin /∂T , which is in any case peaked around f α , we thus expect that a factor λ eff close to 1 is needed in Eq. ( 3) to compensate for the small shift of χ (1) 2;1 . This reasoning is reinforced by several facts: first, in MCT [16] one recovers Eq. ( 3) with λ eff = 1 consistently with the divergence of the singular response which rejects the crossover to f/f α → 0, see also Ref. [START_REF]) 2;1 gradually decreases, which could explain the possible slight increase of κ and λ eff at low T -but this T evolution of κ[END_REF]. Second in the inset of Fig. 4, one sees that, in the limit f/f α 1, X (1) 2;1 (defined below) becomes T independent, which is characteristic of an ideal gas [START_REF] Dejardin | [END_REF]. Third, similar considerations have already been elaborated in Refs. [21,22] for the low frequency behavior of χ (1 or 3) 3 .

IV. DISCUSSION

A. Physical picture

Having qualitatively accounted for the limits of Eq. ( 3) close to T g , we interpret our data as reflecting the fact that E st is a new control parameter of the glass transition since it shifts T g by an amount 3κE 2 st , as illustrated in Fig. 3. To give a physical picture, we recall that translation and rotation are well coupled in glycerol [32], and that dynamical correlations do not correspond to any simple spatial order [START_REF]Dynamical Heterogeneities in Glasses, Colloids and Granular Media[END_REF]3]. As E st tends to rotate the molecules and align their polarization vectors, it also alters their positions because molecules have a given shape and rotating them conflicts with steric constraints: E st thus corresponds to a new constraint unfavorable to dynamical correlations, i.e., it slows down the dynamics, hence the increase of T g with E st . More generally, a static field gives rise to subtle effects strongly dependent of the kind of frustrated system which is considered: in the very specific case of spin glasses, the critical temperature may decrease with the static (magnetic) field (although this is strongly debated [33][34][35][36]), while in other systems reentrant behaviors at very high fields may happen [37]. Some further work is in progress to study the behavior of δT g (E st ) at higher values of the static field.

B. Relating nonlinear cubic responses and T χ T 1. Frequency dependence

The fact that our data of Fig. 2 follow Eq. ( 3) is the first experimental evidence that (some) nonlinear responses can be tightly related to T χ T ≡ T ∂(χ Lin / χ 1 )/∂T where χ 1 = χ Lin (0)χ Lin (∞). This is qualitatively important since two lines of thought have related the same number N corr of dynamically correlated molecules either to the cubic susceptibility χ (3) 3 (see Ref. [17]) or to Max f |T χ T | (see Refs. [18,19]). Up to now, the only experimental indication of the compatibility of these two approaches is that the T dependence of the dimensionless version X (3) 3 of χ (3) 3 is close to that of Max f |T χ T |, see Fig. 4 and Refs. [21,22,38,39]. The fact that the frequency dependence of a cubic susceptibility [χ (1) 2;1 ] can be mapped to that of T χ T on a wide frequency interval is a new result which strongly reinforces the compatibility of the two theoretical approaches.

Temperature dependence

To get the T dependence of N corr precisely [17], the cubic susceptibilities χ (k) n , with n = 3 or "2; 1" and k = 1 or 3, must be put in their dimensionless form

X (k) n = χ (k) n {k B T /[ 0 ( χ 1 ) 2 a 3 ]
} where a 3 is the molecular volume and k B the Boltzmann constant. As evoked above, in the trivial case of an ideal gas the various |X (k) n (f/f α )| are T independent [29], while in a supercooled liquid, for n = 3, they are proportional to N corr (T ) according to Refs. [17,21,22,38]. In Fig. 4 all the observables Z(T ) are normalized at 202 K: (i) the thick line corresponds to the Max f |T χ T | which was found in Refs. [21,38] to be very close to the T dependence of the maximum over the frequency of |X (3) 3 |; (ii) the thin dashed line corresponds to the best estimate of N corr (T < 204.7 K), drawn from the refined analysis of Ref. [22] comparing |X (3) 3 | and |X (1) 3 |; (iii) the points show the T dependence of |X (1) 2;1 | drawn [40] from the inset of Fig. 4. Figure 4 shows that the various |X (k) n | have mainly the same T dependence. This result goes beyond the prediction of Ref. [17], where the effect of E st was not studied, and had been anticipated only within MCT [16]. Lines are fits of X (1) 2;1 with the toy model of Ref. [START_REF] Ladieu | [END_REF] by choosing either N corr = 16 and δ = 0.60; or N corr = 10 and δ distributed onto

[0; 1] as ∝ exp [-(δ -δ 1 ) 2 /(2δ 2
2 )] with δ 1 = 0.30 and δ 1 = 0.29, see the Appendix.

C. Constraints put onto phenomenological models

Box model

Even though our results deepens the consistency of using cubic susceptibilities to get the T dependence of N corr , a problem is still pending. Indeed, a phenomenological model, called the "box model" [42][43][44] has been able to account fairly well (see Ref. [START_REF]The agreement is quantitative for X[END_REF]) for the nonlinear experiments published up to now despite it assumes that the nonlinear behavior comes from the power P dissipated by the electric field in each DH: this yields, finally, χ (k) 3 to be independent of N corr . Our results show that the dissipated power P is not the main source of χ (1) 2;1 . Indeed, as our experiments are in the stationary regime, varying E st does not change the power dissipated by the dipoles, and P only varies due to the small dc conductivity σ coming from the small amount of mobile ions in the liquid. It turns out, see the Appendix, that σ is smaller by a factor 300 than 0 ω α χ Lin (f α ), which is the corresponding quantity when measuring χ (1) 3 (f α ). As the box model fairly well accounts for χ (1) 3 (f > f α ), we predict within this model, for f f α , values of |X (1) 2;1 | which are orders of magnitude smaller than |X (1) 3 |-see the Appendix for a detailed calculation. This is contradicted by Fig. 5, where |X (1) 2;1 | has the same order of magnitude as |X (1) 3 /4|-the factor of 4 is due to the precise definitions of χ (1) 3 and of χ (1) 2;1 given in Sec. II A. Thus, for the first time, the box model is unable to account for a nonlinear experiment, namely χ (1) 2;1 . This is an important step forward clarifying the interpretation of nonlinear responses.

"Toy" model

For completeness we now briefly evoke another phenomenological model, namely the "toy" model of Refs. [45][START_REF] Böhmer | Broadband Dielectric Spectroscopy[END_REF][START_REF] Ladieu | [END_REF], where each dynamical heterogeneity (DH) is modeled as a point evolving in a two well potential of (dimensionless) asymmetry 0 δ 1. A key assumption, added in Ref. [START_REF] Ladieu | [END_REF], is that due to the absence of structural order, the net dipole moment of a given DH is proportional to √ N corr . This yields X (k) 3 ∝ N corr and, with some reasonable values of the free parameters (mostly N corr and δ), an acceptable agreement [START_REF] Ladieu | [END_REF] for the frequency dependence of X (k) 3 at 204.5 K. For f f α , we show in the Appendix that this model predicts the same T dependence for |X (1) 2;1 | as for |X (k) 3 | since all these cubic susceptibilities are found to be proportional to N corr (T ), consistently with Fig. 4. Besides, we show in Fig. 5 that the "toy" model can account for the frequency behavior of X (1) 2;1 : at 202 K we obtain a good fit with δ = 0.60 and N corr = 16, i.e., for parameters reasonably close to those used in Ref. [START_REF] Ladieu | [END_REF] at 204.5 K. The unphysical oscillation close to 0.1f α comes from the heuristic combination of the DH's and of the trivial contribution: this disappears when assuming that δ is distributed among the various DH, but this is at the cost of a supplementary parameter, see Fig. 5. Note that the toy model does not recover exactly Eq. ( 3) with λ eff = 1 since, as shown in the Appendix, χ (1) 2;1 contains two other-quite small-terms in addition to a term ∝∂χ Lin /∂T . Recovering Eq. ( 3) thus requires a deeper theoretical framework.

V. PERSPECTIVES AND CONCLUSION

Finally, our study opens the way for new experiments: due to the very small value of P, one could draw the N corr (T ) at high T 's from χ (1) 2;1 , avoiding the heating problems which limit the reliability of χ (k) 3 when f α is large enough [48]. Below T g , by switching on E st to 100 MV/m-the recipe to get this kind of high fields was found in Ref. [START_REF] Bauer | Fields of 67 MV/m were used successfully[END_REF]-δT g could be as large as 3 K, which will increase τ α by typically a decade. One could thus trigger aging phenomena which could be studied from the early stages since switching on E st is much faster than varying T . Having access to aging on a much larger number of decades in time will probably deepen our understanding of this phenomenon which is the hallmark of glasses.

To conclude, by monitoring a new kind of nonlinear susceptibility, we have shown that a static electric field E st shifts T g in glycerol. We found evidence of a quantitative relation between the two theoretical paths relating N corr to macroscopic observables, clarified the interpretation of nonlinear susceptibilities close to T g , and proposed experiments which should give hints about the physics of glasses.

In the box model, each dynamical heterogeneity (DH) is assumed to be heated by a quantity δT due to the applied field, see Refs. [42,43].

(1) In the case of a pure ac field E ac , one obtains [42] δT ac = 1 2

0 χ 1 E 2 ac c DH (ωτ DH ) 2 1 + (ωτ DH ) 2 1 2 0 χ 1 E 2 ac c DH , (A1)
where c DH is the volumic specific heat of the DH's, while all other quantities are defined in the main text. In the last equality of Eq. (A1), we took into account that nonlinear effects are dominated by the DH's such that (ωτ DH ) 1. This happens for two reasons: (i) these kind of DH's are those for which δT ac reaches its maximum value; (ii) these DH's are the most numerous ones, as soon as ω 1.12ω α , with ω α = 2πf α , see Ref. [44]. We neglect here the part of the heating oscillating at 2ω since it yields a subdominant contribution to the sought χ (1) 3;Box mod , see Ref. [43].

(2) In the case of an experiment measuring χ (1) 2;1 , the applied field is the sum of a static field and of an ac field. For the dipoles, on general grounds [48], the dissipated power is proportional to (∂P /∂t) 2 . As, to the leading order, P ∝ χ Lin (ω = 0)E st + |χ Lin (ω)|E ac cos{ωtarctan[χ Lin (ω)/χ Lin (ω)]}, one finds that the power dissipated by the dipoles is only proportional to E 2 ac and does not contain E st . Therefore, when varying E st and keeping fixed E ac , as we do in our experiments, the power dissipated by the dipoles does not change at all. This is at odds with what happens in an experiment measuring χ (1) 3 where the power dissipated by the dipoles is, according to the box model, the source of the cubic susceptibility. This is the main reason why we shall find below that within the box model one predicts |χ (1) 2;1;Box mod | |χ (1) 3;Box mod |. In a χ (1) 2;1 experiment, the dissipation in the liquid comes from the presence of a small quantity of free ions which yields a contribution χ ions ∝ 1/ω: this contribution becomes visible only at very low frequencies, i.e., below an angular frequency ω ions = tω α with t 2.5 × 10 -2 , the actual value of t depending of the purity of the supercooled liquid (glycerol here). It is convenient to write χ ions = s χ 1 ω ions /ω where s 4 × 10 -2 is a very small number, the actual value of which again depends on the purity of the liquid. Below ω ions , the fact that χ ∝ 1/ω amounts to a conductivity σ (ω ω ions ) = 0 χ ions ω which is independent on frequency. The dissipated power is thus

P = σ (E st + E ac ) 2
= σ E 2 st + 2σ E st E ac cos(ωt) + σ E 2 ac cos 2 (ωt), (A2) The right-hand side of Eq. (A2) contains three terms. The term proportional to E 2 ac has to be discarded since E ac is kept fixed in a χ (1) 2;1 experiment (this is the same reason as above for the power dissipated by the dipoles). The term proportional to E st E ac oscillates at 1ω and gives a nonlinear polarization at 1ω by modulating the ions conductivity term ∝σ E st . However, this term is not the leading one since it involves a factor (st) 2 1, contrary to the term that we calculate now and which is proportional to st [see below Eq. (A8)]. Indeed, focusing on the first term of Eq. (A2), we get, with σ = 0 χ 1 ω α ts, that the dominant contribution to the first harmonics nonlinear polarization comes from the term P st corresponding to the first term of the right-hand side of Eq. (A2):

P st = σ E 2 st = 0 χ 1 ω α tsE 2 st . (A3)
We shall assume that the ions are dispersed homogeneously within the various DH's (the opposite case is treated in Sec. A 3). By solving the power balance equation that one writes for each DH [42,43], one gets

δT st = τ DH P st /c DH . ( A4
)
Combining Eqs. (A1), (A3), and (A4), one can compare the heating produced by a static field and that produced by an ac field of equal amplitude

δT st δT ac E ac =E st = 2ω α τ DH ts < 5ts 5 × 10 -3 , (A5)
where the first inequality is obtained by replacing τ DH by its maximal value equals to 2.5τ α in glycerol [44]; and where the last inequality is obtained by using the maximal experimental values of the two small numbers t and s defined above. Even with these two successive overestimates, we find that δT st is more than two orders of magnitude smaller than δT ac .

Comparison of the nonlinear polarizations produced either by E ac or by E st

We now compare the nonlinear polarization δP produced by the two kinds of heatings calculated in the previous subsection. For each DH one has a nonlinear polarization δP DH given by [42]

δP DH = 0 ∂χ Lin ∂T δT , ( A6
)
where the heating δT depends on the value of the relaxation time τ DH of each DH, and is given either by Eq. (A1) or by Eq. (A4).

To obtain δP , one has to integrate Eq. (A6) over all the possible values of τ DH with a weight G(τ DH )dτ DH . For the two reasons given above just after Eq. (A1), as soon as ω 1.12ω α , the dominant terms for the calculation of the nonlinear effects are given by the DH's such that ωτ DH 1. We thus obtain

δP st δP ac τ 1/ω ∂χ Lin ∂T δT st G(τ DH )dτ DH τ 1/ω ∂χ Lin ∂T δT ac G(τ DH )dτ DH . ( A7
)
As the condition ωτ DH 1 ensures that both the real and the imaginary parts of each integrand in Eq. (A7) do not change sign over the integration domain, we shall overestimate δP st /δP ac by using the maximal value of δT st and the minimal value of δT ac , i.e., δT ac (ωτ DH = 1). We thus obtain by using Eqs. (A1) and (A5)

δP st δP ac < 10ts < 10 -2 (A8)
and therefore, taking into account the fact that the prefactors of the two terms of δP (t) containing χ (1) 3 and χ (1) 2;1 differ by a factor of 4, we conclude that the box model predicts that, for ω 1.12ω α , the values of |χ (1) 2;1 | are more than four hundred times smaller than those of |χ (1) 3 |. This is contradicted by Fig. 5 of the main paper.

Homogeneous heatings

Contrary to the previous paragraphs where each DH has its own temperature, we shall consider here the supercooled liquid as being homogeneously heated by the dissipated electrical power. Here the term "homogeneous" refers to the fact that we assume that the value of dielectric susceptibility is the same in every point of the liquid, hence the volumic dissipated power is the same everywhere in the liquid. For a pure ac field, we have thoroughly studied this question of homogeneous heatings in Ref. [48]: the heat propagation equation was solved within the samples, assuming that the electrodes were thick enough to be ideal heat sinks (i.e., with a zero heating). It was found that below 204 K, the nonlinear susceptibility due to homogeneous heatings were, at any frequency, more than 20 times smaller than the third harmonics cubic susceptibility χ (3) 3 . In the range 204-225 K, the heating contribution to χ (3) 3 increases but remains marginal as long as one does not consider angular frequencies much larger than ω α , see Ref. [48].

For the reasons already explained just after Eq. (A2) above, only two terms matter for the dissipated power in a χ (1) 2;1 experiment: the first term is P st = σ E 2 st which is constant in time; the second term is P 1ω = 2σ E st E ac cos(ωt): it oscillates at 1ω. As in the previous section, P st gives a first harmonics nonlinear polarization ∝ts 1 much larger than the nonlinear polarization ∝(ts) 2 coming from P 1ω . The fact that P 1ω is negligible is reinforced by the fact that it yields a homogeneous heating δT hom 1ω which is confined to the glycerol samples due to the skin effect which prevents ac heatings from extending over long distances [48]. This contrasts with the homogeneous heating δT hom st coming from P st : as it is constant in time, it is not exponentially damped at long distances. Therefore, the corresponding dc heat flux crosses the-unavoidable-insulating piece of thickness L ins which ensures the absence of electrical contact between each electrode and the experimental cell. In our case, this insulating piece is a kind of polymeric resin which conducts heat with a thermal conductivity κ ins , yielding

δT hom st σ E 2 st Se L ins κ ins S = 0 χ 1 ω α E 2 st tseL ins /κ ins , (A9)
where e and S are, respectively, the thickness and the surface of the sample (equal to that of the insulating piece).

We shall now contrast this expression of δT hom st with that of δT hom ac , where δT hom ac is the homogeneous heating in an experiment with a pure ac field. Focusing, for simplicity to the case ω = ω α , we have [48], with χ (ω α ) 0.3 χ 1 :

δT hom ac 0.3 0 χ 1 ω α E 2 ac e 2 /κ gly , (A10)
where κ gly is the thermal conductivity of the glycerol. As in our case one has κ ins ≈ κ gly we obtain

δT hom st δT hom ac tsL ins 0.3e 0.4, (A11)
where the last inequality was obtained by putting the numbers corresponding to our experiment, e.g., L ins = 3 mm and e = 16 μm. We thus conclude that, despite its static character which enlarges its value, δT hom st is weaker than the homogeneous heating δT hom ac obtained in a pure ac field experiment. Since, as long as T 225 K, one finds [48] that δT hom ac (ω α ) yields a nonlinear term more than four times smaller than |χ (3) 3 (ω α )|, and as |χ (3) 3 (ω α )| 0.1|χ (1) 3 (ω α )|, we get with |χ (1) 3 (ω α )| ≈ |χ (1) 2;1 (ω α )| that δT hom st yields a first harmonics nonlinear term which is at least two orders of magnitude smaller than |χ (1) 2;1 (ω α )|. Extending these arguments for all the frequencies and temperatures studied in the main paper, we find that δT st yields a first harmonics nonlinear polarization smaller than the ±4% accuracy corresponding to our experiments measuring χ (1) 2;1 (ω). Note that it would be easy to design an experimental setup where L ins would be as small as e: in that case one would have δT hom st δT hom ac . This shows that these χ (1) 2;1 experiments could be done, in the stationary regime, up to very high temperatures without having any heating problem. This strikingly contrasts with the pure ac field nonlinear experiments, in the steady state regime, where heating becomes an issue as soon as f α is above (typically) 10 kHz, see Ref. [48]. One could thus imagine first to "calibrate" the χ (1) 2;1 temperature dependence over that of χ (1 or 3) 3 in the temperature range where heating effects are negligible, and then to push the χ (1) 2;1 experiments to very high temperatures, taking advantage of the smallness of heatings with a static field. As stated in the main paper, this would allow to draw the temperature evolution of N corr in an unprecedentedly large temperature interval.

APPENDIX B: DERIVATION OF THE TOY MODEL'S PREDICTION FOR χ (1)

2;1

In this section we shall study the "toy model" for the dielectric response close to the glass transition. This toy model was studied in Ref. [START_REF] Ladieu | [END_REF] where its predictions for the cubic susceptibilities χ (1 or 3) 3 in a pure ac field experiment were computed in detail. This model expands on a previous work of Ref. [45] where each DH is represented as a point evolving in an asymmetric double well (ADW), independently of the state of all other DH's. In this model, the energy asymmetry is the difference between the minimal energy in the two wells, the deeper well being labeled by the value of the angle θ 1 between the net dipole moment μ of the considered DH, and the direction of the applied external field. For simplicity, it is assumed that the angle between the dipole moment and the external applied field is θ 2 = θ 1 + π in the well of high energy. The height of the barrier between the two wells, as well as /k B T with k B the Boltzmann's constant, sets the relaxation time τ ∝ 1/ cosh[ /(2k B T )] of the corresponding DH. The distribution of the τ s is assumed to be given by G(τ ) where G(τ ) is the suitable function allowing to recover the experimental frequency dependence of the macroscopic linear dielectric susceptibility. The decisive improvement made in Ref. [START_REF] Ladieu | [END_REF] was to assume that the net dipole moment of any DH is given by μ = μ molec √ N corr , with μ molec the molecular dipole moment: this expresses the fact that no structural order exists in a given DH. This directly gives [START_REF] Ladieu | [END_REF] that the macroscopic linear dielectric susceptibility χ Lin is independent of N corr , while the macroscopic cubic susceptibilities χ (1 or 3) 3 in a pure ac field experiment are directly proportional to N corr .

We shall now compute the prediction of this toy model for χ (1) 2;1 . Let us first consider a set of N identical ADW; i.e., a set of ADW's sharing the same values for all microscopic parameters of the model. When the formula for χ (1) 2;1 in this subset of DH's will have been derived, we shall just have to integrate over the various values of θ 1 , and of τ with the function G(τ ).

We first remind the reader the basic equations for a pure ac field E ac . Denoting n 1 (n 2 ) the number of ADW's in state θ 1 (θ 2 = θ 1 + π ), the polarization P of the considered set of ADW's is, with n = n 1n 2 :

P = (n 1 -n 2 )μ molec √ N corr cos θ 1 NN corr a 3 = M n N where M = μ molec cos θ 1 √ N corr a 3 . (B1)
Combining the two master equations for n 1 (t) and n 2 (t), with n 2 (t) = Nn 1 (t), one gets the equation for the dynamics of P : Adding a constant field E st is equivalent to changing the asymmetry by an amount 2μ molec √ N corr E st cos θ 1 . This amounts to changing the relaxation time τ by an amount ∂τ and the reduced asymmetry δ by an amount ∂δ, where the notation ∂ for differentiation was preferred to the usual δ to avoid the weird expression δδ. By using the exponent "bare" to label (when necessary) the quantities when E st = 0, we obtain ∝ (E ac ) q ∝ e q . This expansion was made iteratively taking advantage that in experiments e 1, which allows to neglect all P bare q >q terms when computing P bare q . As to compute χ (1) 2;1 we are only interested in terms ∝E ac E 2 st , we have only to expand P bare 1 (t) to second order in E st to obtain all the terms contributing to χ (1) 2;1 . In Ref. [START_REF] Ladieu | [END_REF] it was found that

∂δ = d 1 G + d 2 G 2 with d 1 = 1 cosh 2 bare 2k B T ; d 2 = -d 1 δ and G = μ molec √ N corr cos θ 1 k B T E st ∂τ = τ bare (T 1 G + T 2 G 2 ) with T 1 = -δ bare and T 2 = - 1 2 1 cosh 2 bare 2k B T -(δ bare ) 2 . ( B3
P bare 1 (t) = M(1 -δ 2 ) 1 + (ωτ ) 2 F cos [ωt -arctan(ωτ )] , (B4)
and when E st = 0 we thus have

P 1 -P bare 1 = ∂P bare 1 ∂τ ∂τ + ∂P bare 1 ∂δ ∂δ + 1 2
∂ 2 P bare 1 ∂τ 2 ∂τ 2 + 1 2
∂ 2 P bare 1 ∂δ 2 ∂δ 2 + ∂ 2 P bare 1 ∂τ ∂δ ∂τ ∂δ, ( B5
)
which means that we have to compute the five terms of the right-hand side of Eq. (B5) to get all the terms contributing to χ (1) 2;1 . This can be done readily from Eqs. (B3) and (B4). Below we write only the relevant terms, i.e., we keep only the terms proportional to G 2 ∝ E 2 st in each of the five terms (the terms proportional to G indeed vanish when averaging over θ 1 , see below). We obtain

∂P bare 1 ∂τ ∂τ = τ bare ∂P bare 1 ∂τ T 2 G 2 , ∂P bare 1 ∂δ ∂δ = -2δ 1 -δ 2 P bare 1 d 2 G 2 , 1 2
∂ 2 P bare 1 ∂τ 2 ∂τ 2 = 1 2 ∂ 2 P bare 1 ∂τ 2 (τ bare T 1 G) 2 , (B6) 1 2
∂ 2 P bare 1 ∂δ 2 ∂δ 2 = -1 1 -δ 2 P bare 1 d 2 1 G 2 , ∂ 2 P bare 1 ∂τ ∂δ ∂τ ∂δ = -2δ 1 -δ 2 τ bare ∂P bare 1 ∂τ T 1 d 1 G 2 .
As the previous calculations have been done for the DH's having the same relaxation time τ , the expression of the complex linear dielectric susceptibility χ bare Lin at E st = 0 is drawn from Eq. (B4) by using χ bare Lin = P bare 1 /( 0 E ac ). We thus obtain a relation holding both for the real and the imaginary parts of χ (1) 2;1 :

χ (1) 2;1 = G 2 3E 2 st 1 2 ∂ 2 χ bare Lin ∂τ 2 (τ bare T 1 ) 2 + τ bare ∂χ bare Lin ∂τ -2δ 1 -δ 2 T 1 d 1 +T 2 + χ bare Lin -2δd 2 -d 2 1 1 -δ 2 . (B7)
As Eq. (B7) holds for a given value of τ and of θ 1 [since θ 1 is both in G 2 and in χ bare Lin see Eqs. (B2)-(B4)], we have to average over these two quantities.

(1) Averaging over θ 1 : As the wells are assumed to be infinitely deep (see Ref. [START_REF] Ladieu | [END_REF]), we have to neglect the effect of E st on the distribution of θ 1 and to consider that θ 1 is isotropically distributed whatever E st . Therefore the average over θ 1 is straightforward, just as in Ref. [START_REF] Ladieu | [END_REF]. By introducing F bare Lin = χ bare Lin / χ 1 we obtain from Eq. (B7) and cos 2 θ 1 θ 1 = 1/3:

χ (1) 2;1 θ 1 = 0 ( χ 1 ) 2 a 3 k B T N corr 3(1 -δ 2 ) 1 2 ∂ 2 F bare Lin ∂τ 2 (τ bare T 1 ) 2 + τ bare ∂F bare Lin ∂τ -2δ 1 -δ 2 T 1 d 1 + T 2 + F bare Lin -2δd 2 -d 2 1 1 -δ 2 . ( B8
)
(2) Averaging over τ and N corr : Then, we have to average over the τ s with the weight G(τ )dτ already evoked above. Finally, we have to average over the various values of N corr which is likely to be distributed within the various DH's. We note [N corr ] av the corresponding average value of N corr : as in Ref. [START_REF] Ladieu | [END_REF] this yields the macroscopic χ (1) 2;1 to be proportional to [N corr ] av . More precisely, moving to the dimensionless version X (1) 2;1 of χ (1) 2;1 with the definition X (1) 2;1 ≡ χ (1) 2;1 k B T /[ 0 ( χ 1 ) 2 a 3 ] we obtain X (1) 2;1 =

[N corr ] av 3(1 -δ 2 ) J (ωτ α ), (B9)
where for brevity we have denoted J (ωτ α ) the integral with weight G(τ )dτ of the expression written between large brackets in the right-hand side of Eq. (B8). Note that due to time temperature superposition (TTS), we expect that the function J (ωτ α ) does not explicitly depends on the temperature. As in the toy model it is reasonable to choose that δ is a constant in temperature (see Ref. [START_REF] Ladieu | [END_REF]), we conclude that this model predicts that the temperature evolution of X (1) 2;1 (f/f α ) is directly that of [N corr (T )] av . This is what is explained in the main paper, with the simplified notation N corr (T ) instead of the more precise one [N corr (T )] av .

Another important feature of Eq. (B7) is that the second term of the second factor of the right-hand side is proportional to ∂χ bare Lin /∂T . Indeed, we have where we have used the TTS property ensuring that all the various τ s have the same relative variation in temperature. Therefore, the fact that the toy model is able to fit our data of χ (1) 2;1 is not so surprising since Eqs. (B7) and (B10) show that χ (1) 2;1 is related to ∂χ bare Lin /∂T . Of course these two quantities are proportional only when the first and the third terms of the second factor of the right-hand side of Eq. (B7) are negligible. It turns out that, numerically, these two terms play quite a small role. But their mere existence emphasizes the phenomenological nature of the toy model. As stated in the main paper, the ability of this model to fit our data does not change its phenomenological nature and that it cannot explain, from first principles, the existence of a deep relationship between χ (1) 2;1 and ∂χ bare Lin /∂T .

Two last points, mentioned in the main paper, are noteworthy.

(1) To take into account the finite lifetime τ ex of the DH's, we make a linear combination of the expression of χ (1) 2;1 derived above and of χ (1) 2;1;triv corresponding to the "trivial" case of a ideal gas of dipoles. The expression of χ (1) 2;1;triv is drawn from Ref. [29], and a global multiplicative factor close to 0.8 is allowed onto χ (1) 2;1;triv , just as in Ref. [START_REF] Ladieu | [END_REF]. For the same reasons as in Ref. [START_REF] Ladieu | [END_REF], χ (1) 2;1;triv is weighted by the factor 1exp (-f ex /f ) where f ex is a phenomenological parameter of the order of 1/τ ex . Of course, the expression drawn from Eq. (B7) is weighted by the complementary factor exp (-f ex /f ). As the best choice to fit our experiments turns out to be f ex 0.15f α (just as in Ref. [START_REF] Ladieu | [END_REF]), the influence of χ (1) 2;1;triv on the overall fit is really small as soon as f 0.5f α . Of course, χ (1) 2;1;triv plays a dominant role for f 0.1f α . (2) One can average over δ, which is likely to be distributed among the various DH's. As in Ref. [START_REF] Ladieu | [END_REF], we have chosen either a single value of δ, and the best value for fitting our data of χ (1) 2;1 is δ = 0.6 just as in Ref. [START_REF] Ladieu | [END_REF]; or we have taken a Gaussian expression with parameters δ 1 and δ 2 . Note that, as δ is restricted over the [0; 1] interval, δ 1 and δ 2 are in general neither the average of δ nor its standard deviation. Very surprisingly we found that the set of best values of δ 1 and δ 2 to fit our experiments are the same as those selected in Ref. [START_REF] Ladieu | [END_REF] when fitting the data of χ (1 or 3)

FIG. 5 .

 5 FIG. 5. (Color online) At 202 K, modulus (left axis) and phase (right axis) of measured dimensionless susceptibilities (symbols).Lines are fits of X(1) 2;1 with the toy model of Ref.[START_REF] Ladieu | [END_REF] by choosing either N corr = 16 and δ = 0.60; or N corr = 10 and δ distributed onto [0; 1] as ∝ exp [-(δδ 1 ) 2 /(2δ2 2 )] with δ 1 = 0.30 and δ 1 = 0.29, see the Appendix.

  δ sinh e + cosh e) = M (δ cosh e + sinh e) where e(t) ≡ F cos(ωt), F = μ molec √ N corr cos(θ 1 ) k B T E ac and M = μ molec cos θ 1 √ N corr a 3 , δ = tanh 2k B T . (B2)

)

  In Ref.[START_REF] Ladieu | [END_REF] P (t) was expanded in series of powers of the field E ac up to third order P bare (t) = P bare

	P bare 2	(t) + P bare 3	(t) where P bare q	0	+ P bare 1	(t) +

ACKNOWLEDGMENTS

We thank A. Dehaoui for experimental help during his stay, and C. Alba-Simionesco, G. Biroli, J.-P. Bouchaud, P.M. Déjardin, and G. Tarjus for long-lasting help. The support of the european grant "NPRGGLASS" as well as the support of the Institut des Systemes Complexes de la Region Ile de France are warmly acknowledged. 3 . The two continuous lines of Fig. 5 of the main letter correspond to the two kinds of fits (δ = 0.6 single valued and [N corr ] av = 16; or [N corr ] av = 10 and δ distributed with the above-mentioned Gaussian distribution).

APPENDIX A: OVERESTIMATE OF THE BOX

MODEL PREDICTION FOR χ (1) 2;1

In this section we shall compute nonlinear effects in the framework of the box model. As our experiments are in the stationary regime (see Sec. II B above), we do not compute the power dissipated by the switching on of the static field, but only consider the stationary dissipated power. We compare the case of a pure ac field E ac , where one obtains χ (1) 3 , with the case of a static field E st (plus a fixed amplitude of E ac ) where one obtains χ (1) 2;1 by varying E st . χ (1) 3;Box mod and χ (1) 2;1;Box mod will represent the predictions of the box model for χ (1) 3 and χ (1) 2;1 , respectively. As we shall find that χ (1) 2;1;Box mod is by orders of magnitude smaller than χ (1) 3;Box mod , it is sufficient to give an overestimate of dominant terms instead of performing a full calculation.