J Le Cochec

F Ladieu
email: ladieu@drecam.cea.fr

Onset of the nonlinear dielectric response of glasses in the two-level system model

Keywords: numbers: 61.43.Fs, 77.22.Ch, 72.20.Ht

We have calculated the real part χ ′ of the nonlinear dielectric susceptibility of amorphous insulators in the kHz range, by using the two-level system model and a nonperturbative numerical quantum approach. At low temperature T , it is first shown that the standard two-level model should lead to a decrease of χ ′ when the measuring field E is raised, since raising E increases the population of the upper level and induces Rabi oscillations cancelling the ones induced from the ground level. This predicted E-induced decrease of χ ′ is at odds with experiments. However, a good agreement with low-frequency experimental nonlinear data is achieved if, in our fully quantum simulations, interactions between defects are taken into account by a new relaxation rate whose efficiency increases as √ E, as was proposed recently by Burin et al. (Phys. Rev. Lett. 86, 5616 (2001)). In this approach, the behavior of χ ′ at low T is mainly explained by the efficiency of this new relaxation channel. This new relaxation rate could be further tested since it is shown that it should lead: i) to a completely new nonlinear behavior for samples whose thickness is ≃ 10 nm; ii) to a decrease of nonequilibrium effects when E is increased.

Amorphous materials exhibit universal anomalous properties at low temperature. In 1971, Zeller and Pohl [1] discovered below 1 K a quasilinear behavior of the specific heat in a number of glasses contrasting with the Debye law of crystalline materials. Anderson, Halperin, Varma [2] and Phillips [3] proposed an explaination based upon the existence of localized two-level systems (TLS). Their origin may be due to the tunneling of atoms or groups of atoms between two equilibrium positions separated by a narrow energy barrier featuring asymmetric two-well potentials. They are assumed randomly distributed in energy splittings and tunneling barriers as a consequence of the structural disorder of these materials. This model has proven to be successful to understand most salient experimental features.

The standard TLS model assumes defects do not interact with one another. However, defects are strongly coupled to their environment and can emit or absorb phonons. It leads to an indirect interaction between nearest neighbors via the phonon field [4]. Certain recent failures to explain nonequilibrium data (at a few kHz) [5] underscore the likely involvement of these interactions below 100 mK. However, these nonequilibrium effects are small corrections of the kHz stationnary response, and, up to recently, examples of stationnary susceptibilities strongly affected by interactions were very rare : in the kHz regime, it was argued that the ultra-low-T (T ≃ 1 mK) plateau of the dielectric constant in the linear regime, strongly different from the expected logarithmic increase, resulted from interactions [6]. Very recently, such a conclusion was drawn from internal friction experiments [7].

In this work, we show that including interactions in the TLS model with a recently proposed mechanism [8] strongly affects the nonlinear stationnary dielectric susceptibility χ ′ of amorphous insulators at a few kHz. A very complete set of such data was published a few years ago by Rogge et al. [9], twenty years after the pioneering work of Frossati et al. [10]. In the linear regime, χ ′ decreases when T decreases, reaches its minimum at T rev and then increases below T rev (before reaching the above-mentioned ultra-low-T plateau χ ′ plat ). According to the standard TLS model, the χ ′ decrease above T rev is due to the progressive freezing of the diagonal (or relaxational) part of the response, while the χ ′ increase below T rev comes from the induced off-diagonal (or resonant) part of the susceptibility : this effect enlarges as T decreases as do all quantum effects. However, due to the quantum nature of χ ′ below T rev , one expects χ ′ to be strongly depressed by a strong measuring electric field E at a given T . This can be guessed from the quantum saturation phenomenon which is very general in two level systems [START_REF] Cohen-Tannoudji | Atom-Photon Interactions[END_REF]. Indeed, increasing E decreases the population difference between the two energy levels : as the Rabi oscillations produced by E on the upper level are in phase opposition with respect to the ones produced on the ground level, the quantum response, once averaged on many independent TLS's, tends to zero when E is increased. Strikingly, Rogge et al. experiments show the opposite trend : χ ′ (T < T rev ) increases when E is increased.

As it is carefully explained in Ref. [9], this behavior does not result from heating of the sample by E. To give a supplementary argument with respect to Ref. [9], let us note that if E lin is the upper field below which the dielectric susceptibility is measured as being field independent, one expects that the heating of the sample, for a given E > ∼ E lin , is more important when T decreases. A heating effect is thus expected to strech the χ ′ (T ) curve of an amount increasing as T decreases, i.e. one expects

∂χ ′ ∂T E > ∼ E lin < ∂χ ′ ∂T E≤E lin , (1)
to hold at low T , i.e. mainly at T ≤ T rev. . As can be seen, e.g. on the Fig. 3 of Ref. [9], the trend of the data is exactly the opposite of Eq. (1). Finally, since heating effects can be ruled out, the fact that below T rev , χ ′ (E > ∼ E lin ) does not behave as expected from the quantum saturation phenomenon seems extremely intriguing in the framework of the standard TLS model.

However, this was not pointed out since the non linear effects in the TLS model were, up to now, only calculated by using the adiabatic approximation [11]. Such an approximation states that TLS's are at every moment at thermal equilibrium, i.e., it disregards any coherence effects. It predicts an increase of χ ′ with E, i.e. it qualitatively accounts for the experimental behavior. However, in the specific case of the real part of the susceptibility, the consistency of the adiabatic approximation is questionnable [12]. Indeed, as it is very clearly stated in Ref. [11], this approximation does not hold for TLS's whose Tunneling energy ∆ 0 is too small, and yet it finds that the nonlinear part of χ ′ is dominated by the smallest ∆ 0 values (see after Eq. (3.30) in Ref. [11]). More precisely [11], with p 0 ≃ 1 D the TLS dipole, even for the lowest electric fields E ≃ 1 kV/m of frequency ω ≃ 1 kHz, the adiabatic approximation fails when ∆ 0 < ∼ √ hωp 0 E ≃ 3 µK, while it is well known, from instationnarity experiments [5], that smaller Tunneling energies exist in glasses. Besides, the second puzzling point is that, according to the authors themselves [11], the reason of the increase of χ ′ with E in the adiabatic approximation is physically obscure, which leaves unsolved the question of the expected "quantum saturation effect" above mentionned. Finally, several predictions of Ref. [11] are somehow contradicted by experiments [9] : instead of the predicted T rev ∝ E γ with γ > 1, the measured data yield γ < ∼ 1/2; below T rev , at a given E, the predicted peaked behavior of ∂χ ′ /∂T is not observed; at very low T , the observed E dependence of χ ′ plat contradicts the predictions. This work goes beyond the adiabatic approximation, even though, due to the few simplifying assumptions that we have made (see Eq. ( 2)), we do not intend to yield a fully "from first principle calculation". The key point is that phase coherence is not discarded here since non linear effects are treated by a fully quantum non perturbative method. In the first part, we show that the standard TLS model cannot explain the low-frequency experimental data below 100 mK since it yields, at low T , the above-mentioned quantum saturation phenomenon. In a second part, interactions between defects are added by using an interaction mechanism proposed very recently by Burin et al. [8], and a successful agreement is obtained with experiments. Finally, we briefly discuss experimental predictions implied by Burin et al.'s interaction mechanism.

I. STANDARD TWO-LEVEL SYSTEM MODEL

A. Bloch equations of TLS

Dynamics of a unique isolated TLS

Consider a TLS that is sitting in a double-well potential and assume this defect has a dipole moment p 0 . Its energy splitting ǫ is related [15] to the asymmetry energy ∆ and to the tunneling matrix element ∆ 0 , describing transitions between the wells, by ǫ = ∆ 2 + ∆ 2 0 . Due to finite ∆ 0 , the eigen states extend over both sides of the TLS, and the position operator r is no longer diagonal in this eigen basis. As a result, when an external electric field E is applied to p 0 , the coupling Hamiltonian qE.r is not diagonal in the eigen basis [5] (upon which all the operators of this work are expressed), yielding a total Hamiltonian :

