Nonlinear response theory for Markov processes II: Fifth-order response functions

TL;DR

This paper extends nonlinear response theory for Markov processes to fifth-order, analyzing higher-order susceptibilities in models of dipole reorientations and trap models, revealing characteristic temperature regimes and the importance of variable coupling.

Contribution

It introduces fifth-order response calculations for stochastic models, providing new insights into nonlinear susceptibilities and their temperature dependence in glassy systems.

Findings

Humps observed in higher-order susceptibilities for most models.

Two characteristic temperature regimes for the asymmetric double well potential.

Susceptibility behavior depends on the variable coupled to the field.

Abstract

The nonlinear response of stochastic models obeying a master equation is calculated up to fifth-order in the external field thus extending the third-order results obtained earlier (G. Diezemann, Phys. Rev. E{\bf 85}, 051502 (2012)). For sinusoidal fields the $5\om$-component of the susceptibility is computed for the model of dipole reorientations in an asymmetric double well potential and for a trap model with a Gaussian density of states. For most realizations of the models a hump is found in the higher-order susceptibilities. In particular, for the asymmetric double well potential model there are two characteristic temperature regimes showing the occurence of such a hump as compared to a single characteristic regime in case of the third-order response. In case of the trap model the results strongly depend on the variable coupled to the field. As for the third-order response, the low-frequency limit of the susceptibility plays a crucial role with respect to the occurence of a hump. The findings are discussed in light of recent experimental results obtained for supercooled liquids. The differences found for the third-order and the fifth-order response indicate that nonlinear response functions might serve as a powerful tool to discriminate among the large number of existing models for glassy relaxation.

Equations

$${\hat{X}}_3(\omega ,T)=\left|{\chi }_3^{(3)}(\omega ,T)\right|\frac{k_{B}T}{(\Delta {\chi }_1{)}^2{a}^3}$$

$${\chi }^{(1)}(t)$$

$${\chi }^{(3)}(t)$$

$${\chi }^{(5)}(t)$$

$$X_3(\omega ,T)=\frac{T}{(\Delta {\chi }_1{)}^2}\left|{\chi }_3^{(3)}(\omega ,T)\right|\text{mbox}and{X}_5(\omega ,T)=\frac{T^2}{(\Delta {\chi }_1{)}^3}\left|{\chi }_5^{(5)}(\omega ,T)\right|$$

$${\dot{G}}_{kl}(t,t_0)=-\sum _n{W}_{nk}(t)G_{kl}(t,t_0)+\sum _n{W}_{kn}(t)G_{nl}(t,t_0)$$

$$⟨F(t){⟩}_{(H)}=\sum _{kl}{F}_k{G}_{kl}^{(H)}(t,t_0)p_k(t_0)$$

$$W_{kl}^{(H)}(t)=W_{kl}(t)e^{\beta H[\gamma M_k-\mu M_l]}$$

$${\mathbf{G}}^{(H)}=\mathbf{G}+\mathbf{G}\otimes \mathcal{V}\otimes {\mathbf{G}}^{(H)}$$

$${\mathbf{G}}^{(1)}$$

$${\mathbf{G}}^{(3)}$$

$${\mathbf{G}}^{(5)}$$

$$X_3(\omega ,T)=\frac{\Delta {\chi }_3}{\beta (\Delta {\chi }_1{)}^2}|S_3^{(3)}(\omega \tau )|$$

$$T_3=\Delta /\ln [(\sqrt{3}+1)/(\sqrt{3}-1)]\simeq \Delta /1.317.$$

$${\chi }_5^{(5)}(\omega )=\Delta {\chi }_5\times S_5(\omega \tau )\text{mbox}with\Delta {\chi }_5=\left(\frac{M^6}{112}\right){\beta }^5(1-{\delta }^2)$$

