Sharp Polynomial Velocity Decay Bounds for Multidimensional Periodic Schrödinger Operators


Abstract

We investigate periodic Schrödinger operators in arbitrary dimensions in the large coupling regime. Our results establish that both the Lieb–Robinson velocity and the asymptotic velocity decay at an inverse polynomial rate in the coupling, with the precise exponent determined by the period of the underlying potential. In particular, we obtain sharp polynomial decay rates that capture the precise dependence on the periodic structure.

1 Introduction↩︎

Periodic quantum systems, modeled by Schrödinger operators with periodic potentials, are typically associated with ballistic transport, where wave packets spread linearly in time due to the underlying translational symmetry. These operators rigorously model quantum dynamics in crystalline and structured media, capturing the effects of translational symmetry on spectral and transport properties. They provide a natural setting for analyzing the relationship between spectral types and dynamical behavior. Surprisingly, recent work has revealed that while periodic systems do exhibit ballistic transport, the associated transport speed can be asymptotically small. This raises new questions about the subtle interplay between periodicity, spectral structure, and propagation speed.

The mathematical study of periodic Schrödinger operators has a long and rich history, tracing back to the foundational works of Bloch and Floquet in the early 20th century [1], [2]. These contributions established that solutions to the time-independent Schrödinger equation with periodic potentials can be expressed in terms of Bloch waves, leading to the band theory of solids and the identification of spectral bands and gaps [3]. From a mathematical perspective, periodic Schrödinger operators became central objects in spectral theory, particularly following the mid-20th century, with key developments such as the characterization of absolutely continuous spectra and the formulation of Floquet–Bloch theory, which reduces the spectral analysis to a direct integral over quasi-momentum space [4], [5]. In recent decades, attention has expanded beyond static spectral properties to include dynamical behavior, such as quantum transport and wavepacket spreading, with periodic systems often exhibiting ballistic dynamics and dipersive spreading under suitable conditions [6][13]. These studies have deepened our understanding of the mathematical structure underlying quantum systems and have had significant impact on theoretical and applied physics, particularly in the design and analysis of photonic crystals and metamaterials [14], [15]. We refer the reader to the survey [16] for a comprehensive overview and additional references and to [17] for additional details regarding ballistic motion. More recently, interdisciplinary approaches combining spectral theory with techniques from algebra and combinatorics have led to breakthroughs beyond what could be achieved through spectral methods alone. For example, some studies have focused on the irreducibility of associated algebraic varieties and its spectral consequences [18][23].

While periodic Schrödinger operators are classically associated with ballistic transport and absolutely continuous spectrum, recent results have revealed a more nuanced picture. In particular, it has been shown that although wave packets in periodic systems do spread linearly in time, the rate of this spreading (the group velocity) can be made arbitrarily small. This phenomenon, first rigorously demonstrated in one-dimensional settings [24], shows that the group velocity can become vanishingly small as the variations in a non-degenerate periodic potential are made increasingly large. A similar phenomenon has also been observed in quantum walks in periodic fields, as shown in the recent work [25], where it was proved that the propagation velocity remains positive yet can be exponentially suppressed by carefully tuning the periodic structure. Moreover, existing proofs often rely heavily on fairly delicate one-dimensional techniques, particularly in the analysis of Lieb-Robinson-type velocity bounds, so the analysis of this phenomenon in higher dimensions was not understood and would indeed require substantial novel ideas and insights. This raises the question of whether more general or structurally simpler arguments can reveal the same effect and offer broader insight into the relationship between spectral structure and transport speed in periodic media.

In this work, we address the higher-dimensional setting by employing a novel synthesis of complex analysis and perturbation theory to eschew the delicate one-dimensional arguments used previously. Specifically, we show that for non-degenerate periodic potentials, the propagation velocity can be made arbitrarily small by increasing the amplitude of the potential, and we give a sharp rate of decay. This approach not only broadens the dimensional scope of the phenomenon but also gives sharp estimates on the velocities and furthermore even yields a simpler and more transparent proof in the one-dimensional setting.

The remainder of the paper is organized as follows. In Section 2, we recall some terminology and formulate our main results precisely. We use Rayleigh–Schrödinger perturbation series in Section 3 to derive some useful estimates on eigenvalues, which we then employ in Section 4 to prove the main results.

Acknowledgements↩︎

The authors are grateful to Günter Stolz for helpful discussions. W.L. thanks the Department of Mathematics at UC Berkeley for its hospitality where part of this work were done during his visits in Fall 2024 and Fall 2025.

2 Setting and Main Results↩︎

We consider \(V: {\mathbb{Z}}^d \to {\mathbb{R}}\) and the operator \(H = \Delta + V\) on \(\ell^2({\mathbb{Z}}^d)\) given by \[[H\psi](n) = V(n) \psi(n) + \sum_{|n-m|_1=1} \psi(m), \quad n \in {\mathbb{Z}}^d,\] for which the function \(V:{\mathbb{Z}}^d \to {\mathbb{R}}\) is periodic. Given \(p = (p_1,\ldots,p_d) \in {\mathbb{N}}^d\), we denote the period lattice by \(p{\mathbb{Z}}^d := p_1 {\mathbb{Z}}\oplus \cdots \oplus p_d {\mathbb{Z}}\). We say that \(V\) is \(p\)-periodic1 if \(V(n+m)=V(n)\) for all \(n \in {\mathbb{Z}}^d\) and \(m \in p {\mathbb{Z}}^d\). We define the fundamental cell by putting \([\ell] = \{0,1,\ldots, \ell -1\}\) for \(\ell \in {\mathbb{N}}\) and \[\begin{align} W = [p_1] \times [p_2] \times \cdots\times [p_d], \end{align}\] and we say that \(V\) is non-degenerate if \(V\) is injective on \(W\). The separation of \(V\) is defined by \[\mathop{\mathrm{sep}}V = \min\{|V(n) - V(m) |: n,m \in W, \;n \neq m \}.\] One has \(\mathop{\mathrm{sep}}(\mu V) = |\mu| \mathop{\mathrm{sep}}(V)\) for each \(\mu \in {\mathbb{R}}\) and we see that \(V\) is non-degenerate if and only if \(\mathop{\mathrm{sep}}(V) > 0\).

Our main results concern the unitary evolution \(e^{-{{\mathbf{i}}}tH}\), which is known to propagate ballistically in the sense that the evolution of the position operator, \(X\), satisfies: \[\tfrac{1}{t}\underbrace{e^{{{\mathbf{i}}}tH} X e^{-{{\mathbf{i}}}tH}}_{=:X(t)} \xrightarrow{\;\mathrm{s} \;} G,\] with \(G\) a self-adjoint operator having trivial kernel. This was proved first for continuum Schrödinger operators [6], and later generalized to many other systems, including discrete Schrödinger operators, quantum walks, and more general periodic graph operators [7][9], [11], [12].

To study transport in the strong coupling regime, we introduce the scaled periodic Schrödinger \[H_\mu := \Delta + \mu V,\] where the coupling parameter \(\mu>0\) is taken to be large. Note that the introduction of a large parameter \(\mu\) amplifies the variations in the potential landscape. Physically, this corresponds to placing increasingly high barriers between sites of different potential values. As \(\mu\) increases, tunneling between distinct potential sites is suppressed. Consequently, the propagation of the quantum state slows down, and in the limit \(\mu \to \infty\), one expects a strong suppression of transport.