H = 1 2 ǫ 0 0 -ǫ + ∆ ǫ ∆0 ǫ ∆0 ǫ -∆ ǫ p 0 .E cos ωt,
or H = -s.Ω, with s = h 2 Σ where Σ are the three Pauli matrices and Ω is an external effective field (Ω components are given below, note Ω y = 0), which shows an effective spin operator s is associated to the TLS. The systematic use of "spin" language comes from the fact that the three Pauli matrices, combined with the identity matrix, form a general basis for TLS's. Whatever its physical nature, any operator can be expressed as a linear combination of these four matrices, e.g., the density operator ρ can be written : ρ = (1/2)I + (1/h)S.Σ, where S is the quantum mean value of the spin operator s . This shows that S x and S y describe the coherence effects contained in the off-diagonal terms of ρ, while S z is proportional to the population difference between the levels (the occupation probabilities are given by the diagonal terms of ρ).

The movement of p 0 and thereafter the dielectric response of the material stem from the dynamics of S. For a perfectly isolated TLS (note that this implies that T = 0) the evolution of S in the external field Ω is only a precession around the external field Ω, as can be seen from the Schrödinger equation which leads [5] to ∂S/∂t = S × Ω.

Dynamics of an ensemble of non-isolated TLS's

At finite T , the dynamics of the TLS must include the relaxation toward its equilibrium value since each TLS interacts with its environment (phonons or neighboring defects). Since these interactions occur randomly for a given TLS, the dynamical equation must deal with ensemble averaged properties S, i.e. with quantities averaged over many similar TLS's. This evolution is given by the Bloch equations, namely

∂ S ∂t = S × Ω + S-< S > relax τ relax , (2)
where the last term states that the relaxation of S toward the environement equilibrium values < S > relax must be added to the quantum dynamics (see Appendix A). In Eq. ( 2) it is assumed that the relaxation of a given S component, say Sx , occurs with a well defined time constant, say τ x . In the important case of short time scales, one needs to go beyond this approximation since echo signals do not generally decay as a simple exponential ( [13], [14]). This subtle effect is irrelevant here since, as already stated e.g. in Ref. [5], we are only interested in the long time range solution of Eq. ( 2), namely χ ′ (1kHz), i.e. we focus on the particular case ωτ 2 ≪ 1 (see below). Similarly the relaxation term of Eq. ( 2) might become more complicated in the case of very strong fields [START_REF] Geva | [END_REF], leading, e.g., to a Sy /τ x,y term in the relaxation of Sx (see Appendix B). However this should not be the case here since we only focus on the onset of the non linear regime (p 0 E will not much exceed k B T ). As a result, the relaxation terms can be derived quite simply, as we show now. i) Phonon induced relaxation. Let us first focus on phonon field relaxation. The occupation probabilities are altered by the emission or the absorption of phonons, yielding [15] a relaxation of Sz , with the relaxation time τ 1 = κ 1 /(ǫ∆ 2 0 ) tanh ǫ 2kB T , where κ 1 is a sample-dependent constant. Since phonon processes occur randomly and independently for various TLS's, they break the phase coherence of the ensemble of (noninteracting) TLS's, yielding a relaxation time 2τ 1 for Sx and Sy . What are the thermodynamic values < Sx,y,z > to which Sx,y,z relax ? By second order expansion of dynamical correlation functions, it was shown [START_REF] Abragam | The Principles of Nuclear Magnetism[END_REF] that this relaxation occurs towards the so-called "instantaneous equilibrium values", namely, < Sx,y,z (t) >= T r(< ρ (t) > Sx,y,z ) where < ρ(t) >= exp(-H(t)/(k B T ))/T r(exp(-H(t)/(k B T )) is the "instantaneous" thermodynamical density operator and k B is Boltzmann's constant. For this result to be valid, several conditions must be fullfilled, among which the most stringent one is, by far : |p 0 .E| τ c ≤ h where τ c is the correlation time of the random electrical field acting on a given TLS due to its small interactions with its neighbors (see next paragraph ii) ). Finally, these phonon processes yield in the Bloch equations a term ( Sz (t)-< Sz (t) >)/τ 1 for the population relaxation, and ( Sx,y (t)-< Sx,y (t) >)/(2τ 1 ) for the relaxation of the coherence terms. Does τ 1 depend on time ? On one hand, under the above stated assumption |p 0 .E| τ c ≤ h, it was argued [START_REF] Abragam | The Principles of Nuclear Magnetism[END_REF] that τ 1 does not depend on time (see also Ref. [17]). On the other hand, one may argue [11] that, since the applied electric field modulates the asymetry energy ∆, one should use τ

1 (t) = κ 1 /(ǫ ef f ∆ 2 0 ) tanh ǫ ef f 2kB T , where ∆ ef f = ∆ + p 0 E cos ωt and ǫ ef f = ∆ 2 ef f + ∆ 2 0
arise from the diagonalisation of the total t dependent Hamiltonian H. The use of τ 1 (t) is natural within the frame of the adiabatic approximation [11] where the system is assumed to be at thermal equilibrium at every instant. In our fully quantum approach, the question is much more difficult. In the particular case of the low frequency real part of the susceptibility, however, one can easily explain why using either τ 1 or τ 1 (t) lead to very similar results. Indeed, τ 1 (t) and τ 1 mainly differ only for the TLS's whose gap lie in the range ǫ ≤ p 0 E. But, as it will be shown in the insets of Fig. 1 and Fig. 3, the gaps of the TLS's driven in the nonlinear regime by a given E extend on a much larger domain (see sections I.B) and II.B)): this is one of the main results of our fully quantum approach. Thus the possible time dependence of τ 1 is not expected to change the results. This was carefully checked by performing all the calculations reported here twice, once using τ 1 , once using τ 1 (t): the resulting differences between both assumptions turned out in any case to be totally negligible. Hence, throughout the paper τ 1 is considered as time independent, by simplicity.

With the above relations, we get for the diagonal elements < ρ 1,1 (t) > and < ρ 2,2 (t) > :

< ρ 1,1 (t) >= 1 2 + Ω z 2 Ω 2 x + Ω 2 z tanh h Ω 2 x + Ω 2 z 2k B T ,
and

< ρ 2,2 (t) >= 1 2 - Ω z 2 Ω 2 x + Ω 2 z tanh h Ω 2 x + Ω 2 z 2k B T .
For its off-diagonal elements, it is found :

< ρ 1,2 (t) >=< ρ 2,1 (t) >= Ω x 2 Ω 2 x + Ω 2 z tanh h Ω 2 x + Ω 2 z 2k B T ,
where

Ω x (t) = -2 ∆ 0 ǫ p 0 .E h cos ωt, Ω z (t) = - ǫ h -2 ∆ ǫ p 0 .E h cos ωt.
Finally, one finds for the phonon field contribution:

< Sx >= hΩ x 2 Ω 2 x + Ω 2 z tanh h Ω 2 x + Ω 2 z 2k B T , < Sz >= hΩ z 2 Ω 2 x + Ω 2 z tanh h Ω 2 x + Ω 2 z 2k B T ,
and < Sy >= 0.

ii) "Spin-spin" induced relaxation Let us now turn to "spin-spin" interactions : for a given TLS, the effects of thermal transitions of its neighboring TLS's can be modeled as a small (fluctuating in time) electric field, i.e., as small fluctuating terms δH(t) ≪ ǫ, k B T . The latter inequality ensures that the relaxation of the population of the levels (involving Sz ) will not be sensitive to δH(t). It is shown in the Appendix A that, for a given TLS, the oscillations of S x,y (t) are no longer regular but progressively deformed by the random δH(t) terms: due to the absence of correlations between the δH(t) values seen by various TLS's, ensemble averaging leads, by cancelation of phases of many TLS's [21], to a relaxation of Sx,y to zero (while S x,y remains finite for any given TLS). This happens on a short characteristic time scale τ 2 ≪ 0.1ω -1 and yields a supplementary Sx,y /τ 2 for the relaxation of the coherence terms.

The temperature dependence of τ 2 is not clear at present : in echo experiments [START_REF] Piche | [END_REF], [29], both τ 2 ∝ T -1 as well as τ 2 ∝ T -2 were reported [30]. This might come both from the fact that accounting for the detailed shape of echo signals requires a very subtle theory (see e.g. [13]) and from the fact that several mechanisms contributes to τ 2 . Indeed, the pioneering work [22] of Black et al. predicted a τ 2 ∝ T -2 dependence but very recent calculations [23] based upon the mechanism used in part II found that τ 2 ∝ T -1 could be justified at low T . Since this new mechanism will be used in the last section, we use throughout this work τ 2 = κ 2 /T , where κ 2 is a sample dependent constant. In order to try to take into account the various mechanisms which might contributes to τ 2 , the parameter κ 2 will be widely varied, as can be seen in Fig. 2. Last, owing to the smallness of the p 0 E values considered here, we neglect any E effect on τ 2 as explained in Appendix B.

