Exponential number of equilibria and depinning threshold for a directed polymer in a random potential

Introduction

Various aspects of the behaviour of a directed polymer, i.e. an elastic line, in a quenched random potential keep attracting permanent research efforts of both physicists and mathematicians for more than three decades. Among other applications, it was at the center of attention as a model for vortex lines in superconductors, leading to important developments in the physics of pinning (see Refs. for reviews). Its connection to the Kardar-Parisi-Zhang growth (see Ref. for review of earlier works) led to a recent outburst of interest, and it was shown that the probability density of the free energy for a long polymer converges to the famous Tracy-Widom distribution , extending the result for the ground state energy .

In this article we address a somewhat different aspect and consider the problem of counting the total number of equilibria for a directed polymer (DP), in general harmonically confined, immersed in a random potential. Those are defined as the stationary points (minima, maxima, or saddles) of an energy functional (see below). From a broader perspective, describing the statistical structure of the stationary points of random landscapes and fields of various types is a rich problem of intrinsic current interest in various areas of pure and applied mathematics . It also keeps attracting steady interest in the theoretical physics community, and this over more than fifty years , with recent applications to statistical physics , neural networks and complex dynamics , string theory and cosmology .

Note, however, that all the previous works considered only the case of zero internal dimension, equivalent to dealing with a single-particle embedded in a random potential of arbitrary dimension. In such a setting the counting of equilibria (a.k.a, stationary points or force-free configurations) can be placed in a framework of the standard Random Matrix Theory, see Refs. . In contrast we aim here to address the counting problem of manifolds of finite internal dimension. As was already anticipated in Ref. , in the latter case the problem turns out to be intimately related to properties of random Schrödinger operators appearing in the problems of Anderson localization. To the best of our knowledge, this aspect of the counting problem was never investigated before. Its treatment calls for a quite different technique and requires understanding of less studied properties of random Schrödinger operators, such as the modulus of its determinant and generalized Lyapunov exponents. We develop the corresponding approaches, mainly for the 1D case, in the present article.

Model and main results

We consider the following energy functional

where $u(\tau),\,\tau\in[0,L]$ describes the polymer configuration trajectory and $\kappa\geqslant 0$ is the elastic energy coefficient (cf. Fig. 1). Unless stated otherwise, in the main text of the paper we assume the fixed ends configuration $u(0)=u(L)=0$ for simplicity, other types of boundary conditions are briefly discussed in the A. The random potential $V(u(\tau),\tau)$ is chosen to be Gaussian with zero mean and with a translationally-invariant covariance

where we assume the symmetric function $R(u)$ to be at least four times differentiable at $u=0$. To have a better defined problem, the polymer is considered to be confined inside a harmonic well of curvature $m^{2}\geqslant 0$, called the mass parameter, which flattens the line beyond an infrared length, defined as

The limit $m\to 0^{+}$ is of special interest, as the system becomes critical, with a non trivial roughness exponent in the $L/L_{m}\gg 1$ limit .

2 The discrete model

Although most calculations will be performed within the continuous model, technically it will be often more easy to start from a discrete version of the model, passing to the continuous limit in the end of calculations. In this way we replace the continuous variable $\tau$ by a discrete lattice index $i=1,\ldots,K$ with $K=L/a$ (for simplicity we choose units such that the lattice spacing is $a=1$). The energy of the polymer and the correlations of the random potential in such setting are given by

3 Main results

In the present article, we provide some exact results for this problem. Namely, for model (1) we show that the mean total number of equilibria grows exponentially at large $L$, and that the rate $r$ is given in terms of the two length scales $L_{m}$ and $L_{c}$ defined above as

where $g(x)$ is a function calculated below. Although this result is derived for the continuous model (1) and Gaussian disorder, we argue that the function $g(x)$ is universal for a broader class of models, in the limit of weak disorder $L_{c}\gg a$ and small mass $L_{m}\gg a$ (where $a$ is a UV cutoff, such as the lattice spacing). We obtain the asymptotic behaviours

where the constant $c_{0}=1/(8\pi)$ has been calculated analytically (see E written by David Saykin), and we also obtained numerically both $c_{0}\simeq 0.04$ and

for the rate in the zero mass limit $m=0$ (i.e. $L_{m}=\infty$) : the rate of growth of the number of equilibria increases with the disorder strength as $r\sim R^{\prime\prime\prime\prime}(0)^{1/3}$. In the other limit of large confinement, the rate is exponentially suppressed r\sim\exp\big{[}-8m^{3}\sqrt{\kappa}/\big{(}3R^{\prime\prime\prime\prime}(0)\big{)}\big{]}, which shows that as long as the elastic line is shorter than the exponentially large scale, L\lesssim(L_{c}^{3}/L_{m}^{2})\,\exp\big{[}8m^{3}\sqrt{\kappa}/\big{(}3R^{\prime\prime\prime\prime}(0)\big{)}\big{]}, it is typically in a unique equilibrium configuration. This reflects exponentially rare metastable states of a strongly confined polymer induced by the rare events in the Gaussian tail of the disorder.

We can ask the question about the universality of our result Eq. (7). First it is immediate that (7) can be applied to the discrete model (5) in the limit $L_{m},\,L_{c}\gg a$ with parameters $\kappa$ and the same $R(u)$ function (in units such that $a=1$). Our result (7) is based on the asymptotic analysis of the mean number of equilibria, expressed in terms of a ratio of determinants, see Eqs. (21) and (27) below. Below, the analysis is performed in the continuum limit, which leads to evaluating some functional determinants with the Gelfand-Yaglom method. In fact, although we will not pursue it here, one can extend the Gelfand-Yaglom method to calculate the ratio of discrete determinants (21), see A. It is a particular model since it has uncorrelated, Gaussian disorder, and quadratic elastic energy. First it is clear from our derivation that the precise shape of the correlator function $R(u)$ (within the Gaussian class) does not matter for fixed values of the second and fourth derivative at zero. We expect furthermore that the scaling form (7), up to two non-universal scales, $L_{c}$ and $L_{m}$ (which, in some cases, can be independently measured), extends to a broader class of elastic line models (e.g. with short-range correlations in the $\tau$ direction, and with non Gaussian disorder). A remaining question for future work is to which extent the function $g(x)$ is universal.

Furthermore, the method of calculation developed here is interesting in itself as it reveals connections to other stochastic problems such as, multifractality-like scaling of moments for the solution of the initial value problem associated with the 1D Schrödinger operator with a white noise random potential, reflecting anomalous fluctuations of finite-size Lyapunov exponents, thermally activated dynamics of a particle near depinning, and probably more. For example the problem of calculating large deviation function addressed here is closely related to recent work in chaos theory , to be explored.

We discuss the question of counting of stable equilibria in Section 4.3.

We obtain the precise limiting behaviours of the large deviation function :

Finally, we further extend the upper bound for the depinning threshold $f_{c}$, to elastic manifolds of internal dimensionality $d$, relating it in Eq. (55) to the problem of evaluating the mean modulus of the spectral determinant of the Laplace operators with a random potential, as arises in Anderson localization problems.

4 Outline

5 Guideline

This brief description of the content of the article leads us to propose two possible ways to read the article :

Elastic line in disordered medium : Sections 3, 4, 6 and 11.

1D Schrödinger equation with a random potential and generalized Lyapunov exponent : Sections 7, 8, 9 and 10.

Counting of equilibria

As we shall see, the problem of counting equilibria for the energy function describing the lattice version of our model can be treated in a well-established mathematical framework Developing the corresponding formalism directly for the continuum model may represent an interesting mathematical problem.. An equilibrium configuration is found as a solution of the system of $K$ stationarity conditions which can be conveniently written as

where the Hessian is a $K\times K$ matrix given explicitly by

since we have assumed differentiability ($R(u)$ being obviously an even function). More generally, the property (b) is an important consequence of translational invariance and the Gaussian character of the random function $V_{i}(u)$ . Moreover, after taking the average the mod-Hessian factor is obviously independent of ${\bf u}$, and the average of each of the $K$ $\delta-$ factors can be done independently over the distribution of the Gaussian variable $V_{i}^{\prime}(u)$ with the variance $\langle[V_{i}^{\prime}(u)]^{2}\rangle=-R^{\prime\prime}(0)$, which gives

The Gaussian integral yields the constant Jacobian factor $J(m^{2})=|\det(m^{2}\,\mathbf{1}_{K}-\kappa\,\Delta)|^{-1}$ finally implying that

where the averaging goes over the set of independent and identically distributed (i.i.d.) mean-zero Gaussian random variables

with the covariance structure $\left\langle U_{i}U_{j}\right\rangle=2D\,\delta_{ij}$, where the parameter