One measure of the rate of spreading is quantified by the asymptotic velocity. For a state \(\psi \in D(X)\), we denote \[\begin{align} \label{eq:vAsyDef} v_{\mathrm{asy}}(H,\psi) := \limsup_{t \to \infty} \frac{1}{t} \|X e^{-{\mathbf{i}}t H} \psi\|, \end{align}\tag{1}\] where \(X:D(X) \subseteq \ell^2({\mathbb{Z}}^d) \to (\ell^2({\mathbb{Z}}^d))^d\) denotes the (vector-valued) position operator \[\begin{align} X\psi = (X_1\psi,\ldots, X_d\psi), \quad [X_i \psi](x) = x_i \psi(x). \end{align}\] We then define the asymptotic velocity of \(H\) to be \[\label{def:v-asy} v_{\mathrm{asy}}(H)= \sup_{\substack{\psi \in D(X) \\ \|\psi\|=1}} v_{\mathrm{asy}}(H,\psi).\tag{2}\] We use here the refined version 2 , involving the supremum over normalized states in the domain of \(X\), which was first considered in an ongoing work on the velocity of unitary models (quantum walks and CMV matrices) [26].

Here we note that when \(H\) is periodic, the limit of the quantity on the right-hand-side of 1 exists so one can replace \(\limsup\) by \(\lim\) in the setting in which we work.

Our main result shows a sharp rate of decay for \(v_{\mathrm{asy}}\) in the regime of large coupling constant.

Theorem 1. Suppose \(V:{\mathbb{Z}}^d \to {\mathbb{R}}\) is \(p\)-periodic and non-degenerate with \(p_0 := \min\{p_i : i=1,2,\ldots,d\}\). Then, as \(\mu \to \infty\), one has \[\label{eq:vAsyBounds} v_{\mathrm{asy}}(H_\mu)= C\mu^{-p_0+1} + O(\mu^{-p_0}),\tag{3}\] where \(C>0\) is a constant depending on \(V\) and \(p\). Moreover, \[\label{eq:vAsydeltaBounds} v_{\mathrm{asy}}(H_\mu,\delta_0) = c\mu^{-p_0+1} + O(\mu^{-p_0})\tag{4}\] for a constant \(c>0\) depending on \(V\) and \(p\).

In fact, our argument gives more information than is stated in the theorem in the sense that the asymptotic velocity in each coordinate direction is quantified explicitly as a function of the period in that direction (cf.@eq:eq:GinormExpansion ). Broadly speaking, transport is suppressed the most in directions with larger periods. For example, in a two dimensional system, if \(p_1=1\) and \(p_2>1\), then the potential is constant along each line parallel to the horizontal axis, creating a propagating channel along which the wave packet freely travels, leading to the lack of decay of \(v_{\mathrm{asy}}\) seen above. However, the transport in the vertical axis is suppressed in the sense that \(\lim X_2(t)/t\) has order \(\mu^{-p_2+1}\).

A related estimate was obtained in the one-dimensional setting in [24], where an upper bound of the form \(v_{\mathrm{asy}}(H_\mu) \leq C \mu^{-p+1}\) was established (where \(p\) is the period of the potential). Theorem 1 provides not only a bound above and below but a sharp polynomial (in \(\mu\)) asymptotic statement in arbitrary dimension.

One can also quantify the notion of propagation (transport) in quantum systems using a single-body analogue of the Lieb-Robinson velocity [27], of the form \[\label{eq:genericLRbound} |\langle \delta_n, e^{-{\mathbf{i}}tH_\mu} \delta_m \rangle| \lesssim e^{-\rho_0 (|n - m|_1 - v_{{\mathrm{LR}}} |t|)}, \quad \text{for all } n,m \in {\mathbb{Z}}^d, \;t \geq 0,\tag{5}\] for suitable constants \(\rho_0, v_{{\mathrm{LR}}} > 0\). Here and in the following, \(\{\delta_n : n\in{\mathbb{Z}}^d\}\) denotes the canonical basis of \(\ell^2({\mathbb{Z}}^d)\). Inequality (5 ) implies that \(v_{{\mathrm{LR}}}\) serves as an effective upper bound on the speed of information propagation, up to exponentially decaying corrections. The physical interpretation is as follows [28]: starting from the localized initial state \(\delta_m\), the probability that the distance traveled exceeds \(v|t|\) at time \(t\), for any velocity \(v > 0\), is given by \[\begin{align} \sum_{n:\, |n - m|_1 > v|t|} |\langle \delta_n, e^{-{\mathbf{i}}tH_\mu} \delta_m \rangle|^2 &\lesssim& e^{-2\rho |t| (v - v_{{\mathrm{LR}}})}\sum_{n:\, |n - m|_1 > v|t|} e^{-2\rho (|n - m|_1 - v |t|)} \notag\\ &\lesssim& e^{-2\rho |t| (v - v_{{\mathrm{LR}}})}. \end{align}\] Thus, for any \(v > v_{{\mathrm{LR}}}\), the probability of observing propagation beyond the linear bound \(|n - m|_1 = v|t|\) decays exponentially in time. This defines a so-called Lieb-Robinson light cone with effective velocity \(v_{{\mathrm{LR}}}\), see e.g., [29][32].

Our second main result is a bound on the Lieb-Robinson velocity \(v_{\mathrm{LR}}\) decaying at the sharp rate \(\mu^{-p_0+1}\) as \(\mu \to \infty\), with \(p_0:=\min\{p_1,\ldots,p_d\}\).

Theorem 2. Suppose \(V:{\mathbb{Z}}^d \to {\mathbb{R}}\) is \(p\)-periodic and non-degenerate with \(p_0 := \min\{p_i: i=1,2,\ldots,d\}\). Then, for any \(\rho_0>0\), there are constants \(C,C_1,\mu_0>0\) depending on \(d\), \(p\), \(\mathop{\mathrm{sep}}V\), and \(\rho_0\) such that 5 holds with \(v_{\mathrm{LR}}= C_1\mu^{-p_0+1}\). That is, for all \(n,m \in {\mathbb{Z}}^d\) and all \(\mu \geq \mu_0\), one has \[\label{eq:mainBound} | \langle \delta_n, e^{-{{\mathbf{i}}}tH_\mu} \delta_m \rangle| \leq C e^{-\rho_0(|n-m|_1-C_1 \mu^{-p_0+1}|t|)}\tag{6}\] for all \(t \geq 0\).

Furthermore, this is optimal in the sense that if a bound of the form 5 holds for all \(\mu\) large, then \(v_{\mathrm{LR}}\geq c\mu^{-p_0+1}\) for a constant \(c>0\).

Theorem 2 shows that the Lieb–Robinson velocity decays at the sharp rate \(\mu^{-p_{0}+1}\) in the strong coupling regime. In comparison with earlier work in [24], Theorem 2 identifies the precise rate of decay and applies in higher dimensions, both of which appear to lie beyond the scope of [24]. Thus, our result both broadens the scope and sharpens the known velocity bounds. We emphasize that Theorem 2 yields a polynomial rate of decay in the coupling constant rather than an exponential one. The threshold \(\mu_{0}\) can be chosen independently of the period of the potential, and in fact one may take \(\mu_{0}=1/\varepsilon_{0}\), where \(\varepsilon_{0}=\varepsilon_0(d,\mathop{\mathrm{sep}}V,\rho_0)\) given in Lemma 1 below.

We also note that both the Lieb–Robinson velocity \(v_{\mathrm{LR}}\) and the asymptotic velocity \(v_{\mathrm{asy}}\) exhibit the same scaling behavior with respect to the coupling constant, namely \(\mu^{-p_{0}+1}\). The same phenomenon has been observed in the context of quantum walks [26], where these two velocities also scale in the same way. This raises the natural question of whether there exist settings in which the two notions of velocity exhibit different scaling behavior.