$$X_5(\omega ,T)=\frac{|{\chi }_5^{(5)}(\omega ,T)|}{{\beta }^2(\Delta {\chi }_1{)}^3}=\frac{\Delta {\chi }_5}{{\beta }^2(\Delta {\chi }_1{)}^3}\left|S_5(\omega \tau )\right|$$

$$S_5(0)=(2-15\cdot {\delta }^2+15\cdot {\delta }^4)/15$$

$$T_{5;a/b}=\Delta /\ln [(1+z_{a/b})/(1-z_{a/b})]\text{mbox}with{z}_{a/b}=\sqrt{\frac{15\mp \sqrt{105}}{30}}$$

$$T_{5;a}\simeq \Delta /3.145\text{mbox}with{T}_{5;b}\simeq \Delta /0.842$$

$$W({ϵ}^{\prime }|ϵ)=\rho ({ϵ}^{\prime })\kappa (ϵ)\text{mbox}withtheescaperate\kappa (ϵ)={\kappa }_{\infty }{e}^{\beta ϵ}$$

$$⟨M(ϵ)⟩=0;⟨M(ϵ)M({ϵ}^{\prime })⟩=\delta (ϵ-{ϵ}^{\prime })e^{-n\beta ϵ}$$

$${\chi }_3^{(\alpha )}(\omega )=\Delta {\chi }_3\times S_3^{(\alpha )}(\omega \tau )=\frac{M^4}{20}{\beta }^3\left(1-{\delta }^2\right)\times S_3^{(\alpha )}(\omega \tau )$$

$${\chi }_3^{(\alpha )}(0)=\left(\frac{M^4}{20\alpha }\right){\beta }^3\left(3{\delta }^2-1\right)\left(1-{\delta }^2\right)$$

$${\chi }_3^{(\alpha )}(0)\sim {\kappa }_4(M)={\mu }_4-3{\mu }_2^2=2M^4\left(3{\delta }^2-1\right)\left(1-{\delta }^2\right).$$

$$S_{2;1}(x)={\delta }^2\frac{(3-x^2)+ix(6+2x^2+x^4)}{(1+x^2{)}^3}-\frac{2+ix}{2(1+x^2{)}^2}$$

$$S_{2;1}(0)=3{\delta }^2-2\text{mbox}and{S}_{2;1}(\infty )=0$$

$$S_5(x)$$

$$N(x)=(1+x^2)(1+4x^2)(1+9x^2)(1+16x^2)(1+25x^2)$$

$${\chi }_5^{(5)}(0)\sim {\kappa }_6(M)$$

$${\chi }_5^{(5)}(\omega )=\frac{⟨\kappa ⟩}{32}{\beta }^5\left\{{\chi }_k(\omega )+{\chi }_{kl}(\omega )+{\chi }_{klm}(\omega )\right\}$$

$${\chi }_k(\omega )$$

$${\chi }_{kl}(\omega )$$

Full text

PACS: 64.70.P-, 64.70.Q-, 61.20.Lc, 05.40.-a

I. Introduction

The nonlinear response of supercooled liquids and glasses to electrical fields has been investigated intensively during the last decade both, experimentally and theoretically [1, 2, 3]. Motivated by a relation between the maximum value of the modulus of the cubic susceptibility and the volume of correlated domains[4] the data have been mostly analyzed accordingly and the number of correlated particles in a typical domain, $N_{\rm corr}$, have been extracted. The values obtained are compatible with what is to be expected if one assumes a growing length scale to be responsible for the heterogeneous slow dynamics of glass-forming liquids[5, 6, 7, 8, 9, 10, 11, 12].

Experimentally, a number of different techniques have been applied to obtain dynamical nonlinear susceptibilities[3]. In particular, and most important for the present work, it has been found that the scaled modulus of the 3$\omega$-response $\chi_{3}^{(3)}(\omega,T)$,

[TABLE]