measures the strength of the disorder in the problem, and is directly related to the Larkin length at $m=0$ as defined in Eq. (4).

The continuous version of the problem is related to the pair of Schrödinger operators

where $U(\tau)$ is the Gaussian white-noise potential with mean zero and covariance

Taking the appropriate continuous limit makes it rather apparent that the mean number of equilibria should be given in this case by a similar mean modulus of the ratio of two functional determinants for the operators (25) acting on functions vanishing at the two boundaries (Dirichlet boundary conditions) so that (21) takes the form

The calculation of these determinants will be performed in Section 7. In addition we will generalize the method to count equilibria at a given value of the energy in Section 11.

Counting of equilibria in the presence of a force – Depinning

Before evaluating the main objects of our interest, namely the ratios (21) or (27), which we postpone to Section 7, we discuss the effect of a uniform force field and show that the calculation of the rate $r$ is modified in a simple way. This will allow us to establish the connection to depinning threshold.

the no-crossing rule (or Middleton theorems) , known to hold for interface depinning, implies that in any given sample the last equilibrium which disappears upon increasing $f$ is a stable equilibrium,

the sample-dependent threshold force at which this happens has fluctuations decaying to zero at large $L$ (see e.g. Ref. ).

In this article, instead, we consider the ”annealed” rates

This is the force such that the mean total number of equilibria drops below unity in the interval of width $\sim w$, in the $u$-space. Note that the effect of the mass $m$ can be neglected as long as $w\ll\sqrt{|R^{\prime\prime}(0)|}/m^{2}$, and that in $d=0$ the number of equilibria is at most of order twice the number of metastable states. The square-root logarithmic dependence in the width $w$ is thus not surprising for the model of a particle as it originates from rare large barriers in the tail of the Gaussian-distributed potential, and is compatible with calculations of the depinning threshold in Ref. on related models. Due to this effect there is no true finite depinning threshold force for the particle for $w\to+\infty$ (except for a bounded disorder). We now turn to the elastic line, which does admit a well-defined depinning threshold force in the thermodynamic limit, even for the Gaussian disorder.

2.2 Elastic line (d=1𝑑1d=1) – discrete model

Let us now generalise the previous calculation to the DP model. Introducing the restriction

in the integrals in Eqs. (15) and (20), we need to calculate

We have introduced the two length scales (3) and

The formula (41) is valid for an arbitrary number of monomers $K$ and general boundary conditions, i.e. for all three types studied in A. Note that for $w=+\infty$ the result (40) is independent of ${\bf f}$ : this is because one can shift the Gaussian integration measure on ${\bf u}$ in (39) by ${\bf u}_{0}$, using that there exist (for $m>0$) a unique solution ${\bf u}_{0}$ to the equation $(m^{2}\mathbf{1}_{K}-\kappa\Delta){\bf u}_{0}={\bf f}$, which represents the new equilibrium position in the absence of disorder, displaced by the force.

2.3 Elastic line (d=1𝑑1d=1) – continuous model

Note that the dimension of $w^{2}$ is now $[u]^{2}[L]$. The extension of (41) is

where $-\partial_{\tau}^{2}\psi_{n}(\tau)=q_{n}^{2}\psi_{n}(\tau)$. The determinant in the denominator is easily obtained for various boundary conditions, Dirichlet, Neumann or periodic (cf. A) :

where the $\gamma$-independent prefactor is fixed by zeta-regularisation. We stress that the leading behaviour of the determinant is independent of the boundary conditions in the large $L$ limit

The modulus of the determinant in the denominator of (44,45) is

when $L\gg L_{m},L_{w}$. Finally we deduce

where $\Lambda(1)=\lim_{L\to\infty}(1/L)\ln\left\langle\left|\det\left(L_{m}^{-2}-\partial_{\tau}^{2}+U(\tau)\right)\right|\right\rangle$ will be determined in Sections 7 and 8. In this section we only need to known that taking $L_{m}\to\infty$ yields the value $\Lambda(1)=C/L_{c}$, Eq. (10), where $C$ is a dimensionless number of order unity and $L_{c}$ the Larkin length.

The form (49) is now appropriate to consider first the limit $m^{2}\to 0$ at fixed $w$, and second the limit $w\to\infty$ (which clearly do not commute). The first limit leads to $r_{w}(f):=\lim_{m^{2}\to 0}r_{w,m}(f)$, with

as displayed in the text. It is interesting to note that it has a similar order of magnitude as the Larkin-Ovchinnikov (LO) formula

with $R^{\prime\prime\prime\prime}(0)=|R^{\prime\prime}(0)|/v_{p}^{2}$, which, however is only an order of magnitude estimate. As discussed in the text, our result (52) is an exact upper bound for the true $f_{c}$ in the continuous model, or the discrete one in the limit $L_{c}\gg a$.

It is useful to recall that a similar robustness of the depinning threshold $f_{c}$ with respect to boundary conditions was observed in Ref. (and previous works cited there). There $f_{c}$ was studied for an elastic line on a cylinder of width $W=\alpha\,v_{p}(L/L_{c})^{\zeta}$, where $\zeta$ is the roughness exponent at depinning (and $m=0$) and $\alpha$ a dimensionless constant. The latter is measured from the roughness of the last metastable configuration encountered as $f$ is increased towards $f_{c}$. The value of $f_{c}$ was found independent of the aspect ratio $\alpha$ of the cylinder when both $L$ and $W$ become large. Only finite size corrections, which are subdominant, depend on the aspect ratio and other details. These subdominant sample to sample fluctuations of the depinning threshold force were also studied in Ref. .

2.4 Interface model on arbitrary graph and dimension d𝑑d

where $\Omega$ is the volume of the system and $U_{i}$ a Gaussian random potential with correlator \left\langle U_{i}U_{j}\right\rangle=\big{(}R^{\prime\prime\prime\prime}(0)/\kappa^{2}\big{)}\,\delta_{ij}. This formula further generalizes to an arbitrary graph.

where $\Theta(M)=1$ if all eigenvalues of the matrix $M$ are positive, and otherwise. Because $V^{\prime}_{i}(u_{i})$ and $V^{\prime\prime}_{i}(u_{i})$ are independent, Eq. (18), we have the same simplification as for the calculation of the total number of equilibria :

where we have also made use of the fact that $\partial_{i}\partial_{j}\mathcal{H}$ does not depend explicitly on ${\bf u}$, but only through $V^{\prime\prime}_{i}(u_{i})$, which are i.i.d. Gaussian random variables. Finally we can perform the same sequence of manipulations for a constant force

Numerical calculation of the depinning threshold

We have computed $f_{c}$ numerically for a discrete elastic chain of $K$ monomers using the algorithm developed in Ref. . The energy of the chain is given by Eq. (5) with $m^{2}=0$ and $\kappa=1$. The DP lives on a torus $(\tau,u)\in[0,L]\times[0,2\pi/v_{p}]$ (i.e. we choose periodic boundary conditions in the two directions). A convenient model for the disordered force is the harmonic model

where $\xi_{\alpha,i}$ are independent real Gaussian numbers of zero mean, $\left\langle\xi_{\alpha,i}\xi_{\beta,j}\right\rangle=\sigma^{2}\,\delta_{\alpha,\beta}\delta_{i,j}$. This leads to the correlations (6) for the periodic correlation function

This model provides a simple practical way to ensure that the disordered strength is translational invariant. Setting $v_{p}=1$ for convenience, the variance of the disorder is $-R^{\prime\prime}(0)=R^{\prime\prime\prime\prime}(0)=\sigma^{2}$. The presence of a stationary state is detected for a given force $f$. The equations of motion are (for $\kappa=1$)

The range of disordered strength considered in the simulation corresponds to $L_{c}/a=\sigma^{-2/3}$ between $1/\sqrt{2}$ to $2$. Thus it is quite remarkable that we do not observe any significant deviation from the behaviour $f_{c}\propto\sigma^{4/3}$ (Fig. 2), which is expected to hold in the continuum limit ($L_{c}\gg a$).

Equilibria configurations: spatial correlations in the annealed measure

It is interesting to ask what is the typical spatial configuration ${\bf u}$ of an elastic line at a force-free stationary point (equilibrium) chosen at random. In a fixed environment, it involves a quenched measure, which is difficult to study. It is simpler, but still quite instructive, to ask the same question over the set of all stationary points in all environments, i.e. to define the annealed measure

where $\rho({\bf u})$ was defined in (16). We also introduce the annealed averaging of any function of the monomers positions :

where $K$ is the number of monomers and $\left\langle\cdots\right\rangle$ denotes averaging over the disorder. The normalization in (63) obviously ensures that $\left\langle 1\right\rangle_{a}=1$. We thus calculate the generating function in presence of a source ${\bf j}$ as