iii) Final form of the Bloch equations

Inserting the above relaxation terms in Eq.( 2), the three Bloch equations can be written as follows:

d Sx dt -Ω z (t) Sy + Sx -< Sx > 2τ 1 + Sx τ 2 = 0, (3a) d Sy dt -Ω x (t) Sz + Ω z (t) Sx + Sy 2τ 1 + Sy τ 2 = 0, (3b)
d Sz dt + Ω x (t) Sy + Sz -< Sz > τ 1 = 0, (3c)
where all the S/τ terms come from the relaxation processes, while all the Ω S terms arise from the quantum dynamics, i.e. from the fact that H and s do not commute.

Equations (3a) and (3b) also write

d S+ dt + iΩ z (t) S+ + S+ τ * 2 = iΩ x (t) Sz + < Sx > 2τ 1 , (4)
with

S+ = Sx + i Sy ,
and τ * 2 = 2τ1τ2 2τ1+τ2 . Let us note that τ * 2 appears due to the existence in eqs.(3a)-(3b) of the two terms Sx,y /(2τ 1 ). Even if they are required by consistency (see above and Ref. [START_REF]τ2 becomes infinite and the relaxation of both Sx and Sy[END_REF]), these two terms do not exist in the pionneering works accounting either for the small instationnarities [5] or for echo experiments [START_REF] Piche | [END_REF], [29], [30]. In fact these two terms play a negligible role in the nonlinear susceptibility. To show this, let us first note that as long as τ 1 > τ 2 , one gets τ * 2 ≃ τ 2 , i.e. the Eqs. (3a)-(3c) amount to the simpler Bloch equations used before (especially in pulse echo experiments). The key point is that, in the (∆, ∆ 0 ) plane, this domain where τ 1 > τ 2 is quite large : it is shown in the inset of Fig. 1 and in Ref.

[25] that this domain contains, at least, all the TLS's such that ǫ ≤ e 1,2 = (κ 1 T /κ 2 ) 1/3 . As shown in the inset of Fig. 1, e 1,2 ≃ 0.2 K is much larger than the p 0 E values studied in this work. This indicates that the TLS's standing out of the τ 1 > τ 2 domain should not be affected by E, i.e. they should be in the linear regime (see Ref. [START_REF]TLS's out of the τ1 > τ2 domain are in the linear regime, one can use Ref[END_REF]). To summarize, nonlinear effects should come mainly from the τ 1 > τ 2 region where the two terms Sx,y /(2τ 1 ) are negligible. This will be analytically demonstrated in section B)2).

Non perturbative resolution of the Bloch equations

The Bloch equations cannot be solved analytically and even their numerical resolution is so far a great challenge. However, in the audio-frequency range, some approximations can be made which strongly simplify the calculations. As τ * 2 is much shorter than the typical time (∼ 0.1 ω ) to modify the populations, Sz may be considered constant [START_REF] Cohen-Tannoudji | Atom-Photon Interactions[END_REF] in the right hand-side of Eq. ( 4). The coherence terms follow adiabatically the population evolution. They reach at every moment the stationary state corresponding to the "frozen" occupation numbers.

Therefore, Eq. ( 4) can be solved independently of Eq. ( 3). The stationary solution of Eq. ( 4) is

S+ = iΩ x Sz + < Sx > /2τ 1 iΩ z + 1/τ * 2 , (5)
which inserted into Eq. ( 3) leads to a differential equation for Sz :

d Sz dt + Ω 2 x /τ * 2 Ω 2 x + 1/τ * 2 2 Sz + Sz -< Sz > τ 1 = Ω x Ω z Ω 2 x + 1/τ * 2 2 < Sx > 2τ 1 , (6)
Sz (t) in Eq. ( 6) is expanded into its Fourier series to get its stationary state. The expansion is limited to a finite number of harmonics. This number, of the order of 10, is found a posteriori when a stable and accurate result is obtained. So the differential equation is equivalent to a linear system whose solutions are the harmonics Sn z . The inverse Fourier-transform gives the periodic evolution of Sz (t). The coherence terms Sx and Sy are deduced from Eq. ( 5) where Sz (t), the solution of Eq. ( 6), is inserted. Finally, the first harmonics S1

x of Sx (t) is sought, to be included into the dielectric susceptibility (see Eq. ( 7) below).

Indeed, the susceptibility [5] of a single TLS reads

χ = -2 |p 0 | |E| cos θ ∆ ǫ S1 z h + ∆ 0 ǫ S1 x h , (7)
and it must be averaged over the distribution of TLS's [5] and over the dipole-orientation angle θ to yield the total susceptibility of the sample :

χ = P ∆max 0 d∆ ∆0max ∆0min d∆ 0 ∆ 0 1 -1 d (cos θ) χ (∆, ∆ 0 , θ) . (8)
In the remainder of this article, we concentrate on the real part χ ′ of χ which is linked to the capacitance of the sample, i.e., to its dielectric constant ǫ r by :

ǫ r -1 = χ ′ ǫ 0 .
B. The quantum saturation effect: the quantum part of χ ′ (T ) is depressed by a E increase

Numerical results

We have used the standard values for amorphous-SiO 2 : p 0 = 1 D, P = 3 × 10 44 Jm -3 , κ 1 = 10 -8 sK 3 (all the energies in τ 1 taken in K), ∆ 0min = 10 -6 K, ∆ 0max = 10 K, ∆ max = 10 K. As explained above, we took τ 2 = κ 2 /T , where κ 2 was ranged from 3.10 -11 sK to 10 -7 sK, allowing to check our fundamental assumption ωτ 2 ≪ 1 provided T ≥ 0.5 mK. Last, the numerical relative accuracy of our simulations was, in any case, better than 10 -3 : this was checked very carefully, both by increasing the number of harmonics when solving Eq. ( 6) and by letting the successive integration procedures converge to better than 10 -4 . For each set of parameters E, T, κ 1 , κ 2 , ∆ 0,min at least 4 × 10 4 couples of (∆, ∆ 0 ) were computed.

The simulations are displayed on Fig. 1. The resonant response (low temperature) is strongly depressed by the drive level, while the relaxation contribution (high temperature) is little affected. This is at odds with the experiments [9] where increasing E leads to an increase of both the resonant response and of its slope |∂ǫ ′ r /∂T | below T rev . Let us note that the curve labeled "linear response" was obtained independently by a standard series expansion of the Bloch equations keeping only, as in Ref. [5], the terms proportionnal to E : as E is made very small, the nonlinear calculations very precisely converge towards the linear regime.

However, the extreme sensitiveness of the resonance to the external field is very striking. It decreases rapidly while |p 0 .E| << k B T . The low-temperature phase-coherent upturn is destroyed by its environment (the external field), although the perturbation is much smaller than any thermodynamical quantity, which suggests that this effect has a quantum origin. This is further confirmed by the inset of Fig. 2 showing the influence of T and τ 2 on δχ ′ (E, T ) = 1 -χ ′ (E, T )/χ ′ (0, T ) : for a given E, the smaller T , the larger δχ ′ , which is expected since quantum effects generally increase as T decreases. Similarly, δχ ′ is larger when κ 2 is made smaller, i.e., when quantum coherence is made more "fragile". Finally, the dimensionless δχ ′ appears to depend not only on E, T, κ 2 but also on κ 1 , and it is shown in the main part of Fig. 2 that all these dependencies are a universal function of a dimensionless scale η, namely, :

δχ ′ = 0.1 × √ η if η < ∼ 1 0.1 × ln(η) if η ≫ 1 with η = p 0 E k B T ( κ 1 T 2 κ 2 ) α , (9)
where α ≃ 0.45 ± .05 and ln(η) might be replaced by a power law of η with an exponent lower than 0.1. This universal δχ ′ (η) dependence holds only when the relaxational part of χ ′ can be totally neglected, i.e., well below T rev ≃ 50 mK : in Fig. 2, only data corresponding to T ≤ 10 mK have been plotted. For these low T , δχ ′ (η) remains universal even when (κ 1 , κ 2 , E) are varied over several decades. The factor κ 1 /(T 2 κ 2 ) in η becomes very large at low T , yielding nonlinear effects even for very small E: this expresses that the lower T , the smaller the onset field of the nonlinear regime, as already seen on Fig. 1. Let us mention that the data of Fig. 2 correspond to the particular case θ = 0.

ε onset e 1,2 =(κ 1 T/κ 2 ) 1/3 τ 1 > τ 2 τ 1 < τ 2
FIG. 1: Inset: At T = 10 mK, κ1 = 10 -8 sK 3 , and κ2 = 10 -8 sK, the domain of TLS's such that τ2 < τ1 is quite large and contains all the gaps smaller than e1,2 = (κ1T /κ2) 1/3 -see [25]-. Even at p0E = 5.12mK this domain is larger than the one of the TLS's driven in the nonlinear regime defined by ǫ ≤ ǫonset ≃ 70 mK (see Eq. ( 11c) ). Note that ǫonset ≫ p0E (p0E is the small black area very near the origin): this explains that the nonlinear effects are visible even at very low fields, as shown in the main Figure . Main Figure : Dielectric susceptibility of amorphous-SiO2 at 1 kHz vs temperature simulated at various fields -the value of p0E in Kelvin labels each curve-within the standard two-level system model with the following set of parameters: p0 = 1 D, κ1 = 10 -8 sK 3 , κ2 = 10 -9 sK, ∆0,min = 10 -6 K, ∆max = ∆0max = 10 K, P = 3 × 10 44 Jm -3 . The low-temperature response vanishes rapidly as the electric field is increased due to the quantum saturation phenomenon. The linear response was obtained by an independent perturbative method.

[p 0 E / (k B T)](T -2 κ 1 /κ 2 ) α 1 - χ'(E,T)/χ'(E=0,T) 0 0.2 0.4 1.E-08 1.E-06 1.E-04 1.E-02 A A, κ 2 = 0.03nsK A,T=4 mK A, κ 1 = 0.1nsK 3 p 0 E [K] [K] 1 - χ'(E,T)/χ'(E=0,T) FIG. 2: Inset : δχ ′ = 1 -χ ′ (E, T )/χ ′ (0, T ) plotted versus p0E (in kelvins).
Curve A corresponds to p0 = 1 D, κ1 = 10 -8 sK 3 , κ2 = 10 -9 sK, ∆0,min = 10 -6 K, ∆max = ∆0max = 10 K, P = 3 × 10 44 Jm -3 and T = 2 mK. The other three curves show the effect upon quantum saturation of the parameter which was changed with respect to A : increasing T , as well as decreasing κ1, decreases δχ ′ ; while decreasing κ2 increases δχ ′ , as expected due to the quantum nature of δχ ′ . Main figure : The various influences of the simulation parameters can be reduced to a universal function of the dimensionless variable η

= p 0 E k B T ( κ 1 T 2 κ 2
) α with α = 0.45 ± .05 numerically. The dashed line shows that δχ ′ ∝ √ η when η < ∼ 1. The various parameters were ranged over several decades : 10 -10 sK 3 ≤ κ1 ≤ 10 -8 sK 3 ; 3 × 10 -11 sK≤ κ2 ≤ 10 -7 sK; 10 -6 K≤ ∆0,min ≤ 10 -4 K; 10 -8 K≤ p0E ≤ 3 mK. The data of this figure correspond to the particular case θ = 0.

Physical interpretation

To further understand the universal δχ ′ (η) and demonstrate its quantum origin, let us briefly go into the structure of the Bloch equations. By using the identity Ω x < Sz >= Ω z < Sx >, Eq. ( 6) can be written :

d Sz dt + Sz τ z = < Sz > τ z,1 ,    1 τz,1 = 1 τ1 1 + 1 2 (Ωxτ * 2 ) 2 1+(Ωz τ * 2 ) 2 1 τz = 1 τ1 1 + τ1 τ * 2 (Ωxτ * 2 ) 2 1+(Ωz τ * 2 ) 2 , (10)
In Eq. ( 10), one gets at E → 0 : τ z = τ z,1 = τ 1 . As argued in section I)A)2), the nonlinear behavior should come from the TLS's such that τ 1 > τ 2 : in this case we see indeed from Eq. ( 10) that increasing E decreases τ z much more than τ z,1 . This strongly depresses the off diagonal suceptibility, as we show now.