where $\Delta\chi_{1}$ denotes the static linear response, $k_{B}$ the Boltzmann constant and $a^{3}$ the molecular volume, exhibits a hump-like structure. It is the height of this hump that has been related to $N_{\rm corr}$ for supercooled liquids[1, 13, 14]. Such a hump has not only been observed in molecular glass forming liquids, but also in a particular alcohol[13], in plastic crystals[15, 16] and a ionic liquid[17].

There have been some model calculations that aim at a quantitative description of the data obtained. The so-called box model was found to account for some of the features observed[18, 19]. A modification of a similar model shows results consistent with the experimental ones without the assumption of any spatial correlations[20]. A toy model considers a number of dipoles reorienting in an asymmetric double well potential and includes the so-called trivial reorientations at low frequencies[21].

The nonlinear response theory of dipole reorientations has been worked out for various models already some time ago[22, 23, 24, 25] and has also been extended to include long-range dipolar interactions[26]. In ref.[27], denoted as I in the following, I have computed the third-order nonlinear response for Markov processes using time-dependent perturbation theory for the propagator of the master equation (ME) governing the stochastic dynamics. In particular, I considered two simple stochastic models and computed the modulus of the cubic response for these. One model considers dipole reorientations in an asymmetric double well potential, described as a two-state model (ADWP model)[28, 29]. This model shows a hump in the reduced modulus in a narrow temperature range located around a characteristic temperature $T_{3}$ at which the low-frequency response vanishes. For temperatures below $T_{3}$, the height of the hump increases with temperature at variance with experimental findings. In the range $T>T_{3}$, the height decreases with increasing $T$. The other model I considered is the trap model with a Gaussian density of traps[30, 31, 32, 33, 34]. Depending on the particular choice of the variable that couples to the field a hump is obtained in most cases the height of which shows different temperature dependences. For specific variables it increases, for others it decreases or is essentially temperature-independent. The conclusion from these calculations is that it is indeed possible to obtain a hump in the modulus of the third-order response from such mean-field models without assuming the existence of glassy correlations.

Quite recently, the fifth-order susceptibility ($\chi_{5}$) of glycerol and propylene carbonate has been measured and also in this case a hump at a frequency somewhat below the inverse $\alpha$-relaxation time has been observed. The relation between the relative moduli has been interpreted in terms of the existence of compact dynamically correlated regions[35]. In addition, it is argued that a number of models for slow dynamics cannot account for the observed features, including the box model. This is very interesting because, as noted above, the box model has been shown to be able to describe the third-order response[18, 19] and it also has been used successfully to analyze nonlinear hole-burning data[36, 37]. Therefore, it appears possible that measurements of higher-order nonlinear response functions might provide additional information that allows to discriminate among different models for the slow relaxation in glass-forming liquids.

In the present paper, I will compute the fifth-order response function $\chi_{5}^{(5)}(\omega)$ for the same models considered in I. In particular, I will discuss the effect of sinusoidal fields of the form $H(t)=H_{0}\cos{(\omega t)}$. For this oscillating field the various response functions for times long compared to the initial transients can be written as:

[TABLE]

where $c.c.$ denotes the complex conjugate.

In the following sections, I will mainly discuss the scaled moduli of the maximum frequency-component,

[TABLE]

The prefactors $T/(\Delta\chi_{1})^{2}$ and $T^{2}/(\Delta\chi_{1})^{3}$ eliminate the trivial temperature dependences stemming from $\Delta\chi_{1}\sim T^{-1}$, and $\chi_{3}^{(3)}\sim T^{-3}$, $\chi_{5}^{(5)}\sim T^{-5}$. This allows to compare the different moduli directly. As outlined in the following section, the calculations are performed employing perturbation theory for the propagator of a master equation for Markov processes.

II. Fifth-order response functions

Here, I briefly discuss the calculation of the response of a dynamical system with a kinetics described by a master equation[38] using the same notation as in I[27].