Hence the annealed measure $\mathcal{P}_{a}({\bf u})$ over ${\bf u}$ is centered Gaussian and its two point correlation function is

It is then immediate to provide a few extensions. For instance

Relation with one-dimensional Anderson localization

The main outcome of Section 3 was the representation of the number of equilibria of the elastic line in terms of the determinant of a discrete random Schrödinger operator, Eq. (21), or a random differential operator (27). Such determinants are common in the problem of one-dimensional localization of a wave by a random medium. A famous example is the Herbert-Jones-Thouless relation (see also the appendix of Ref. for a discussion of the continuous case). This remark will allow us to make a precise connection with the study of the one-dimensional Schrödinger equation for a random (time-independent) potential.

in terms of the solutions $y(\tau)$ of the initial value (Cauchy) problem for the same operators : $(m^{2}/\kappa+H)\,y(\tau)=0$, i.e.

Cauchy problems related to products of random $2\times 2$ matrices of various types as well as their continuum limits were under active consideration recently, with many important analytical insights, see Refs. and references therein.

At this point it is appropriate to note that the spectral problem defined by (70) with $y(0)=y(L)=0$ is the classical problem of one-dimensional localization for the Schrödinger equation with a Gaussian white-noise potential, studied extensively since the seminal work , see Refs. for further details. To make contact to notations used for that problem in the literature we introduce

which plays the role of the energy in the Schrödinger equation (70). As is well-known the main qualitative feature of the Cauchy problem (70,71) is the exponential growth of its solution from chosen initial conditions in every realization of the disorder . It is conventional to define the localization length as the inverse of the Lyapunov exponent

where $\gamma_{n}L$ is the cumulant of order $n$ of $\ln|y(L)|$ at large $L$.

defined for $E<0$. For the study of the localization properties it is natural to consider the so-called semiclassical regime of large positive energy, $E\gg D^{2/3}$. In that regime the fluctuations can be considered as Gaussian with subdominant higher cumulants , $\Lambda(q)\simeq\gamma_{1}\,(q+q^{2}/2)$, where the Lyapunov exponent can be computed perturbatively , $\gamma_{1}\simeq D/(4E)$, thus

The property $\gamma_{2}\simeq\gamma_{1}$ is known as “single parameter scaling” (see also the discussion in Ref. ). The GLE is plotted in Fig. 3 for $q=1$. This regime would correspond to $m^{2}<0$ with $|L_{m}|\ll L_{c}$ while here, in the elastic line counting problem, we are interested in the opposite case $m^{2}\geqslant 0$ of negative effective energies $E<0$. The explicit evaluation of (72) remains an outstanding and non-perturbative problem where non-Gaussian fluctuations dominate (for early works see e.g. Refs. where it was analysed for integer value of $q$ ; see also Section 10).

2 Stochastic Riccati equation

A recursive method was developed in Ref. allowing to obtain integral representations for the $\gamma_{n}$’s in terms of multiple integrals, but they are quite complicated and not convenient to obtain limiting behaviours. An alternative way of calculating these quantities, suitable for the $E\to-\infty$ limit (large mass limit for the DP problem), was proposed in Ref. for a different model. We adjust this approach here: the transformation

relates (70) to the stochastic Riccati equation Note that this equation also describes the thermally activated motion at temperature $T$ of a particle near depinning with $-E=f-f_{c}$ and $D=T$

For $E<0$ and $U(\tau)=0$, $z=\sqrt{|E|}$ is a fixed point. The “noise” $U(\tau)$ thus generates fluctuations around this point.

We present first a perturbative approach to the determination of the GLE, valid for small disorder $D$ and large negative energy $E$. This method provides the main $q$-dependence of $\Lambda(q)$ in this regime, however it will turn insufficient in order to obtain the GLE for the specific value $q=1$, which will be shown to involve non perturbative contributions. We follow the method introduced in Ref. (section 6 of this reference) in a different situation : since in the limit $E\to-\infty$ the process $z(x)$ is most of the time trapped near $z=\sqrt{-E}$, this suggests to linearize the “force”, $-\mathscr{U}^{\prime}(z)=-E-z^{2}\simeq-2\sqrt{|E|}(z-\sqrt{|E|})$, leading to the Ornstein-Uhlenbeck process. The method developed below is a systematic perturbative expansion around the Ornstein-Uhlenbeck process. This can be most conveniently achieved at the level of the stochastic differential equation (SDE) (80) (this is more straightforward that from the Fokker-Planck equation (FPE)). For convenience, we write

and rescale the coordinate and the process as

Eq. (80) leads to the SDE in terms of dimensionless variables

where $\eta(u)$ is a normalized Gaussian white noise, $\left\langle\eta(u)\eta(u^{\prime})\right\rangle=\delta(u-u^{\prime})$, and

is the small perturbative parameter. We now expand the process in powers of $\epsilon$ as $\zeta=\zeta_{0}+\zeta_{1}+\zeta_{2}+\cdots$. The different terms can be found recursively

etc. Making use of the Gaussian nature of the Ornstein-Uhlenbeck process and of

We can check the method on the Lyapunov exponent: it leads to the expansion

which is in perfect correspondence with the analytic part of the expansion obtained from the exact result. The complex Lyapunov exponent for $D=1$ is (see also Ref. )

The Lyapunov exponent exhibits non-analytic contributions, \sim\exp\big{[}-8|E|^{3/2}/(3D)\big{]}, which are associated to the possibility of rare excursions of the process $z(x)$ to $\pm\infty$ related to the exponentially small probability current (this problem was not present in the case studied in Ref. by the same method).

where $\langle XY\rangle_{c}=\langle XY\rangle-\langle X\rangle\langle Y\rangle$. The limit $u\to\infty$ ensures that the correlator is computed in the stationary regime. Using the expressions (86,87,88), some lengthy algebra gives

The leading term to the third cumulant is more easy to compute

This approximate expression of the GLE is compared with numerics of Section 8 in Fig. 4 : the agreement is excellent. Quite remarkably, Eq. (98) shows that the rate $r$ vanishes up to third order in $s\ll 1$ (large mass limit). We conjecture that this remains true at all orders in $D$, i.e.

which is confirmed by numerics (see Figs. 4 and 6, and discussion below).

2.2 The need of a non perturbative analysis and a first estimate

For $E<0$, the potential (81) is not confining, thus the noise is not only responsible for small fluctuations around $+\sqrt{|E|}$, characterized by (98), but can also produce large excursions of the process at $\pm\infty$: if the process overcomes the potential barrier at $-\sqrt{|E|}$, it is rapidly driven towards $-\infty$, reinjected at $+\infty$, from which it eventually goes back to $+\sqrt{|E|}$. The rare jumps are separated by time intervals exponentially distributed . The probability rate for a jump is exponentially small :

We stress an important point : as noticed above, the Lyapunov exponent $\gamma_{1}$ provides a non-analytic contribution to the GLE $\Lambda(q)=\gamma_{1}\,q+\gamma_{2}\,q^{2}/2+\cdots$. This contribution, \sim\exp\big{[}-8|E|^{3/2}/(3D)\big{]}, cf. Eq. (93), is however much smaller than (101), which is therefore fully controlled by fluctuations. In the next section, we confirm the behaviour (101) by a more detailed analysis and provide the pre-exponential function.

Generalized Lyapunov exponent and spectral analysis

The study of the non-analytic contribution to the GLE requires powerful methods. In this section we use the relation with a spectral problem to develop some accurate numerical methods. The diffusion (80) is characterized by the (“backward”) generator $\mathscr{G}$ and its adjoint $\mathscr{G}^{\dagger}$ (“forward generator”)

where the initial condition is $z_{0}=\infty$ for Dirichlet boundary condition $y(0)=0$. We introduce the biorthogonal set of right and left eigenvectors

of the operator $\mathscr{O}_{q}$ . To make connection with the study of elastic line, we can rewrite the number of equilibria as

The FPE $\partial_{\tau}\mathcal{P}_{\tau}(z)=\mathscr{G}^{\dagger}\mathcal{P}_{\tau}(z)$ can be discretized in space as follows : we write $z=n\,b$ and introduce the transition rate from site $m$ to site $n=m\pm 1$

where $\mathscr{U}_{n}=\mathscr{U}(n\,b)$. The continuous time random walk is thus described by the master equation

(in the limit $b\to 0$ we recover the continuous diffusion for $D=1$ ; cf. Ref. for example). The boundaries must be discussed in detail. In order to mimic the absorption at $z=-\infty$ with reinjection at $z=+\infty$, we consider a finite lattice $n\in\{-M,-M+1,\cdots,+M\}$ and choose the following transition rates connecting the two boundaries