Let us first derive, from Eq. ( 10), the critical value E * such that 1/τ z becomes larger than 1/τ z,1 : focusing on the gaps ǫ lying within the τ 1 > τ 2 domain, i.e. in the domain where τ * 2 ≃ τ 2 , E * is determined by the condition

τ 1 τ 2 Ω 2 x ≃ 1 + Ω 2 z τ 2 2
, yielding :

   p0E * ǫ = hT kB √ κ1κ2 if k B ǫτ 2 ≤ h (11a) p0E * ǫ = ǫ κ2 κ1 if k B ǫτ 2 ≥ h (11b) ,
where all the energies are expressed in kelvins. With the standard values κ 1 = 10 -8 sK 3 and κ 2 = 10 -9 sK, we see that p 0 E * is much smaller than ǫ. Indeed, for T = 10 mK we get p 0 E * /ǫ = 2.10 -5 for the smallest gaps following Eq. (11a), and, for example, p 0 E * /ǫ ≤ 3.10 -3 for the gaps ǫ ≃ k B T which follow Eq. (11b). Solving Eq. (11b) with respect to ǫ, for a given E, leads to a characteristic gap

ǫ onset = p 0 E κ 1 κ 2 1/4 , (11c)
where all the energies are in Kelvins. For the highest p 0 E ≃ 5.12 mK studied here, we get ǫ onset ≃ 70 mK. As shown in the inset of Fig. 1, ǫ onset is both much larger than p 0 E and corresponds to a domain smaller than the one defined by our assumption τ 1 > τ 2 .

To show that E * in Eq. (11b) is indeed the critical field for a given TLS, at which the kind of nonlinearities of Figs. 1-2 onsets, let us now compare χ ′ (E ≪ E * ) and χ ′ (E * ).

i) If E ≪ E * , we get from Eq. ( 10) τ z ≃ τ z,1 ≃ τ 1 . Solving Eq. ( 10) is straightforward and leads for the n th harmonics of Sz (t) :

Sn z = < Sn z > 1 + n 2 ω 2 τ 2 1 , (12)
where < Sn z > is the n th harmonics of < Sz (t) >. Remembering that the region of interest is ǫ < ǫ onset , it can be checked that ωτ 1 ≫ 1 for basically all the considered TLS's. This yields, from Eq. (12), Sz (t) ≃< S0 z >. Furthermore, since p 0 E ≪ ǫ due to Eqs. (11), we get < Sz (t) >≃< S0 z >, which, once combined with the identity

Ω x < Sz (t) >= Ω z < Sx (t) >, yields Sz (t) ≃ Ω z < Sx (t) > /Ω x .
Once reported into Eq. ( 5), this yields :

Sx (t) ≃ < Sx (t) > 1 + Ω 2 z τ 2 2 (Ω 2 z τ 2 2 ), (13)
where in the last factor the fact that Ω 2 z τ 2 2 ≫ τ 2 /(2τ 1 ), which holds for any reasonable set of (κ 1 , κ 2 ), was used to drop the term τ 2 /(2τ 1 ).

ii) For E = E * , we get from Eq. ( 10), τ z,1 ≃ τ 1 and τ 1 /2 ≤ τ z (t) ≤ τ 1 . The fact that τ z is now smaller than τ z,1 is responsible for the onset of nonlinear effects. This can be seen by setting τ z = τ 1 /2 throughout the electrical period. With this simplification, one gets, with a derivation similar to the one yielding Eq. ( 13) :

Sx (t) ≃ < Sx (t) > 1 + Ω 2 z τ 2 2 ( 1 2 Ω 2 z τ 2 2 ), (14)
The off-diagonal part of the response in phase with E is χ′ x ∝ S1

x /E : it is read directly from Eqs. ( 13)-( 14), remembering that < Sx >∝ E cos ωt. This yields χ′

x (E = E * ) ≃ 1 2 χ′