Writing $G_{kl}(t,t_{0})$ for the conditional probability to find the system in state $k$ at time $t$ provided it was in state $l$ at time $t_{0}$, the ME has the form

[TABLE]

where the rates for a transition from state $k$ to state $l$ are given by $W_{lk}(t)$. The response of the system to an external field applied at time $t_{0}$ and measured by an observable $F(t)$,

[TABLE]

is determined by the one-time probabilities $p_{k}(t)$ obeying the same ME and given by $p_{k}(t)=\sum_{l}G_{kl}(t,t_{0})p_{l}(t_{0})$. Here, $G^{(H)}_{kl}(t,t_{0})$ denotes the propagator in the presence of the field and for the transition rates I use

[TABLE]

where $\gamma$ and $\mu$ can be chosen arbitrarily[28, 39, 40] and $\beta=T^{-1}$ with the Boltzmann constant set to unity. However, usually one considers models that obey detailed balance and therefore relates the parameters via $\gamma=1-\mu$.

A perturbation expansion is achieved via the expansion of the transition rates $W_{kl}^{(H)}(t)$ and a concomittant expansion of ${\bf G}^{(H)}(t,t_{0})$ in terms of the corresponding ’field-free’ propagator ${\bf G}(t,t_{0})$ using the decomposition ${\cal W}^{(H)}(t)={\cal W}(t)+{\cal V}(t)$ with ${\cal V}(t)=\sum_{n=1}^{\infty}{\cal V}^{(n)}(t)$, cf. I.

The perturbation expansion for the propagator starts from the Dyson-like equation

[TABLE]

where I omitted the time arguments and the convolution is abbreviated by ${\bf G}\otimes{\cal V}\otimes{\bf G}^{(H)}\equiv\int_{t_{0}}^{t}\!dt^{\prime}{\bf G}(t,t^{\prime}){\cal V}(t^{\prime}){\bf G}^{(H)}(t^{\prime},t_{0})$.

Using the series expansion for the transition rates, one easily finds for the relevant terms:

[TABLE]

Here, $\left\{{\bf G}\otimes{\cal V}^{(x)}\otimes{\bf G}\otimes{\cal V}^{(y)}\otimes{\bf G}\otimes{\cal V}^{(z)}\otimes{\bf G}\right\}$ is a shorthand notation for the sum of all permutations of $x$, $y$ and $z$.

In the next step, one uses the expression for the matrix elements of ${\bf G}^{(n)}(t,t_{0})$, $G^{(n)}_{kl}(t,t_{0})$, in eq.(5) in order to compute the nth-order response. For the calculation of the fifth-order response one then proceeds in the same way as in I[27]. I will not dwell on the general results further, but discuss the response functions obtained for specific models in the following.

III. The ADWP-model

As in refs.[27, 29], two dipole orientations characterized by polar angles $\theta_{1}=\theta$ and $\theta_{2}=\theta+\pi$ are assumed to describe the minima of the double well potential of the system. The flips of the dipoles take place with rates $W_{12}=We^{-\beta\Delta/2}$ and $W_{21}=We^{+\beta\Delta/2}$, where $\Delta$ denotes the asymmetry and the bare rate for $\Delta=0$ is given by $W$. If no field is applied, the solution of the ME reads as $G_{kl}(t)=p_{k}^{\rm eq}\left(1-e^{-t/\tau}\right)+\delta_{kl}e^{-t/\tau}$ with the relaxation rate $\tau^{-1}=2W\cosh{\!(\beta\Delta/2)}$ and the equilibrium populations $p_{k}^{\rm eq}=\tau\cdot W_{kl}$. The moment coupling to the field is $M_{k}=M\cos{\!(\theta_{k})}$, i.e. $M_{1}=M\cos{\!(\theta)}$ and $M_{2}=-M\cos{\!(\theta)}$ where $M$ denotes the static molecular dipole moment. Assuming isotropic distributions of the moments leads to an average over the angle $\theta$ and one has $\langle\cos^{n}{\!(\theta_{k})}\rangle=(n+1)^{-1}$ for $n$ even and $\langle\cos^{n}{\!(\theta_{k})}\rangle=0$ otherwise.