These equations define the $(2M+1)\times(2M+1)$ matrix \big{(}\mathscr{G}^{\dagger}\big{)}_{n,m}. Adding $(q\,n\,b)\,\delta_{n,m}$ we obtain the matrix \big{(}\mathscr{O}_{q}\big{)}_{n,m} and perform exact diagonalization.

that $r=\Lambda(1)-k$ is non-analytic in the small parameter $D/k^{3}$,

its main exponential behaviour is $\sim\exp[-4k^{3}/(3D)]$, corresponding to (101).

We find the value $r=\Lambda(1)\simeq 0.59\,D^{1/3}$ for $E=0$ (i.e. $m^{2}=0$).

For $q=0$, when $\Lambda(0)=0$, we have shown that the two coefficients are equal, $A_{+}(0)=A_{-}(0)$, cf. Eq. (322). We have verified numerically that this property remains true for $q>0$, see an example of solution in Fig. 15. This key observation has led us to propose the following method for the determination of the lowest eigenvalue $\mathscr{E}_{0}(q)=-\Lambda(q)$ : we solve (115) and choose $\lambda$ large enough to get the power law behaviours $\varphi(z)\sim|z|^{-2-q}$. Then, decreasing progressively $\lambda$, the eigenvalue is given when the two coefficients match exactly :

We have compared the numerical values obtained in this way with the one deduced by direct diagonalization of the matrix \big{(}\mathscr{O}_{q}\big{)}_{n,m} (previous Subsection): we have obtained a perfect agreement for the smallest values of the parameter $k$ (Fig. 6), when diagonalization is reliable (cf Fig. 5). The method is sufficiently accurate to make accessible relatively large $k$ (up to $k=2.5$) leading to the precision on the eigenvalue, and thus the rate $r=\Lambda(1)-k$, up to $\sim 10^{-11}$.

We obtain the value of the GLE in the limit $E\to 0^{-}$ (i.e. $m^{2}\to 0^{+}$) :

where we have obtained numerically $a_{1}\simeq 0.47$ (cf. Fig. 3). Correspondingly, we get

In the limit $E\to-\infty$ (i.e. $m^{2}\to+\infty$), we confirm once again the main exponential behaviour (101). Moreover, the method allows to extract unambiguously the behaviour of the pre-exponential function : cf. Fig. 8. The results of the numerical calculation are compatible with

This asymptotic behaviour corresponds to the limiting behaviour (8) for $x\to 0$. In E, written by D. Saykin, it is demonstrated that the above coefficient $0.08$ is actually $1/(4\pi)$.

3 Higher eigenvalues

The other eigenvalues of the operator $\mathscr{O}_{q}$, $\mathscr{E}_{n}(q)$ for $n>0$, are also of interest as they control the relaxation towards equilibrium. We have obtained it by a direct diagonalization of the discretized operator \big{(}\mathscr{O}_{q}\big{)}_{n,m} : we have plotted the three first energies in Fig. 7 (we have also checked that the diagonalization of the discretized operator \big{(}\mathscr{H}_{q}\big{)}_{n,m} gives the same result). Interestingly, we see that the two lowest energy levels become exponentially close in the large $k/D^{1/3}$ limit (cf. Fig. 7). The higher excited eigenvalues are well separated.

A more precise method is to study neatly the differential equation (115) : the first node in the solution $\varphi(z)$ appears when $\lambda\leqslant-\mathscr{E}_{1}(q)$. This allows to determine the non-analytic contribution to $\mathscr{E}_{1}(q)$ for relatively large $|E|^{3/2}/D$. Such a precise analysis of $\mathscr{E}_{1}(1)$ shows that the non-analytic contribution to $\mathscr{E}_{1}(1)$ is half of the correction to $\mathscr{E}_{0}(1)$ :

compare with Eq. (121). See Fig. 8 : the left plot shows the $1/|E|$ behaviour of the pre-exponential factor and the right plot exhibits the ratio $1/2$ in the large $|E|/D^{2/3}$ limit (we attribute the deviation of the last point to some numerical error). In Section 9, it will be shown that such non-analytic contribution to $\mathscr{E}_{1}(1)$ can indeed be obtained by a judicious extension of the classical WKB analysis. Furthermore, we will demonstrate that the dimensional factor is $1/(8\pi)\simeq 0.04$.

We can connect this spectral analysis with the study of equilibria for the elastic line : Fig. 7 shows that only the two lowest energies are positive and lead to exponentially growing contributions in Eq. (109) :

where $r_{1}=-\mathscr{E}_{1}(1)-\sqrt{|E|}<r$ and $c_{2}=\mathscr{E}_{2}(1)+\sqrt{|E|}>0$. Moreover, in the limit $|E|^{3/2}\gg D$ ($L_{m}\ll L_{c}$), Eqs. (121,122) show that $r_{1}\simeq r/2$.

In this section we study the eigenvalue $\mathscr{E}_{1}(q)$, which controls finite size effect (elastic line of finite length $L$), Eq. (123), and also provides a lower bound for the rate $r$. The calculation confirms the $D\to 0$ behaviour obtained numerically in the previous section, Eq. (122).

As we have seen in the article, the results seem to give a strong evidence in favour of very fast, hence non-perturbative vanishing of the $\Lambda(1)=-|E|^{1/2}\mathcal{E}_{0}(1)$ in this regime. Although $\Lambda(1)$ is given by an eigenvalue which does not belong to the spectrum of $\mathcal{H}_{q}$, we still can analyse semiclassically the ground state energy $\mathcal{E}_{1}(q)$ of the latter Hamiltonian, which provides a lower bound for $\Lambda(q)=-\mathscr{E}_{0}(q)>-\mathscr{E}_{1}(q)$ (moreover the numerics has showed that the non-analytic corrections to $\mathscr{E}_{0}(1)$ an $\mathscr{E}_{1}(1)$ only differ by a factor $1/2$ in the $s\to 0$ limit ; cf. Section D). The only way the energy $\mathscr{E}_{1}(q)$ of the ground state (located semiclassically in the right well around $\zeta\approx 1$) can get such a non-perturbative shift is by a tunnelling admixture from the lowest-level eigenstate located semiclassically in the higher (left) potential well.

We now develop an improved WKB procedure in order to analyse the ground state of the Schrödinger equation

as $s\to 0$. The potential $V_{q}(\zeta)$ has two deep minima around $\zeta=\pm 1$ separated by a big barrier around $\zeta=0$ (Fig. 9). The potential is well approximated by two quadratic wells near these two points

Before starting the presentation of the WKB method, it is useful to have in mind the correspondence with the standard formulae in quantum mechanics. This is done by identifying $s=\hbar^{2}/2m$ : we can set $\hbar=1$ and $m=1/(2s)$. We see that the classical oscillator frequency in each well in our problem is given by $\omega=2$. The spectra related to the two harmonic wells (127,128) are :

which shows that the two spectra are in correspondence for integer $q$, apart for the $q+1$ lowest levels (Fig. 9).

2 Weakly asymmetric double well (|q+1|≪1much-less-than𝑞11|q+1|\ll 1)

𝑞11|q+1|\ll 1) For $q=-1$ the double well potential is symmetric and it is possible to use known results to express the energy of the ground state . The analysis of the bottom of the spectrum can be mapped onto a simple two level problem

3 Strongly asymmetric double well (|q+1|≳1greater-than-or-equivalent-to𝑞11|q+1|\gtrsim 1)

𝑞11|q+1|\gtrsim 1) In the “strongly” asymmetric potential limit, i.e. for $|q+1|\gtrsim 1$, it is not anymore possible to restrict the problem to a two level problem. We have to develop a different strategy involving two different approximation schemes : in the neighbourhood of each potential well ($\zeta\sim\pm 1$), the potential is replaced by parabolas, which allows to express the exact solution of the approximated Schrödinger equation locally. In between, inside the potential barrier, we write the approximate WKB solution of the (exact) Schrödinger equation and match the three expressions. This strategy is borrowed from Ref. .

1\zeta\sim+1 In the neighbourhood of the right well, the Schrödinger equation (126) takes the form

It will be convenient to introduce the variable $\xi=(\zeta-1)/\sqrt{s}$ and parametrize the energy as

where $\epsilon$ is a small (negative) shift to the ground state energy of the right well. Extracting the Gaussian function $\Psi(\zeta)=y(\xi)\,\exp[-(1/2)\xi^{2}]$ we obtain that the function $y(\xi)$ obeys the Hermite equation $y^{\prime\prime}(\xi)-2\xi\,y^{\prime}(\xi)+2\epsilon\,y(\xi)=0$, whose solution can be expressed in terms of the Hermite function with integral representation (for $\epsilon<0$)

This integral representation is suitable to extract the asymptotic behaviour