x (E ≪ E * ), where the factor 1/2 comes from the above relation τ z = 1 2 τ z,1 , which was a simplification of the case E = E * . The comparison of Eqs. ( 13)-( 14) is thus only semi-quantitative, but it yields the main two features of the quantum saturation phenomenon : first χ′ (E * ) < χ′ (E ≪ E * ), second this effect comes from the off-diagonal part of the susceptibility, i.e., it is purely quantum (the diagonal susceptibility χ′ z ∝ S1 z /E is much smaller than χ′ x due to the fact that ωτ 1 ≫ 1 below T rev ). We have here an example of quantum decoherence [33]. It is not surprising that these effects were missed by the adiabatic approximation mentionned in the introduction since, in this approach, τ 2 has disappeared, yielding for the nonlinear onset [9] no other possibility than |p 0 .E| ∼ k B T , as expected for a system at equilibrium. Moreover we have shown that the quantum saturation depends on the precise coupling of the three Bloch equations, i.e. of the fact that τ z evolves faster with E than τ z,1 : this is out of reach for the adiabatic approximation since it contains only one differential equation [11] instead of Eqs. (3a)-(3c). Finally, the results of Figs 1-2 do not depend on the precise microscopic mechanism involved in τ 2 , but only on the fact, well established by echo experiments, that, for a vast subclass of TLS's one has τ 2 ≪ τ 1 : this is the main reason of the E-induced depression of χ ′ of Figs 1-2.

Effect of the density of states

More can be learnt from Eqs. (11), and more precisely from Eq. (11b) which holds for the vast majority of the TLS's responsible for the nonlinear behavior. First, let us note that the onset field E * increases as κ 2 /κ 1 : this suggests that the depression of χ ′ , when E is increased, depends on E κ 1 /κ 2 , which, remembering that κ 1 /κ 2 is the square of a temperature, leads to the dimensionless scale p 0 E/(k B T ) κ 1 /(T 2 κ 2 ) as the natural parameter for the quantum saturation phenomenon. This dimensionless scale matches exactly the definition of η in Eq. ( 9).

Second, from the above discussion of Eqs. ( 13)-( 14), the TLS's such that ǫ ≤ ǫ onset are already in the saturation regime, while the gaps larger than ǫ onset are hardly altered by E. It is thus natural to consider the number of TLS's such that ǫ ≤ ǫ onset as an estimate of the amplitude of the quantum saturation phenomenon 1 -χ ′ (E, T )/χ ′ (0, T ), stating :

1 -χ ′ (E, T )/χ ′ (0, T ) ∝ ǫonset ǫmin P (ǫ)dǫ ∝ √ E ∝ √ η, (15)
where the last equality was obtained by using the above-stated relationship E ∝ η; while the second equality uses both Eq. (11c) and the fact that the energetic density of states P (ǫ) is a constant due to the standard distribution P (∆, ∆ 0 ) = P /∆ 0 . Equation (15) yields exactly Eq. (9) derived from the numerical simulations. This argument enables to state that the small corrections to the standard P /∆ 0 previously proposed only yield small changes to the behavior of Figs. 1-2 : this is true, e.g., for the slight depression of the density of states at small gaps derived by Burin [31], as well as for P /∆ 1+y 0 with |y| ≪ 1 proposed in Ref. [32]. To summarize this section I), solving the Bloch equations leads to the quantum saturation effect, i.e., to a strong decrease of the off-diagonal part of χ ′ when E is raised. This effect holds for a very large set of κ 1 and κ 2 -the main parameters of the model-, and it mainly comes from the TLS's such that ǫ ≤ ǫ onset < e 1,2 . For an ensemble of TLS's with a P /∆ 0 density of states, quantum saturation goes as E 0.5 , and such an exponent justifies a posteriori the nonperturbative character of the method used here. Last, the quantum saturation phenomenon onsets for fields E * ≪ k B T /p 0 , as seen from Eq. ( 9). It is thus non-negligible since the field is, in most experiments, decreased well below k B T /p 0 . However, in the literature, the trend of the data is systematically the opposite of the one of Figs. 12. Since -see Appendix B-more general Bloch equations, corresponding to larger E, should not qualitatively change the results of Figs. 12, we conclude that the standard TLS model cannot account for the basic features of the nonlinear experimental data in the kHz range.

II. ADDING INTERACTIONS

A. Burin et al 's mechanism

At this step, at least one drive-dependent parameter must be added into the model to explain the large discrepancy with the experimental data. Moreover, it must enhance the relaxation process at low temperature, since coherence is broken by the external field as shown in Figs. 12.

Recently, Burin et al. [8] proposed an additional field-induced relaxation mechanism. They show that the resonant dipole-dipole coupling, which is so small in glasses, can be strongly increased by a low-frequency electric field. Indeed, thermal excitations, which are at zero-field localized on each TLS, tend to delocalize by hopping to resonant nearest neighbors. This is due to the fact that resonant hopping demands both TLS's to have very close values of both ∆ and ∆ 0 : as the electrical field modulates the TLS parameter ∆, the probability of finding, for a given TLS, a resonant TLS, increases from a negligible value at very low E, to a non-negligible value above a threshold of the external field. Let us note that this mechanism transports energy; hence it can be treated as a new relaxation mode.

The frequency must be small for the electric field to have time to modulate the coupling parameters. This is of no consequence here, since our crucial assumption ωτ 2 ≪ 1, leading to Eq. ( 5), already restricts our work to the low frequency case. Another assumption is that the external field amplitude is smaller than the characteristic splitting energy ∼ k B T , in order to treat the field as a weak perturbation. The typical values of the frequency and |p 0 .E| are respectively 100 Hz and 1 mK but may be softened as a rigourous determination is out of reach.

When the electric field increases, so does the probability of finding a resonant neighbor close enough to yield tunneling with not too small a probability : the one-particle excitation will relax more rapidly at high E towards another site. One can show the relaxation rate is proportionnal to the square root of the drive level [8]. To include this new enregy relaxation channel, we set in Eqs.(3a)-(3c) τ -1

1 = τ -1 1,ph. + τ -1
1B where τ 1,ph is the phonon field induced relaxation mechanism used throughout section I) and where

τ 1B = B |p 0 .E| , ( 16
)
with the constant B = 10 -5 sK 1/2 for physically reasonnable parameters [8]. As a result, increasing E at any given T leads to an increase of the susceptibility χ ′ : this shows that Burin et al.'s mechanism is strong enough to overcome the decrease due to the "quantum saturation phenomenon". However, the agreement between the set of calculated curves (unreported) and the data is very poor since the net increase of χ ′ (T ) when E is increased is stronger at high T than at low T . This is due to the fact that, since relaxation dominates the total response, the most influent TLS's are such that ǫ ≤ k B T : their number enlarges with T and so does their supplementary relaxational response due to the new relaxation channel τ 1B .

To interpolate between Fig. 1 and Eq. ( 16) which appear as extreme cases, one might use the very general argument that interaction effects should disappear at high T , e.g., above 100 mK -see Ref. [20]-. This demands that the chosen τ 1B (T ) becomes infinite, i.e., negligible, at high T . Such a general requirement can be of course modeled by different laws but all the ones we tried gave the same kind of behavior for the susceptibility. This is why we report on the calculations using a simple law, namely,

τ 1B (T ) = τ 1B 1 -e -TB /T with T B = 15 mK, ( 17
)
where τ 1B is given by Eq. ( 16) and the thermally activated behavior models a dipole-dipole coupling constant of T B = 15 mK : the energy scale T B can be deduced from Fig. 3 of Rogge et al.'s data [9] on a-SiO x since χ ′ becomes T -independent below 15 mK even for E values ten times larger than the range of the linear regime. Of course, this T B scale can be adjusted empirically since the T where χ ′ becomes T -independent depends on the material. As the coupling constant goes as g/ |rr ′ | 3 and as [5], for a-SiO 2 , g ∼ 10 Knm 3 , we get a mean distance λ B between interacting dipoles of nearly 10 nm.

B. Numerical results

The modified-model predictions using Eq. ( 17) are displayed on Fig. 3. The values of |p 0 .E| have been limited to 10 mK because of the restrictions on both the Bloch equations and the field-induced mechanism. A trend completely different from the one of Fig. 1 is obtained at low temperature since an increase of the response is observed when the drive level increases.

By computing separately (unreported) in Eq. ( 8) the two terms of the right hand side of Eq. ( 7), we checked that χ ′