The linear response for the model is given by $\chi_{1}(\omega)=\Delta\chi_{1}[1-i\omega\tau]^{-1}$ with $\Delta\chi_{1}=\beta(M^{2}/3)\left(1-\delta^{2}\right)$ and $\delta=\tanh{\!(\beta\Delta/2)}$ and its properties have been discussed in I[27].

Also the results for the third-order response have been presented there already and I repeat them in Appendix A for completeness. One important finding is that the scaled modulus, eq.(3)

[TABLE]

with $\Delta\chi_{3}=\beta^{3}(M^{4}/20)(1-\delta^{2})$ shows a hump in a finite temperature regime (cf. Fig.3 of I), determined by the vanishing of the static nonlinear susceptibility $\chi_{3}^{(3)}(\omega\to 0,T)$, eq.(A.2) This determines the characteristic temperature

[TABLE]

The consequences of the existence of $T_{3}$ and the resulting hump in $X_{3}$ has been discussed in detail in I.

Another third-order response that has been investigated experimentally and that shows a very similar behavior as $X_{3}(\omega)$ is the response quadratic to a dc field and linearly to an ac field, $\chi_{2;1}(\omega)$[41]. Without going into details regarding the behavior of this function here, I just mention that it can be written in the form $\chi_{2;1}(\omega)=\Delta\chi_{3}\cdot S_{2;1}(\omega\tau)$ and the scaled modulus is given by $X_{2;1}(\omega)=(T\Delta\chi_{3})/(\Delta\chi_{1})^{2}\left|S_{2;1}(\omega\tau)\right|$ with the spectral function given in eq.(A.3). Due to the low-frequency limit of $S_{2;1}(\omega)$, $S_{2;1}(0)=(3\delta^{2}-1)$, one finds that $X_{2;1}(0)$ vanishes at the same characteristic temperature as $X_{3}$, $T_{2;1}=T_{3}$. Of course, this behavior is expected, as it is solely determined by the zero frequency limit of the respective spectral functions.

For the fifth-order response, one finds the following expression:

[TABLE]

where I again have performed the isotropic average. The spectral function $S_{5}(x)\equiv S_{5}^{(5)}(x)$ is given in Appendix A, eq.(Appendix A: Nonlinear response functions for the ADWP model). In Fig.1, the fifth-order response $\chi_{5}^{(5)}(\omega)$ is shown for different values of the asymmetry $\Delta$.

This can be compared directly to $\chi_{3}^{(3)}(\omega)$ shown in Fig.2 of I. The overall behavior is very similar. The corresponding scaled modulus is given by, cf. eq.(3):

[TABLE]

As in case of the third-order response functions, the modulus $X_{5}(\omega,T)$ shows a hump in the temperature regime where the low-frequency limit of the susceptibility vanishes, i.e. where $S_{5}(0)=0$. As this quantity is given by

[TABLE]

one finds two characteristic temperatures:

[TABLE]

which yields ($z_{a}\simeq 0.398$, $z_{b}\simeq 0.917$)

[TABLE]

In Fig.2 $X_{5}(\omega)$ is shown in the two relevant temperature regimes centered around $T_{5;a/b}$ on a logarithmic scale. As in case of the third-order response, a hump is observed in a rather narrow temperature range around $T_{5;a/b}$.

Note that in the temperature range around $T_{3}$ only trivial behavior is observed in $X_{5}$, cf. the black line in Fig.2. This means that the experimental finding that both $X_{3}$ and $X_{5}$ show a hump in the investigated temperature regime with decreasing hump height with increasing temperature[35] cannot be reproduced by the ADWP model in its form used here. This fact is further exemplified in Fig.3, where the relative height of the hump, $X_{k}(\omega_{\rm max})/X_{k}(0)$ for $k=3$, $5$ is shown.