and, splitting the integral in (134) as $\int_{0}^{\infty}=\int_{-\infty}^{+\infty}-\int_{-\infty}^{0}$ :

decays exponentially for $\xi=(\zeta-1)/\sqrt{s}\to+\infty$ and grows exponentially for $\xi=(\zeta-1)/\sqrt{s}\to-\infty$. Keeping the variable $\xi$, we write :

3.2 Step 2: wave function for ζ∼−1similar-to𝜁1\zeta\sim-1

In the neighbourhood of the left harmonic well, the Schrödinger equation (126) takes the form

which now decays in the negative direction and grows in the positive direction :

3.3 Step 3: WKB solution inside the barrier

Inside the potential barrier, the solution of (126) is well approximated by the WKB wave function

3.4 Step 4: matching

In order to match the WKB solution (142) with (137) and (140), it is convenient to introduce a specific notation for the action : we define

The total action associated to the trajectory going from $\zeta=-1$ to the turning point is denoted

thus $S_{q}=\Phi_{L}(\zeta)+\Phi_{R}(\zeta)$, obviously.

In the vicinity of the right harmonic well at $\zeta=+1$, the asymptotic form (138) for $\xi\to-\infty$ should match with the WKB solution (142), which provides a first constraint on the two coefficients $A$ and $B$. We must therefore consider $\zeta$ sufficiently far from the turning point, i.e. $(\zeta_{0}-\zeta)\gg\sqrt{s}$, so that the asymptotic behaviour (138) holds. $\zeta$ must however be sufficiently close to the turning point so that the parabolic approximation (for the potential) is still justified. This last observation allows us to simplify the action $\Phi_{R}(\zeta)$ : the parameter $t_{*}=(1-\zeta)/(1-\zeta_{0})\gg 1$ will be treated as a large parameter. Using $V_{q}(\zeta_{0})=\mathcal{E}$ we have $p(\zeta)=(1/s)\sqrt{(1-\zeta)^{2}-(1-\zeta_{0})^{2}}$. Introducing the variable $t=(1-\zeta^{\prime})/(1-\zeta_{0})$ we find in this approximation:

and finally using $t_{*}\sqrt{t_{*}^{2}-1}\approx t_{*}^{2}-1/2$ and $\mbox{arccosh}(t_{*})\approx\ln{(2t_{*})}$, we get

We now find convenient to express $\Phi_{R}(\zeta)$ in terms of the variable \xi=(\zeta-1)/\sqrt{s}=-\big{[}(1-\zeta_{0})/\sqrt{s}\big{]}t_{*} introduced above :

from which we write the WKB wave function as

where we made use that $\epsilon\ll 1$. Matching with (138) for $\xi\to-\infty$ gives

where we have used $\Gamma(-\epsilon)\simeq-1/\epsilon$. We get the first condition

We proceed the same way in the neighbourhood of the left harmonic well and match the asymptotic behaviour (141) for $\xi\to+\infty$ with the WKB approximation (142). In the region where the two expressions of the wave function match, we can use the parabolic approximation for the potential :

(we recall that we can neglect $\epsilon$ there). Introducing the variable $t=(\zeta^{\prime}+1)/\sqrt{s(2q+1)}$, this rewrites

(142) and (141) match in the regime where $t_{*}=(\zeta+1)/\sqrt{s(2q+1)}\gg 1$, thus, using $t_{*}\sqrt{1+t_{*}^{2}}\approx t_{*}^{2}+1/2$ and $\mbox{arcsinh}(t_{*})\approx\ln{(2t_{*})}$, we find

At this stage it is convenient to use the variable $\xi=(\zeta+1)/\sqrt{s}=\sqrt{(2q+1)}t_{*}$ introduced above :

3.5 Ground state energy

Comparing (152) and (159) finally provides the expression of the shift of the ground state energy

Making use of (144), we can write the action (146) as

In F, we show that it presents the limiting behaviour

We conclude that the ground state energy of the Schrödinger operator $\mathcal{H}_{q}$ is given by

which presents a different pre-exponential dependence, compared to (132) obtained for weakly asymmetric double-well. Going back to the initial notation

It is however important to remember that (164,165) are not the full result, but only the non-analytic contribution to the ground state energy (which is not expected to be the dominant correction). As discussed at length in the main text and in Section 8, the ground state energy is also shifted by analytic contributions (in $s$) related to the non-harmonicity of the potential. As explained above, the analytic terms are also given by (98) :

For $q=1$, we have demonstrated that the first analytic contributions vanish (up to $\mathcal{O}(s^{2})$) and observed numerically that this is true at all orders in $s$. Hence the final result for $q=1$ is expected to be

Note that not only the power law of the pre-exponential term perfectly agrees with the numerical result, cf. Eq. (122), but moreover the dimensionless factor $1/(8\pi)\simeq 0.0397887$ coincides with the value extracted numerically in Subsection 8.3.

The generalized Lyapunov exponent for even integer argument

We recall that the starting point is the Schrödinger equation

Let us start from the simplest non trivial case $q=2$ which requires to analyse $\left\langle y(\tau)^{2}\right\rangle$. Differentiating $y^{2}$ leads to consider the set of three coupled SDEs

which are interpreted in the Stratonovich convention, as usual in physics when Gaussian white noises arise in order to model regular physical noises . The SDEs can be rewritten in the Itô convention as

The largest eigenvalue of the matrix $M_{2}$ controls the exponential growth of $\left\langle y(\tau)^{2}\right\rangle$ and therefore coincides with $\widehat{\Lambda}(2)=\Lambda(2)$. The characteristic polynomial is $P_{2}(\lambda)=\det(M_{2}-\lambda\,\mathbf{1}_{3})=-\lambda^{3}-4E\,\lambda+4D$, whose appropriate root is (for $E=-k^{2}<0$)

We thus have recovered a result of . For $D=0$ we have $\Lambda(2)=2k$ as it should, and for $E=0$, one gets $\Lambda(2)=(4D)^{1/3}$. One can also check that the $s\to 0$ expansion coincides with (98) for $q=2$, which corresponds to the limit $E\to-\infty$. For large positive energy $E\gg D^{2/3}$, we have $\Lambda(q)\simeq(q+q^{2}/2)\,\gamma_{1}\simeq(q+q^{2}/2)\,D/(4E)$ (cf. Section 7.1) : we can check that the expansion of (175) for $E\to+\infty$ gives $\Lambda(2)\simeq D/E$, as it should. The GLE is plotted in Fig. 10 for $q=1$ and $q=2$ : in order to match the asymptotic behaviours, $\Lambda(q)\simeq q\sqrt{-E}$ for $E\to-\infty$, we choose to plot $\Lambda(q)/q$.

for $q=1$, the analytic part of $\Lambda(q)-kq$ vanishes and $\Lambda(1)$ is a purely non-analytic function of $s=D/|E|^{3/2}$ (Section 8).

For $q=2n$, the non-analytic contributions to $\Lambda(q)$ vanish and the GLE is analytic in $s=D/|E|^{3/2}$.

In particular, this shows that the replica trick must be used with caution as it is not possible to deduce $\langle|y(\tau)|^{q}\rangle$ for arbitrary real $q$ from $\langle y(\tau)^{2n}\rangle$.

Let us now discuss systematically the calculation of $\Lambda(2n)$. For $q=2n$, one has to consider the $q+1$ equations

for $m=0,\,1,\cdots,\,q$. The stochastic differential equation can be rewritten in the Itô convention :

The idea that the case of integer $q$ leads to consider a closed system of equations for correlators of $y$ and $y^{\prime}$ has been used earlier (the matrix $M_{q}$ was obtained in Refs. ).

2.2 Perturbative analysis (D/k3≪1)D/k^{3}\ll 1)

As a check, let us recover the perturbative result (98). For $E=-k^{2}$ and in the absence of disorder, the matrix $M_{q}$ has the spectrum

The largest eigenvalue $\widehat{\Lambda}(q)=+qk$ is associated with the right and left eigenvectors

The result is in perfect agreement with (98). We emphasize that only for even integer $q=2n$ does the systematic perturbative expansion provide the exact result, without any additional non-analytic contribution in $D/k^{3}$ : according to the notation of Section 7 we can write

The determination of $\widehat{\Lambda}(q)$ reduces to analyse the largest eigenvalue of a $(q+1)\times(q+1)$ matrix $M_{q}$, which is easy to implement. For $E=0$, we can find the analytic expression of $\widehat{\Lambda}(q)$ up to $q=6$ (see Table 1).

This simple method also allows us to study the large $q$ behaviour. We obtain numerically (Fig. 11)