x behaves qualitatively as in section I) and that the new trend of Fig. 3 is due to the diagonal part χ ′ z . To explain this new behavior, one first note that τ 1B (T ) is now the upper bound of τ 1 , even for the numerous TLS's whose small ∆ 0 value lead, in section I), to a very large τ 1 . With ωτ 1B (T < ∼ T B ) < ∼ 1, the 1/(ω 2 τ 2 1 ) cutoff of Sz seen on Eq. ( 12) has now disappeared, i.e. the d Sz /dt term in Eq. ( 6) can be dropped, yielding : 8) and the same parameters as in the main part of Fig. 1. The calculations were done within a modified TLS model where excitations are no longer localized but can experience field-induced hops to neighboring sites, which is modeled by an additive relaxation channel (see the definition of τ1B(T ) in Eq. ( 17)). The data show a linear behavior at low enough drive levels (the p0E values label the curves), an evolution of Trev with E compatible with experiments and a substantial decrease of the T dependence of χ ′ at low T ; Inset : For p0E ≪ kBT , in the (∆, ∆0) plane, < S1 z > is not negligible only within the ǫ < 2T domain. Even for p0E = 0.8 mK, the hatched area where τz < τz,1 has a non negligible size with respect to this ǫ < 2T domain : this yields a supplementary T -dependent contribution to the diagonal susceptibility χ ′ z which overcomes the E-induced depression of χ ′ x seen on Fig. 1, and yields the E-enhanced χ ′ trend seen on the main part of the figure.

Sz (t) ≃ τ z τ z,1 Sz (t) , (18)
where τ z , τ z,1 are defined in Eq. (10). At E → 0, one has τ z ≃ τ z,1 ≃ τ 1 , yielding with Eq. ( 18), Sz (t) ≃ Sz (t) . With the additionnal remark that S1 z (ǫ < 2T ) ≃ hp 0 E/(4k B T ) while S1 z (ǫ > 2T ) ≃ 0, one gets, with the standard P /∆ 0 density of states, that χ ′ z (T ) ∝ + ln T : this is the trend seen above T rev . To explain the behavior below T rev , the key point is that for quite a large domain in the (∆, ∆ 0 ) one has τ z /τ z,1 < 1 : since this factor is T dependent, it will modify the T dependence just above derived for χ ′ z from Eq. ( 18). Focusing on the ǫ ≤ 2T gaps, we get, below T rev , τ 1 ≃ τ 1,B (T ) and

τ 1 ≫ τ 2 : the condition τ z /τ z,1 < 1 amounts to τ 1 τ 2 Ω 2 x ≥ Ω 2 z τ 2 2 , i.e. : ǫ ≤ 2p 0 E τ 1B (T )/τ 2 sin φ + cos φ with φ = arctan ∆ 0 ∆ , (19)
The τ z /τ z,1 < 1 condition is shown, as a hatched domain, in the inset of Fig. 3. Even for the lowest E studied here, it is not negligible with respect to the ǫ < 2T area. Since in the hatched domain one has

τ z /τ z,1 ≃ τ 2 Ω 2 z /(τ 1 Ω 2
x ), this factor remains T dependent even below T B when τ 1B (T ) has reached its maximum value : this is due to the fact that τ 2 remains T dependent even at very low T . With S1 z (ǫ < 2T ) ≃ hp 0 E/(4k B T ), integration of Eq. ( 18) within the hatched area yields a contribution δχ ′ z ∝ E 3/4 /T 1/2 . Thus: i) this term increases as T decreases; ii) δχ ′ z increases with E, i.e. it can overcome the E-induced depression of χ ′ x . Disregarding the slight difference -see [START_REF]The E 3/4 exponent is not exact since[END_REF]between the δχ ′

x ∝ -E 1/2 seen for the quantum saturation phenomenon and the δχ ′ z ∝ +E 3/4 , the linear regime of Fig. 3, up to p 0 E = 0.32mK can be seen as resulting from the compensation of both effects. At higher E, the δχ ′ z increase dominates over the E-induced depression of χ ′ x , yielding a net increase of χ ′ with E. Note that χ ′ z (E) becomes T independent when T ≤ p 0 E/k B : in this case, indeed, S1 z (ǫ < 2T ) is no longer T dependent. This yields the substantial decrease of the T dependence of χ ′ seen for the two highest E values on Fig. 3.

Last, the off-diagonal susceptibility χ ′ x ∝ S1

x mainly behaves as in section I), i.e. we recover the quantum saturation phenomenon yielding, when E is raised, both a decrease of χ ′

x and of the slope |∂χ ′ x /∂T |. With respect to section I) the quantum saturation effect is somehow weakened, which can be understood since, for a given E, the number of TLS's lying within the ǫ ≤ ǫ onset domain of Fig. 1 is larger than the corresponding one in the hatched area of Fig.

Finally, the variations of χ ′

x with T remain smaller than the ones of χ ′ z , excepted in the case where T < p 0 E/k B : the small T dependence of χ ′ (T < T B ; p 0 E > ∼ 5 mK) is thus the only case where χ ′ x dominates the T behavior of χ ′ on Fig. 3.

To summarize, the bigger the electric field, the smaller the field-induced relaxation time (see Eqs. ( 16)-( 17)), which enhances the relaxationnal part of the response, leading to a net increase of χ ′ with E at a given T . At a given E, when T decreases below T B , the χ ′ (T ) increase is due to the fact that τ 2 is still T dependent : this is, of course, out of reach for the adiabatic approximation where τ 2 has disappeared. Finally, inserting Burin et al.'s new relaxation rate in Bloch equations allow to account for the main trend of the nonlinear data : however, in this approach, the so-called "resonant" regime below T rev is not a coherent one but, mainly, a field-enhanced relaxation regime.

C. Comparison with experiments

On Fig. 3, one observes a pseudolinear regime up to p 0 E ≃ 0.05k B T where the dielectric response is quasiindependent on the external field. This value of the electrical field agrees with the experimental linear regime, which, depending on the materials, extends up to p 0 E/(k B T ) in the range [0.02; 0.12] (see Figs. 345of Ref. [9]). We checked that this pseudolinear regime comes from the form of τ 1,B ∝ E -β where β takes the highly nontrivial value 1/2. Setting lower values for β, such as β = 0.1, yields the quantum saturation phenomenon to dominate, leading to the same trends as in Fig. 1, at odds with experiments. Setting β = 1 leads to the tendency of Fig. 3 but with a linear regime reduced to p 0 E/(k B T ) < 0.01. The second key point is the trend of the reversion temperature T rev with E : using Eq. ( 17), i.e., β = 1/2, leads T rev to increase by a factor three when E = 30 × E rev , where E rev is the electrical field such that the nonlinearities onset at T rev . This is in good agreement with Fig. 3 of Ref. [9]. On the contrary, using β = 1 leads T rev to increase much faster with E : T rev (E = 30 × E rev ) = 30 × T rev (E = 0). Finally, the key role of β = 1/2 is somehow reminiscent of Eq. ( 9) where δχ ′ ∝ √ E, even if an analytical argument supporting this idea is still lacking.

With respect to experimental data, a failure, at this step of the discussion, is the ratio between the two slopes ∂ǫ ′ r /∂ ln T below and above the reversion temperature. In Fig. 3 this ratio is near -1.7:1 instead of -1:1 in most experiments. Furthermore, the low-temperature linear -susceptibility data tend to a T -independent plateau while they do not in our simulations. At very low temperature, interactions are likely to be so strong that the independent TLS model does not apply anymore, even with a renormalized relaxation time such as that of Eq. (17). A transition toward a dipole-glass was invoked to explain the behavior of the samples whose χ ′ no longer depends on T below a few mK. In this picture, dipole orientation is progressively frozen, which would lead to a plateau of the susceptibility [6], [START_REF] Wurger | From coherent tunneling to relaxation, dissipative quantum dynamics of interacting defects[END_REF] : by continuity, this would weaken the slope ratio near -1:1. Since the TLS model should not apply at very low T , it is not surprising that the plateau of the susceptibility measured in the nonlinear regime is not well accounted for by Fig. 3. Indeed, Fig. 3 does not show a completely T -independent plateau but only a substantial reduction of the T -dependence of χ ′ at low T : as stated in II)B), this is due to χ ′

x which still exhibits a small T dependence, even when χ ′ z has turned into its T independent regime. However, if, on Fig. 3, the susceptibility is frozen below a given T , one gets plateaus for χ ′ whose heights depend on E, as in experiments. Finally, pushing β toward 1 strengthens the tendency of χ ′ to become T independent at low T (unreported), even if β ≃ 1 leads to the above-mentioned discrepencies with respect to experimental data. Let us note that some materials (see Rogge et al. [9]) do not yield any sign of such a glass transition even at T = 0.6 mK.