Remark 3. It can be seen from the proofs of the main results that the argument is fundamentally graph theoretic in nature, so it generalizes readily to the case of periodic Schrödinger operators on \({\mathbb{Z}}^d\)-periodic graphs. If \(\mathcal{G} = (\mathcal{V}, \mathcal{E})\) is \({\mathbb{Z}}^d\)-periodic and \(V:\mathcal{V} \to {\mathbb{R}}\) is \(p\)-periodic and non-degenerate, then the same proofs as in the case \(\mathcal{G} = {\mathbb{Z}}^d\) produce similar results for the operator \(\Delta + V\) (with \(\Delta\) the graph Laplacian) with \(p_0\) equal to the minimal combinatorial distance separating two distinct vertices in the same \(p{\mathbb{Z}}^d\)-orbit.

Let us reiterate the strength of the results: we improve the rate of decay from [24], obtain a sharp rate of decay, and generalize to higher dimensions, all with a conceptually simpler argument that can even be applied to \({\mathbb{Z}}^d\)-periodic graphs other than the square lattice \({\mathbb{Z}}^d\) itself.

3 Asymptotics of Eigenvalues and Eigenvectors↩︎

We will need to derive some asymptotic statements for eigenvalues and eigenvectors of Floquet matrices. Throughout the discussion, we fix a nonempty finite set2 \(S\) and consider \({\mathbb{C}}^{S \times S}\), the collection of \(S \times S\) matrices with complex entries. Let \(\{e_n : n \in S\}\) be the standard basis of \({\mathbb{C}}^S\). In the discussion below, one can picture \(B\) as the Floquet operator of a Laplacian on a suitable fundamental domain and \(D\) as the associated potential. We note in particular that it is not needed to assume the operators are Hermitian, so this can be applied to complex potentials and complex quasimomenta.

Theorem 4. Assume \(S\) is a nonempty set \(B,D \in {\mathbb{C}}^{S\times S}\) are matrices such that:

  • \(D = \mathrm{diag}(D_n)\) with pairwise distinct entries,

  • \(B_{nn}=0\) for all \(n \in S\).

Let \(A_\varepsilon = D + \varepsilon B\). For \(\varepsilon \in {\mathbb{R}}\) sufficiently small, the eigenvalues and eigenvectors have convergent expansions of the form \[\begin{align} \eta_n(\varepsilon) = \sum_{r=0}^\infty \eta_n^{(r)} \varepsilon^r, \qquad u_n(\varepsilon) = \sum_{r=0}^\infty u_n^{(r)}\varepsilon^r, \end{align}\] respectively, where \(\langle e_n, u_n^{(r)} \rangle = 0\) for all \(r \geq 1\), \[\begin{align} \label{eq:etaUr610} \eta_n^{(0)} = D_n, \quad u_n^{(0)} = e_n, \end{align}\tag{7}\] and \[\begin{align} \label{eq:etaUr611} \eta_n^{(1)} = 0, \quad (u_n^{(1)})_m = \frac{B_{mn}}{D_n-D_m}, \quad m \neq n. \end{align}\tag{8}\] Furthermore, for all \(r \geq 2\): \[\begin{align} \tag{9} \eta_n^{(r)} &= \sum_{m \neq n} B_{nm}(u_n^{(r-1)})_m \\ \tag{10} \left(u_n^{(r)} \right)_m &= \frac{1}{D_n-D_m} \left( \sum_{\ell \neq n} B_{m\ell} (u_n^{(r-1)})_\ell - \sum_{s=2}^{r-1} \eta_n^{(s)} (u_n^{(r-s)})_m\right), \quad m \neq n. \end{align}\]

Remark 5. Note that the sum in 9 starts from \(s=2\) and in particular is absent for \(r=2\).

Proof of Theorem 4. This follows from the Rayleigh–Schrödinger expansion. To keep the paper more self-contained, we give the details here. By definition, \(A_\varepsilon\) is an analytic function of \(\varepsilon\). Since \(D\) has simple spectrum, the existence of an expansion follows from perturbation theory (compare e.g.[33]). We thus fix \(n \in S\), \(\eta_n(0) = D_n\), and \(u_n(0) = e_n\), and consider their analytic continuations (note that this immediately gives us 7 ). Since any analytic vector-valued function has analytic components, we can always normalize by the \(n\)th component to obtain an expansion for \(u_n\) such that \((u_n(\varepsilon))_n \equiv 1\), whence \((u_n^{(r)})_n =0\) for all \(r \geq 1\).

Expanding \(A_\varepsilon u_n(\varepsilon) = \eta_n(\varepsilon)u_n(\varepsilon)\) and collecting terms gives \[\label{eq:seriesCollected} Du_n^{(0)} + \sum_{r=1}^\infty (Du_n^{(r)} + Bu_n^{(r-1)})\varepsilon^r = \eta_n^{(0)}u_n^{(0)} + \sum_{r=1}^\infty \sum_{s=0}^r \eta_n^{(s)} u_n^{(r-s)} \varepsilon^r.\tag{11}\] Since \(D\) is diagonal and \(\langle e_n, u_n^{(r)} \rangle = 0\) for \(r \geq 1\), we have \[\langle e_n, Du_n^{(r)} \rangle = 0\] for all \(r \geq 1\). Using this, project both sides of 11 onto \(e_n\) to get: \[D_n + \sum_{r=1}^\infty \langle e_n, B u_n^{(r-1)} \rangle \varepsilon^r = \eta_n^{(0)} + \sum_{r=1}^\infty \eta_n^{(r)} \varepsilon^r,\] giving3 \[\eta_n^{(r)} = \langle e_n, B u_n^{(r-1)} \rangle = \sum_{m \neq n} B_{nm}(u_n^{(r-1)})_m,\] which proves 9 . Note that \(\eta_n^{(1)} = 0\) follows directly.4

Projecting 11 onto \(e_m\) with \(m \neq n\) yields \[D_m (u_n^{(r)})_m + (Bu_n^{(r-1)})_m = \sum_{s=0}^{r} \eta_n^{(s)} (u_n^{(r-s)})_m,\] which, with the help of \(u_n^{(0)} = e_n\), \(\eta_n^{(0)}=D_n\), becomes: \[(D_m-D_n) (u_n^{(r)})_m + (Bu_n^{(r-1)})_m = \sum_{s=1}^{r-1} \eta_n^{(s)} (u_n^{(r-s)})_m,\] which proves 10 (and the second identity in 8 ) after rearranging and recalling that \(\eta_n^{(1)}=0\). ◻

Let us compute a few terms in these expansions by hand and then establish a general pattern. As already noted, we have \[\eta_n^{(0)}=D_n, \quad u_n^{(0)} = e_n, \quad \eta_n^{(1)}=0, \quad (u_n^{(1)})_m = \frac{B_{mn}}{D_n - D_m}, \quad m \neq n.\] In general, we see from Theorem 4 that \(u_n^{(r)}\) is determined by \(\{\eta_n^{(s)} : s<r\}\) and by the vectors \(\{u_n^{(s)} : s<r\}\), (as well as the requirement \(\langle e_n, u_n^{(r)} \rangle = 0\) for \(r\geq 1\)) whereas \(\eta_n^{(r)}\) is determined by \(u_n^{(r-1)}\).