Finally, we note that the main behaviour $\Lambda(q)\sim D^{1/3}q^{4/3}$ was obtained for integer $q$ in the conference’s proceedings . It is also in agreement with the numerical calculations of $L(q)=\Lambda(q)/q\sim q^{\alpha-1}$ of Zillmer & Pikovsky , who obtained the exponents $\alpha\simeq 1.28$ and $\alpha\simeq 1.38$ for two different values of the energy.

3 Large deviations for the wave function amplitude

As it is clear from its definition (75), the generalized Lyapunov exponent is the generating function of the cumulants of the logarithm of the wave function $\ln|y(x)|$ [more precisely, $y(x)$ is the solution of the Cauchy initial value problem (71,70)]. As such, the GLE can be related to the large deviation function $\Phi(u)$ controlling the distribution

$\Lambda(q)$ and $\Phi(u)$ are related by a Legendre transform

We can therefore relate the behaviour (183) for $q\to+\infty$ to the asymptotic behaviour of the large deviation function

Here we show how our method can be extended to study the mean number of equilibria at fixed values of the energy, defined, for our discrete model of an elastic line, as

It is in fact convenient to split the total energy (5) into an elastic part and a disorder part as follows

($\Delta$ is here the discrete Laplacian) and define the mean number of equilibria at fixed elastic energy ${H_{e}}$ and disorder energy ${H_{d}}$ as

The easiest observable to study is the Laplace transform

using that the total energy is $H=H_{e}+H_{d}$. Note that the elastic energy is always positive, while the disorder part can be of either sign (hence the $s_{d}$ dependence really involves a double-sided Laplace transform).

We now use that $V_{i}(u_{i})$ and $V_{i}^{\prime\prime}(u_{i})$ are correlated Gaussian variables, described by the covariance matrix

and independent from $V^{\prime}_{i}(u_{i})$. We also use, for each monomer, the following property of Gaussian variables (the ”complete the square” trick)

for any function $F(x)$, where the second average runs over the marginal distribution of $V_{i}^{\prime\prime}(u_{i})$. Equivalently one can write for any constant $V$

where $\chi$ is a Gaussian random variable of variance $R(0)$ and of mean $V/R(0)$, independent of $V^{\prime\prime}(u_{i})$. So imposing the value of $V_{i}(u_{i})$ amounts to shift the random potential by an independent Gaussian random variable. This leads to the following decoupling

As a result, the Laplace transforms $\widetilde{N}(s_{e},s_{d})$ and $\widetilde{N}(s,s)$, defined in (191) are expressed very simply in terms of the determinants studied in this paper

The formula (21) can be similarly written for the continuum model, for which the elastic and disorder energies are defined respectively as

The joint Laplace transform (defined in the same way as for the discrete model) now reads

where we recall that $\langle U(\tau)U(\tau^{\prime})\rangle=2D\,\delta(\tau-\tau^{\prime})$ with $D=R^{\prime\prime\prime\prime}(0)/(2\kappa^{2})=1/(2L_{c}^{3})$. The above product form shows the statistical independence of the elastic and the disorder energy at a force-free point.

Let us now study some implications of this formula. We first review some of our previous results needed here

The determinants in (204) not containing $U(\tau)$ can be read from Eqs. (46) for various boundary conditions, all having the same leading behaviour in the large $L$ limit $\det(\gamma-\partial_{\tau}^{2})\sim\exp(\sqrt{\gamma}L)$.

The following ratio of determinants was studied for large $L$

and its rate of growth $r(m^{2})$ as a function of $m^{2}$ was shown to be [see Eq. (8)]

where $C\simeq 0.46$ and we recall that $L_{m}=\sqrt{\kappa}/m$ and $L_{c}=\kappa^{2/3}R^{\prime\prime\prime\prime}(0)^{-1/3}$ [the first limiting behaviour was given in Eq. (120), recalling the correspondence of notations $E=-m^{2}/\kappa$ and $(2D)^{1/3}=1/L_{c}$].

In view of (i) and (ii) there are thus two main applications, studied respectively in the two subsections below :

2 Large deviations and rate functions in the large L𝐿L limit

We will denote $H_{e}=h_{e}L$, $H_{d}=h_{d}L$ and $H=hL$, where $h_{e}$, $h_{d}$ and $h$ are respectively the elastic, disorder and total energy densities. We expect in the large $L$ limit that

Hence, from the definition of the Laplace transforms (191)

On the other hand from (i) and (ii) above we can write, in the limit of large $L$

One can use a saddle point to estimate the leading large $L$ behaviour of (209) and (210), and we find that the rate functions $r(m^{2};h_{e},h_{d})$ and $r(m^{2};h)$ are related to $\Gamma(s_{e},s_{d})$ and $\Gamma(s,s)$ by Legendre transforms. More precisely one has

Below we will also study the rates associated to fixing the elastic energy alone, and the disorder energy alone, namely

Note the two important general observations, which will be confirmed below by explicit calculations:

The maximum over $h_{e}$ of $r_{e}(m^{2};h_{e})$ occurs at the field $h_{e}^{*}$ which corresponds to $s_{e}=0$ and its value is $\Gamma(0,0)=r(m^{2})$. It is easy to check that this property holds from the derivative conditions. The field $h_{e}^{*}$ then corresponds to the typical (i.e. the most probable) value of $h_{e}$. Same property holds, respectively, for the rate $r_{d}(m^{2};h_{d})$ (with typical value $h_{d}^{*}$ occurring at $s_{d}=0$).

The form obtained in (212) satisfies that

which reflects the independence of the elastic and potential energy noted above. Hence we have

To proceed further, it is convenient to introduce dimensionless variables. Three main dimensions are involved here : the energy $[\mathcal{H}]=E$, the length $[\tau]=L$ and the field $[u]$. We deduce the dimensions of the main quantities : the correlator $[R(u)]=E^{2}/L$, the mass $[m]=E^{1/2}L^{-1/2}[u]^{-1}$ and the elastic constant $[\kappa]=E\,L[u]^{-2}$. Because we will be interested in the $m\to 0$ limit, we choose to rescale all observables with respect to the length scale $L_{c}$. We define the rescaled energy density and its conjugate variable

We recall that $\widetilde{\Lambda}(\mu)$ is a monotonously increasing function of $\mu$, see Fig. 3 and Fig. 12, with the three limiting behaviours obtained in the previous sections

We now study respectively the rates at fixed elastic and disorder energy, and finally at fixed total energy.

As a warm up simple exercise, we consider first the number of equilibria constrained by the elastic energy

Here we only compute the large deviation function $r_{e}(m^{2};h_{e})$ (a more precise calculation of the distribution of elastic energy will be presented in Subsection 11.3).

which also corresponds with its annealed average value as indicated by $\langle\cdots\rangle_{a}$. From (234), we deduce the variance of the annealed distribution of the elastic energy density as

2.2 Mean number of equilibria constrained by the disorder energy

In a second stage, we consider the number of equilibria constrained by the disorder energy $H_{d}$ :

Although it is not possible to obtain a simple analytical form, as for the elastic energy, one can establish various general features and limiting behaviours.

The minimizer of Eq. (239) is given by $\sigma_{d}=\sigma_{*}$ solution of

As discussed above for the elastic energy, from general considerations the maximum of $r_{d}(m^{2};h_{d})$ over $h_{d}$ should equal $r(m^{2})$ and occur at the typical field $h_{d}^{*}$ corresponding to the argmin value $\sigma_{*}=0$. We thus immediately obtain the typical, most probable value of the (dimensionless) disorder energy density as

Because $\widetilde{\Lambda}^{\prime}(\mu)>0\ \forall\mu$ (cf. Fig. 12) we see that the typical disorder energy is negative

We recall that $\mu=(L_{c}/L_{m})^{2}\propto m^{2}$. In the limit of zero mass, $m\to 0$ (absence of external confinement), we can use the value obtained from the numerics $\widetilde{\Lambda}^{\prime}(0)=a_{1}\simeq 0.47$ (cf. Fig. 12) to get a good estimate for the typical disorder energy. In the other limit of strong confinement ($L_{m}\ll L_{c}$), using (229), we deduce that the typical disorder energy takes the form

i.e. is twice the elastic energy (236), with the opposite sign. Note that our numerics indicate that $\widetilde{\Lambda}^{\prime}(\mu)$ reaches its maximum at $\mu\approx 0.4$ (when the curvature $\widetilde{\Lambda}^{\prime\prime}(\mu)$ changes in sign). Intuition about ground states would suggest that the lower the mass, the larger the available space for the elastic line to explore better locations in the random potential, hence the lower the typical disorder energy. However here we are dealing with all stationary points, which seems, from our numerics, to behave differently (with the minimum $h_{d}^{*}$ occurring at a non zero value of $\mu$).

