On the finite-size Lyapunov exponent for the Schroedinger operator with skew-shift potential

TL;DR

This paper investigates the finite-size Lyapunov exponent for a quantum system with a skew-shift potential, providing numerical verification of conditions related to localization phenomena in such systems.

Contribution

It offers the first numerical verification of finite-size Lyapunov exponent conditions for skew-shift potentials, advancing understanding of localization in these models.

Findings

Numerical verification of finite-size Lyapunov exponent conditions.

Support for localization conjectures in skew-shift potential models.

Extension of finite-size criteria to non-golden mean frequencies.

Abstract

It is known that a one-dimensional quantum particle is localized when subjected to an arbitrarily weak random potential. It is conjectured that localization also occurs for an arbitrarily weak potential generated from the nonlinear skew-shift dynamics: $v_n=2\cos\left(\binom{n}{2}\omega +ny+x\right)$ with $\omega$ an irrational number. Recently, Han, Schlag, and the second author derived a finite-size criterion in the case when $\omega$ is the golden mean, which allows to derive the positivity of the infinite-volume Lyapunov exponent from three conditions imposed at a fixed, finite scale. Here we numerically verify the two conditions among these that are amenable to computer calculations.

Tables

Table 1. Table 1. Finite-size Lyapunov exponents calculated by numerical integration over a regular 2001 × 2001 2001 2001 2001\times 2001 grid of ( x , y ) ∈ [ 0 , 1 ] 2 𝑥 𝑦 superscript 0 1 2 (x,y)\in[0,1]^{2} .

	$E=E_0$	$E=E_1$
$L_{N_0}(1/2,E)$	0.08071	0.46561
$L_{2N_0}(1/2,E)$	0.08070	0.46559
$\frac{L_{N_0}(1/2,E)-L_{2N_0}(1/2,E)}{L_{N_0}(1/2,E)}$

Equations

$${\psi }_{n+1}+{\psi }_{n-1}+\lambda v_n{\psi }_n=E{\psi }_n,$$

$$(H\psi {)}_n={\psi }_{n+1}+{\psi }_{n-1}+\lambda v_n{\psi }_n.$$

$$v_n=2\cos \left(\left(\genfrac{}{}{0ex}{}{n}{2}\right)\omega +ny+x\right)$$

$$\displaystyle\left(\begin{array}[]{c}\psi_{n+1}\\ \psi_{n}\end{array}\right)=$$

$$\text{where}{M}_n(x,y;\lambda ,E):=$$

$$\text{and}{A}_j(x,y;\lambda ,E):=$$

$$L(\lambda ,E):=\underset{n\to \infty }{\lim }\frac{1}{n}{\int }_{{\mathbb{T}}^2}\log \Vert M_n(x,y;\lambda ,E)\Vert \mathrm{d}x\mathrm{d}y,$$

$$\frac{1}{n}\log \Vert M_n(x,y;\lambda ,E)\Vert \stackrel{n\to \infty }{⟶}L(\lambda ,E)$$

$$L(\lambda ,E)>0.$$

$$L_n(\lambda ,E):=\frac{1}{n}{\int }_{{\mathbb{T}}^2}\log \Vert M_n(x,y;\lambda ,E)\Vert \mathrm{d}x\mathrm{d}y$$

$$u_n(x,y;\lambda ,E):=\frac{1}{n}\log \Vert M_n(x,y;\lambda ,E)\Vert .$$

$${\mathcal{B}}_n(\lambda ,E):=\left\{(x,y)\in {\mathbb{T}}^2:\left|u_n(x,y;\lambda ,E)-L_n(\lambda ,E)\right|>\frac{L_n(\lambda ,E)}{10}\right\}.$$

$$L(\lambda ,E)\ge \frac{1}{2}{L}_{N_0}(\lambda ,E)\ge 10^{-3}.$$

$$E_0=0,\text{ and }{E}_1=-\mathrm{2.49512326.}$$