D. New predictions

Let us move briefly to the physical predictions implied by Burin et al.'s mechanism. Remembering that the inequality ωτ 2 ≪ 1 allowed the key simplification for the derivation of χ ′ (E, T ) -see Eq. ( 5)-we restrict ourselves to the kHz range where this condition is fulfilled. Two main predictions can be done : i) τ 1B (T ) will be suppressed in samples whose thickness h is smaller than the distance λ B separating the quasisimilar TLS's required by Burin et al.'s mechanism. Indeed, at distances larger than h, dipolar interactions within the dielectric will be suppressed by the screening effect of the numerous electrons of the electrodes. Thus, if h < ∼ λ B , one should observe a non linear behavior such as the one calculated in section I) -see Fig. 1-, where the quantum saturation of the levels only remains. In other words, ranging h from a fraction of λ B to a few λ B in a series of samples and studying χ ′ (E, T ) should lead to a gradual transition from Fig. 1 to Fig. 3 if Burin et al 's mechanism is relevant, while it should not affect the non linear behavior in the standard TLS model. Note that such an experiment looks feasible due to the quite large value of λ B ≃ 10 nm, -see II)A)-. This is due to the fact that Burin et al.'s mechanism requires the two interacting TLS's to have both very close values of ∆ and very close values of ∆ 0 : these conditions are stringent enough to make λ B much larger than the distance between a given TLS and its nearest neighbor.

ii) The net relaxation frequency τ -1

1 + τ -1
1B of a given TLS increases as E increases. Thus, nonequilibrium data should be of smaller amplitude when E is raised. Indeed, they are currently interpreted as resulting from the very large τ 1 existing in any glass due to the subclass of TLS's whose energy barrier is so high that ∆ 0 is very small. These very "slow" TLS's have an extremely delayed response to any change of the external constraints, such as the d.c. electrical, or strain, field imposed to the sample: these TLS's yield an excess of states at low energy with respect to the equilibrium density of states, the latter having a small depression at low energies due to TLS-TLS interactions. To our knowledge, the influence of E on nonequilibrium phenomena has been reported only once, in Rogge et al.'s work devoted to nonequilibrium phenomena on a mylar sample [41]. Applying a relative strain field F to the sample leads to a sudden jump of the dielectric capacity C, measured at 5 kHz, followed by a logarithmic relaxation. At T = 11 mK, i.e., well below T rev , and with F = 2.7 × 10 -6 , the initial relative jump is dC/C = 13 × 10 -7 if the measuring field is E = 5 × 10 4 V/m (see Fig. 1 of Ref. [41]), while it decreases to dC/C = 4.5 × 10 -7 if the measuring field is E = 8.5 × 10 4 V/m (see Fig. 2 of Ref. [41]). Let us note that, with p 0 = 1 D and a relative dielectric constant of 5, E = 5 × 10 4 V/m amounts to an energy of 10 mK, of the order of T : in terms of our Fig. 3 this means that one stands just above the pseudolinear regime, i.e., in a regime where our calculations, as well as Burin's mechanism, should apply. Even if this was not investigated systematically, this single experimental datum favors the idea that nonequilibrium effects should be of smaller amplitude when E is increased, due to the interaction-induced reduction of the diagonal relaxation time.

III. CONCLUSIONS

In conclusion, we have simulated the nonlinear dielectric susceptibility of amorphous materials by using the TLS model. Phase coherence effects have been taken into account, which is the main difference with the adiabatic approximation. In the kHz range, the standard TLS model yields a nonlinear behavior at odds with experiments due to the field induced depression of the quantum response. However, it was possible to fit in many details the experimental low-temperature field-induced rising response by adding a new relaxation mechanism based upon the existence of interactions below 100 mK. In this approach, the low temperature response mainly loses its quantum origin at low frequency. Our work stresses the necessity to inject interactions into the TLS model to get satisfactory predictions.

where the electric field E rand is random in time for the considered TLS, and, at a given instant t, varies randomly for various TLS's. Note that Eq. (A1) is expressed in the eigen basis of the TLS.

Defining the density operator ρ(t) by :

ρ(t) = 1 2 + z x + iy x -iy 1 2 -z , (A2)
it is clear that x, y, z are, respectively, the quantum mean values of the three spin operators ( Sx , Sy , Sz are the corresponding symbols once the ensemble average over many similar TLS's is made). By using ih ρ = Hρ -ρH, where the dot stands for time derivation, the dynamics of x, y, z follows :

   ż = -Ω 1 y , hΩ 1 = -2 ∆0 ǫ p 0 E rand (A3a) ẋ = -Ω 0 y , hΩ 0 = ǫ + 2 ∆ ǫ p 0 E rand (A3b) ẏ = Ω 0 x + Ω 1 z (A3c)
.

To characterize the random fluctuations in time of E rand we model its autocorrelation function by <

E rand (t)E rand (t + t ′ ) > t = u 2 p 2 0 τc [θ(t ′ + τ c ) -θ(t ′ -τ c )]
where θ(t) stands for the Heaviside step function, τ c is the characteristic time scale of the fluctuations and u/ √ τ c the typical scale of the fluctuating part of the Hamiltonian H. This means that E rand (t) is drawn at random once every τ c and can be considered constant over time intervals [nτ c , (n + 1)τ c ], where n is an integer. Within each of these intervals, E rand (t) takes the constant value E n . This allows to solve exactly the equation for ÿ obtained from Eqs. (A3) : ÿ + (Ω 2 0,n + Ω 2 1,n )y = 0. This yields :

y(nτ c + t) = y(nτ c ) cos Ω n t + ẏ(nτ c ) Ω n sin Ω n t, (A4)
where Ω n = Ω 2 0,n + Ω 2 1,n with Ω 0,n and Ω 1,n defined as in Eqs. (A3) by setting E rand (nτ c + t) = E n . Inserting Eq. (A4) into Eq. (A3a) and Eq. (A3b), with the notation X n = X(nτ c ) for any quantity X, we get :

     x n+1 = x n - Ω0,nyn Ωn s n - Ω0,n ẏn Ω 2 n (1 -c n ) (A5) z n+1 = z n - Ω1,nyn Ωn s n - Ω1,n ẏn Ω 2 n (1 -c n ) (A6) ẏn+1 = Ω 0,n+1 x n+1 + Ω 1,n+1 z n+1 (A7)
, where s n = sin Ω n τ c , c n = cos Ω n τ c . The four equations (A4)-(A7) allow to deduce x, y, z at step (n + 1) provided the corresponding quantities are known at step n. Choosing the initial conditions x 1 , y 1 , z 1 , yields ẏ1 = Ω 0 x 1 + Ω 1 z 1 which allows to initiate the recurrence. Finally, let us note that choosing the initial quantum state as |Φ 1 >= a 1 |+ > + 1 -|a 1 | 2 exp(iϕ 1 )|->, where |+ >, |-> are the eigen states of the TLS, amounts to setting :

x 1 = |a 1 | 1 -|a 1 | 2 cos ϕ 1 , y 1 = |a 1 | 1 -|a 1 | 2 sin ϕ 1 , z 1 = |a 1 | 2 -1/2.
Figure 4 shows the dynamics of a TLS defined by ∆ = 1 K, ∆ 0 = 0.01 K evolving from the initial state a 1 = 1/2; ϕ 1 = π/2, i.e., from x 1 = 0; y 1 = √ 3/4; z 1 = -1/4. The random field characteristics were set to u/ √ τ c = 0.1 K and τ c = h/(4ǫ), i.e., τ c was chosen four times lower than the Bohr period. Without 'noise', y(t) exhibits the well-known regular Bohr oscillations (short-dashed line on Fig. 4). The effect of 'noise' is to deform these ocillations (continuous line on Fig. 4) by an amount increasing with time : as a result the periodicity of y(t) gradually disappears. This is illustrated in the inset of Fig. 4 showing the exponential decrease in time of the absolute value |C y | of the autocorrelation of y, defined by C y (t) =< δy(t ′ )δy(t ′ + t) > t ′ /λ 2 with δy(t) = y(t)-< y > and λ 2 =< (δy) 2 >.

Since the value y n depends on the set of values E n drawn for the considered TLS from n = 1, ensemble averaging (over many TLS's with the same ∆, ∆ 0 ) will lead to a cancelation of y due to the absence of correlations between the noise series seen by different TLS's. This cancelation happens on a time scale τ 2 which should be of the order of the one of C y shown in the inset of Fig. 4. This cancelation of y after ensemble averaging amounts to a supplementary relaxation term Sy /τ 2 in the Bloch equation describing Sy dynamics.

The dynamics of x(t) (unreported on Fig. 4) is similar to the one of y, yielding a corresponding relaxation term Sx /τ 2 . This contrasts totally with the dynamics of z(t), depicted on Fig. 4: provided the amount of noise δH(t) is much smaller than the gap ǫ, z(t) stands very close to its initial value z 1 , even at large times. In fact small fluctuations exist, with an autocorrelation decrease similar to the one of C y , but the key point is that |z(t)/ < z > -1| ≪ 1. Hence the 'noise' does not yield any supplementary relaxation term in the Bloch equation governing the population dynamics Sz . ). z, the quantum mean value of Sz, is basically constant (solid line with square symbols), i.e., mostly unchanged by the random electric field. On the contrary, y, the mean quantum value of Sy, is strongly affected by random electric field : the periodic Bohr oscillations (short dashed line) seen in the absence of random electric field, are progressively distorded when random electric field is present. Inset: As a result, Cy, the normalized autocorrelation function of y(t), decreases exponentially with time.

Interaction effects with a finite measuring field E.

When the measuring E is no longer zero, the whole dynamics should be recalculated, with the supplementary dipolar Hamiltonian corresponding to E. However, the fact that the measuring frequency ω is much lower than 1/τ 2 greatly simplifies the problem. Indeed, if ω were zero, taking into account of E would strictly amount to replace ∆ by ∆ + p 0 .E : with this new definition of ∆, all the previous calculations apply, yielding the same relaxation terms in the Bloch equations. We will assume that this holds true for finite ω, due to the fact that for the kHz frequencies considered here, the experimental values of τ 2 ensure ωτ 2 ≪ 1, even at the lowest T studied in the body of the paper.

B. Validity of Bloch equations.

The three Bloch equations Eqs. (3a)-(3c) are valid in the quasilinear response [34]. When the electric field becomes strong enough, the relaxation terms form a nondiagonal matrix, e.g. a Sz /τ x,z term might come into play in the first Bloch equation, and the corresponding Bloch equations are usually named in the litterature Generalised Bloch Equations (G.B.E.). However, up to our knowledge, these generalized relaxation terms have been calculated only in the case of transverse fields in the rotating wave approximation [START_REF] Geva | [END_REF]. This is at odds with our physical situation : i) The transverse field case amounts to ∆ = 0, which, by far, is not the case considered here (remember that, due to the 1/∆ 0 density of states, for most TLS's one has ∆ ≥ ∆ 0 ) ; ii) The measuring field E ∼ cosωt is an oscillating one, not a rotating one ∼ exp iωt and the rotating wave approximation would be valid only close to the resonance ω ≃ ǫ/h, a condition totally irrealistic here due to the extreme smallness of hω = 2 × 10 -7 K.

However, even if they do not apply in our case, one can use the G.B.E. derived in the rotating wave approximation for transverse fields to guess qualitatively what could be the influence of the off-diagonal relaxation terms. Two points are worth mentioning : i) One can easily check that the GBE still yield qualitatively the quantum saturation phenomenon, even if the off-diagonal relaxation terms are responsible for quantitative modifications. In particular, it was shown, in the limit of infinite E, that the GBE reduce to the standard Bloch equations with τ 2 = 2τ 1 and that one gets a vanishing susceptibility.