One might argue that the experimental results for the third-order susceptibilities can be interpreted in terms of the ADWP model if one assumes that the relevant temperature range is above $T_{3}$. On the basis of the calculated fifth-order response, one would additionally have to assume that the relevant temperature range is above $T_{5;b}$.

IV. The Gaussian trap model

As in I, another model for glassy relaxation that I will consider is the well known trap model with a Gaussian density of states[30, 31, 32]. In brief, one considers the metastable states of a glass-forming liquid to be characterized by a free energy $\epsilon$ and assumes transitions among different values of $\epsilon$ to be given by

[TABLE]

This means that after an activated escape the destination trap is chosen at random, i.e. according to the density of states $\rho(\epsilon^{\prime})$. The propagator $G(\epsilon,t|\epsilon^{\prime})$ is obtained from the solution of the ME and all properties of the system can be obtained from this. As in I, it is assumed that the field couples to a variable $M(\epsilon)$, for which a Gaussian factorization is assumed to hold[42]:

[TABLE]

Here, $n$ is a model parameter. Of course, this choice is quite arbitrary and can be relaxed. However, eq.(15) has the advantage that the scaled linear response becomes independent of $n$, i.e. on the particularly chosen variable, cf. the discussion in I. As for the four point correlations needed for $X_{3}$, a Gaussian factorization property is assumed to hold for the six-point functions $\langle M(\epsilon_{1})M(\epsilon_{2})M(\epsilon_{3})M(\epsilon_{4})M(\epsilon_{5})M(\epsilon_{6})\rangle$ relevant for the computation of $X_{5}$. For more details about the actual calculation, I refer to Appendix B.

In Fig.4a, the results for $\chi_{3}^{(3)}(\omega)$ and $\chi_{5}^{(3)}(\omega)$ are shown for $n=0$ and two different temperatures.

It is evident, that the overall behavior is very similar. Fig.4b shows the moduli $X_{k}(\omega)$ as a function of temperature. Both, $X_{3}$ and $X_{5}$ show a hump, the height of which increases with increasing temperature, cf. the inset of Fig.4b. Thus, for $n=0$, i.e. energy independent variables, a similar behavior for the nonlinear responses is observed. Compared to the experimental observations, in particular the temperature dependence of the height of the humps is different.

As in I, the next step consists in varying $n$, which means to consider different dynamical variables. For $n=1$, the results for $X_{3}$ and $X_{5}$ differ because for $X_{3}$ it has been observed that a hump only exists above a certain threshold temperature and this is basically temperature-independent. This does not hold for $X_{5}$, as shown in Fig.5.

Here, the maximum height decreases as a function of temperature and it becomes very large for low temperatures. Therefore, for this choice, $n=1$, the humps observed in the $X_{k}$ behave differently at low temperatures and this only changes gradually with increasing temperature.

In I, $X_{3}$ also for $n=-1$ was considered, in which case a hump is observed that becomes smaller for increasing temperature. Without showing the results here, it is noticed that a similar behavior is found for $X_{5}$ with the relative height of the hump for $X_{5}$ being larger than for the third-order response. However, it was already mentioned in I that the case $n=-1$ is special in the sense that the mean relaxation time is temperature independent. Thus, I will not discuss this case further.

Another interesting choice is given by $n=-4$ because in this case the mean relaxation time $\langle\tau^{(n)}\rangle=\kappa_{\infty}^{-1}e^{(n+1)(n+3)/2\beta^{2}\sigma^{2}}$ concides with $\tau_{\rm eq}$, the corresponding quantity for $n=0$. In Fig.6, the moduli $X_{k}(\omega)$ (a) and the height of the humps (b) are plotted.

For this case, it is evident, that both nonlinear response functions show a behavior exhibiting a hump with a temperature dependent height that is decreasing with increasing temperature.