$$H_{N}(x,y;\lambda)=\left(\begin{array}[]{cccccc}\lambda v_{1}(x,y)&1&0&0&\ldots&0\\ 1&\lambda v_{2}(x,y)&1&0&\ldots&0\\ 0&1&\lambda v_{3}(x,y)&1&\ldots&0\\ 0&0&1&\lambda v_{4}(x,y)&\ddots&0\\ \vdots&\vdots&\vdots&\ddots&\ddots&1\\ 0&\ldots&\ldots&0&1&\lambda v_{N}(x,y)\\ \end{array}\right)$$

$$\Gamma <5.708\times 10^{-4}.$$

$$\mathrm{dist}(E_1,\mathrm{spec}H)<10^{-8}.$$

$$\underset{E\in \{E_0,E_1\}}{\min }{L}_{N_0}(\lambda ,E)\ge$$

$$\underset{E\in \{E_0,E_1\}}{\max }\frac{L_{N_0}(\lambda ,E)-L_{2N_0}(\lambda ,E)}{L_{N_0}(\lambda ,E)}\le$$

$$Q_N(P)=\frac{1}{P}\sum _{i=1}^P{u}_N(\lambda ,E;x_i,y_i)$$

$$\mathbb{P}(|Q_N(P)-L_N|>ϵ\times 10^{-2})\le \frac{1}{300{ϵ}^2}.$$

$$\mathrm{Var}\left(\frac{1}{N}\log\|M_{N}(\lambda,E;x,y)\|\right)\leq\max_{j}\log\left\|\left(\begin{array}[]{cc}E-\lambda v_{j}&-1\\ 1&0\end{array}\right)\right\|.$$

$$\left\|\left(\begin{array}[]{cc}a&-1\\ 1&0\end{array}\right)\right\|=\frac{\sqrt{2+a^{2}+|a|\sqrt{4+a^{2}}}}{\sqrt{2}},$$

$$\mathrm{Var}\left(\frac{1}{N}\log \Vert M_N(\lambda ,E;x,y)\Vert \right)\le \log \frac{\sqrt{2+(|E_1|+1{)}^2+(|E_1|+1)\sqrt{4+(|E_1|+1{)}^2}}}{\sqrt{2}}.$$

$$\mathrm{Var}(Q_N(P)-L_N)=\frac{\mathrm{Var}(\frac{1}{N}\log \Vert M_N(\lambda ,E;x,y)\Vert )}{P}<\frac{4/3}{P}<\frac{1}{3}\times 10^{-6}.$$

Full text

1. Introduction

A one-dimensional quantum particle living on $\mathbb{Z}$ with energy $E\in\mathbb{R}$ is described by the discrete Schrödinger equation

[TABLE]

where $\psi=(\psi_{n})_{n\in\mathbb{Z}}$ is a sequence in $\ell^{2}(\mathbb{Z};\mathbb{C})$. The real-valued potential sequence $v=(v_{n})_{n\in\mathbb{Z}}$ represents the environment that the particle is subjected to. (The “coupling constant” $\lambda>0$ is factored out for convenience.) As first famously realized by Anderson in 1958 [2], the decoherence introduced by a random (meaning independent and identically distributed) sequence of potentials can drastically affect the spectral and dynamical properties of the quantum particle. Physically, one observes a sudden onset of insulating behavior in the presence of a random environment (“Anderson localization”). Mathematically, it is known that for arbitrarily small $\lambda>0$, the one-dimensional Schrödinger operator

[TABLE]

has pure point spectrum with exponentially decaying eigenfunctions [6, 9, 16].

A natural follow-up question is then: How random does the environment have to localize the quantum particle? Alternative “quasi-random” environments are generated by sampling a nice function along the orbit of an ergodic dynamical system. This question is interesting from a purely mathematical ergodic theory perspective, but it also has practical implications, since computer simulations are mostly based on appropriate pseudo-random number sequences. For example, one can consider $v_{n}=2\cos(n\alpha+\theta)$ generated from sampling cosine along an irrational circle rotation; this is the well-known Harper (or Almost-Mathieu) model . It turns out that these linear underlying dynamics only produce localization for sufficiently strong potentials, namely only for $\lambda>1$ [15]. One would thus like to consider dynamics which are slightly more quasi-random than the shift.