ii) In the GBE, the off-diagonal relaxation times become infinite (i.e. negligible) when τ c → 0, where τ c is the correlation time of the random field created, on a given TLS, by its neighbors. In the same spirit [37], in the GBE, τ 2 is affected by a multiplicative factor 1 + Ξ 2 τ 2 c where Ξ = |p 0 .E| /h is the Rabi frequency. The order of magnitude of τ c in glasses was measured only once by Devaud and Prieur [38] who found τ c ≃ 10 -8 s at T = 70 mK with an expected τ c ∼ 1/T temperature dependence. The E dependence in the relaxation times can be neglected if Ξτ c ≤ 1.

Aware of these limits, we guess the standard Bloch equations can give a fair approximation as long as |p 0 .E| does not exceed 0.1 -1 mK at low temperature. As an additional remark, the validity domain of our calculations extends as τ c decreases.

To summarize, the GBE do not suppress the quantum saturation phenomenon, on the contrary, they are intended to quantitatively account for the various measurable quantities in the saturation regime (such as linewidths, etc...). The problem of the strong depression of χ ′ when p 0 E is increased from extremely small values up to 10 -4 -10 -3 K is thus unavoidable and is at odds with Rogge et al.'s experiments [9] which were carried out on various glasses and showed absolutely no sign of field induced depression of χ ′ (T < T rev ), despite the fact that p 0 E was varied from 0.05 mK to 50 mK : the fact that the domain p 0 E ≤ 1 mK was experimentally investigated is of special importance since, as stated above, in this domain, at least, the Bloch equations used here should be valid.

τ z < τ z, 1 2TFIG. 3 :

 13 FIG. 3: Main Figure :Simulation of a-SiO2 susceptibility at 1 kHz vs temperature with Eq. (8) and the same parameters as in the main part of Fig.1. The calculations were done within a modified TLS model where excitations are no longer localized but can experience field-induced hops to neighboring sites, which is modeled by an additive relaxation channel (see the definition of τ1B(T ) in Eq. (17)). The data show a linear behavior at low enough drive levels (the p0E values label the curves), an evolution of Trev with E compatible with experiments and a substantial decrease of the T dependence of χ ′ at low T ; Inset : For p0E ≪ kBT , in the (∆, ∆0) plane, < S1 z > is not negligible only within the ǫ < 2T domain. Even for p0E = 0.8 mK, the hatched area where τz < τz,1 has a non negligible size with respect to this ǫ < 2T domain : this yields a supplementary T -dependent contribution to the diagonal susceptibility χ ′ z which overcomes the E-induced depression of χ ′ x seen on Fig.1, and yields the E-enhanced χ ′ trend seen on the main part of the figure.

FIG. 4 :

 4 FIG.4: Dynamics of a TLS (∆ = 1 K, ∆0 = 0.01 K) submitted to a random electric field (u/ √ τc = 0.1 K, τc is the quarter of the Bohr period h/ ∆ 2 + ∆ 2 0 ). z, the quantum mean value of Sz, is basically constant (solid line with square symbols), i.e., mostly unchanged by the random electric field. On the contrary, y, the mean quantum value of Sy, is strongly affected by random electric field : the periodic Bohr oscillations (short dashed line) seen in the absence of random electric field, are progressively distorded when random electric field is present. Inset: As a result, Cy, the normalized autocorrelation function of y(t), decreases exponentially with time.

Acknowledgments

Many thanks to P. Pari, P. Forget, P. Ailloud (CNRS/LPS) and P. Trouslard (CEA/INSTN/LVdG) for help in experiments motivating this theoretical work. Scientific dicussions with J. Joffrin and J.-Y. Prieur (CNRS, Orsay) turned out to be crucial for this work, as well as the indication of Ref. [27] by Pf. O. H. Rousseau (University of Perpignan). Useful discussions with D. Boutard, M. Ocio, M. Rotter, are also acknowledged.

IV. APPENDIX A. Phase decoherence induced by small TLS interactions

In this Appendix, we aim at giving some physical insight into the relaxation term introduced in the dynamics of an ensemble of TLS's due to their small mutual interactions. Expanding on the assumption that these interactions are much smaller than the other relevant energy scales (such as T or the gap ǫ), the basic idea [START_REF] Abragam | The Principles of Nuclear Magnetism[END_REF] is to model these interactions by a small random electric field acting on each TLS. This idea is not new [22], [START_REF] Abragam | The Principles of Nuclear Magnetism[END_REF], and numerical results are presented here only to help understand theoretical results.

Interactions effects when the measuring field E = 0

Consider first the case where the measuring field E = 0. Modeling mutual interactions between TLS's by a random electric field leads, for a given TLS, to a total Hamiltonian given, by :