Our numerics also indicates that $\widetilde{\Lambda}^{\prime\prime}(\mu)$ is always above $\sim-0.17$ (Fig. 12). Thus the combination $\eta+\widetilde{\Lambda}^{\prime\prime}(\mu)>0$ always remain positive, as needed. This combination enters the variance of the disorder energy, from the annealed distribution, as one finds, from (249),

We compare in Fig. 13 these limiting behaviours with the result of a numerical resolution of Eqs. (241,242) : the agreement is excellent.

An interesting feature is that the decay of the rate at large positive energy is not controlled by the mass, as it was the case for the elastic energy, cf. Eq. (234): instead it has a parabolic shape controlled by $R(0)$ :

We have used the numerical data (Sections 8.1 and 8.2) to compute the rate $r_{d}(m^{2};h_{d})$ in order to check our analysis : the agreement is good, as one can see in Fig. 13.

2.3 Mean number of equilibria constrained by the total energy

We now study the rate (225) at fixed total energy $H$. We recall that the total energy can have both signs, contrary to the elastic energy $H_{e}$, which is positive. Here also we can obtain the asymptotic behaviours easily. We have now to consider

hence, not surprisingly, the typical total energy is the sum $h^{*}=h_{e}^{*}+h_{d}^{*}$ of the typical disorder and elastic energies obtained in the previous sections. In the original units, the typical total energy density reads

From (229) we see that the last factor $(1-4\sqrt{\mu}\,\widetilde{\Lambda}^{\prime}(\mu))$ changes sign as the mass increases and varies from $1$ (positive typical energy, elastic energy dominates) for $m\to 0$, corresponding to weak confinement, to $-1$ (negative typical energy, disorder energy dominates) for $m\to+\infty$, corresponding to strong confinement.

The expansion of (252) for $\sigma_{*}\to 0$ allows to study the typical fluctuations, in the same way as for the disorder energy. As expected, the variance is given by the sum of variances of elastic and disorder energy computed above :

We now derive limiting behaviours for the rate. Asymptotic analysis of (252) using (229) leads to (the analysis is quite similar to the two previous sections)

Correspondingly, using (253) we deduce the following behaviours for the large deviation rate function

Using the data of the numerical calculation (Sections 8.1 and 8.2), we have solved numerically (252) and computed the rate (253) : the result is plotted in Fig. 14. We see that the agreement with the limiting behaviours discussed in the text is excellent.

Let us discuss some main features of this result. In the limit of large total negative energy density we thus obtain the dominant term [corresponding to $-(\eta/2)\sigma_{*}^{2}$ in Eq. (253)]

which is identical to the leading behaviour for large negative disorder energy from (250). Since the elastic energy is positive, it is reasonable to expect that the two tails coincide to leading order. Notice by comparing (258) and (249) that they differ however by the prefactor of the next (subleading) order.

Consider now the limit of large positive total energy density. From (258) we obtain

This is the same behaviour as (234), obtained for the elastic energy. Again it is quite reasonable that these two tails coincide, since the disorder energy rate function was found to decay much faster at large positive disorder energy density than the elastic energy one. Hence the elastic energy dominates in that regime. Note however that the next (subleading) order term, which is $\mathcal{O}(1)$, is different for the total and the elastic energy rate functions.

As was already discussed in the context of the elastic energy, the decay of the large deviation function with $h$ is controlled by the mass : when the mass vanishes, $r(0;h)$ is a monotonous function which saturates

3 Mean number of equilibria at fixed value of the elastic energy for any L𝐿L

Let us study the mean number of equilibria at fixed value of the elastic energy $H_{e}$. As discussed in C, the annealed distribution of elastic energy $H_{e}$ for our model, exactly coincides with the one of the Larkin model. The present calculation, however, which treats carefully the boundary conditions, has not, to our knowledge, been reported previously, even in the context of the Larkin model.

Note that since the elastic energy is positive, these are bona-fide Laplace transforms. We can now use the expressions of these determinants, Eq. (46).

Let us focus on periodic boundary conditions, which leads to the simplest formula, and indicate the results for others. Then we have, with $L_{m}=\sqrt{\kappa}/m$

The natural unit of ${H_{e}}$ is ${H_{e}}_{m}={|R^{\prime\prime}(0)|}/{m^{2}}$. The quantity of most interest is the annealed (i.e. over samples) probability distribution, $\mathcal{P}_{L,m}(H_{e})$, that an equilibrium chosen at random has elastic energy $H_{e}$, takes the form (for all classes of boundary conditions)

For the periodic case, we have to compute the inverse Laplace transform

Let us also indicate the result for Dir/Dir

The above formulae can be inverted in principle to obtain $\mathcal{P}_{L,m}(H_{e})$ for any $L$. For periodic boundary condition, using

From this we conclude that for $L/L_{m}\gg 1$ the leading (i.e. typical) fluctuations of the elastic energy are Gaussian

where $\chi$ is a Gaussian random variable of unit variance. The form (271) is furthermore consistent with the large deviation rate function (in the notations of the previous section) associated to $H_{e}=h_{e}L$, that is

which indeed coincides with the Legendre transform w.r.t $s_{e}$ of $\Gamma(s_{e},s_{d}=0)$ in Eq. (212), cf. Eq. (234).

3.2 Case m=0𝑚0m=0

In this final paragraph, we consider the distribution of the elastic energy in the zero mass limit. One of our aim is here to clarify the origin of the saturation of the rate $r(m^{2};h)$, Eq. (261), or equivalently of the rate $r_{e}(m^{2};h_{e})$, as the tail of the distribution of the total energy is controlled by the distribution of the elastic energy. We consider the case of Neumann/Dirichlet boundary conditions, which describes the case where the line is attached by one end, the other end being free. We choose these boundary conditions for the simplicity of the determinant, Eq. (46). From (265), we see that the distribution then takes the form

This expression shows that the typical elastic energy density grows with the length

in terms of the dimensionless variable introduced above.

In G, we show that the limiting behaviours of the inverse Laplace transform (274) are

If the distribution is rewritten under the form \mathcal{P}_{L,0}(H_{e})\sim\exp\{L\,\big{[}r_{e}(m^{2};h_{e})-r(m^{2})\big{]}\} for $L\to\infty$, we have

The first behaviour is in exact correspondence with Eq. (261) or Eq. (234) for $m=0$. The second behaviour ensures the decay of $\mathcal{P}_{L,0}(H_{e})$ for large $H_{e}$, as it should. This makes clear that the saturation of the rate Eq. (261) only reflects a subtlety concerning the order of the two limits $L\to\infty$ and $h_{e}\to\infty$ when $m=0$.

Conclusion

By extending the Kac-Rice approach to manifolds of finite internal dimension, we have shown that the mean number of equilibria of an elastic line in a random potential in presence of a parabolic confinement, grows exponentially with its length, and developed a theory to calculate the rate of growth $r$. The mean number equilibria is shown to be related to the expectation value of the modulus of the determinant of the Laplace operator in presence of diagonal disorder, the same operator which features in the Anderson localization problem. The relation is valid for elastic interface of arbitrary internal dimension $d$. In $d=1$ (elastic line) using the Gelfand-Yaglom theorem, the rate of growth can be related with the large deviation function of the Lyapunov exponent fluctuations associated to a 1D random Schrödinger problem, but in the negative energy range. This large deviation function (generalized Lyapunov exponent GLE) is given by the lowest eigenvalue of an associated Fokker-Planck operator which is analysed in detail by several complementary techniques.

From these methods, we found that the rate $r$ is described by a universal function of the disorder strength for which we obtained analytical and numerical results. The disorder strength is naturally parameterized using the Larkin length $L_{c}$ and the dimensionless control parameter is the ratio of $L_{c}$ to the length $L_{m}$ imposed by the parabolic confining potential. We extract analytically the asymptotic behaviours of this scaling function for small and large value of the argument. For strong confinement, the rate $r$ is small and given by a non-perturbative (instanton, Lifshitz tail-like) contribution to GLE. For weak confinement, the rate $r$ is found to be proportional to the inverse Larkin length of the pinning theory. We have also discussed the question of counting of stable equilibria. Finally, we have shown how to extend the method to calculate the asymptotic number of equilibria at fixed energy : we have obtained large deviation rate functions controlling the number of equilibria constrained by either the total, the elastic or the disorder energy. These large deviation functions control the distribution of the (total, elastic or disorder) energy density of the line at equilibrium. Some connections with the Larkin model have been discussed.