A standing conjecture in this direction concerns the case when the potential is generated from the nonlinear skew-shift dynamics $T:\mathbb{T}^{2}\to\mathbb{T}^{2}$, $T(x,y)=(x+y,y+\omega)$, namely it is of the form

[TABLE]

with $\omega$ irrational (say Diophantine). The key difference between (3) compared to $v_{n}=2\cos(n\alpha+\theta)$ is the appearance of the nonlinear quadratic term $n^{2}\omega$. The conjecture states that the associated Schrödinger operator $H$ defined by (2) is Anderson localized for arbitrarily small $\lambda>0$ everywhere in the spectrum. Partial results in this vein are due to Bourgain [3] and Bourgain-Goldstein-Schlag [5]. Note that the conjecture says in particular that the skew-shift dynamics is appreciably more random-like than the circle rotation where $v_{n}=2\cos(n\alpha+\theta)$. (Recall that the latter is only localized for $\lambda>1$.) The observation that the skew-shift is more quasi-random than the shift has been made in another context by Rudnick-Sarnak-Zaharescu [21] and others [13, 19, 20] (concerning the spacing distribution) and also recently in [1] (concerning eigenvalues of large Hermitian matrices).

A crucial ingredient for localization on which we will focus is the positivity of the Lyapunov exponent of the associated cocycle, which is defined as follows. From now on, $v_{n}=v_{n}(x,y)$ is given by (3) with $\omega$ irrational. The second-order difference equation (1) can be solved by using transfer matrices:

[TABLE]

The Lyapunov exponent is defined as

[TABLE]

where the limit exists by subadditivity (Fekete’s lemma). We remark that the Furstenberg-Kesten theorem [23] implies that

[TABLE]

for Lebesgue-almost every initial condition of the skew-shift $(x,y)\in\mathbb{T}^{2}$, as long as $\omega$ is irrational. This leads us to the following relaxed version of the conjecture from above.

Conjecture 1**.**

Let $\omega$ be irrational. For every $\lambda>0$ and every $E\in\mathbb{R}$, it holds that

[TABLE]

Some limited progress on Conjecture 1 was made in [3, 5, 11, 12, 17, 18]. We collect some remarks concerning this conjecture.

Remark 2*.*

1. (i)

The lower bound on the Lyapunov exponent is nontrivial only when $E$ lies in the spectrum of the Schrödinger operator $H$ defined by (2). A concurrent conjecture, to which we will return later, says that the skew-shift model has no gaps in the spectrum. 2. (ii)

In the case of i.i.d. random $v=(v_{n})_{n\in\mathbb{Z}}$ it is known that the Lyapunov exponent is positive by Furstenberg’s theorem [8]. 3. (iii)

Herman’s subharmonicity trick [14] implies $L(\lambda,E)\geq\log\lambda$ and thus establishes a lower bound on the Lyapunov exponent if and only if $\lambda>1$. The interesting regime for Conjecture 1 is therefore $\lambda\in(0,1]$.

2. A finite-size criterion

Our paper is a follow-up to an approach on Conjecture 1 which was recently initiated in [11]. That paper was based on the methods developed in [5], specifically large deviation estimates for Lyapunov exponents and an inductive multi-scale machine based on the Avalanche Principle, and [11] rendered these methods effective, i.e., explicit constants were obtained for every step of the argument. As a result, one obtains finite-size criteria for the validity of Conjecture 1. Namely, deriving the conjecture for a fixed choice of parameters $\lambda,E$ and $\omega=\frac{\sqrt{5}-1}{2}$ (the golden mean) is reduced to verifying $3$ numerical conditions on the Lyapunov exponent at a fixed initial scale, called $N_{0}$ below. If true, these conditions can be fed into the effective inductive machine from [11] to obtain $L(\lambda,E)>0$.

For definiteness, we focus on the following finite-size criterion obtained in [11] (see Theorem 1.4 there). We define the Lyapunov exponent at scale $n\geq 1$ by

[TABLE]

and its non-averaged analog by

[TABLE]

We also define the bad set $\mathcal{B}_{n}$, where $u_{n}$ deviates by more than $10\%$ from its average:

[TABLE]

Theorem 3** ([11]).**

Let $\omega=\frac{\sqrt{5}-1}{2}$ and let $\lambda\in[1/2,1]$. Let $N_{0}=30,000$. Assume that for some energy $E\in[-2-2\lambda,2+2\lambda]$ the following hold:

1. (i)

$L_{N_{0}}(\lambda,E)\geq 2\times 10^{-3}$, 2. (ii)

$\frac{L_{N_{0}}(\lambda,E)-L_{2N_{0}}(\lambda,E)}{L_{N_{0}}(\lambda,E)}\leq\frac{1}{8}$, 3. (iii)

$\max(|\mathcal{B}_{N_{0}}(\lambda,E)|,|\mathcal{B}_{2N_{0}}(\lambda,E)|)\leq N_{0}^{-165}$.

Then:

[TABLE]

We see that Theorem 3 reduces the proof of Conjecture 1 for specific parameters $\lambda,\omega,E$ to a numerical calculation at the initial scales $N_{0}=30,000$ and $2N_{0}$. We remark that the methods in [11] are flexible and can be used to obtain variant finite-size criteria; some further examples are stated in [11]. We work with Theorem 3 because the initial scale $N_{0}=30,000$, while large, is amenable to numerical verification which is our goal here.

Our main contribution is to conduct a detailed numerical study of conditions (i)-(iii). The results are obtained by running MATLAB code on Harvard’s research computing cluster Odyssey. We mention that numerics for the skew-shift are a delicate matter because of the highly oscillatory matrix elements of $A_{j}$ defined in (4), particularly the quadratic $n^{2}\omega$ term in (3) combined with the irrational nature of $\omega$. These difficulties were partially overcome in an earlier numerical study [7] which went up to the scale $N=200$.

The remainder of the paper is organized as follows. In Section 3.1, we conduct some preliminary investigations into the spectrum of the Schrödinger operator by studying its finite-size approximations and employing some standard mathematical bounds controlling the error of these approximations. In Section 3.2, we verify conditions (i) and (ii) of Theorem 3 numerically for specific parameter choices—see Table 1 below. In Section 3.3, we investigate condition (iii) of Theorem 3. Given the large power $165$, condition (iii) cannot reasonably be verified numerically. Nonetheless, we include some graphs that hopefully shed some light on this condition.

We recall Remark 1.2 in [11] about condition (iii): Condition (iii) is a large-deviation estimate at the initial scales $N_{0}$ and $2N_{0}$ and hence plays a different role from conditions (i) and (ii). For instance, analogous large deviations estimates are analytically known to hold a priori, independently of the positivity of the Lyapunov exponent, for the Harper model and others [4, 5, 10]. At any rate, an analytical proof of condition (iii) is warranted, and the importance of this problem is elevated further in light of our numerical verification of conditions (i) and (ii) in the present work.

From now on, we fix the parameters $\omega=\frac{\sqrt{5}-1}{2}$ and $\lambda=1/2$.

3. Numerical results

3.1. Preliminaries on the spectrum of $H$

As mentioned after Conjecture 1, the only energies $E$ for which the positivity of the Lyapunov exponent $L(\lambda,E)$ is nontrivial are those lying in the spectrum of the infinite-volume Schrödinger operator $H(x,y)$ defined in (2) with $v_{n}=v_{n}(x,y)$ given by (3). General facts from ergodic theory and spectral theory imply that $\mathrm{spec}\,H(x,y)\subset[-2-2\lambda,2+2\lambda]$ and it does is independent of $(x,y)$ if a zero-measure subset of $\mathbb{T}^{2}$ is ignored.

Apart from these general facts, however, not much is known about $\mathrm{spec}\,H(x,y)$. Hence, it is a priori unclear on which energies $E$ we should focus our numerical investigation. Before we explain how we choose which energies to study (given our choice of other parameters $\omega=\frac{\sqrt{5}-1}{2}$ and $\lambda=1/2$), we note that this issue is not as concerning as it may seem at first sight, in light of another standing conjecture that the skew-shift Schrödinger operator $H$ given by (2) has no gaps in the spectrum. This concurrent conjecture is also supported by numerics [7]; see (8) below for a numerical upper bound on the spectral gap.