As discussed in detail in I for $X_{3}$, it is the low-frequency limit (cf. eq.(B.17)) that determines the existence of a hump-like structure also for $X_{5}$. These results obtained for $X_{5}$ show that different scenarios are possible ranging from a similar to a very different behavior of the different $X_{k}$.

V. Conclusions

In the present paper I have extended the calculation of nonlinear response functions for Markov processes presented in I[27] to the fifth-order response. In particular, I have considered the same examples of simple stochastic models, namely the ADWP model and the trap model with a Gaussian density of states.

For the ADWP model, I find that in the temperature range where $X_{3}$ shows a significant hump with a height that decreases with temperature the fifth-order response $X_{5}$ shows either trivial or the opposite behavior. Only for temperatures that are also higher than the higher characteristic temperature for $X_{5}$, one finds a decreasing height of the corresponding hump. In that temperature regime, however, the hump in $X_{3}$ will be very small or vanishes. It is therefore unlikely that the ADWP model in the present form could to be applicable in a straightforward way for the interpretation of the dielectric relaxation in the investigated systems.

In case of the trap model, a strong dependence of the behavior of $X_{5}$ on the choice of the dynamical variable coupling to the field is observed. This situation is similar to the corresponding one for $X_{3}$. For various values of the parameter $n$ the overall behavior of the hump is found to be similar for $X_{3}$ and $X_{5}$. If one chooses $n=1$, however, one even finds a different temperature dependence of the corresponding heights, cf. Fig.5b. A choice that might resemble the experimental findings partially would be $n=-4$ although no comparison is attempted in the present paper.

A note regarding the nature of the models considered appears in order. In the present paper (and also in I) only (mean-field) models for relaxation with a well-defined stochastic dynamics have been considered. In particular, no assumptions about the relation between the relaxation times and thermodynamic quantities (such as configurational entropy) have been made. Furthermore, I did not not assume any distributions of relaxation times. In phenomenological models, often a relation like the one assumed in the Adam-Gibbs model is assumed to hold[43]. This is a perfectly valid approach to obtain a meaningful parameterization of experimental data. However, from a theoretical point of view it appears more sound to start from given dynamical rules and calculate experimentally relevant quantities on that basis.

Calculations such as the ones performed in the preceeding sections can be further applied to other nonlinear sequences of external fields as they have been considered recently[41, 44]. Furthermore, in the framework of models similar to the ones considered here, it is possible to compare the response of the system after a temperature-jump with the effects of strong external fields. In a forthcoming publication, the analysis of higher-order nonlinear response functions will be performed for stochastic models obeying a Langevin equation instead of a master equation and the results will be compared.

The conclusion to be drawn from the calculations presented in the present paper is that the models considered here can only be used to describe the higher-order nonlinear response functions of supercooled liquids as observed experimentally for rather limited sets of parameters. This means that these response functions might be of great value when it comes to discriminate among various models for the dynamics.

Acknowledgment

I thank Roland Böhmer, Gerald Hinze, Francois Ladieu and Jeppe Dyre for fruitful discussions.

Appendix A: Nonlinear response functions for the ADWP model

$\chi_{3}^{(\alpha)}(\omega)$:

Here, I repeat the results for the third-order response for the two-state ADWP model given in I[27]. One has for $\alpha=1$, $3$:

[TABLE]

with $S_{3}^{(\alpha)}(x)$ only depending on the product $x=\omega\tau$ and given explicitly in eq.(15) in I[27].

When compared to the model of Brownian rotational diffusion, the following can be observed. For $\Delta=0$, $\chi_{3}^{(\alpha)}(\omega)$ is very similar to the corresponding expression for the model of rotational Brownian motion. For finite $\Delta$, however, the third-order response for the ADWP-model shows a characteristic temperature dependence, that is absent in the latter model.