The next coefficients of \(\eta\) and \(u\) can then be computed as \[\label{eq:eta2expression} \eta_n^{(2)} = \sum_{m \neq n} B_{nm} (u_n^{(1)})_m = \sum_{m \neq n} \frac{B_{nm}B_{mn}}{D_n-D_m}.\tag{12}\] and \[\begin{align} \nonumber\left(u_n^{(2)} \right)_m & = \frac{1}{D_n - D_m}\sum_{m_1 \neq n} B_{mm_1}(u_n^{(1)})_{m_1} \\ & = \sum_{m_1 \neq n} \frac{ B_{mm_1} B_{m_1n}}{(D_n - D_m)(D_n - D_{m_1})}, \end{align}\] which allows us to compute the third-order coefficients: \[\label{eq:eta3expression} \eta_n^{(3)} = \sum_{m \neq n} B_{nm} (u_n^{(2)})_m = \sum_{\substack{m \neq n \\ m_1 \neq n} } \frac{ B_{nm} B_{mm_1} B_{m_1n}}{(D_n - D_m)(D_n - D_{m_1})}\tag{13}\] and \[\begin{align} \nonumber\left(u_n^{(3)} \right)_m & = \frac{1}{D_n - D_m}\sum_{m_1 \neq n} B_{mm_1}(u_n^{(2)})_{m_1} - \frac{1}{D_n-D_m}\left( \sum_{m_1 \neq n} \frac{B_{nm_1}B_{m_1n}}{D_n-D_{m_1}} \right)\left( \frac{B_{mn}}{D_n - D_m} \right) \\ & = \sum_{m_1,m_2 \neq n} \frac{ B_{mm_1} B_{m_1m_2} B_{m_2n} }{(D_n - D_{m})(D_n - D_{m_1})(D_n - D_{m_2})} - \sum_{m_1 \neq n} \frac{B_{mn} B_{nm_1} B_{m_1n}}{(D_n-D_m)^2(D_n - D_{m_1})}. \end{align}\]

We can give this a useful interpretation by thinking of \(B\) as an operator on a (directed) graph \(\Gamma\) with vertex set \(S\) and an edge from \(n\) to \(m\) whenever \(B_{nm} \neq 0\). We then define a path of length \(\ell\) from \(m\) to \(n\) (notation: \(\gamma:m\to n\)) in \(\Gamma\) to be a finite sequence \(\gamma = (n_0,\ldots,n_\ell)\) where \(n_0=m\), \(n_\ell=n\), and \((n_{s-1},n_s)\) is an edge for every \(1 \le s \le \ell\). Write \(|\gamma|:=\ell\) for the length of \(\gamma\). We say that \(\gamma\) is irreducible if \(n_s \neq n\) for all \(s=1,2,\ldots,\ell-1\). For a path \(\gamma=(n_0,\ldots, n_\ell)\), we denote \(B_\gamma = B_{n_0,n_1} \cdots B_{n_{\ell-1},n_\ell}\). We call \(\gamma\) a loop if \(n_0 = n_\ell\).

With this interpretation, we see that \(\eta_n^{(r)}\) has the form \[\label{eq:etaloopExpansion} \eta_n^{(r)} = \sum_{\substack{\gamma:n \to n \\ |\gamma|=r}} f(\gamma,D)B_\gamma,\tag{14}\] where the notation indicates that the sum runs over loops of length \(r\) from \(n\) to itself, and the coefficient \(f(\gamma,D)\) depends on the loop \(\gamma\) and the values assumed by \(D\) on vertices in the loop. Likewise, \(u_n^{(r)}\) has the form \[\label{eq:uloopExpansion} (u_n^{(r)})_m = \sum_{\substack{\gamma:m \to n \\ |\gamma|=r}} g(\gamma,D)B_\gamma,\tag{15}\] where again the notation indicates that the sum runs over paths from \(m\) to \(n\) of length \(r\) and once again the coefficient depends on the path \(\gamma\) and the values assumed by \(D\) on the path.

Corollary 6. With assumptions as above, \(\eta_n^{(r)}\) has an expansion of the form 14 and \(u_n^{(r)}\) has an expansion of the form 15 for every \(n \in S\) and every \(r \geq 1\).

For an irreducible loop \(\gamma:n \to n\), \(\gamma = (n_0,\ldots,n_r)\), we have \[\label{eq:fgamma:irrLoop} f(\gamma,D) = \prod_{j=1}^{r-1} \frac{1}{D_n-D_{n_j}}.\qquad{(1)}\] Likewise, for \(m\neq n\) and an irreducible path \(\gamma:m \to n\), \(\gamma = (n_0,\ldots,n_r)\), we have \[\label{eq:fgamma:irrPath} g(\gamma,D) = \prod_{j=0}^{r-1} \frac{1}{D_n-D_{n_j}}.\qquad{(2)}\]

Proof. The proof is by induction on \(r\). The desired statement (including ?? and ?? ) is already established for \(r =1,2,3\) in the discussion above. Inductively, assume the statement of the theorem holds up to order \(r-1\).

If \(\gamma = (n_0,\ldots,n_\ell)\) and \(\gamma' = (m_0,\ldots,m_{\ell'})\) are paths with \(n_\ell = m_0\), we write \(\gamma \circledast \gamma' = (n_0,\ldots,n_\ell,m_1,\ldots,m_{\ell'})\) for their amalgamation. Notice that \[B_{\gamma \circledast \gamma'} = B_\gamma B_{\gamma'}.\] We then have (with the help of the inductive hypothesis) \[\begin{align} \eta_n^{(r)} = \sum_{m} B_{nm}(u_n^{(r-1)})_m & = \sum_m B_{nm} \sum_{\substack{\gamma:m \to n \\ |\gamma| = r-1}} g(\gamma,D)B_\gamma \\ & = \sum_m \sum_{\substack{\gamma:m \to n \\ |\gamma| = r-1}} g(\gamma,D) B_{(n,m) \circledast \gamma}. \end{align}\] For \(\gamma\) a path of length \(r-1\) from \(m\) to \(n\), we see that \((n,m) \circledast \gamma\) is a loop from \(n\) to \(n\) of length \(r\), which is irreducible if and only if \(\gamma\) is irreducible. Thus, the inductive step for \(\eta\) and ?? is concluded.

Similarly, we have \[\begin{align} \left(u_n^{(r)} \right)_m & = \frac{1}{D_n - D_m} \left( \sum_{m_1 }B_{mm_1}(u_n^{(r-1)})_{m_1} - \sum_{s=2}^{r-1} \eta_n^{(s)} (u_n^{(r-s)})_m\right) \\ & = \frac{1}{D_n - D_m} \left( \sum_{m_1} B_{mm_1} \sum_{\substack{\gamma:m_1\to n \\ |\gamma|=r-1}} g(\gamma,D) B_\gamma - \sum_{s=2}^{r-1} \left[ \sum_{\substack{\gamma:n\to n \\ |\gamma|=s}} f(\gamma,D)B_\gamma\right] \left[\sum_{\substack{\gamma':m\to n \\ |\gamma'|=r-s}} g(\gamma',D)B_{\gamma'}\right]\right) \\ & = \sum_{m_1} \sum_{\substack{\gamma:m_1\to n \\ |\gamma|=r-1}} \frac{g(\gamma,D)}{D_n - D_m} B_{(m,m_1) \circledast \gamma} - \sum_{s=2}^{r-1} \sum_{\substack{\gamma:n\to n \\ |\gamma|=s}} \sum_{\substack{\gamma':m\to n \\ |\gamma'|=r-s}} \frac{f(\gamma,D) g(\gamma',D)}{D_n - D_m}B_{\gamma'\circledast \gamma}. \end{align}\] In the first term, we observe that \((m,m_1)\circledast \gamma\) is a path of length \(r\) from \(m\) to \(n\) (which is irreducible if and only if \(\gamma\) is irreducible), while, in the second term, \(\gamma'\circledast\gamma\) is a path from \(m\) to \(n\) of length \(r\) (and \(\gamma' \circledast\gamma\) is never irreducible), concluding the inductive step for \(u\) and ?? . ◻

Remark 7. These ideas and expansions are useful and have been employed in other contexts. As this paper was being completed, the preprint [34] was posted, which, independent of us, employed a similar expansion in the study of flat bands.

4 Proof of Main Results↩︎