We choose to focus our investigations on the energies

[TABLE]

The energy $E_{0}=0$ is a natural choice by symmetry considerations. Indeed, the spectrum of $H$ is symmetric under reflection at [math]; note also the reflection symmetry in Figure 1 for the finite-size analog of this fact.

The energy $E_{1}$ is chosen as follows. While $\mathrm{spec}\,H$ is not directly accessible by numerics, we can consider the finite-size approximation to $H$, the $N\times N$ Hamiltonian matrix $H_{N}(x,y;\lambda)$ defined by

[TABLE]

Standard spectral estimates can be used to approximate $\mathrm{spec}\,H$ by $\mathrm{spec}\,H_{N}$, see for instance Corollary 2 in [7]. In [7], these estimates are combined with the eigenvalue and eigenvector results for $H_{100}$ to derive the following upper bound on the largest spectral gap in $\mathrm{spec}\,H$:

[TABLE]

We set $N=N_{0}=30,000$ and diagonalize $H_{N}(x,y;1/2)$ over a $200\times 200$ regular grid of pairs $(x,y)\in[0,1]^{2}$. The resulting empirical spectral measure is shown as a histogram in Figure 1. We choose $E_{1}$ as the location of the bin with the tallest peak (using higher precision than the one used to produce the histogram). In other words, $E_{1}$ is close to an eigenvalue of $H_{N_{0}}$ for the maximal number of pairs $(x,y)\in[0,1]^{2}$ and in this sense is the number of maximal likelihood on scale $N_{0}$ to lie in $\mathrm{spec}\,H$. (We remark that we determine $E_{1}$ to 16-digit accuracy as $E_{1}=-2.495123260049612$ and conduct all computations at this accuracy. The approximation (7) is chosen in the main text for readabiity.)

Moreover, the spectral estimate from Corollary 2 in [7] can be used to show that $E_{1}$ lies inside the spectrum to that level of accuracy, namely,

[TABLE]

This estimate is relevant insofar as we want to check $L(\lambda,E)$ for $E\in\mathrm{spec}\,H$ and by the above considerations $E_{1}$ is our best guess for an element of $\mathrm{spec}\,H$ from spectral information at scale $N_{0}$.

3.2. Numerical verification of conditions (i) and (ii)

We recall that $\omega=\frac{\sqrt{5}-1}{2}$, $\lambda=1/2$, and $N_{0}=30,000$. Our main results are the numerical approximations of the finite-size Lyapunov exponents $L_{N_{0}}(1/2,E)$ and $L_{2N_{0}}(1/2,E)$ for $E\in\{E_{0},E_{1}\}$ given in (7). The results are summarized in Table 1. From these results, we see that

[TABLE]

Hence, conditions (i) and (ii) of Theorem 3 hold numerically for $E_{0}$ and $E_{1}$ (in fact, comfortably so).

Recall that the Lyapunov exponents $L_{N_{0}}(1/2,E)$ and $L_{2N_{0}}(1/2,E)$ are defined as integrals over $\mathbb{T}^{2}$ in (5). We make some remarks concerning numerical integration errors. We would like to emphasize that there is no viable deterministic bound on the numerical integration error because the integrand $\log\|M_{N_{0}}(\lambda,E;x,y)\|$ can have a large derivative in $x$ and especially $y$. Instead, we independently verified the numbers in Table 1 with a Monte-Carlo integration using $P=2001^{2}$ pseudorandom integration points, e.g., we still obtain $L_{N_{0}}(1/2,E)\approx 0.08071$ in that case. For Monte-Carlo integration the integration error for $L_{N_{0}}(1/2,E)$ with $E\in\{E_{0},E_{1}\}$ can be estimated as follows.

The random numerical integration can be studied through the random variable

[TABLE]

where we assume that the $P$ random variables $(x_{1},y_{1}),\ldots,(x_{P},y_{P})$ are independently drawn from Lebesgue measure on $[0,1]^{2}$.

Lemma 4**.**