From the point of view of Anderson localization, we have also discovered several interesting properties of the generalized Lyapunov exponent (GLE) $\Lambda(q)$. We have shown a crucial difference between the GLE for $q=1$, which is a non-analytic function of the disorder strength $D$ (for $E\to-\infty$) while $\Lambda(q)$ is analytic when $q$ is an even integer. It would be interesting to see if the difference persists for larger odd integer arguments. We were able to obtain few exact results for even integer $q$. The limit of large positive argument $q\to+\infty$ was analysed, which is related to the large deviations of the wave function (solution of the Cauchy initial value problem) for large values. The analysis of the limit $q\to-\infty$ of the GLE has remained an open question, which can be related to the large deviations of the wave function for small values.

The connection discussed here with the generalized Lyapunov exponents of a localization problem opens a bridge between pinning theory and localization theory that should inspire further works.

Acknowledgements

CT acknowledges stimulating discussions and suggestions from Philippe Bougerol, Alain Comtet, Aurélien Grabsch, Jean-Marc Luck, Nicolas Pavloff, Yves Tourigny and Denis Ullmo. We thank David Saykin for sharing his results and writing E. We thank the PCMI Summer School 2017, where some of this work was performed. The research at King’s College London was supported by EPSRC grant EP/N009436/1 The many faces of random characteristic polynomials. This research was also supported by ANR grant ANR-17-CE30-0027-01 RaMaTraF. We are grateful to an anonymous referee for numerous useful remarks.

Appendix A Boundary conditions and determinants

In this section, we discuss several formulae for the determinant $\det(H-E)$ of the Schrödinger operator to justify and make more precise the formula of the main text and of B below.

Consider the discrete Hamiltonian which appears in the text, $H_{i,j}=-\Delta_{i,j}+U_{i}\delta_{i,j}$ and the associated Schrödinger equation

for $i\in\{1,\cdots,K\}$ ; boundary conditions set the values for $\psi_{0}$ and $\psi_{K+1}$.

There are three natural boundary conditions for the elastic line problem studied in this paper. We detail them here in the discrete and continuum settings.

In the continuum limit it corresponds to the usual Laplacian with Dirichlet boundary conditions $\psi(0)=\psi(L)=0$.

A.1.2 Neumann boundary conditions

In the continuum limit it corresponds to usual Laplacian with Neumann boundary conditions $\psi^{\prime}(0)=\psi^{\prime}(L)=0$.

A.1.3 Periodic boundary conditions

In the continuum limit this corresponds to usual Laplacian for $\psi(0)=\psi(L)$ and $\psi^{\prime}(0)=\psi^{\prime}(L)$.

A.2 Determinants

We now discuss separately the formulae for the determinant $\det(H-E)$ for the three types of boundary conditions.

We denote by $y_{i}(E)$ the solution of the initial value problem with

The $K$ eigenvalues $\{E_{\alpha}\}_{\alpha=1,\cdots,K}$ of the $K\times K$ Hamiltonian $H$ are the solutions of the quantization condition

We see by recursion that $y_{2}(E)=2+U_{1}-E$, $y_{3}(E)=(2+U_{2}-E)(2+U_{1}-E)-1$, etc. By induction, it is straightforward to prove that $y_{j+1}(E)$ is a polynomial of degree $j$ in $E$ with higher degree term $(-E)^{j}$. Using these remarks we can write

This can now be used for the discrete polymer model, formula (21) in the text, with the correspondence $E=-m^{2}/\kappa$. As a simple illustration consider the free case with $U_{n}=-2$. We deduce $y_{n}=\sin qn/\sin q$ where $E=-2\cos q$. As a consequence the Dirichlet determinant is $\det(H-E)=\sin q(K+1)/\sin q$, corresponding to the spectrum $q_{n}=n\pi/(K+1)$ with $n=1,\cdots,K$.

Denoting $y(\tau)$ the solution of $(H-E)y(\tau)=0$ with initial conditions

we see that Eq. (284) is obviously the discrete version of the Gelfand-Yaglom formula

up to a $E$ independent multiplicative factor which depends on the ultraviolet regularization (the factor chosen here and formulae given below correspond to zeta regularization ). Let us illustrate the formula in the free case where $H=-\partial_{\tau}^{2}$ : setting $\gamma=-E$ for convenience for the following, we find $y(\tau)=\sinh(\sqrt{\gamma}\tau)/\sqrt{\gamma}$ thus $\det(\gamma-\partial_{\tau}^{2})=2\sinh(\sqrt{\gamma}L)/\sqrt{\gamma}$. Rewriting the hyperbolic function as an infinite product

We recognize the eigenvalues of the Laplacian. In this case we recall the eigenfunctions, $-\partial_{\tau}^{2}\psi_{n}(\tau)=q_{n}^{2}\psi_{n}(\tau)$,

A.2.2 Free boundary conditions (Neumann)

$\psi_{0}=\psi_{1}$ and $\psi_{K}=\psi_{K+1}$. We now denote by $x_{i}(E)$ the solution of the initial value problem with

A.2.3 Periodic boundary conditions

For the sake of completeness, let us add a phase and now write the periodic boundary conditions as

We recall also the related eigenfunctions for completeness

For references on functional determinants, cf. Refs. and the review .

Appendix B Pinning of the line for fixed boundary conditions

We provide some additional information for the analysis of § 4.2.3 : we show how the discussion can be extended to the case of a line with fixed endpoints. We can also perform the calculation for an elastic line pinned at its two ends (Dirichlet boundary conditions). although it is not the standard one for the depinning problem, where the elastic line is usually free to move, it is a usual setting in the related sandpile problem . The Fourier coefficients are

In the limit $L_{m}/L\to\infty$ (vanishing mass) and $L/L_{w}\to\infty$, we find

Appendix C The Larkin model

The Larkin model is a simplified model for pinning. Let us recall it here for an elastic line in the continuum (discrete versions, and extensions to higher $d$ are immediate, for reviews see Refs. ). In the Larkin model the non-linear random potential term in Eq. (1) is replaced simply by a linear random force term

which is centered Gaussian with correlator $\left\langle\phi(\tau)\phi(\tau^{\prime})\right\rangle=-R^{\prime\prime}(0)\,\delta(\tau-\tau^{\prime})$. It amounts to expand $V(u,\tau)=V(0,\tau)+V^{\prime}(0,\tau)\,u+\cdots$ and discard higher order terms. Being quadratic, there is a unique energy minimum (i.e. a single equilibrium)

which is distributed as a Gaussian field with correlator

where in the second equation we assumed $L$ large enough to ignore boundary conditions: this integral behaves as $L_{m}^{2\zeta_{L}}$ where $\zeta_{L}=3/2$ is the Larkin roughness exponent.

It is shown that at zero temperature the Larkin model, a zero dimensional version of the elastic line model (1), has the same correlation functions for the field $u(\tau)$, to all orders in perturbation theory in the disorder, as the original model (1) (the so-called dimensional reduction phenomenon). However it lacks the (non-perturbative) feature of multiple equilibria. Indeed, one notes that it depends on a single parameter, $R^{\prime\prime}(0)$, and lacks the physics of pinning arising from a non vanishing $R^{\prime\prime\prime\prime}(0)$. However, it does provide a good approximate description of the correlation functions of the true model (in its ground state) at scales smaller than $L_{c}$, or for $m\gg m_{c}$ (such that $L_{m_{c}}=L_{c}$), although even in this regime it lacks the non-perturbative corrections originating from rare events, such as the one obtained in this paper.

The elastic energy at the minimum is easily obtained as

and its Laplace transform is simply obtained by integration over the Gaussian field $\phi(\tau)$ as

which is exactly the same factor which occurs in (265). Thus the prediction of the Larkin model (which has a single equilibrium) coincides exactly with the annealed probability distribution of the elastic energy of the true model over all equilibria. As an example, by averaging (311) over $\phi(\tau)$ one recovers its average (or most probable) value as

In the Larkin model as written in (311) the disorder energy is simply $H_{d}=-2H_{e}$ and total energy $H=-H_{e}$, so $H_{d}$ and $H_{e}$ are not independent as in the true model (in Section 11.2, the relation $H_{d}=-2H_{e}$ was shown to hold only for the typical values in the strong confinement regime, $L_{m}\ll L_{c}$). Trying to improve on the model (as adding the $V(0,\tau)$ term from the expansion of the potential) does not lead to a consistent description of $H_{d}$. That no such simple approximation exist for $H_{d}$ and $H$ is corroborated by the results of Section 11.2.2 where the typical disorder energy density is obtained in terms of the (quite non-trivial) GLE of a 1D Anderson localization problem.

We make several remarks on the spectrum of the forward generator

when the potential is such that the diffusion is characterized by a non-equilibrium stationary state (NESS) on the full real axis. This situation requires that $\mathscr{U}(z)\to\pm\infty$ for $z\to\pm\infty$, and the drift $|\mathscr{U}^{\prime}(z)|$ grows sufficiently fast at infinity so that the particle is driven in a finite time from $+\infty$ to $-\infty$ :