4.1 Floquet Theory↩︎

Let us briefly review a few aspects of Floquet theory that we need. For more background, see [16] and for proofs of relevant statements that apply in the present setting, see [18].

Recall that \(W = [p_1] \times \cdots \times [p_d] \cong {\mathbb{Z}}^d / p {\mathbb{Z}}^d\) is the fundamental domain, where \([m]=\{0,1,\ldots,m-1\}\). Put \[P =\# W = p_1p_2 \cdots p_d.\] Given \(z = (z_1, z_2, \ldots, z_d) \in ({\mathbb{C}}^*)^d\) (where \({\mathbb{C}}^* = {\mathbb{C}}\setminus \{0\}\)), we write \(H(z)\) for the restriction of \(H\) to the space \[\mathscr{H}_p (z) := \{\psi \in {\mathbb{C}}^{{\mathbb{Z}}^d} : \psi(n+p_j e_j) = z_j \psi(n)\}.\]

The matrix \(H(z)\) is self-adjoint and hence has real eigenvalues for any \(z \in (\partial {\mathbb{D}})^d\). For such \(z\), the spectral band functions \(\lambda_j(z)\), \(1 \le j \le P\), are then defined by listing the eigenvalues of \(H(z)\) in order with multiplicity: \[\mathop{\mathrm{spec}}H(z) = \{\lambda_j(z) : 1 \le j \le P \}, \quad \lambda_1(z) \leq \lambda_2(z) \leq \cdots \leq \lambda_P(z).\] For each \(z\) for which \(H(z)\) is diagonalizable, there is an invertible similarity transformation \(Q(z)\) such that \[[Q(z)]^{-1} H(z)Q(z) = \mathop{\mathrm{diag}}(\lambda_j(z)).\]

In this work, we are interested in the behavior of \(H_\mu = \Delta+ \mu V\) as \(|\mu| \to \infty\). Rescaling by \(\mu\), it suffices to consider \(A_\varepsilon := \varepsilon \Delta + V\) as \(|\varepsilon| \to 0\). Let us write \(A_\varepsilon(z)\) for the corresponding Floquet matrices.

We also define for \(\rho>0\) the set \[\mathop{\mathrm{{\mathcal{A}}}}(\rho) := \{ z \in ({\mathbb{C}}^*)^d : |\log|z_j||<\rho \text{ for all } 1\le j \le d\}.\]

Lemma 1. Assume \(V:{\mathbb{Z}}^d \to {\mathbb{R}}\) is \(p\)-periodic and non-degenerate, let \(\rho_0>0\) be given, and put \[\varepsilon_0 = \varepsilon_0(d,\mathop{\mathrm{sep}}V,\rho_0) := \frac{\mathop{\mathrm{sep}}V}{8d(1+\cosh (2\rho_0))}.\]

For each \(|\varepsilon| < 2 \varepsilon_0\), the spectrum of \(A_\varepsilon(z)\) is simple for every \(z \in \mathop{\mathrm{{\mathcal{A}}}}(2\rho_0)\), and each \(\eta_n(\varepsilon,z)\) is an analytic function of \((\varepsilon,z) \in (-2\varepsilon_0, 2\varepsilon_0) \times \mathop{\mathrm{{\mathcal{A}}}}(2\rho_0)\) for each \(n\).

Proof. If \(|\varepsilon| < 2 \varepsilon_0\), simplicity follows from the Gershgorin circle theorem, so the desired analyticity follows from eigenvalue perturbation theory [33]. ◻

Choosing \(|\varepsilon|\) small enough that Lemma 1 applies and \(z \in \mathop{\mathrm{{\mathcal{A}}}}(2\rho_0)\), we define \(\eta_n(\varepsilon,z)\) to be the Floquet eigenvalue of \(A_\varepsilon(z)\) analytically continued from \(\eta_n(0,z)=V(n)\), and let \(Q_\varepsilon(z)\) denote the matrices diagonalizing \(A_\varepsilon(z)\): \[\label{eq:QepsilonConjugatesAepsilon} [Q_\varepsilon(z)]^{-1} A_\varepsilon(z)Q_\varepsilon(z) = \mathop{\mathrm{diag}}(\eta_n(\varepsilon,z)).\tag{16}\] We choose \(Q_\varepsilon(z)\) by selecting the columns to be the analytic eigenvectors as in Theorem 4.

With notation as above, we define \[L^2({\mathbb{T}}^d, {\mathbb{C}}^W; \tfrac{{\mathrm{d}}\theta}{|{\mathbb{T}}^d|}) = \left\{ f:{\mathbb{T}}^d \to {\mathbb{C}}^W \;\vert \; \|f\| := \left[ \int_{{\mathbb{T}}^d} \|f(\theta)\|_{{\mathbb{C}}^W}^2 \, \frac{{\mathrm{d}}\theta}{|{\mathbb{T}}^d|} \right]^{1/2} < \infty \right\},\] where \({\mathbb{T}}:={\mathbb{R}}/(2\pi {\mathbb{Z}})\) and we write \({\mathrm{d}}\theta = {\mathrm{d}}\theta_1 \, \cdots {\mathrm{d}}\theta_d\) for the standard Lebesgue measure on \({\mathbb{T}}^d\). The Floquet transform \({\mathscr{F}}: \ell^2({\mathbb{Z}}^d) \to L^2({\mathbb{T}}^d, {\mathbb{C}}^W; \frac{{\mathrm{d}}\theta}{|{\mathbb{T}}^d|})\) is then given by \[\delta_{n+x \odot p} \mapsto e^{-{{\mathbf{i}}} \langle x, \cdot \rangle}e_n, \quad n \in W, \;x \in {\mathbb{Z}}^d,\] where we define \[x \odot p = (x_1p_1,\ldots,x_dp_d), \quad \langle x, \theta \rangle = \sum_{j=1}^d x_j\theta_j.\] By direct computations, one can see that \({\mathscr{F}}\) is unitary and conjugates \(A_\varepsilon\) to a decomposable operator, viz.: \[({\mathscr{F}}A_\varepsilon {\mathscr{F}}^{-1} g)(\theta) = A_\varepsilon(e^{{{\mathbf{i}}}\theta})g(\theta), \quad g \in L^2({\mathbb{T}}^d,{\mathbb{C}}^W; \tfrac{{\mathrm{d}}\theta}{|{\mathbb{T}}^d|}),\] where \(e^{{\mathbf{i}}\theta} =(e^{{\mathbf{i}}\theta_1}, \ldots, e^{{\mathbf{i}}\theta_d})\). For \(x,y \in {\mathbb{Z}}^d\) and \(n,m \in W\) given, we note by unitarity of \({\mathscr{F}}\) that \[\begin{align} \nonumber \langle \delta_{m+y \odot p}, e^{-{{\mathbf{i}}}tA_\varepsilon} \delta_{n+x \odot p} \rangle & = \langle {\mathscr{F}}\delta_{m+y \odot p}, {\mathscr{F}}e^{-{{\mathbf{i}}}t A_\varepsilon }{\mathscr{F}}^{-1}{\mathscr{F}}\delta_{n+x \odot p} \rangle \\ \label{eq:FloquetBlockRepn} & = \int_{{\mathbb{T}}^d} \langle e_m, e^{-{{\mathbf{i}}} t A_\varepsilon(e^{{{\mathbf{i}}}\theta})}e_n \rangle e^{-{{\mathbf{i}}} \langle (x-y),\theta \rangle} \, \frac{{\mathrm{d}}\theta}{|{\mathbb{T}}^d|}. \end{align}\tag{17}\] In our discussion of bounds on the Lieb–Robinson velocity, the desired estimates are most conveniently formulated in terms of bounds on the blocks of the unitary propagator with respect to the period lattice. Concretely (with the period \(p\) fixed), we denote \[A(x,y) = \Pi(W+y\odot p) \, A \, \Pi(W+x \odot p)^*, \quad x,y \in {\mathbb{Z}}^d,\] where \(\Pi(S)\) denotes the canonical projection \(\ell^2({\mathbb{Z}}^d) \to \ell^2(S)\). We note that each such block is represented as a matrix in \({\mathbb{C}}^{W \times W}\). With the help of 17 , we obtain the following representation of the blocks of \(e^{-{{\mathbf{i}}}t A_\varepsilon}\): \[\label{eq:duhamel} e^{- {{\mathbf{i}}} t A_\varepsilon}(x,y) = \int_{{\mathbb{T}}^d} e^{-{{\mathbf{i}}}t A_\varepsilon(e^{{{\mathbf{i}}}\theta})} e^{-{{\mathbf{i}}} \langle (x-y),\theta \rangle} \, \frac{{\mathrm{d}}\theta}{|{\mathbb{T}}^d|}.\tag{18}\]

Remark 8. Using the triangle inequality and well-known inequalities for matrix norms, we can pass between estimates on blocks and estimates on matrix elements of the propagator \(e^{-{\mathbf{i}}t A_\varepsilon}\) at the expense of some explicit \(p\)-dependent constants. Concretely, for \(n \in W + x\odot p\) and \(m \in W + y \odot p\), we have \[p_0|x-y|_1- C \leq |n-m|_1 \leq p_{\max} |x-y|_1 + C.\] where \(C = 2\sum p_j\) and \(p_{\max} = \max\{p_j\}\).

Likewise, with \(\|\cdot\|\) the standard Euclidean matrix norm and \(\|\cdot \|_\infty\) denoting the maximum modulus of an entry, we have the beloved bound \[\|B\|_\infty \leq \|B\| \leq P \|B\|_\infty.\] Taken together, these observations enable one to pass between estimates on blocks and matrix elements.

4.2 Eigenvalue Perturbation↩︎

By the nondegeneracy assumption and Lemma 1, the Floquet eigenvalues of \(A_\varepsilon\) have an expansion \[\eta_n(\varepsilon,z) = \sum_{r=0}^\infty \eta_n^{(r)}(z) \varepsilon^r\] for small \(\varepsilon\). Let us begin by discussing the low-order expansion coefficients.

Lemma 2. Assume \(V\) is \(p\)-periodic and non-degenerate, and put \(p_0 = \min \{p_1,\ldots,p_d\}\). We have:

  1. For \(r=0,1,2,\ldots,p_j-1\), \(\eta_n^{(r)}(z)\) is real and independent of \(z_j\). In particular, if \(r<p_0\), then \(\eta_n^{(r)}(z)\) is real and independent of \(z\).

  2. For \(r= p_j\), \(\eta_n^{(r)}(z)\) is of the form \[\label{etajp950:cosineform} \eta_n^{(r)}(z) = \sum_{j:p_j=p_0} c_{j,n}(z_{j}+z_{j}^{-1}) + \Phi_n(z)\tag{19}\] where the \(c_{j,n}\)’s are nonzero real constants independent of \(z\) and \(\Phi_n\) is a Laurent polynomial with real coefficients depending solely on the variables \(\{z_i : p_i < p_j\}\). In particular, if \(p_j=p_0\), then \(\Phi_n\) is a real constant.

If \(p_0=1\), then \(\Phi_n=0\) and \(c_{j,n} = 1\) for each \(j\) such that \(p_j = p_0\).

Proof. The statements in the case \(p_0 \geq 2\) follow from Corollary 6 while the statements for \(p_0=1\) follow from computations with the help of the Feynman–Hellmann theorem. ◻

This can then be put to use to get suitable estimates.

Lemma 3. Assume \(V:{\mathbb{Z}}^d \to {\mathbb{R}}\) is \(p\)-periodic and non-degenerate, let \(\rho_0>0\) be given, define \(\varepsilon_0\) as in Lemma 1, and let \(p_0 := \min\{ p_1,\ldots,p_d\}\).

  1. There is a constant \(C_1 = C_1(d, p, \mathop{\mathrm{sep}}V, \rho_0) > 0\) such that \[|\mathop{\mathrm{Im}}\eta_n(\varepsilon,z)| \leq C_1 \varepsilon^{ p_0}\] for all \(|\varepsilon| \leq \varepsilon_0\) and \(z \in \overline{\mathop{\mathrm{{\mathcal{A}}}}(\rho_0)}\), and \(n \in W\).

  2. There is a constant \(C_2 = C_2(d, p, \mathop{\mathrm{sep}}V, \rho_0) > 0\) such that \[\max\{\|I - Q_\varepsilon(z)\|,\|I-[Q_\varepsilon(z)]^{-1}\|\} \leq C_2\varepsilon\] for all \(|\varepsilon| \leq \varepsilon_0\) and \(z \in \overline{\mathop{\mathrm{{\mathcal{A}}}}(\rho_0)}\).

Proof. [item:imaginaryPart] On account of Lemma 2, we deduce that \[\eta_n(\varepsilon,z) = V(n) + \sum_{r=1}^\infty \eta_n^{(r)}(z)\varepsilon^r,\] where the \(\eta_n^{(r)}\) are analytic and each is real and constant for \(r < p_0\). The desired bound follows presently. The dependence of the constant follows from the form of \(\eta_n^{(r)}\) in Theorem 4 and induction.

[item:QepsBound] Noting that \(Q_\varepsilon(z)\) is analytic as a function of \((\varepsilon,z)\) with \(Q_0 = I\), this follows directly. The dependence of the constant follows from the form of \(u_n^{(r)}\) in Theorem 4 and induction. ◻

4.3 Proofs of Main Results↩︎

We now put everything together to prove the main results.

Proof of Theorem 1. From the definitions, we see that \[\label{eq:vAsyRelationGi} v_{\mathrm{asy}}(A_\varepsilon,\psi) = \|(G_1\psi,\ldots,G_d\psi)\|,\tag{20}\] where \(G_i\) denotes the \(i\)th component of the group velocity, which in turn is given by: \[\begin{align} \nonumber G_i&:= \operatorname{s-lim}_{t \to \infty} \frac{1}{t} X_i(t) = \operatorname{s-lim}_{t \to \infty} \frac{1}{t} e^{{\mathbf{i}}t A_\varepsilon}X_i e^{-{\mathbf{i}}t A_\varepsilon} \\ & \qquad = \frac{p_i}{|{\mathbb{T}}^d|} {\mathscr{F}}^{-1} \left[ \int^\oplus_{{\mathbb{T}}^d} \sum_n \frac{\partial \eta_n(\varepsilon,e^{{\mathbf{i}}\theta})}{\partial \theta_i} \Pi_n(\varepsilon, e^{{\mathbf{i}}\theta}) \, {\mathrm{d}}\theta \right] {\mathscr{F}} \label{eq:vasy} \end{align}\tag{21}\] in the strong sense, where we recall \(X_i\) denotes the position operator \(X_i \psi(x) = x_i \psi(x)\), \(\eta_n(\varepsilon,e^{{\mathbf{i}}\theta})\) denotes the eigenvalue of \(A_\varepsilon(e^{{\mathbf{i}}\theta})\) that is analytically continued from \(V(n)\), and \(\Pi_n(\varepsilon, e^{{\mathbf{i}}\theta})\) denotes projection onto the corresponding eigenspace; see e.g.[12] for a proof. In particular, \(G_i\) is, up to a unitary conjugation, given by a generalized multiplication operator and therefore \[\label{eq:GiNorm} \|G_i\| = p_i \max\left\{ \left| \frac{\partial \eta_n}{\partial \theta_i}(\varepsilon, e^{{\mathbf{i}}\theta}) \right| : \theta \in {\mathbb{T}}^d, \;n \in W \right\}.\tag{22}\]

For each \(i\), we note from Lemma 2 that one has \[\label{eq:eta95nDerivExpansion} \frac{\partial \eta_n}{\partial \theta_i} = 2c_{i,n} \sin( \theta_i) \varepsilon^{p_i} + O(\varepsilon^{p_i + 1}) ,\tag{23}\] with \(c_{i,n} \neq 0\). Putting together 22 and 23 , \[\label{eq:GinormExpansion} \|G_i\| = \widetilde{c}_i \varepsilon^{p_i} + O(\varepsilon^{p_i+1}), \quad i = 1,2,\ldots, d,\tag{24}\] where \(\widetilde{c}_i = 2p_i \max_n |c_{i,n}| \neq 0\).

Putting together 24 and 20 (and taking the maximum over \(i\)) gives \(v_{\mathrm{asy}}(A_\varepsilon) = c\varepsilon^{p_0} + O(\varepsilon^{p_0+1})\). The result for \(v_{\mathrm{asy}}(H_\mu)\) in 3 follows after rescaling \(H_\mu = \mu A_{1/\mu}\). The result for \(v_{\mathrm{asy}}(H_\mu, \delta_0)\) follows from the computations above and the observations \({\mathscr{F}}\delta_0 = e_0\) and \(\Pi_n = e_ne_n^* + O(\varepsilon)\), viz: \[\begin{align} \|G_i\delta_0\|^2 & = p_i^2 \left\| \left[ \int^\oplus_{{\mathbb{T}}^d} \sum_n \frac{\partial \eta_n(\varepsilon,e^{{\mathbf{i}}\theta})}{\partial \theta_i} \Pi_n(\varepsilon, e^{{\mathbf{i}}\theta}) \, \frac{{\mathrm{d}}\theta}{|{\mathbb{T}}^d|} \right] e_0 \right\|^2 \\ & = p_i^2 \int_{{\mathbb{T}}^d} \sum_n\left\| \frac{\partial \eta_n(\varepsilon,e^{{\mathbf{i}}\theta})}{\partial \theta_i} \Pi_n(\varepsilon, e^{{\mathbf{i}}\theta}) e_0 \, \right\|^2 \frac{{\mathrm{d}}\theta}{|{\mathbb{T}}^d|} \\ & = p_i^2 \int_{{\mathbb{T}}^d} \left| \frac{\partial \eta_0(\varepsilon,e^{{\mathbf{i}}\theta})}{\partial \theta_i} \right|^2(1 + O(\varepsilon))\, \frac{{\mathrm{d}}\theta}{|{\mathbb{T}}^d|} \\ & = p_i^2 \int_{{\mathbb{T}}} \left| 2c_{i,0} \sin\theta_i \varepsilon^{p_i} + O(\varepsilon^{p_i+1})\right|^2(1 + O(\varepsilon))\, \frac{{\mathrm{d}}\theta_i}{|{\mathbb{T}}|} \\ & = 4p_i^2c_{i,0}^2(\tfrac{1}{2}) \varepsilon^{2p_i} + O(\varepsilon^{2p_i+1}), \end{align}\] which shows \(v_{\mathrm{asy}}(A_\varepsilon,\delta_0) = c \varepsilon^{p_0} + O(\varepsilon^{p_0+1})\) after taking the maximum over \(i\). Rescaling to \(H_\mu = \mu A_{1/\mu}\) as before concludes the argument. ◻

Proof of Theorem 2. Fix \(\rho_0>0\), define \(\varepsilon_0\) as in Lemma 1, and let us consider \(A_\varepsilon = \varepsilon \Delta+V\) with \(|\varepsilon| \leq \varepsilon_0\). In order to bound (block) matrix elements, we use 18 and use the substitution \(z_j = \exp({\mathbf{i}}\theta_j)\), \({\mathrm{d}}z_j/z_j = {\mathbf{i}}\, {\mathrm{d}}\theta_j\) to transform to an iterated contour integral: \[\begin{align} \nonumber e^{-{{\mathbf{i}}}tA_\varepsilon}(x,y) & = \int_{{\mathbb{T}}^d} e^{-{{\mathbf{i}}} t A_\varepsilon(e^{ {{\mathbf{i}}} \theta})} e^{-{{\mathbf{i}}}\langle (x-y), \theta \rangle} \, \frac{{\mathrm{d}}\theta}{|{\mathbb{T}}^d|} \\ \label{eq:unitaryGroupIteratedContour} & = \frac{1}{{\mathbf{i}}^d|{\mathbb{T}}^d|} \oint_{|z_1|=1} \cdots \oint_{|z_d|=1} e^{-{\mathbf{i}}t A_\varepsilon(z)} z^{-(x-y+\mathbf{1})}\, {\mathrm{d}}z, \end{align}\tag{25}\] where we used the multi-index notation \(z^n = z_1^{n_1}\cdots z_d^{n_d}\) for \(z \in ({\mathbb{C}}^*)^d\), \(n \in {\mathbb{Z}}^d\), and \(\mathbf{1}=(1,1,\ldots,1)\).

Upper Bound. The first step is to deform the contours as follows: for each \(1 \le j \le d\), define \[\sigma_j = \mathop{\mathrm{sgn}}(x_j-y_j) = \begin{cases} -1 & x_j < y_j, \\ 0 & x_j = y_j, \\ 1 & x_j > y_j. \end{cases}\] We then shall pass to the contours \(|z_j|=e^{\sigma_j \rho_0}\). Since the integrand is holomorphic in \(\mathop{\mathrm{{\mathcal{A}}}}(2\rho_0)\), we have \[\begin{align} & \frac{1}{{\mathbf{i}}^d|{\mathbb{T}}^d|} \oint_{|z_1|=1} \cdots \oint_{|z_d|=1} e^{-{\mathbf{i}}t A_\varepsilon(z)} z^{-(x-y+\mathbf{1})}\, {\mathrm{d}}z \\ & \qquad = \frac{1}{{\mathbf{i}}^d|{\mathbb{T}}^d|} \oint_{|z_1| = e^{\sigma_1 \rho_0}} \cdots \oint_{|z_d|=e^{\sigma_d \rho_0}} e^{-{\mathbf{i}}t A_\varepsilon(z)} z^{-(x-y+\mathbf{1})}\, {\mathrm{d}}z. \end{align}\] Estimating, diagonalizing, and using Lemma 3.[item:QepsBound] to estimate \(\|Q_\varepsilon^{\pm 1}(z)\|\), we deduce for \(t \geq 0\): \[\begin{align} & \|e^{-{{\mathbf{i}}}tA_\varepsilon}(x,y)\| \\ & \qquad \leq \frac{1}{|{\mathbb{T}}^d|} \oint_{|z_1| = e^{\sigma_1 \rho_0}} \cdots \oint_{|z_d| = e^{\sigma_d \rho_0}} \|Q(z)\|\| \mathop{\mathrm{diag}}(\exp(-{{\mathbf{i}}}t \eta_n(\varepsilon, z)))\| \times \cdots \\ &\cdots \times\| [Q(z)]^{-1}\| |z|^{-(x-y+\mathbf{1})}\, |{\mathrm{d}}z|, \\ & \qquad\leq (1+C_2\varepsilon+O(\varepsilon^2))^2 \exp\left(t \max_{n,z} |\mathop{\mathrm{Im}}\eta_n(\varepsilon,z)| \right)\prod_{j=1}^d e^{-\rho_0|x_j-y_j|}. \end{align}\] Using Lemma 3.[item:imaginaryPart], we deduce \[\begin{align} \|e^{-{{\mathbf{i}}}tA_\varepsilon}(x,y)\|\leq C e^{C_1|t|\varepsilon^{p_0}} e^{-\rho_0|x-y|_1}, \end{align}\] which concludes the proof of the upper bound after rescaling to \(H_\mu = \mu A_{1/\mu}\) and using Remark 8 to pass from blocks to matrix elements.

Lower Bound. The lower bound follows readily from 4 and a straightforward generalization of [24] to higher dimensions (the formulation in [24] is for \(d=1\), but the proof applies in higher dimensions with cosmetic changes). ◻

References↩︎

[1]
F. Bloch. Über die Quantenmechanik der Elektronen in Kristallgittern. Zeitschrift für Physik, 52:555–600, 1929.
[2]
G. Floquet. Sur les équations différentielles linéaires à coefficients périodiques. Annales scientifiques de l’École Normale Supérieure, 12:47–88, 1883.
[3]
N. W. Ashcroft and N. D. Mermin. Solid State Physics. Brooks Cole, 1976.
[4]
M. Reed and B. Simon. Methods of Modern Mathematical Physics IV: Analysis of Operators. Academic Press, 1978.
[5]
M. S. P. Eastham. The Spectral Theory of Periodic Differential Equations. Scottish Academic Press, 1973.
[6]
J. Asch and A. Knauf. Motion in periodic potentials. Nonlinearity, 11(1):175–200, 1998.
[7]
A. Ahlbrecht, H. Vogts, A. H. Werner, and R. F. Werner. Asymptotic evolution of quantum walks with random coin. J. Math. Phys., 52(4):042201, 36, 2011.
[8]
A. Boutet de Monvel and M. Sabri. Ballistic transport in periodic and random media. In From complex analysis to operator theory—a panorama, volume 291 of Oper. Theory Adv. Appl., pages 163–216. Birkhäuser/Springer, Cham, 2023.
[9]
D. Damanik, J. Fillman, and D. C. Ong. Spreading estimates for quantum walks on the integer lattice via power-law bounds on transfer matrices. J. Math. Pures Appl. (9), 105(3):293–341, 2016.
[10]
D. Damanik, J. Fillman, and G. Young. Optimal dispersion for discrete periodic schrödinger operators. arxiv:2505.14475.
[11]
D. Damanik, M. Lukic, and W. Yessen. Quantum dynamics of periodic and limit-periodic Jacobi and block Jacobi matrices with applications to some quantum many body problems. Commun. Math. Phys., 337(3):1535–1561, 2015.
[12]
J. Fillman. Ballistic transport for periodic Jacobi operators on \(\mathbb {Z}^d\). In From operator theory to orthogonal polynomials, combinatorics, and number theory—a volume in honor of Lance Littlejohn’s 70th birthday, volume 285 of Oper. Theory Adv. Appl., pages 57–68. Birkhäuser/Springer, Cham, 2021.
[13]
A. Sagiv, R. Kassem, and M. I. Weinstein. Dispersive decay estimates for periodic jacobi operators on the half-line. arxiv:2505.14498.
[14]
P. Kuchment. The mathematics of photonic crystals. In Mathematical modeling in optical science, pages 207–272. SIAM, Philadelphia, PA, 2001.
[15]
J. D. Joannopoulos, S. G. Johnson, J. N. Winn, and R. D. Meade. Photonic crystals: molding the flow of light. Princeton University Press, Princeton, NJ, 2nd edition, 2008.
[16]
P. Kuchment. An overview of periodic elliptic operators. Bull. Amer. Math. Soc. (N.S.), 53(3):343–414, 2016.
[17]
D. Damanik, T. Malinovitch, and G. Young. What is ballistic transport? in press.
[18]
J. Fillman, W. Liu, and R. Matos. Irreducibility of the Bloch variety for finite-range Schrödinger operators. J. Funct. Anal., 283(10):Paper No. 109670, 22, 2022.
[19]
J. Fillman, W. Liu, and R. Matos. Algebraic properties of the Fermi variety for periodic graph operators. Journal of Functional Analysis, 286(4):110286, 2024.
[20]
L. Fisher, W. Li, and S. P. Shipman. Reducible Fermi surface for multi-layer quantum graphs including stacked graphene. Comm. Math. Phys., 385(3):1499–1534, 2021.
[21]
W. Li and S. P. Shipman. Irreducibility of the Fermi surface for planar periodic graph operators. Lett. Math. Phys., 110(9):2543–2572, 2020.
[22]
W. Liu. Irreducibility of the Fermi variety for discrete periodic Schrödinger operators and embedded eigenvalues. Geom. Funct. Anal., 32(1):1–30, 2022.
[23]
S. P. Shipman. Reducible Fermi surfaces for non-symmetric bilayer quantum-graph operators. J. Spectr. Theory, 10(1):33–72, 2020.
[24]
H. Abdul-Rahman, M. Darras, C. Fischbacher, and G. Stolz. Slow propagation velocities in Schrödinger operators with large periodic potential. Ann. Henri Poincaré, 202X.
[25]
H. Abdul-Rahman and G. Stolz. Exponentially decaying velocity bounds of quantum walks in periodic fields. Communications in Mathematical Physics, 403:1297–1327, 2023.
[26]
H. Abdul-Rahman, M. Cedzich, G. Stolz, and A. H. Werner. Exponential suppression of transport in periodic electric quantum walks and skew-shift cmv matrices.
[27]
E. H. Lieb and D. W. Robinson. The finite group velocity of quantum spin systems. Communications in Mathematical Physics, 28(3):251–257, 1972.
[28]
M. Aizenman and S. Warzel. Absolutely continuous spectrum implies ballistic transport for quantum particles in a random potential on tree graphs. Journal of Mathematical Physics, 53(9):095205, 06 2012.
[29]
J. Arbunich, J. Faupin, F. Pusateri, and I. M. Sigal. Maximal speed of quantum propagation for the hartree equation. Communications in Partial Differential Equations, 48(15):1–34, 2023.
[30]
S. Breteaux, J. Faupin, M. Lemm, D. H. Ou Yang, I. M. Sigal, and J. Zhang. Light cones for open quantum systems. arXiv:2303.08921.
[31]
C. Cedzich, A. Joye, A. H. Werner, and R. F. Werner. Exponential tail estimates for quantum lattice dynamics, 2024. arxiv:2408.02108.
[32]
M. C. Tran, A. Y. Guo, C. L. Baldwin, A. Ehrenberg, A. V. Gorshkov, and A. Lucas. ieb-Robinson light cone for power-law interactions. Physical Review Letters, 127(16):160401, 2021.
[33]
T. Kato. Perturbation Theory for Linear Operators. Classics in Mathematics. Springer, Berlin, Heidelberg, 1995.
[34]
M. Faust and I. Kachkovskiy. Absence of flat bands for discrete periodic graph operators with generic potentials, 2025.

  1. One could say that \(V\) is periodic if there exists some full-rank lattice \(\Gamma \subseteq {\mathbb{Z}}^d\) such that \(V(\cdot-\gamma)=V\) for every \(\gamma \in \Gamma\). However, it is not hard to check that any such \(\Gamma\) contains a subgroup of the form \(p {\mathbb{Z}}^d\) for suitable \(p\) and hence no generality is lost.↩︎

  2. It is convenient to use a general index set, such as the fundamental cell in the main application↩︎

  3. Here let us remark that, formally, the “\(m \neq n\)” on the summation is not required, since \(B_{nn}=0\) by assumption. However, it is a useful reminder, so we leave it in the notation.↩︎

  4. Let us point out this also follows from the Feynman–Hellmann theorem directly.↩︎