Let $N\geq 1$, $P=2001^{2}$ and $0<\epsilon<1$. With respect to uniform probability measure, it holds that

[TABLE]

For example, if we take $\epsilon=7$ and use the numerical result $Q_{N_{0}}(P)\geq 0.08\times 10^{-2}$, then we find that $\mathbb{P}(L_{N_{0}}<10^{-2})<0.05\%$. The bound is also informative for $\epsilon=0.1$ in which case it implies $L_{N_{0}}>0.0797$ with probability $>75\%$. These bounds and their analogs for $L_{2N_{0}}$ indicate the correctness of the numerical data in Table 1.

Proof of Lemma 4.

By submultiplicativity of the norm, $\|AB\|\leq\|A\|\|B\|$, we have

[TABLE]

Since

[TABLE]

and $|v_{j}|\leq 2$, it follows that

[TABLE]

Using $|E_{1}|+1\leq 3.5$, the right-hand side is bounded by $4/3$ and so, by independence,

[TABLE]

The lemma now follows from applying Chebyshev’s inequality with deviation $\epsilon\times 10^{-2}$. ∎

3.3. Numerical investigation of condition (iii)

As mentioned in the discussion after Theorem 3, condition (iii) is a large-deviation estimate at the initial scale and plays a different role from conditions (i) and (ii).

In this section, we present two different kinds of graphs of the Lyapunov exponent $L_{N}(1/2,E;x,y)$ along slices of $\mathbb{T}^{2}$ of constant $x$- and $y$-values for comparison purposes.

• •

Figure 2 shows $L_{N}(1/2,0;x,y)$ along normal $x$- and $y$-slices. In this case, there is no visible bad set $\mathcal{B}_{N}(1/2,0)$.

• •

Figure 3 shows $L_{N}(1/2,E_{*};x,y)$ along slices of constant $x$- and $y$-values which are designed to contain a point in the bad set $\mathcal{B}_{N}(1/2,E_{*})$ where the Lyapunov exponent dips close to zero. We find such “bad slices” by taking the energy $E_{*}$ to be an eigenvalue of the corresponding finite-size operator $H_{N}(0.5,0.5)$ and looking for a pair of slices that intersect at a nearby point along which the Lyapunov exponent is minimal.

From both figures, we see that the $x$- and $y$-dependence of $L_{N}(1/2,E;x,y)$ are starkly different. The fact that the $y$-dependence is more oscillatory is expected in view of the expression $ny+x$ in (4). The main difference between Figures 2 and 3, is that in Figure 3 the Lyapunov exponent gets close to [math] at a few special points along both slices. These special points contribute to the bad set $\mathcal{B}_{N}$ defined in (6).

It turns out that the method we use to construct these “bad slices” by hand (described in the second bullet point above) only works for relatively small $N$ and this is the reason why we take $N=100$ in Figure 3. This fact can be considered good news in view of condition (iii). Going a step further, we might assume that the only possible reason for the occurrence of a bad point $(x,y)\in\mathcal{B}_{N}(\lambda,E)$ is that $E$ is an eigenvalue of the finite-size operator $H_{N}(x,y)$. With this assumption, we can estimate the size of the bad set using the data gathered in Section 3.1. Indeed, it suffices to count what is the maximum number of times that an energy $E$ occurs as an eigenvalue of $H_{N}(x,y)$ for $(x,y)$ ranging over the $2001\times 2001$ toroidal grid. The resulting count, when rounding energies to 8 decimal digits, is $31$. This suggests a rough estimate on $|\mathcal{B}_{N_{0}}(1/2,0)|$ of $31/(2001)^{2}<8\times 10^{-5}$, which is still quite far from what is required for condition (iii).

In light of these results, we believe that finding an analytical proof of condition (iii) makes for an important open problem.

Supplemental material

The MATLAB code used to perform the computations presented in this paper and the complete numerical output (to 16-digit accuracy) are available through an online repository [22].

Acknowledgments

The authors thank Wilhelm Schlag for encouragement and advice. The computations in this paper were run on the Odyssey cluster supported by the FAS Division of Science, Research Computing Group at Harvard University.