The static nonlinear susceptibilites are determined by the limiting values of the spectral function, $S_{3}^{(\alpha)}(0)=(3\delta^{2}-1)/\alpha$, and thus is given by:

[TABLE]

I note, that $\chi_{3}^{(\alpha)}(0)$ is determined by the fourth-order cumulant, given here in terms of the central moments $\mu_{n}=\langle(M-\langle M\rangle)^{n}\rangle$,

[TABLE]

$\chi_{2;1}(\omega)$:

Here, for completeness and the convenience of the reader I give the spectral function for the response that is quadratic in a dc field and linear in an ac field with frequency $\omega$, a situation that has been considered in detail in ref.[41]. Again denoting $x=\omega\tau$, one finds in a calculation that is competely equivalent to the one performed in the case of an ac field:

[TABLE]

This spectral function has limiting values

[TABLE]

$\chi_{5}^{(5)}(\omega)$:

Here, I give the result for the spectral function determining the frequency-dependence of the $5\omega$-component of the fifth-order response to an ac field. The calculation is performed as outlined in Section II, using $\mu=0$ and $\gamma=1$, but the results are independent of this particular choice for a two-state model. One finds:

[TABLE]

with the denominator given by

[TABLE]

Similar to the situation for the third-order response, one has

[TABLE]

Appendix B: Fifth-order response functions for a trap model

In the present Appendix, I will give the expression for $\chi_{5}^{(5)}(\omega)$ for the Gaussian trap model also considered in I. The calculation is performed as outlined in Section II and also in App.B of I. As the results for the third-order response only weakly depend on the values of $\mu$ and $\gamma$ in eq.(6), I restrict the calculation of the fifth-order response to the case $\mu=1$. This means that the coupling of the field takes place only via the initial state of the transition. In view of the properties of the trap model this is a meaningful assumption because the field-free transition rates given by eq.(14) depend on the energy of the initial state of the transition but not on the energy of the destination state. Additionally, I again use the Gaussian factorization approximation for the variable $M(\epsilon)$ as has been employed in eq.(22) of I but now applied to products of the form $\langle M(\epsilon_{1})M(\epsilon_{2})M(\epsilon_{3})M(\epsilon_{4})M(\epsilon_{5})M(\epsilon_{6})\rangle$, yielding 15 independent terms if one assumes $\langle M(\epsilon)\rangle=0$ and $\langle M(\epsilon_{1})M(\epsilon_{2})\rangle=\delta(\epsilon_{1}-\epsilon_{2})e^{-n\beta\epsilon_{1}}$, cf. eq.(15).

Using a discrete notation (as in Appendix B of I), i.e. $G_{kl}(t)=G(\epsilon_{k},t|\epsilon_{l},0)$, $p^{eq}_{k}=p^{eq}(\epsilon_{k})$, $\rho_{k}=\rho(\epsilon_{k})$ and $\kappa_{k}=\kappa(\epsilon_{k})$, one finds the following expression for the fifth-order response:

[TABLE]

where the individual terms are given by:

[TABLE]

[TABLE]

[TABLE]

where I used the following definitions

[TABLE]

The spectral functions are defined by:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

with the individual spectral densities given by:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Here, $\omega_{ab...}=\omega_{a}+\omega_{b}+...$ and all frequencies sum up to $5\omega$, i.e. $\sum_{n}\omega_{n}=5\omega$. From these functions all others can be obtained in a straightforward manner by replacing the corresponding indices, e.g. $S_{kk}(\omega_{1},\omega_{2})=S_{k(l=k)}(\omega_{1},\omega_{2})$ etc..

Using the above results, it is straightforward to compute the zero-frequency limit of the fifth-order response which is given by:

[TABLE]

with

[TABLE]

the high-temperature limit of which coincides with $\overline{\langle M^{2}\rangle^{x}}$ according to eq.(B.5) because of $p_{k}^{eq}\to\rho_{k}$ for $T\to\infty$.