Transfer matrices and excitations with matrix product states

It is well known (and reiterated in section 4.1) that static correlation functions with respect to a uniform MPS ground state are obtained using the regular MPS transfer matrix (MPS-TM) [23], which is given by

$$\begin{eqnarray}&&\mathcal{T}_{A}^{{}}=\mathop{\sum }\limits_{s}{{\bar{A}}^{s}}\;\otimes \;{{A}^{s}},\end{eqnarray} \tag{ 3 }$$

with ${{\bar{A}}^{s}}$ being the complex conjugate of A^s. To simplify notation, we will generally omit the subscript A denoting the MPS matrix if it is not necessary.

For continuum results we define the (generator of a) uniform cMPS-TM

$$\begin{eqnarray}&&\mathcal{P}_{Q,R}^{{}}=\bar{Q}\;\otimes \;{\mathbb{1}}+{\mathbb{1}}\;\otimes \;Q+\bar{R}\;\otimes \;R,\end{eqnarray} \tag{ 4 }$$

where again $\bar{Q}$ and $\bar{R}$ denote the complex conjugates of Q and R, respectively. We will again omit subscripts whenever they are not necessary. The relation to the lattice TM is given by

$$\begin{eqnarray}&&\mathcal{P}_{{}}^{{}}=\mathop{{\rm lim} }\limits_{\epsilon \to 0}\frac{1}{\epsilon }{\rm log} \mathcal{T}\end{eqnarray} \tag{ 5 }$$

with being the lattice spacing of an underlying lattice discretization and $\mathcal{T}$ being the TM of the corresponding MPS defined on the discretized lattice.

Finally, for two-dimensional systems studied using PEPS, we work in the setting of infinitely long cylinders. By blocking the PEPS tensors, A^s_udlr, on a ring along the (finite) transversal y-direction of the cylinder, we can then interpret this contracted object as a uniform MPS along the (infinite) longitudinal x-direction of the cylinder, and we define the longitudinal TM as in (3). For a square lattice geometry, this MPS has bond dimension ${{D}^{{{N}_{y}}}}$ and physical dimension ${{d}^{{{N}_{y}}}}$, where N_y is the number of sites along the circumference of the cylinder. Equivalently, an elementary TM $\mathbb{E}$ can be constructed from the individual PEPS tensors, A^s_udlr, and the TM along a ring is obtained by contracting these elementary TMs along a ring. A graphical representation of the obtained TM is given in figure 1.

Zoom In Zoom Out Reset image size

If a uniform MPS defined by a set of matrices and A^s₁ is invariant under a local unitary symmetry operation u, one can show [32] that

$$\begin{eqnarray}&&A_{2}^{s}=\mathop{\sum }\limits_{k}\langle s|u|k\rangle A_{1}^{k}={{{\rm e}}^{{\rm i}}}\theta {{V}^{\dagger }}A_{1}^{s}V,\end{eqnarray} \tag{ 6 }$$

where A₂^s defines the transformed state. Here, V is a unitary gauge transformation on the auxiliary space and ${{{\rm e}}^{{\rm i}\theta }}$ is the dominant eigenvalue with magnitude one of the mixed TM

$$\begin{eqnarray}&&\mathcal{T}_{{{A}_{1}}}^{{{A}_{2}}}=\mathop{\sum }\limits_{s}\bar{A}_{2}^{s}\;\otimes \;A_{1}^{s}.\end{eqnarray} \tag{ 7 }$$

In fact, the MPS is invariant under the local symmetry, u, if and only if the spectral radius of the mixed TM $\rho (\mathcal{T}_{{{A}_{1}}}^{{{A}_{2}}})$ is one (i.e., the fidelity per lattice site is one). In the case of higher-dimensional lattice systems, a relation similar to (6) holds for PEPS [32].

For phases with a spontaneously broken symmetry, the ground state is degenerate and the variationally best ground-state approximations within the manifold of MPS of bond dimension D are minimally entangled states that exhibit maximal symmetry, breaking. Such states can be transformed into each other by applying the symmetry the operations of the broken symmetry.

In the ground state of a one-dimensional quantum system, continuous symmetries for which the order parameter does not commute with the Hamiltonian cannot be spontaneously broken [33, 34]. Nevertheless, close to or within a gapless phase with a continuous symmetry, it is sometimes energetically favorable for a variational cMPS approximation of the ground state to break this symmetry and to approximate an excited state with very small excitation energy and much smaller entanglement. This pseudo-symmetry-breaking is purely an effect of the finite bond dimension and also gives rise to a pseudo-order parameter [35, 36]. The symmetry is restored in the limit $D\to \infty$.

In a phase with broken symmetry on a lattice, let A^s₁ and A^s₂ be MPS approximations of two ground states with the same variational energy but different order parameters and maximally broken symmetry. The orthogonality of these states requires that the fidelity per lattice site must be strictly smaller than one (i.e., the spectral radius of the mixed TM, $\rho (\mathcal{T}_{{{A}_{1}}}^{{{A}_{2}}})\lt 1$).

Equivalently, for continuum systems we define the (generator of the) mixed cMPS-TM as

$$\begin{eqnarray}&&\mathcal{P}_{{{Q}_{1}},{{R}_{1}}}^{{{Q}_{2}},{{R}_{2}}}={{\bar{Q}}_{2}}\;\otimes \;{\mathbb{1}}+{\mathbb{1}}\;\otimes \;{{Q}_{1}}+{{\bar{R}}_{2}}\;\otimes \;{{R}_{1}}\end{eqnarray} \tag{ 8 }$$

where ${{Q}_{1}},{{R}_{1}}$ and ${{Q}_{2}},{{R}_{2}}$ are two different cMPS representations. Similar to the lattice case, if ${{Q}_{1}},{{R}_{1}}$ and ${{Q}_{2}},{{R}_{2}}$ describe two equally good ground-state approximations with the same variational energy, but with different order parameters and maximum symmetry-breaking, the spectrum $\mathcal{P}_{{{Q}_{1}},{{R}_{1}}}^{{{Q}_{2}},{{R}_{2}}}$ has strictly negative real parts.

The degeneracy of the ground state in phases with broken symmetries gives rise to topologically nontrivial excitations (kinks or domain walls), which typically correspond to the elementary excitations of the model. The mixed cMPS-TM of type (7) or (8) of these symmetry broken ground states plays a crucial role in obtaining a variational approximation for such excitations, whereas the regular cMPS-TM of type (3) or (4) is the central object for topologically trivial excitations [36, 37].

We will first focus on three prototypical one-dimensional lattice models. We start with the spin-1/2 XY model in an external magnetic field

$$\begin{eqnarray}&&{{H}_{{\rm XY}}}=-\mathop{\sum }\limits_{j}(1+\gamma )S_{j}^{x}S_{j+1}^{x}+(1-\gamma )S_{j}^{y}S_{j+1}^{y}+gS_{j}^{z},\end{eqnarray} \tag{ 10 }$$

which can be solved exactly [38–40]. Here, $S_{j}^{\alpha }$ denote spin-1/2 operators defined on site j. We consider the gapped ferromagnetic regime, $0\lt \gamma \lt 1$ and $0\lt g\lt 1$, where the system is in a symmetry-broken phase and the ground state is twofold degenerate with local order parameter ${{m}_{x}}=\langle S_{j}^{x}\rangle$. Here the elementary excitations are domain-wall-like and therefore well approximated by a topologically nontrivial MPS ansatz. Specifically, we consider the incommensurate phase ${{\gamma }^{2}}+{{g}^{2}}\lt 1$, where correlations oscillate with arbitrary wave vectors.

As a second example we consider the spin-1/2 XXZ model in an external magnetic field

$$\begin{eqnarray}&&{{H}_{{\rm XXZ}}}=-\mathop{\sum }\limits_{j}S_{j}^{x}S_{j+1}^{x}+S_{j}^{y}S_{j+1}^{y}+\Delta S_{j}^{z}S_{j+1}^{z}+hS_{j}^{z}.\end{eqnarray} \tag{ 11 }$$

This model is solvable as well, and the ground state and elementary excitations in the thermodynamic limit can be obtained via Bethe ansatz [41–43]. Here we consider the antiferromagnetic gapless incommensurate phase specified by $-1\lt \Delta \lt 0$ and $0\lt |h|\lt 1-\Delta$, where there are gapless excitations at multiples of the Fermi momentum ${{k}_{{\rm F}}}=(\frac{1}{2}-{{m}_{z}})\pi$ with ${{m}_{z}}=\langle S_{j}^{z}\rangle$ as the ground-state magnetization. In this phase there is no spontaneous symmetry-breaking, but due to criticality the finite D MPS ground-state approximation breaks the continuous rotational symmetry in the XY plane (see section 2.2). This makes it possible to use a topologically nontrivial variational MPS ansatz for the excitations.

As a third example we study the (extended) Fermi-Hubbard model

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} {{H}_{{\rm HUB}}} & = & -\mathop{\sum }\limits_{\sigma j}c_{\sigma j}^{\dagger }{{c}_{\sigma j+1}}-{{c}_{\sigma j}}c_{\sigma j+1}^{\dagger }+V{{n}_{j}}{{n}_{j+1}} \\ {} & {} & +\mathop{\sum }\limits_{j}U\left( {{n}_{\uparrow j}}-\frac{1}{2} \right)\left( {{n}_{\downarrow j}}-\frac{1}{2} \right)-\mu {{n}_{j}}, \\ \end{array}\end{eqnarray} \tag{ 12 }$$

where $c_{\sigma }^{\dagger },{{c}_{\sigma }}$ denote the creation and annihilation operators of fermions with spin σ, ${{n}_{\sigma }}=c_{\sigma }^{\dagger }{{c}_{\sigma }}$ and $n={{n}_{\uparrow }}+{{n}_{\downarrow }}$. For $V\ne 0$, this model is non integrable. We consider the repulsive regime where $U,V\gt 0$, away from half filling ($\mu \ne 0$), which again corresponds to a gapless incommensurate phase. There is no spontaneous symmetry-breaking in this phase and we consider topologically trivial excitations only.

The last example studied is the Kondo lattice model (KLM) [44]

$$\begin{eqnarray}&&{{H}_{{\rm KLM}}}=-\mathop{\sum }\limits_{\sigma j}c_{\sigma j}^{\dagger }{{c}_{\sigma j+1}}-{{c}_{\sigma j}}c_{\sigma j+1}^{\dagger }-\mu {{n}_{j}}+J{\boldsymbol{S}} _{j}^{c}\cdot {\boldsymbol{S}} _{j}^{d}\end{eqnarray} \tag{ 13 }$$

where $c_{\sigma }^{\dagger },{{c}_{\sigma }}$ denote the creation and annihilation operators of conduction electrons (c) with spin σ, $n=c_{\uparrow }^{\dagger }{{c}_{\uparrow }}+c_{\downarrow }^{\dagger }{{c}_{\downarrow }}$; ${\boldsymbol{S}} _{j}^{c}$ and ${\boldsymbol{S}} _{j}^{d}$ are the spin operators for conduction electrons (c) and localized electrons (d), respectively. We consider the paramagnetic metallic phase away from half filling ($\mu \ne 0$). This model in this phase has also been studied at finite temperatures with TMRG techniques in [46].

In figures 2 and 3 we show results for (1) the XY model at $\gamma =0.3$, g = 0.2, and D = 40, (2) the XXZ model at $\Delta =-0.5$, h = 1, and $D=100$, (3) the extended Hubbard model at U = 5, V = 1, $\mu =2$, and D = 100, and (4) the KLM at J = 2, $\mu =-1,$ and D = 120. On the left we plot the eigenvalues, ${{\lambda }_{j}}={{{\rm e}}^{-{{\varepsilon }_{j}}+{\rm i}{{\phi }_{j}}}}$, of the regular MPS-TM on the complex plane within the unit circle, whereas on the right we plot ${{\varepsilon }_{j}}=-{\rm log} |{{\lambda }_{j}}|$ versus the complex argument ${{\phi }_{j}}$, along with the lowest variational excitation energies obtained from a topologically trivial variational MPS ansatz [37]. We do not plot the dominant eigenvalue, ${{\lambda }_{0}}=1$. In figure 4 we show results for the mixed MPS-TM and topologically nontrivial excitations for the XY model and XXZ model only, with the same parameters as above. If available, we also plot the exact dispersion of the elementary excitations, as well as the lower boundaries of multiparticle continua, for reference.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

From figure 4 it is apparent that the topologically nontrivial variational ansatz captures elementary excitations with high accuracy. The two-particle continuum consists of combinations of two elementary excitations, and is thus partially captured by a topologically trivial ansatz (see figures 2 and 3), which is consistent with the results in [36]. The full continuum can be recovered by using a variational MPS ansatz, including the scattering states of elementary excitations [45]. Consequently, low-lying states with higher odd particle numbers are partially captured by a topologically nontrivial ansatz, whereas low-lying states with even particle numbers are partially captured by a topologically trivial ansatz. In the case where there are no topologically nontrivial excitations, there is no such distinction.

Concerning the eigenvalues of the MPS-TM, we can now make the following remarkable observations. In the plots on the left of figures 2–4, we see that most of the eigenvalues arrange themselves along several lines with an approximately constant complex argument, ${{\phi }_{j}}={{\phi }_{\alpha }}$, where α labels distinct lines. This fact is reflected in the arrangement of the ${{\varepsilon }_{j}}$ in columns in the plots on the right, where low-lying ${{\varepsilon }_{j}}$ correspond to eigenvalues ${{\lambda }_{j}}$ close to the unit circle. Denote the lowest-lying ${{\varepsilon }_{j}}$ for each line, α, as ${{\varepsilon }_{\alpha }}$.

We can also observe that the complex arguments, ${{\phi }_{\alpha }},$ of these ${{\varepsilon }_{\alpha }}$ coincide very precisely with the momenta, ${{k}_{{\rm min} }},$ of the minima in the dispersion of excitation energies. This fact has been exploited in [46] to locate the Fermi momentum, ${{k}_{{\rm F}}},$ in the KLM. While there was no justification given there as to how the phase, $\phi ,$ of the second-largest eigenvalue of the TM is related to low-energy excitations, we contribute several arguments for this fact in section 4. Connections between peaks in the static structure factor and the logarithm of the absolute value and the complex argument of TM eigenvalues have also been exploited in [14, 47] to study the temperature-dependence of static correlation functions.

The values of the lowest-lying ${{\varepsilon }_{\alpha }}$ are related to minima, $E_{{\rm min} }^{\alpha },$ in the dispersion of elementary excitation energies via some characteristic velocities, ${{v}_{\alpha }}$ ¹ , and can serve as a first approximation for this energy if the velocity can be estimated. It appears that these velocities, which determine the energy scale for each line, α, can also vary between different excitation minima, $E_{{\rm min} }^{\alpha }$, within each respective model in the examples shown.

Indeed, for the XY model, the momenta, ${{k}_{{\rm min} }},$ of the minima of the elementary excitations and the three-particle continuum are well reproduced by the eigenvalues of the mixed MPS-TM with the largest magnitude, with deviations of ${\rm O}({{10}^{-3}})$. Consequently the same holds for the minima of the two-particle continuum and the regular MPS-TM. For the elementary excitations we estimate the characteristic velocity, ${{v}_{1}},$ relating the lowest excitation energy, ${{E}_{{\rm min} }}$ and ${{\varepsilon }_{1}}=-{\rm log} |{{\lambda }_{1}}|$, where ${{\lambda }_{1}}$ is the eigenvalue with, second-largest magnitude, as ${{v}_{1}}={{E}_{{\rm min} }}/{{\varepsilon }_{1}}\approx 0.9409$ (see figure 4 inset), where we have extrapolated the value of ${{\varepsilon }_{1}}$ for $D\to \infty$. Towards the end of section 4.3 we show how the value of this velocity can be estimated from assuming a Lorentz-invariant low-energy behavior. There, we obtain an estimate for v₁ that agrees with the value obtained above within 1% accuracy.

For the XXZ model, the momenta of the gapless excitations at multiples of ${{k}_{{\rm F}}}$ are even more precisely reproduced by the arguments, ${{\phi }_{j}},$ of the eigenvalues of the regular and mixed MPS-TM with magnitudes close to one (i.e., ${{\varepsilon }_{j}}$ close to zero), with relative deviations of ${\rm O}({{10}^{-6}})$. Notice that in the limit $D\to \infty$, we expect ${{\varepsilon }_{j}}\to 0$; that is, the spectral radius of the mixed MPS-TM also becomes one and the rotational symmetry in the XY plane is restored (see section 2.2).

For the extended Hubbard model, the star-like structure of the eigenvalue spectrum in figure 3 is very pronounced. In the right plot, the ratios of the variational dispersion minima and the lowest ${{\varepsilon }_{j}}$ at $k={{\phi }_{j}}\approx 0.4\pi$ and $k={{\phi }_{j}}\approx \pi$ might suggest a characteristic velocity, ${{v}_{\alpha }}\gt 1$. However, one would expect to have gapless excitations at these momenta, as well as at $k\approx 0.2\pi$ suggesting that the corresponding variational energies are not converged. It is instructive to either use a larger bond dimension or to enhance by using an ansatz including scattering, which we have not performed here. It is, however, interesting to note for these minima that the accuracy of the variational energies and the low-lying ${{\varepsilon }_{j}}$ appears to be roughly on the same level.

At last, for the KLM we obtain a variational low-energy dispersion that exhibits overall low variation in k, which is typical for the heavy fermion regime. We also observe dents in the slowly varying part coming from gapless excitations at multiples of the Fermi momentum, for which we obtain an estimate of ${{k}_{{\rm F}}}\approx 0.41541(3)\pi$ by extrapolating $D\to \infty$. While in [46], ${{k}_{{\rm F}}}$ is estimated by approximating the finite temperature QTM of the quantum system using TMRG. Subsequently extrapolating $T\to 0$ –which is inherently prone to technical difficulties—we directly work at T = 0 and construct an approximation of the zero temperature QTM from an MPS ground-state approximation. Indeed, our estimate of ${{k}_{{\rm F}}}$ at T = 0 is of the same order as the value in [46] for lowest T; however, there are no claims about whether the value at T = 0 was possible due to technical difficulties.

The above observations are truly remarkable. The mere knowledge of the ground state MPS-TM already yields important information about the excitation spectrum of the underlying Hamiltonian. More specifically, the momenta, ${{k}_{{\rm min} }},$ of the excitation energy minima in the respective particle sectors can be determined accurately, and the corresponding energies can be estimated in a first approximation just from static ground-state properties. An advantage of this TM based approach over just considering static correlation functions is discussed in section 4.2.

We will now turn to continuous (1+1)-dimensional field theories and study the Lieb-Liniger model [48] using cMPS methods. The Hamiltonian is given by

$$\begin{eqnarray}&&{{H}_{{\rm LL}}}=\int _{-\infty }^{\infty }{\rm d}x\;{{\partial }_{x}}{{\psi }^{\dagger }}{{\partial }_{x}}\psi +c{{\psi }^{\dagger }}{{\psi }^{\dagger }}\psi \psi -\mu {{\psi }^{\dagger }}\psi ,\end{eqnarray} \tag{ 14 }$$

with repulsive interaction strength $c\gt 0$ and chemical potential $\mu \gt 0$, where ψ and ${{\psi }^{\dagger }}$ are bosonic field operators. The model depends only on a single parameter, $\gamma =\frac{c}{\rho }$, with $\rho =\langle {{\psi }^{\dagger }}\psi \rangle$ being the ground state particle density, which is critical for all values of γ.

In figures 5 and 6 we show results for the eigenvalues of the regular and mixed cMPS-TM, similar to the lattice case for (a) $\gamma \approx 1.35$ and (b) $\gamma \approx 311.5$ and D = 64. Given the relation in (5), the right column now plots the real part of the eigenvalues, $-{{\sigma }_{j}}$ of the generator $\mathcal{P}_{{}}^{{}},$ versus their imaginary part, which is now interpreted as momentum. In the left column the eigenvalues, ${{\lambda }_{j}},$ of $\mathcal{T}_{{}}^{{}}={\rm exp} (\mathcal{P}_{{}}^{{}})$ are plotted on the complex plane within the unit circle, as in the lattice case. In the continuum setting, momentum space is no longer $2\pi$-periodic, and the definition of $\mathcal{T}_{{}}^{{}}$ is not fully justified, as it can come with any real power, $x\geqslant 0,$ in continuum correlation functions, where as only integer powers appear in lattice correlation functions. Nevertheless, it helps to illustrate that the spectrum of eigenvalues of the TM exhibits a similar structure. The fact that this structure in figures 5 and 6 is less outspoken than for some of the lattice models indicates a larger contribution of microscopic effects for this specific case.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Next we will study the simplest Lorentz-invariant theory available: the free (1+1)-dimensional Klein-Gordon boson described by the Hamiltonian

$$\begin{eqnarray}&&{{H}_{{\rm KG}}}=\frac{1}{2}\int _{-\infty }^{\infty }{\rm d}x\left[ {{\pi }^{2}}+{{\left( {{\partial }_{x}}\phi \right)}^{2}}+{{m}^{2}}{{\phi }^{2}} \right],\end{eqnarray} \tag{ 15 }$$

where we have taken the speed of light to be unity. The field operators, ϕ and $\pi ,$ can be written in terms of the cMPS Fock space operators, ψ and ${{\psi }^{\dagger }},$ as:

$$\begin{eqnarray}&&\phi =\frac{1}{\sqrt{2\nu }}(\psi +{{\psi }^{\dagger }})\qquad \pi =-\frac{{\rm i}}{2}\sqrt{2\nu }(\psi -{{\psi }^{\dagger }}),\end{eqnarray} \tag{ 16 }$$

where an arbitrary scale, $\nu ,$ is introduced [49]. The Hamiltonian (15) needs to be regularized, and this is achieved by adding the term ${{\left( {{\partial }_{x}}\pi \right)}^{2}}{{\nu }^{-2}}$ to the Hamiltonian.

We make the following observations regarding the eigenspectrum of the generator of the cMPS-TM, $\mathcal{P}_{{}}^{{}},$ corresponding to ${{H}_{{\rm KG}}}$, as plotted in figure 7. The eigenvalues, ${{\sigma }_{j}},$ of $\mathcal{P}_{{}}^{{}}$ are all real and negative for all values of m and D. This reflects the fact that the relativistic dispersion relation has a single minimum at momentum zero. The eigenvalue with the largest real part of $\mathcal{P}_{{}}^{{}}$ converges to $-m$ as $D\to \infty$. As is argued in section 5.1, it follows from the Euclidean invariance of the quantum transfer matrix of a relativistic theory that in the limit $D\to \infty ,$ the eigenspectrum of $\mathcal{P}_{{}}^{{}}$ should be the same as that of ${{H}_{{\rm KG}}}$ (up to the sign), corresponding to a characteristic speed, $v=c=1$. The above observation provides numerical support for this by using and extrapolating from finite D data.

Zoom In Zoom Out Reset image size

We can also study the distribution of the eigenvalues of $\mathcal{P}_{{}}^{{}}$ as a function of the bond dimension. For any value of m, the (negative) eigenvalues, ${{\sigma }_{j}},$ become dense in the region $[m,+\infty )$. In the gapped phase ($m\gt 0$), the density of eigenvalues diverges at m (i.e., the ratio of the nth largest and second-largest eigenvalue of $\mathcal{P}_{{}}^{{}}$ converges to unity for low-lying n). This is similar to the density of states in a gapped single particle excitation branch near the minimum of the dispersion relation. As the mass, $m,$ is taken to zero, the theory becomes critical, and the cMPS approximation enters the so-called 'finite entanglement regime' [50]. The eigenvalues of $\mathcal{P}_{{}}^{{}}$ are still all real and converge to zero as $D\to \infty$, but their ratios now converge to values larger than one, and are expected to encode universal data [49]. This is tantamount to the statement that the effect of the finite bond dimension is only to introduce a single scale into the underlying conformal field theory, and implies that universal quantities can be extracted straightforwardly from cMPS data [49].

Let us exemplify this by spelling out the results obtained by taking m = 0.2 in (15) and using the modest range of bond dimensions up to D = 36. Scaling with $1/D$ and extrapolating to $D\to \infty$ yields that the second-largest eigenvalue of $\mathcal{P}_{{}}^{{}}$ tends to 0.201, thus reproducing the mass accurately. The ratio of the third and second eigenvalue of $\mathcal{P}_{{}}^{{}}$ is estimated to converge to 1.040 as $D\to \infty$, and the ratio of the fourth and second is estimated to converge to 1.086. We note that, since the theory is free, the deviation from unity cannot be due to convergence to some bound state just above the lowest branch, and can only be an effect of numerical accuracy. The same ratios for m = 0 converge to approximately 2.51 and 3.1, respectively, and are related to properties of the excitation spectrum of the underlying conformal field theory [49].

The observed connection between the eigenvalues of the TM and the minimum of the dispersion relation opens up a way to infer properties of the dispersion relation of two-dimensional systems, which are notoriously difficult to deal with. To this end, given a two-dimensional translation-invariant PEPS [28] on a square lattice cylinder with periodic boundary conditions in the y direction, we block all PEPS tensors in a ring around the cylinder (i.e., with the same x coordinate along the cylinder). We then consider the quasi-one-dimensional system along the cylinder obtained that way, described by blocked tensors A^s_lr, and its TM, $\mathcal{T}_{{}}^{{}}$.

As the original state was also translational invariant in the y direction, we can label the eigenvectors of $\mathcal{T}_{{}}^{{}}$ by eigenvalues ${{{\rm e}}^{{\rm i}{{k}_{y}}}}$ of the action of the translation operator on the auxiliary degrees of freedom, as given by (6). On a hexagonal lattice cylinder, we additionally block two neighboring PEPS tensors to obtain a quasi square lattice before further blocking all obtained tensors in a ring around the cylinder. Assuming that the observed connection between the leading eigenvalues of the TM and the minimum of the dispersion relation holds for each k_y independently, we obtain the location and relative energy of the minima of the dispersion relation for each possible value of k_y, which yields a cut through the dispersion relation. By closing the periodic boundaries in different ways, we can obtain this information along different symmetry axes, allowing us to reconstruct the overall form of the dispersion relation.

We now apply this strategy to the Affleck–Kennedy–Lieb–Tasaki (AKLT) model [51] on the square and hexagonal lattice. Its ground state is constructed by placing spin-$\frac{1}{2}$ singlets on the edges of the lattice and projecting the spin-$\frac{1}{2}$'s at each vertex onto the sector with maximal spin (${{S}_{{\rm phys}}}=2$ for the square lattice and ${{S}_{{\rm phys}}}=\frac{3}{2}$ for the hexagonal lattice); this construction yields the unique ground state of the ${\rm SU}(2)$ -invariant Hamiltonian $H={{\sum }_{\langle i,j\rangle }}{{\Pi }_{i,j}}$, where ${{\Pi }_{i,j}}$ is the projector onto the subspace with spin $S={{S}_{{\rm phys}}}$ on neighboring sites i and j. This construction corresponds exactly to a PEPS, with bond dimension D = 2. Even though the exact ground state is known, little is known about the excited states, although recently an indirect method was proposed to estimate the gap by means of a tensor network renormalization group method [52].

We can use an iterative eigensolver to exactly determine the low-lying spectrum of $\mathcal{T}_{{}}^{{}}$ on cylinders with sufficiently large circumference. Since the model possesses ${\rm SU}(2)$ symmetry, we can additionally label the eigenvectors of the TM by their spin (which corresponds to the spin of the excitation), thereby aiding the identification of different excitation branches. As the TM of the AKLT model is Hermitian (up to a gauge transformation), its eigenvalues are real, and thus ${{k}_{x}}=0,\pi$. The pairs $({{k}_{x}},{{k}_{y}})$ for all eigenvalues of $\mathcal{T}_{{}}^{{}}$ are therefore arranged along lines in the Brillouin zone; by properly closing the periodic boundaries, we can thus obtain data points along different symmetry axes. The results for the square and hexagonal lattice are shown in figure 8. In both cases, we find an isolated branch of antiferromagnetic spin-1 excitations, with a two-particle continuum starting at about twice the elementary quasi energy gap, in agreement with known results for one-dimensional systems. In particular, for the square lattice we find the minimum of the dispersion at momentum $({{k}_{x}},{{k}_{y}})=(\pi ,\pi )$, whereas for the hexagonal lattice the minimum is found at $({{k}_{x}},{{k}_{y}})=(0,0)$. For both lattices the minima appear on the isolated S = 1 branches.

Zoom In Zoom Out Reset image size

This approach has further been used by some of the authors to study anyon condensation in the toric code model with string tension [53].

Let us first recall how the regular MPS-TM gives access to all static correlation functions in the corresponding MPS. For this we assume a complete eigendecomposition of the TM

$$\begin{eqnarray}&&\mathcal{T}_{{}}^{{}}=\mathop{\sum }\limits_{j=0}^{{{D}^{2}}}{{\lambda }_{j}}|j)(j|,\end{eqnarray} \tag{ 18 }$$

where we have dropped the subscript A for notational simplicity. Here $|j)$ and $(j|$ are the right and left eigenvectors of $\mathcal{T}_{{}}^{{}},$ respectively, with $(i|j)={{\delta }_{ij}}$; we again write ${{\lambda }_{j}}={{{\rm e}}^{-{{\varepsilon }_{j}}+{\rm i}{{\phi }_{j}}}}$ for the eigenvalues. We demand the ground state MPS representation to be injective, such that there is a unique dominant eigenvalue, ${{\lambda }_{0}}=1$ and $|{{\lambda }_{j\gt 0}}|\lt 1,$ for all other eigenvalues. Furthermore, we define the operator, transfer matrix (OTM) for operators, $O,$ acting on n sites as

$$\begin{eqnarray}&&{{\mathcal{J}}_{O}}=\mathop{\sum }\limits_{\frac{{{i}_{1}}\ldots {{i}_{n}}}{{{j}_{1}}\ldots {{j}_{n}}}}O_{{{j}_{1}}\ldots {{j}_{n}}}^{{{i}_{1}}\ldots {{i}_{n}}}\left( {{{\bar{A}}}^{{{i}_{1}}}}\ldots {{{\bar{A}}}^{{{i}_{n}}}} \right)\;\otimes \;\left( {{A}^{{{j}_{1}}}}\ldots {{A}^{{{j}_{n}}}} \right),\end{eqnarray} \tag{ 19 }$$

where $O_{{{j}_{1}}\ldots {{j}_{n}}}^{{{i}_{1}}\ldots {{i}_{n}}}=\langle {{i}_{1}}\ldots {{i}_{n}}|O|{{j}_{1}}\ldots {{j}_{n}}\rangle$.

With the above definitions, it is well known that (17) can be written as

$$\begin{eqnarray}&&{{C}_{AB}}(n+1)=\mathop{\sum }\limits_{j\gt 0}f_{j}^{{\rm AB}}\ \lambda _{j}^{n}=\mathop{\sum }\limits_{j\gt 0}f_{j}^{{\rm AB}}\ {{{\rm e}}^{-{{\varepsilon }_{j}}n}}{{{\rm e}}^{{\rm i}{{\phi }_{j}}n}}\end{eqnarray} \tag{ 20 }$$

where we have defined the form factors

$$\begin{eqnarray}&&f_{j}^{{\rm AB}}=(0|{{\mathcal{J}}_{A}}|j)(j|{{\mathcal{J}}_{B}}|0).\end{eqnarray} \tag{ 21 }$$

As for a finite bond dimension, (20) is a finite sum of exponentials; connected correlation functions for sufficiently large distances, $n,$ must decay as a pure exponential, $\| {{C}_{AB}}(n)\| \sim {\rm exp} (-n/\xi ),$ where the correlation length, $\xi ,$ corresponds to the inverse of the smallest ${{\varepsilon }_{j}}$ with a nonzero form factor.

In contrast, typical correlation functions in gapped phases are expected to decay at large distances, $n,$ as

$$\begin{eqnarray}&&\| {{C}_{AB}}(n)\| \sim {{n}^{-\eta }}{{{\rm e}}^{-\frac{n}{\xi }}},\end{eqnarray} \tag{ 22 }$$

where there is an additional power law contribution to the decay with an exponent, η, which in principle depends on operators A and B. Close to a critical point, this form can be motivated from conformal field theory or from a general renormalization group argument. Approaching the critical point takes the correlation length, $\xi \to \infty$, and a pure power law decay of correlations remains, where the scaling exponent, $\eta ,$ can depend on the choice of operators A and B. Sufficiently deep in a gapped phase, on the other hand, (22) is known as the Ornstein-Zernike form, and η typically depends on the number of spatial dimensions, $d,$ as $\eta =d/2$ [54, 55] (i.e., for a one-dimensional quantum system

$$\begin{eqnarray}&&\| {{C}_{AB}}(n)\| \sim \frac{1}{\sqrt{n}}{{{\rm e}}^{-n/\xi }}\end{eqnarray} \tag{ 23 }$$

for large distances, n).

In the limit $n\to \infty$, an MPS with a finite bond dimension would correspond to $\eta =0$. Nevertheless, for $D\to \infty$ the correct form of the correlation functions should be restored. We proceed by noting that the scaling form in (22) can be obtained from (20) under the following assumptions:

• (i)  
  The leading eigenvalues, ${{\lambda }_{j}},$ arrange themselves on a finite number of lines in the complex plane with constant complex argument ${{\phi }_{\alpha }}$, where we label these lines with index α. Consider also the form factors $f_{j}^{{\rm AB}},$ as defined in (21), along such a line. We define the subset of indices $\left\{ {{j}_{\alpha }}|{{\phi }_{j}}={{\phi }_{\alpha }}\ {\rm and}\ f_{j}^{{\rm AB}}\ne 0 \right\}$ (i.e., indices ${{j}_{\alpha }}$ label all eigenvalues with constant complex argument ${{\phi }_{\alpha }},$ for which the form factors, ${{f}_{j}},$ are non zero).
• (ii)  
  On each of these lines, the ${{\lambda }_{{{j}_{\alpha }}}}$ become sufficiently dense for $D\gg 1,$ and the corresponding ${{\varepsilon }_{{{j}_{\alpha }}}}$ will follow some dispersion. We will then reorder the indices ${{j}_{\alpha }}$, such that the ${{\varepsilon }_{{{j}_{\alpha }}}}$ are in ascending order and ${{\Delta }_{\alpha }}={{\varepsilon }_{{{j}_{\alpha }}=0}}$ is the minimum. In the complex plane this parameterization corresponds to going from the ${{\lambda }_{{{j}_{\alpha }}}}$ closest to the unit circle along line α towards the center. For small ${{j}_{\alpha }}$, ${{\varepsilon }_{{{j}_{\alpha }}}}$ changes smoothly as ${{\varepsilon }_{{{j}_{\alpha }}}}={{\Delta }_{\alpha }}+g\;j_{\alpha }^{\kappa }$ to the leading order, with some constants, $g,\kappa \gt 0$, possibly depending on α.
• (iii)  
  The form factors $f_{{{j}_{\alpha }}}^{{\rm AB}}$ also vary smoothly and follow, to the leading order, some power law $f_{{{j}_{\alpha }}}^{{\rm AB}}\sim j_{\alpha }^{\rho },$ where the exponent ρ depends on the operators A and B and possibly α. For $\rho =0,$ the leading order corresponds to a non zero constant.

With the above assumptions and using the Euler-McLaurin formula to approximate the discrete sum with an integral, the contribution of one line, $\alpha ,$ of eigenvalues in (20) becomes approximately

$$\begin{eqnarray}&&{{{\rm e}}^{{\rm i}{{\phi }_{\alpha }}n}}{{{\rm e}}^{-{{\Delta }_{\alpha }}n}}\int _{0}^{\infty }{\rm d}z\;f_{\alpha }^{{\rm AB}}(z){{{\rm e}}^{-ng{{z}^{\kappa }}}}\end{eqnarray} \tag{ 24 }$$

where we have replaced ${{j}_{\alpha }}$ with a continuous parameter, $z,$ and integrate to $\infty$ for convenience, as ${{{\rm e}}^{-ng{{z}^{\kappa }}}}$ decays sufficiently fast with increasing z, even for moderate n.

Using the saddle point approximation allows one to replace $f_{\alpha }^{{\rm AB}}(z)$ by its dominant behavior near z = 0, and we obtain for n sufficiently large

$$\begin{eqnarray}&&{{{\rm e}}^{{\rm i}{{\phi }_{\alpha }}n}}{{{\rm e}}^{-{{\Delta }_{\alpha }}n}}\int _{0}^{\infty }{\rm d}z\;{{z}^{\rho }}{{{\rm e}}^{-ng{{z}^{\kappa }}}}\sim {{{\rm e}}^{{\rm i}{{\phi }_{\alpha }}n}}{{n}^{-\frac{1+\rho }{\kappa }}}{{{\rm e}}^{-{{\Delta }_{\alpha }}n}}\end{eqnarray} \tag{ 25 }$$

and we have recovered (22) with $\eta =\frac{1+\rho }{\kappa }$. Deep in a gapped phase, we can reasonably expect that $\kappa =2,$ and with $\rho =0,$ we then recover the typical Ornstein-Zernike form for one-dimensional quantum systems with $\eta =1/2$. The correlation length, $\xi ,$ is then given by $\frac{1}{\xi }={{\varepsilon }_{1}}={{\Delta }_{1}},$ with ${{\varepsilon }_{1}}$ being the smallest non zero ${{\varepsilon }_{j}}$.

In figure 9 we show corresponding numerical evidence for the XY model, defined in (10), at $\gamma =0.5$ and g = 1.05 (i.e., in the gapped paramagnetic phase, with bond dimension D = 32). In this phase all the eigenvalues of the TM are real (i.e., ${{\phi }_{j}}=0$). We plot ${{\varepsilon }_{j}}$ along with a quadratic fit, $\varepsilon (j)=\Delta +a\;{{j}^{2}}$, as well as the form factors $f_{j}^{{\rm XX}}$, where X stands for S^x, versus index j for $j\lt 7$. Again, j labels eigenvalues for which the form factors, $f_{j}^{{\rm XX}},$ are non zero, in ascending order. The data confirms the expected values of the exponents as $\kappa =2$ and $\rho =0$, yielding the Ornstein-Zernike behavior of the correlation function

$$\begin{eqnarray}&&{{C}_{{\rm XX}}}(n)\sim \frac{{{{\rm e}}^{-\Delta n}}}{\sqrt{n}}\end{eqnarray} \tag{ 26 }$$

expected for this model in this parameter regime [39].

Zoom In Zoom Out Reset image size

For a more rigorous analysis in the framework of an exact MPS formulation for the ground state of the XY model, see [56].

Whereas the previous subsection indicates why the peculiar structure of the eigenvalue spectrum of the MPS-TM allows one to recover the typical form of static correlation functions for $D\to \infty$, it makes no connection between the branches of eigenvalues appearing in this spectrum and the dispersion relations of the elementary excitations of the model. To make this connection more explicit, we now reiterate the well-known result that the single mode approximation (SMA) produces dispersion relations that are strongly dependent on the static structure factor [57–61]. In particular, we will illustrate why the energies, $E(k),$ of the lowest energy-momentum eigenstates, $|E_{k}^{\alpha }\rangle ,$ of a local translation-invariant Hamiltonian, $H={{\sum }_{n}}{{h}_{n}},$ become very small at momenta k where the TM has eigenvalues ${{\lambda }_{j}}={{{\rm e}}^{-{{\varepsilon }_{j}}+{\rm i}{{\phi }_{j}}}},$ with ${{\varepsilon }_{j}}$ approaching zero and $\pm k={{\phi }_{j}}$.

The generality of this subsection is based on the recent proof that elementary excitations on top of a gapped, strongly correlated ground state tend to be localized or particle-like [62]. This means that there exists a representation of all lowest-lying excited states by acting with the Fourier transform of a quasi local operator, ${{O}^{(\ell )}},$ centered around site n = 0 and having support on sites $n\in [-\ell ,\ell ]$, on the ground state:

$$\begin{eqnarray}&&|{{\phi }_{k}}({{O}^{(\ell )}})\rangle =\frac{1}{\sqrt{V}}\mathop{\sum }\limits_{n}{{{\rm e}}^{{\rm i}kn}}O_{n}^{(\ell )}|{{\psi }_{0}}\rangle ,\end{eqnarray} \tag{ 27 }$$

where $O_{n}^{(\ell )}={{U}_{n}}{{O}^{(\ell )}}U_{n}^{\dagger }$ with U_n being the lattice translation operator over n sites; V is the volume of the system. For the remainder of this subsection, we work in a finite system with periodic boundary conditions to be able to define normalizable energy-momentum states. The unnormalized state in (27) becomes exponentially close to a true, isolated energy-momentum eigenstate, $|E_{k}^{\alpha }\rangle ,$ with increasing linear size, $\ell ,$ of the support of ${{O}^{(\ell )}}$ [62].

We can readily use this set of states in a variational ansatz for low-energy excited states, and we will consider their energy expectation value

$$\begin{eqnarray}&&E(k)=\frac{\langle {{\phi }_{k}}({{O}^{(\ell )}})|H|{{\phi }_{k}}({{O}^{(\ell )}})\rangle }{\langle {{\phi }_{k}}({{O}^{(\ell )}})|{{\phi }_{k}}({{O}^{(\ell )}})\rangle }\end{eqnarray} \tag{ 28 }$$

which is a good approximation for the lowest excitation energies of H at momentum k. The operator, ${{O}^{(\ell )}},$ should have a zero vacuum expectation value, and can be a chosen Hermitian² . Together with the assumption of parity invariance for both H and ${{O}^{(\ell )}}$, it is well known [60] that (28) can be rewritten as

$$\begin{eqnarray}&&E(k)=\frac{1}{2}\frac{\langle {{\psi }_{0}}|[{{O}_{-k}},[H,{{O}_{k}}]]|{{\psi }_{0}}\rangle }{\langle {{\psi }_{0}}|{{O}_{-k}}{{O}_{k}}|{{\psi }_{0}}\rangle }=\frac{1}{2}\frac{F(k)}{S(k)},\end{eqnarray} \tag{ 29 }$$

where ${{O}_{k}}=\frac{1}{\sqrt{V}}{{\sum }_{n}}{{{\rm e}}^{{\rm i}kn}}{{O}_{n}}$; we have omitted the superscript $(\ell )$ for notational simplicity. Here, F(k) is the double commutator expectation value known as the oscillator strength and $S(k)={{S}_{OO}}(k)$ is the static structure factor, which is related to the static correlation function by a Fourier transform,

$$\begin{eqnarray}&&{{C}_{AB}}(n)=\int _{0}^{2\pi }\frac{{\rm d}k}{2\pi }\;{{S}_{AB}}(k){{{\rm e}}^{{\rm i}kn}}.\end{eqnarray} \tag{ 30 }$$

In (29), we have again discarded subscripts denoting the operators, $O,$ for notational simplicity.

For a Hamiltonian consisting of strictly local terms, $h_{j}^{(m)},$ acting on sites $[j,j+m]$ (i.e., with support on $m+1$ sites), F(k) is a finite polynomial in powers of ${{{\rm e}}^{{\rm i}k}}$ and can thus be bounded as

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} |F(k)| & = & |\mathop{\sum }\limits_{nj}{{{\rm e}}^{{\rm i}kn}}\langle [{{O}_{0}},[{{h}_{j}},{{O}_{n}}]]\rangle | \\ {} & {} & \leqslant 4(4\ell +2m+1)(2\ell +m+1)\parallel O{{\parallel }^{2}}\parallel h\parallel . \\ \end{array}\end{eqnarray} \tag{ 31 }$$

However, the static structure factor, $S(k),$ can become very large. This will happen when the momentum, $k,$ corresponds to the period of an oscillating static correlation function, $C(n),$ with a very large correlation length. For these momenta, one can therefore optimize over possible operators, $O,$ such that the resulting energy expectation value is very small, and thus, by virtue of the variational principle, excitations with small energies exist.

Generically, there is a one-to-one correspondence between the momenta, $k,$ where S(k) peaks and the complex arguments, ${{\phi }_{j}},$ of the TM eigenvalues. This can easily be shown in the framework of matrix product states where we use the same notation as in section 4.1. Due to translation invariance, the static structure factor can be written as

$$\begin{eqnarray}&&S(k)=\mathop{\sum }\limits_{n}{{{\rm e}}^{{\rm i}kn}}\langle {{O}_{0}}{{O}_{n}}\rangle .\end{eqnarray} \tag{ 32 }$$

This sum can be split into one part, $|x|\leqslant 2l,$ where operators ${{O}_{0}}$ and O_n overlap and a nother part, $|x|\gt 2l,$ where they commute. The first part can be bounded as

$$\begin{eqnarray}&&\mathop{\sum }\limits_{n=-2\ell }^{2\ell }{{{\rm e}}^{{\rm i}kn}}\langle {{O}_{0}}{{O}_{n}}\rangle ={{D}_{\ell }}(O)\leqslant (4\ell +1)\|O{{\|}^{2}}.\end{eqnarray} \tag{ 33 }$$

The remaining part can be written as

$$\begin{eqnarray}&&\begin{array}{ccccccccccccccc} 2\Re \mathop{\sum }\limits_{n\gt 2\ell }{{{\rm e}}^{{\rm i}kn}}\langle {{O}_{0}}{{O}_{n}}\rangle \\ \quad =2\Re \mathop{\sum }\limits_{n\gt 2\ell }{{{\rm e}}^{{\rm i}kn}}(0|{{\mathcal{J}}_{O}}\mathcal{T}_{{}}^{n-(2\,l+1)}{{\mathcal{J}}_{O}}|0) \\ \quad =2\Re \;{{{\rm e}}^{{\rm i}k(2\ell +1)}}\mathop{\sum }\limits_{n=0}^{\infty }(0|{{\mathcal{J}}_{O}}{{[{{{\rm e}}^{{\rm i}k}}\mathcal{T}_{{}}^{{}}]}^{n}}{{\mathcal{J}}_{O}}|0), \\ \end{array}\end{eqnarray} \tag{ 34 }$$

where ${{\mathcal{J}}_{O}}$ is the OTM defined in (19). To perform the geometric sum in the last line, we define the projector, $Q={\mathbb{1}}-|0)(0|$, which projects out the dominant eigenvector with eigenvalue one (i.e., $\mathcal{E}=Q\mathcal{T}_{{}}^{{}}=\mathcal{T}_{{}}^{{}}Q=Q\mathcal{T}_{{}}^{{}}Q=\mathcal{T}_{{}}^{{}}-|0)(0|$). We can then write $\mathcal{T}_{{}}^{n}=Q{{\mathcal{E}}^{n}}+|0)(0|={{\mathcal{E}}^{n}}Q+|0)(0|$, where the second term will not contribute to S(k) due to the zero vacuum expectation value of O. Assuming that there are no other eigenvalues with magnitude one, which is guaranteed if the MPS is injective, we can now safely perform the geometric sum and obtain

$$\begin{eqnarray}&&S(k)={{D}_{\ell }}(O)+2\Re (0|{{\mathcal{J}}_{O}}\frac{{{{\rm e}}^{{\rm i}k(2\ell +1)}}}{{\mathbb{1}}-{{{\rm e}}^{{\rm i}k}}\mathcal{E}}\ {{\mathcal{J}}_{O}}|0)\end{eqnarray} \tag{ 35 }$$

$$\begin{eqnarray}&&={{D}_{\ell }}(O)+2\Re \mathop{\sum }\limits_{j\gt 0}f_{j}^{OO}\frac{{{{\rm e}}^{{\rm i}k(2\ell +1)}}}{1-{{{\rm e}}^{-{{\varepsilon }_{j}}}}{{e}^{{\rm i}(k+{{\phi }_{j}})}}},\end{eqnarray} \tag{ 36 }$$

where the form factors, $f_{j}^{OO},$ are given by (21). It is apparent that the fraction can become very large for an eigenvalue of $\mathcal{T}_{{}}^{{}}$ with ${{\varepsilon }_{j}}$ close to zero (i.e., magnitude close to one) with k approaching the argument $\pm {{\phi }_{j}}$ ³ . By the above arguments, we can therefore expect low-energy excited states at these momenta, assuming that the form factors, $f_{j}^{OO},$ do not vanish.

In the case of symmetry-breaking and topologically nontrivial excitations, similar arguments lead to the same form of S(k) as in (36), where $\mathcal{E}$ is replaced with the mixed MPS-TM⁴ .

As a final point, we elaborate further on the relation between the static structure factor and the spectrum of the TM. For an injective MPS, it will always be possible to find an operator whose OTM exactly excites one or more of the eigenvalues of the TM on a branch with fixed ${{\phi }_{\alpha }}$ (i.e., $({{j}_{\alpha }}|{{{\mathcal{J}}}_{O}}|0)\ne 0$). This, however, could be an operator with very large support. Conversely, it is possible that for operators with small support, several eigenvalues are excited on branches with different arguments, ${{\phi }_{\alpha }},$ which are close together. For these operators, the static structure factor might then have a maximum at a momentum $k,$ which does not exactly correspond to one of the arguments, ${{\phi }_{\alpha }}$.

An example of this is the transition from commensurate to incommensurate order in the bilinear-biquadratic S = 1 Heisenberg chain [63–65]

$$\begin{eqnarray}&&{{H}_{{\rm BLBQ}}}=\mathop{\sum }\limits_{n}{\rm cos} (\theta ){{{\boldsymbol{S}} }_{n}}\cdot {{{\boldsymbol{S}} }_{n+1}}+{\rm sin} (\theta ){{({{{\boldsymbol{S}} }_{n}}\cdot {{{\boldsymbol{S}} }_{n+1}})}^{2}}.\end{eqnarray} \tag{ 37 }$$

We consider the regime $0\lt \theta \lt 0.25\pi$, where the oscillation period of static correlation functions changes from π to some incommensurate period exactly at the AKLT point, ${{\theta }_{{\rm VBS}}}={\rm arctan} (1/3)\approx 0.1024\pi$. However, by looking at the static structure factor of a simple one-site operator such as S^z, the peak stays at $k=\pi$ until some significantly larger value, $\tilde{\theta }\approx 0.1314\pi$. Based on the relation between the structure factor and the minima of the dispersion relation, we have reason to expect that the dispersion relation starts shifting away from π at this latter value of θ. It appears, however, that the minimum of the dispersion relation starts shifting at $\theta \approx 0.12\pi$, a value in between ${{\theta }_{{\rm VBS}}}$ and $\tilde{\theta }$. In figure 10 we observe from the full spectrum of the TM that this happens when the eigenvalues with the largest magnitudes have finished aligning along the line of constant phase, ${{\phi }_{\alpha }}\approx 0.755\pi$. At this point it appears that the support of the operator, $O,$ generating the excitation has shifted from the—by now very small—branch with ${{\phi }_{\alpha }}=\pi$ to the now fully aligned branch with ${{\phi }_{\alpha }}\approx 0.755\pi$. This value of ${{\phi }_{\alpha }}$ further shifts towards $2\pi /3$ with $\theta \to 0.25\pi$, where the gap then closes at momenta $k=0,\pm 2\pi /3$ [66–68].

Zoom In Zoom Out Reset image size

From this example it is apparent that, especially in the vicinity of such peculiar crossover points, it is worthwhile to look at the full spectrum of the TM. It generically contains more information than a simple static structure factor, as the TM is completely independent of the choice of operator. It illustrates that there is a crossover regime where the MPS-TM has to react by developing additional branches with constant ${{\phi }_{\alpha }},$ and the support of the operator, $O,$ generating the excitation has to shift to this newly developed branch. It appears that this process is completed after the eigenvalues have fully aligned along the newly developed branch. It would be interesting to investigate this process further with the precise knowledge about the operator, $O,$ generating the excitation and how the form factors $({{j}_{\alpha }}|{{{\mathcal{J}}}_{O}}|0)$ develop with θ.

For now we conclude with the observation that there are indeed also situations where the locations of the lowest-lying TM eigenvalues do not precisely coincide with the minima of the dispersion. These are, however, very special cases like the one discussed above, where the TM has to adapt to changing conditions within some crossover regime. There, the peculiar structure of the TM eigenvalue spectrum is, however, also different from normal situations as investigated in section 3, which can serve as an indicator for such exceptional situations.

In the previous two subsections, we have discussed two effects of the peculiar distribution of the eigenvalues of the MPS-TM. First, the clustering of eigenvalues onto lines starting from the origin allows one to recover the typical form of static correlation functions in gapped quantum ground states in the limit $D\to \infty$. Second, this distribution causes peaks in the static structure factor, which can be related to minima in the dispersion relation of excitations using the single mode approximation. In retrospect, the first effect only requires a dense distribution of eigenvalues along a line, without any connection between the location of these lines and the dispersion relation of the excitations of the system. The second argument only requires the existence of a single dominant eigenvalue with a phase corresponding to the momentum of the minimum of the dispersion relation, and does not explain why there needs to be a dense distribution of eigenvalues. This leaves open the question of whether a different structure of the eigenvalue distribution could give rise to similar effects. Put differently, we would like to answer the reverse question: to what extent are static correlation functions and the clustering of the eigenvalues of the MPS-TM determined by the excited states of the Hamiltonian?

For Lorentz-invariant field theories, the Källén-Lehmann representation of two-point correlation functions provides a direct connection to the excitation spectrum of the Hamiltonian. Whereas the Källén-Lehmann representation exists for arbitrary dynamical correlation functions, we here focus on the case of static correlation functions between two scalar operators, A(x) and B(y), where it is given by

$$\begin{eqnarray}&&\begin{array}{ccccccccccccccc} \langle \Psi |A(x)B(y)|\Psi \rangle \\ =\int {\rm d}{{M}^{2}}\rho ({{M}^{2}})\int \frac{{{{\rm d}}^{d}}k}{{{(2\pi )}^{d}}}\frac{{{{\rm e}}^{{\rm i}k(x-y)}}}{2\sqrt{{{M}^{2}}+{{k}^{2}}}}\langle \Psi |A(0)|{{M}^{2}},0\rangle \langle {{M}^{2}},0|B(0)|\Psi \rangle , \\ \end{array}\end{eqnarray} \tag{ 38 }$$

where—contrary to standard field theory notation—x and y denote spatial vectors, k denotes momentum, and d represents the number of spatial dimensions. The first integral is over all possible masses in the theory, and $\rho ({{M}^{2}})$ corresponds to the density of states. If the lowest-lying excitations correspond to single-particle excitations with discrete masses, ${{M}_{\alpha }}$, then $\rho ({{M}^{2}})$ will contain a contribution, ${{\sum }_{\alpha }}\delta ({{M}^{2}}-M_{\alpha }^{2})$. The state $|M_{\alpha }^{2},0\rangle$ corresponds to the presence of such an excitation with mass ${{M}_{\alpha }}$ and momentum zero.

Let us again restrict to the case of d = 1. We are interested in the long-range behavior of correlation functions in a lattice model, which is clearly dictated by the low-energy behavior of the model. If this low-energy behavior can be captured by a Lorentz-invariant theory with masses ${{M}_{\alpha }}$, a remnant of the Källén-Lehmann representation of correlation functions should exist in the lattice model. For the field theory, the minimum of all dispersion relations is at momentum zero. However, in many cases taking the continuum limit of a lattice theory requires that, for example, N sites are blocked, and a single lattice dispersion relation with several minima for momenta ${{\phi }_{\alpha }}=2\pi \alpha /N$ with $\alpha =0,\ldots ,N-1$ gives rise to N independent dispersion relations of the field theory. A prototypical example for N = 2 is the XX-model, which corresponds to the staggered fermion discretization of relativistic Dirac fermions, where the two components of the Dirac spinor are put on even and odd sites, respectively [69].

Equation (38) represents the static correlation function in real space as an integral over momentum space. The corresponding representation for the static structure factor in momentum space can therefore easily be identified. In the Lorentz-invariant case, the single-particle excitations add a contribution to the static correlation function that is given by a constant times the inverse of the dispersion relation of that excitation. For typical operators, these will be the dominant contributions. If the low-energy behavior of a lattice model is Lorentz-invariant, we can thus expect that S(k) should also receive contributions, ${{S}_{\alpha }}(k)$, of the form

$$\begin{eqnarray}&&{{S}_{\alpha }}(k)={{c}_{\alpha }}{{[m_{\alpha }^{2}+v_{\alpha }^{2}{{(k-{{\phi }_{\alpha }})}^{2}}]}^{-1/2}}.\end{eqnarray} \tag{ 39 }$$

Here, ${{c}_{\alpha }}$ is a constant depending on the choice of operators A and B, ${{m}_{\alpha }}$ is the mass of the excitations in units of the inverse lattice spacing, and ${{v}_{\alpha }}$ is the characteristic velocity, which for a proper Lorentz-invariant low-energy behavior should be the same (and thus independent of α) for all excitations.

Since this form is only expected to hold for $k\approx {{\phi }_{\alpha }}$ we can instead choose a $2\pi$-periodic version of the dispersion relation to get contributions of the form

$$\begin{eqnarray}&&{{S}_{\alpha }}(k)={{c}_{\alpha }}{{\{m_{\alpha }^{2}+v_{\alpha }^{2}[2-2{\rm cos} \;(k-{{\phi }_{\alpha }})]\}}^{-1/2}}.\end{eqnarray} \tag{ 40 }$$

We can now transform from the $2\pi$-periodic variable, k, to the complex variable $z={{{\rm e}}^{{\rm i}k}}$, and we change the notation of the static structure factor as $S(k)\to S(z={{{\rm e}}^{{\rm i}k}})$ to write

$$\begin{eqnarray}&&C(n)={{\oint }_{\mathcal{C}}}\frac{{\rm d}z}{2\pi {\rm i}}\;S(z){{z}^{n-1}},\end{eqnarray} \tag{ 41 }$$

where the contour integral is over the unit circle $\mathcal{C}$, which is where S(z) is originally defined (see figure 11).

Zoom In Zoom Out Reset image size

If now we were to construct an analytic continuation of S(z), we expect that every contribution of the form of (40) produces a square root singularity at z = 0 and inverse square root singularities at the points

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} z_{\alpha }^{(\pm )} & = & {{{\rm e}}^{{\rm i}{{\phi }_{\alpha }}}}\left[ 1\pm {{\Delta }_{\alpha }}\left( \sqrt{1+\Delta _{\alpha }^{2}/4}\pm {{\Delta }_{\alpha }}/2 \right) \right] \\ {} & {} & \approx {{{\rm e}}^{{\rm i}{{\phi }_{\alpha }}\pm {{\Delta }_{\alpha }}+{\rm O}(\Delta _{\alpha }^{3})}} \\ \end{array}\end{eqnarray} \tag{ 42 }$$

with ${{\Delta }_{\alpha }}={{m}_{\alpha }}/{{v}_{\alpha }}$. We can thus choose the branch cuts to go from $z_{\alpha }^{(-)}$ to 0 and from $z_{\alpha }^{+}$ to $+\infty$. Assuming that there are no other singularities, we can then deform the integration contour as in figure 11. In the limit $\eta ,\eta ^{\prime} \to 0$, the arc segments ${{\mathcal{A}}_{\alpha }}$ and $\mathcal{A}_{\alpha ,\beta }^{\prime }$ do not contribute, as they correspond to square root and inverse square root singularities, respectively, whereas line segments ${{\mathcal{L}}_{\alpha }}$ and $\mathcal{L}_{\alpha }^{\prime }$ produce equivalent contributions

$$\begin{eqnarray}&&{{{\rm e}}^{{\rm i}{{\phi }_{\alpha }}n}}\int _{0}^{{{{\rm e}}^{-{{\Delta }_{\alpha }}}}}{\rm d}y\;{{y}^{n-1/2}}{{[(y-{{{\rm e}}^{-{{\Delta }_{\alpha }}}})(y-{{{\rm e}}^{+{{\Delta }_{\alpha }}}})]}^{-1/2}}\end{eqnarray} \tag{ 43 }$$

where we have transformed to $y=z{{{\rm e}}^{-{\rm i}{{\phi }_{\alpha }}}}$ and the contributions to the correlation function, C(n), are now written as integrals over $0\leqslant y\leqslant {{{\rm e}}^{-{{\Delta }_{\alpha }}}}$ along the branch cuts.

It is apparent that a finite D MPS approximation of the static correlation function tries to reproduce this continuum form with a discrete sum over the eigenvalues of the TM. These eigenvalues cluster on the branch cuts of the static structure factor, which we have related to the single-particle excitation spectrum by assuming a Lorentz-invariant low-energy behavior.

For large n, where the low-energy contributions dominate and the approximations are valid, we can again use the saddle point approximation around the point $y={{{\rm e}}^{-{{\Delta }_{\alpha }}}}$ and change variables to $y={{{\rm e}}^{-({{\Delta }_{\alpha }}+x)}}$ dominate to obtain for contributions ${{C}_{\alpha }}(n)$ to C(n) up to some factor

$$\begin{eqnarray}&&{{C}_{\alpha }}(n)\sim {{{\rm e}}^{{\rm i}{{\phi }_{\alpha }}n}}{{{\rm e}}^{-{{\Delta }_{\alpha }}n}}\int _{0}^{\infty }{\rm d}x\;\frac{{{{\rm e}}^{-x(n-\frac{1}{2})}}}{\sqrt{x}}(1+{\rm O}(x))\sim \frac{{{{\rm e}}^{{\rm i}{{\phi }_{\alpha }}n}}{{{\rm e}}^{-{{\Delta }_{\alpha }}n}}}{\sqrt{n-1/2}},\end{eqnarray} \tag{ 44 }$$

which reproduces correlations of the Ornstein-Zernike form, with correlation lengths $\xi _{\alpha }^{-1}={{\Delta }_{\alpha }}$. Therefore, the dominant contribution for n large stems from the branch with the smallest ${{\Delta }_{\alpha }}={{m}_{\alpha }}/{{v}_{\alpha }}$.

It is interesting to compare this to the discussion in the previous section. The single-mode approximation for the dispersion relation, E(k), was given in (29) as $E(k)=F(k)/[2S(k)]$, which allowed us to conclude that peaks in the static structure factor, S(k), produce minima in the corresponding dispersion relation. Assuming Lorentz invariance, the Källén-Lehmann representation produces contributions to S(k) of the form ${{S}_{\alpha }}(k)\sim {{c}_{\alpha }}/E(k)$, with ${{c}_{\alpha }}$ being a k-independent constant, which looks like the reverse relation: minima in the dispersion relation give rise to peaks in the static structure factor. However, the single-mode approximation also allows for minima that are not caused by S(k), but rather by a small or vanishing value for the oscillator strength, F(k). It would be interesting if one could show that such excitations are necessarily related to low-energy features which have an intrinsically non relativistic description.

As a final justification for the argumentation in this subsection, we apply the above results to the case of the XY-model in the incommensurate gapped phase investigated in section 3.1. There, it is observed that for the elementary excitations, the smallest excitation energy, ${{E}_{{\rm min} }}$ (i.e., the energy gap), and the eigenvalue of the TM with the second-largest magnitude, ${{\lambda }_{1}}$, are related by some characteristic velocity, v₁, via ${{\varepsilon }_{1}}={{E}_{{\rm min} }}/{{v}_{1}}$, where ${{\varepsilon }_{1}}=-{\rm log} |{{\lambda }_{1}}|$ and ${{\varepsilon }_{1}}={{\xi }^{-1}}$, with ξ being the correlation length, as established in section 4.1. If we assume a Lorentz-invariant low-energy behavior—with some characteristic velocity, v_c–of the dispersion, E(k), around the minimum of the form

$$\begin{eqnarray}&&E(k)\approx \sqrt{E_{{\rm min} }^{2}+v_{c}^{2}{{(k-{{k}_{{\rm min} }})}^{2}}},\end{eqnarray} \tag{ 45 }$$

as also used in (39), we can deduce by virtue of (42) and (44) that ${{\xi }^{-1}}={{\Delta }_{1}}={{E}_{{\rm min} }}/{{v}_{c}}$, where we have interpreted the energy gap, ${{E}_{{\rm min} }}$, as the lowest mass, m₁. On the other hand, with ${{\xi }^{-1}}={{\varepsilon }_{1}}={{E}_{{\rm min} }}/{{v}_{1}}$, we thus see that ${{v}_{1}}={{v}_{c}}$ is exactly the characteristic velocity appearing in (45).

With $E(k)=\sqrt{{{(g-{\rm cos} (k))}^{2}}+{{\gamma }^{2}}{{{\rm sin} }^{2}}(k)}$ known exactly for the XY-model [38–40], and assuming (45) to hold around $k={{k}_{{\rm min} }}$, we can then estimate v₁ as

$$\begin{eqnarray}&&v_{1}^{2}=\frac{1}{2}\frac{{{{\rm d}}^{2}}E{{(k)}^{2}}}{{\rm d}{{k}^{2}}}{{|}_{k={{k}_{{\rm min} }}}}=\frac{{{(1-{{\gamma }^{2}})}^{2}}-{{g}^{2}}}{1-{{\gamma }^{2}}}.\end{eqnarray} \tag{ 46 }$$

For the parameters $\gamma =0.3$ and g = 0.2 considered in section 3.1, this yields an estimate of ${{v}_{1}}=0.9306$, which differs from the value ${{v}_{1}}=0.9409$ obtained from ${{v}_{1}}={{E}_{{\rm min} }}/{{\varepsilon }_{1}}$ by only $\approx 1\%$.

In this section we derive more rigorous statements connecting the decay of static connected correlation functions to the dispersion, E(k), of low-energy excitations of local translation-invariant Hamiltonians in the thermodynamic limit. We generalize the seminal work of Hastings [3], which proved that the inverse of the energy gap, Δ, of a local translation-invariant Hamiltonian times a constant serves as an upper bound for the correlation length, ξ, of connected static correlation functions. This implies that if the gap vanishes, these correlations may (and in most cases will) be long range with a diverging correlation length.

The proof gives a statement relating the smallest overall excitation energy and the largest correlation length in the system, but it does not take into account momentum information. In this work we extend the results in [3] to derive bounds on the decay of momentum-filtered correlation functions and relate them to the dispersion, E(k), of low-energy excitations. Specifically, we show that the inverse of the energy gap, E(k), at a specific momentum, k, times a constant serves as an upper bound for a momentum-resolved correlation length, ${{\xi }_{k}}$. Conversely, the existence of a finite correlation length, ${{\xi }_{k}},$ thus implies an upper bound for the energy, E(k).

The detailed derivation of the bound in the most general setting is given in the appendix; here we, for the sake of simplicity, present the result for one-dimensional lattice systems and operators acting on single sites. We start from the static connected correlation function of operators ${{A}_{i=0}}$ and ${{B}_{j=\ell \gt 0}}$, for which we assume zero vacuum expectation value. We now attempt to extract momentum-space information while retaining real-space information by replacing operator ${{B}_{\ell }}$ at site ℓ with a Gaussian wave packet centered around site ℓ, defined as

$$\begin{eqnarray}&&{{\tilde{B}}_{\ell }}(k)={{N}_{r}}\mathop{\sum }\limits_{n}{{{\rm e}}^{-\frac{{{n}^{2}}}{2r}}}{{{\rm e}}^{{\rm i}kn}}{{B}_{\ell +n}}\end{eqnarray} \tag{ 47 }$$

where N_r is a normalization constant. We define the momentum-filtered correlation function as

$$\begin{eqnarray}&&{{C}_{k}}(\ell )=\langle {{A}_{0}}{{\tilde{B}}_{\ell }}(k)\rangle ,\end{eqnarray} \tag{ 48 }$$

which corresponds to a Fourier transform of the product of the static correlation function and a Gaussian wave packet centered around site ℓ. In momentum space this yields the convolution of the static structure factor, S(k), and another Gaussian wave packet in momentum space.

In the appendix, it is proven that by tuning r as a fraction of ℓ, the momentum-filtered correlation function, ${{C}_{k}}(\ell )$, can—for sufficiently large ℓ–be bounded by

$$\begin{eqnarray}&&|{{C}_{k}}(\ell )|\leqslant {{c}_{1}}\parallel A\parallel \parallel B\parallel {{{\rm e}}^{-{{c}_{2}}\ell }},\end{eqnarray} \tag{ 49 }$$

where c₁ and c₂ are some constants and the inverse of c₂ is an upper bound for the momentum-filtered correlation length ${{\xi }_{k}}$, given by

$$\begin{eqnarray}&&{{\xi }_{k}}\leqslant \frac{1}{{{c}_{2}}}=\frac{1}{\delta }+\frac{{{v}_{{\rm LR}}}}{{{E}^{*}}(k,\delta )}.\end{eqnarray} \tag{ 50 }$$

Here, ${{v}_{{\rm LR}}}$ is the characteristic Lieb-Robinson velocity [1] and ${{E}^{*}}(k,\delta )={{{\rm min} }_{|k-k^{\prime} |\leqslant \delta }}E(k^{\prime} )$ is the minimum of the dispersion E(k) in an interval around k given by δ.

The constant, δ, is introduced in the proof and and can be tuned to obtain the sharpest bound. From (50) it is clear that there is a tradeoff, as increasing δ generally leads to a decrease of ${{E}^{*}}(k,\delta )$. If k, however, corresponds to a minimum in the dispersion E(k), then the function ${{E}^{*}}(k,\delta )$ is largely insensitive to δ in some region around the minimum, and we can choose δ as large as possible within this region. If k, on the other hand, corresponds to a regular point where $\frac{{\rm d}E}{{\rm d}k}\ne 0,$ then there is a direct effect from increasing δ to decreasing E*. An optimal choice of δ is thus dependent on the form of E(k).

Colloquially, the above result means that the decay of the momentum-filtered correlation functions is dictated by the corresponding low-energy states around that given momentum. Equivalently, one can say that a large correlation length for a given momentum, k, implies a small excitation energy, E(k).

In practice, this means that one can deduce an upper bound for the low-energy spectrum of the Hamiltonian by looking at the momentum dependence of momentum-filtered correlation functions, where one would have to optimize over the parameter δ ⁵ .

Consider a one-dimensional local-lattice Hamiltonian, H, with translation invariance and a unique ground state, $|{{\psi }_{0}}\rangle ,$ with ground-state energy ${{E}_{0}}=0$. To avoid issues with the infrared orthogonality catastrophe, we work in a finite system of N sites with periodic boundary conditions for this section. The ground state can be obtained by an evolution in imaginary time β

$$\begin{eqnarray}&&|{{\psi }_{0}}\rangle =\mathop{{\rm lim} }\limits_{\beta \to \infty }\frac{{{{\rm e}}^{-\beta H}}|{{\phi }_{{\rm init}}}\rangle }{\parallel {{{\rm e}}^{-\beta H}}|{{\phi }_{{\rm init}}}\rangle \parallel },\end{eqnarray} \tag{ 51 }$$

with $|{{\phi }_{{\rm init}}}\rangle ,$ some initial state which is nonorthogonal to $|{{\psi }_{0}}\rangle$. Since H is a sum of local terms, h_n, we can break β into small imaginary-time steps, δ, and use a Suzuki-Trotter decomposition of ${{{\rm e}}^{-\delta H}}\approx \prod _{n}^{{}}{{{\rm e}}^{-{{\delta }^{{}}}{{h}_{n}}}}$ [6, 70]. There are several strategies to write this (or alternative) decomposition(s) as a two-dimensional tensor network with translation invariance in the spatial direction [14, 15, 71]. Obviously, all the information of H is thus encoded into this tensor network.

If the state $|{{\phi }_{{\rm init}}}\rangle$ is initially in the form of a translation-invariant MPS, then we can also interpret $|{{\psi }_{0}}\rangle$ as a translation-invariant MPS by grouping contractions along imaginary time. Each column of the tensor network is an MPS matrix A^s, itself being a half-infinite MPO. A^s then exactly represents the ground state up to a Trotter error⁶ . A graphical representation of this construction is given in figure 12.

Zoom In Zoom Out Reset image size

We immediately see that $\langle {{\psi }_{0}}|{{\psi }_{0}}\rangle =Tr_{A}^{N}$ where the MPS-TM $\mathcal{T}_{A}^{{}},$ is defined in (3). For a system with a unique ground state, the boundary conditions at $\beta =0$ are irrelevant, and an equivalent network with periodic boundary conditions in the temporal direction would be obtained for the thermal partition function, ${{Z}_{\beta }}=T{{r}^{-\beta H}},$ in the limit $\beta \to \infty$. The exact MPS-TM for the ground state $|{{\psi }_{0}}\rangle$ thus corresponds to the QTM at zero temperature, defined in [6, 7].

It is important to note that all the information about the QTM—in particular its eigenvalues—is thus contained within the ground state $|{{\psi }_{0}}\rangle$ and its exact MPS representation, A^s. Note, however, that this exact representation with an exponentially diverging bond dimension differs from a finite D approximation, ${{\tilde{A}}^{s}},$ which can be obtained, for instance, from some variational algorithm. For an example of an analytic derivation of such an exact ground state MPS representation, A^s, and the effect of truncating to finite D for the case of the $S=1/2$ XY model, see [56].

In section 5.3 we present a construction how such a finite D approximation—which only retains degrees of freedom relevant for the physical degree of freedom, s–can be obtained from this exact MPS representation, A^s, and thus from the true QTM. To understand how this relates the MPS-TM to the Hamiltonian of the system, we first need to discuss how the latter relates to the exact QTM, which is the topic of the next subsection.

To understand what information about the underlying Hamiltonian can be extracted from the knowledge about the QTM, we first consider the case of relativistic (1+1)-dimensional field theories. The vacuum or ground state can be expressed in terms of a path integral formulation very similar to figure 12, but where both real space and imaginary time are continuous. We can then identify a factor, ${{{\rm e}}^{-\delta H}},$ for some infinitesimally small δ with a narrow horizontal slice of the entire network. Due to relativistic invariance, however, such slices are invariant under Euclidean rotations between real space and imaginary time. This means in particular that a vertical slice–corresponding to a field theory analogue of $\mathcal{T}_{{}}^{{}},$—is also equal to ${{{\rm e}}^{-\delta H}}$. Hence, knowing the spectrum of $\mathcal{T}_{{}}^{{}}$ immediately yields knowledge about the spectrum of H. This in turn implies that all the information about the eigenvalues of H is already contained within the ground state. Of course, Lorentz invariance also strongly restricts the dispersion relations of the theory, such that they are completely determined by a single parameter—the mass of the excitation.

The relevant question is thus how much Euclidean invariance between real space and imaginary time remains in nonrelativistic lattice systems. As a concrete example, we consider the one-dimensional XYZ model on a chain with N sites and periodic boundary conditions

$$\begin{eqnarray}&&{{H}_{{\rm XYZ}}}=-\mathop{\sum }\limits_{j}{{J}^{x}}S_{j}^{x}S_{j+1}^{x}+{{J}^{y}}S_{j}^{y}S_{j+1}^{y}+{{J}^{z}}S_{j}^{z}S_{j+1}^{z}.\end{eqnarray} \tag{ 52 }$$

By invoking a Suzuki-Trotter decomposition [6, 70], the finite temperature partition function, ${{Z}_{\beta }},$ can be mapped onto a classical two-dimensional eight-vertex model [43, 72] of dimension N × M, where the statistical weights, W, of the arrow configurations depend on the model parameters, the inverse temperature β, the Trotter number M, and the type of chosen decomposition. Here we use a generalized Suzuki-Trotter decomposition, ${{{\rm e}}^{-\delta H}}={{{\rm e}}^{-\delta {{\sum }_{n}}{{h}_{n}}}}\approx \prod _{n}^{{}}{{{\rm e}}^{-{{\delta }^{{}}}{{h}_{n}}}}$, with $\delta =\beta /M$. This leads to a real space decomposition introduced by Suzuki [6], where the QTM, $\mathcal{T}_{{}}^{{}},$ is translation invariant on a slanted lattice [7], and we assume M to be a multiple of N. The statistical weights are then given by ${{W}_{abcd}}=\langle cd|{{{\rm e}}^{-\delta {{h}_{n}}}}|ab\rangle$ [43]. The same network is obtained after a 90-degree rotation of the lattice, and considering an equivalent effective model on a chain of M sites and Trotter number N at an effective inverse temperature, $\tilde{\beta }$. The QTM of the original lattice can therefore be written in terms of the effective model Hamiltonian $\tilde{H}$ as $\mathcal{T}_{{}}^{{}}={{{\rm e}}^{-\tilde{\delta }\tilde{H}}}$ with $\tilde{\delta }=\tilde{\beta }/N$ and effective statistical weights $\tilde{W}$ up to a Trotter error in N. The effective model parameters and the effective inverse temperature $\tilde{\beta }$ can be obtained from the vertex weights after rotating (i.e., from ${{\tilde{W}}_{bdac}}={{W}_{abcd}}$). For a graphical representation of this mapping, see figure 13. From this relation the spectrum of the effective Hamiltonian can, in principle, be studied by looking at the QTM generated by the ground state of the original Hamiltonian .

Zoom In Zoom Out Reset image size

Note that depending on the model and parameter regime, the effective parameters may also become complex, resulting in a non-Hermitian effective Hamiltonian, and thus accounting for non-Hermitian transfer matrices. This particular fact was already discussed in the case of systems with Lorentz invariance, as the continuum limit can often only be taken after blocking N sites. In that case, all eigenvalues of the TM would have phases ${{\phi }_{\alpha }}=2\pi n/N$ with $n=0,\ldots ,N-1$, and by defining $\tilde{H}=-{\rm log} (\mathcal{T}_{{}}^{N})$ up to some energy scale, a Hermitian-effective Hamiltonian would be obtained. For systems with incommensurate order, where the eigenvalues can have arbitrary phases that are not fractions of $2\pi$, this is no longer possible.

In this subsection we show how to obtain an MPS approximation, ${{\tilde{A}}^{s}},$ with a finite bond dimension from the exact QTM constructed in the previous subsections. As was shown previously, an exact MPS representation of the ground state can be constructed from imaginary time evolution, where the MPS matrices A^s are given as a semi-infinite MPO an with exponentially diverging virtual dimension (see figure 12). This construction allows one to identify the MPS-TM with the exact QTM at zero temperature.

In the following, we assume that the QTM can be written in terms of an effective local Hamiltonian, $\tilde{H}={{\sum }_{n}}{{\tilde{h}}_{n}}$ as $\mathcal{T}_{{}}^{{}}={{{\rm e}}^{-\tilde{H}}}$, for instance via the construction of the previous subsection. In this representation, the MPS-TM appears to be completely translation invariant in the imaginary time direction, with no special role being played by the physical index s, where the matrices, A^s and ${{\bar{A}}^{s}},$ of the exact MPS representation of ket and bra are connected. The point is of course that, for expectation values of operators, O, different from the identity, there would be an extra insertion at imaginary time $\tau =0,$ corresponding to the OTM, ${{\mathcal{J}}_{O}},$ defined in (19). Here we have set the origin $\tau =0$ of the imaginary time axis at the point where A^s and ${{\bar{A}}^{s}}$ are connected.

We can then interchange the roles of real space and imaginary time and give two distinct new interpretations to the expectation value $\langle \Psi |O(x)|\Psi \rangle$. First, we can interpret the new imaginary time direction as the evolution of a pure state of an infinite one-dimensional virtual system in the x direction according to the MPS-TM $\mathcal{T}_{{}}^{{}},$ or equivalently, the corresponding Hamiltonian $\tilde{H}$. At certain 'times' x, there is an insertion of some impurity, O, at the fixed coordinate $\tau =0$, which destroys the translation invariance of the virtual system. An alternative interpretation is given below while constructing an approximation with a finite bond dimension of the exact ground-state MPS, A^s.

To arrive at a finite virtual dimension D, it is necessary to restrict the exponentially diverging amount of virtual degrees of freedom (DOF) of A^s to a finite subset that is relevant for the physical DOF, s, at $\tau =0$. We thus need to identify the relevant low-energy subspace of the virtual system to describe the evolution of the impurity at position $\tau =0$ as a function of x; the finite D approximation, ${{\tilde{A}}^{s}},$ is then obtained from the exact MPS, A^s, by projecting onto this relevant subspace.

The relevant low-energy subspace for an impurity problem can be obtained by applying real space RG transformations, as was first shown in the seminal work of Wilson using his NRG [73]. However, we here follow the more recent construction using the MERA [16–20]. This approach allows one to identify the relevant DOF as those living at the causal cone of the impurity [19, 20]. For completeness, we repeat this argument for our specific case.

The following construction starts from the assumption that there exists a sequence of real space RG transformations, ${{\mathcal{U}}_{r}},$ that renormalizes $\tilde{H}$ onto its low-energy subspace, where r labels the layers of successively applied RG transformations. Equivalently, ${{\mathcal{U}}_{r}}$ renormalizes $\mathcal{T}_{{}}^{{}}$ onto the subspace spanned by its dominant eigenvectors. For concreteness, we consider a real space RG procedure that coarse grains four neighboring sites into two renormalized sites, which can be realized, for example, by a modified binary MERA. If $\tilde{H}$ is scale invariant, the transformations become independent of the layer index, $r,$ after some amount of initial layers, r^*. At this point, all RG-irrelevant terms have been removed and the renormalized Hamiltonian is a fixed point of the scale-invariant RG transformation, $\mathcal{U}$. For non-scale-invariant Hamiltonians, the series of RG transformations terminates after ${{r}_{{\rm max} }}\approx {\rm log} (\mathcal{R})$ layers, where $\mathcal{R}$ is the dominant length scale of the Hamiltonian. At this point there are no relevant DOF left. These two cases can, for example, be represented by a scale-invariant or a finite-range MERA, respectively [17]. An example for the latter case is shown in figure 14 .

Zoom In Zoom Out Reset image size

The MERA construction allows one to conclude that any perturbation at position $\tau =0$ can only affect the DOF living at its causal cone [20]. This causal cone is shown as red solid lines in figure 14. In particular, let us focus on one side of this causal cone (e.g., the lower half of figure 14). It is apparent that only one index of the RG network protrudes from the boundary of the causal cone for each layer. We can therefore interpret these legs as the sites of an effective lattice system, ${{\mathcal{L}}^{W}},$ defined along the boundary of the causal cone, which we call the Wilson chain. We label the sites along this chain by the layer index, r. Note that the r-th site on this chain is an effective renormalized description of ${{N}_{r}}\approx {{2}^{r}}$ sites of the original lattice, reminiscent of the logarithmic discretization introduced by Wilson, hence the chosen nomenclature [19, 20, 73].

Contracting the RG network outside the causal cone (non shaded region in figure 14) allows, one to renormalize the MPS-TM, $\mathcal{T}_{A}^{{}},$ into a new TM, $\mathcal{T}_{{}}^{W},$ along the Wilson chain, ${{\mathcal{L}}^{W}}$ [19, 20]. It is immediately clear that $\mathcal{T}_{{}}^{W}$ corresponds to the TM, $\mathcal{T}_{\tilde{A}}^{{}},$ of a new MPS with matrices ${{\tilde{A}}^{s}}$, which are obtained by projecting the exact MPS representation, ${{A}^{s}},$ exactly onto this subspace of relevant DOF along the Wilson chain. For non critical systems, the Wilson chain is of finite length $L,$ and we obtain a virtual system with finite dimension $D\approx {{\chi }^{{{r}_{{\rm max} }}}}$, where χ is the bond dimension of the RG network. For critical, scale-invariant systems, however, the Wilson chain is still infinite and any truncation to a finite system introduces some error.

It is of course well known that ground states with finite correlation length ξ can be well approximated by a finite-D MPS [4]. By virtue of the previous sections, the effective Hamiltonian $\tilde{H}$ is gapped, where the gap is exactly equal to ${{\xi }^{-1}}$. It is interesting to contrast this interpretation with the more usual considerations regarding the gap of the actual physical Hamiltonian, $H,$ of the system. For Lorentz-invariant systems, H and $\tilde{H}$ are equal and the correlation length is directly given by the inverse of the gap in units where the speed of light c = 1. Without Lorentz invariance, there is no obvious relation between H and its Euclidean rotation, $\tilde{H}$. The result of [3] allows one to bound ξ (and thus the gap of $\tilde{H}$) in terms of the gap, $\Delta ,$ of H. When the low-energy behavior of H still allows for an effective relativistic description, it can indeed be expected that $\xi \sim v/\Delta ,$ where v replaces the speed of light with some characteristic speed of the system. However, there is also the possibility of an inherently non relativistic low-energy behavior, resulting in, for example, gapless Hamiltonians, $H,$ for which the ground state is an exact MPS [74], so that the corresponding $\tilde{H}$ is gapped.

Finally, we can obtain an alternative interpretation by applying the Jamiołkowski isomorphism to map the pure state defined on the infinite, one-dimensional virtual system of the previous discussion to the density matrix of a half-infinite, one-dimensional virtual system with a boundary at $\tau =0$. By again interchanging the roles of real space and imaginary time, this half-infinite system undergoes dissipative evolution in the new time direction x, corresponding to a Hamiltonian containing the terms of $\tilde{H}$ with support on the region $\tau \gt 0$, and additional Lindblad operators corresponding to the action of $\tilde{H}$ across the boundary and the additional action of an operator, $O,$ inserted at certain 'times', x. Hence, all the Lindblad operators are acting near the boundary, and the truncation of the bond dimension corresponds to selecting the relevant DOF to describe the boundary of the system. It is a virtue of the MERA that it naturally unifies the process of selecting the relevant DOF for boundaries and impurities [20]. The resulting virtual system becomes effectively zero-dimensional, and can be interpreted as providing a holographic description of the physical system [26].

We start by defining ${{\tilde{A}}_{X}}$ as

$$\begin{eqnarray}&&{{\tilde{A}}_{X}}=\frac{1}{2\pi }\int _{-\infty }^{+\infty }\frac{{{{\rm e}}^{{\rm i}Ht}}{{A}_{X}}{{{\rm e}}^{-{\rm i}Ht}}}{-{\rm i}t+\epsilon }{{{\rm e}}^{-\frac{{{t}^{2}}}{2q}}}{\rm d}t\end{eqnarray} \tag{ A.9 }$$

and first show that

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} |\langle {{\Phi }_{k^{\prime} ,\alpha }}|{{\tilde{A}}_{X}}|{{\Psi }_{0}}\rangle | & = & \left| \langle {{\Phi }_{k^{\prime} ,\alpha }}|{{A}_{X}}|{{\Psi }_{0}}\rangle \right|\left| \frac{1}{2\pi }\int _{-\infty }^{+\infty }{\rm d}t\;\frac{{{{\rm e}}^{{\rm i}{{E}_{k^{\prime} ,\alpha }}t-\frac{{{t}^{2}}}{2q}}}}{-it+\epsilon } \right| \\ {} & {} & =\left| \langle {{\Phi }_{k^{\prime} ,\alpha }}|{{A}_{X}}|{{\Psi }_{0}}\rangle \right|\left| \frac{1}{2\pi }\int _{-\infty }^{-{{E}_{k^{\prime} ,\alpha }}}{\rm d}E\int _{-\infty }^{+\infty }{\rm d}t\;{{{\rm e}}^{(-{\rm i}t+\epsilon )E}}{{{\rm e}}^{-\frac{{{t}^{2}}}{2q}}}dt \right| \\ {} & {} & =\left| \langle {{\Phi }_{k^{\prime} ,\alpha }}|{{A}_{X}}|{{\Psi }_{0}}\rangle \right|\left| \int _{-\infty }^{-{{E}_{k^{\prime} ,\alpha }}}{\rm d}E\sqrt{\frac{q}{2\pi }}{{{\rm e}}^{-\frac{q{{E}^{2}}}{2}}} \right| \\ {} & {} & \leqslant {{c}_{{\rm erf}}}{\rm exp} \left( -\frac{qE_{k^{\prime} ,\alpha }^{2}}{2} \right)\left| \langle {{\Phi }_{k^{\prime} ,\alpha }}|{{A}_{X}}|{{\Psi }_{0}}\rangle \right| \\ \end{array}\end{eqnarray} \tag{ A.10 }$$

and

$$\begin{eqnarray*}\begin{array}{ccccccccccccccc} |\langle {{\Psi }_{0}}|{{A}_{X}}-{{\tilde{A}}_{X}}|{{\Phi }_{k^{\prime} ,\alpha }}\rangle | & = & |\langle {{\Psi }_{0}}|{{A}_{X}}|{{\Phi }_{k^{\prime} ,\alpha }}\rangle |\left| 1-\frac{1}{2\pi }\int _{-\infty }^{+\infty }{\rm d}t\;\frac{{{{\rm e}}^{-{\rm i}{{E}_{k^{\prime} ,\alpha }}t-\frac{{{t}^{2}}}{2q}}}}{-it+\epsilon } \right| \\ {} & = & |\langle {{\Psi }_{0}}|{{A}_{X}}|{{\Phi }_{k^{\prime} ,\alpha }}\rangle |\left| 1-\frac{1}{2\pi }\int _{-\infty }^{{{E}_{k^{\prime} ,\alpha }}}{\rm d}E\int _{-\infty }^{+\infty }{\rm d}t\;{{{\rm e}}^{(-{\rm i}t+\epsilon )E}}{{{\rm e}}^{-\frac{{{t}^{2}}}{2q}}} \right| \\ \end{array}\end{eqnarray*}$$

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} {} & = & |\langle {{\Psi }_{0}}|{{A}_{X}}|{{\Phi }_{k^{\prime} ,\alpha }}\rangle |\left| \int _{{{E}_{k^{\prime} ,\alpha }}}^{+\infty }{\rm d}E\sqrt{\frac{q}{2\pi }}{{{\rm e}}^{-\frac{q{{E}^{2}}}{2}}} \right| \\ {} & {} & \leqslant {{c}_{{\rm erf}}}{\rm exp} \left( -\frac{qE_{k^{\prime} ,\alpha }^{2}}{2} \right)|\langle {{\Psi }_{0}}|{{A}_{X}}|{{\Phi }_{k^{\prime} ,\alpha }}\rangle |, \\ \end{array}\end{eqnarray} \tag{ A.11 }$$

where ${{c}_{{\rm erf}}}$ is a constant that allows one to bound the error function by a Gaussian. We can use these inequalities to show that

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} \parallel {{\tilde{A}}_{X}}|{{\Psi }_{0}}\rangle \parallel & = & {{\left( {{\int }_{\mathcal{B}}}{\rm d}k^{\prime} \mathop{\sum }\limits_{\alpha }|\langle {{\Phi }_{k^{\prime} ,\alpha }}|{{\tilde{A}}_{X}}|{{\Psi }_{0}}\rangle {{|}^{2}} \right)}^{1/2}} \\ {} & {} & \leqslant {{c}_{{\rm erf}}}\parallel {{A}_{X}}|{{\Psi }_{0}}\rangle \parallel \leqslant {{c}_{{\rm erf}}}\parallel {{A}_{X}}\parallel \\ \end{array}\end{eqnarray} \tag{ A.12 }$$

and similarly using the triangle inequality

$$\begin{eqnarray}&&\parallel \langle {{\Psi }_{0}}|{{\tilde{A}}_{X}}\parallel \leqslant \parallel \langle {{\Psi }_{0}}|({{\tilde{A}}_{X}}-{{A}_{X}})\parallel +\parallel \langle {{\Psi }_{0}}|{{A}_{X}}\parallel \leqslant ({{c}_{{\rm erf}}}+1)\parallel {{A}_{X}}\parallel .\end{eqnarray} \tag{ A.13 }$$

We now write ${{A}_{X}}{{\tilde{B}}_{Y}}(k)$ as $[{{\tilde{A}}_{X}},{{\tilde{B}}_{Y}}(k)]+({{A}_{X}}-{{\tilde{A}}_{X}}){{\tilde{B}}_{Y}}(k)+{{\tilde{B}}_{Y}}(k){{\tilde{A}}_{X}}$ and again use the triangle inequality to bound $|C(k)|$ as

$$\begin{eqnarray}&&|C(k)|\leqslant |\tilde{C}(k)|+|\langle {{\Psi }_{0}}|({{A}_{X}}-{{\tilde{A}}_{X}}){{\tilde{B}}_{Y}}(k)|{{\Psi }_{0}}\rangle |+|\langle {{\Psi }_{0}}|{{\tilde{B}}_{Y}}(k){{\tilde{A}}_{X}}|{{\Psi }_{0}}\rangle |,\end{eqnarray} \tag{ A.14 }$$

where we have defined a new correlator, $\tilde{C}(k)=\langle {{\Psi }_{0}}|[{{\tilde{A}}_{X}},{{\tilde{B}}_{Y}}(k)]|{{\Psi }_{0}}\rangle$.

For both the second and third term on the right-hand side of (A.14), we introduce a resolution of the identity, which we separate into two parts: one coming from momentum sectors with momentum k' satisfying $\parallel k^{\prime} -k\parallel \leqslant \delta$ (for the second term) or $\parallel k^{\prime} +k\parallel \leqslant \delta$ (for the third term), and one coming from the rest. For the latter contribution, we use the Cauchy-Schwarz inequality to write, for example, the third term

$$\begin{eqnarray*}&&\begin{array}{ccccccccccccccc} \left| {{\int }_{\parallel k^{\prime} +k\parallel \gt \delta }}{\rm d}k^{\prime} \mathop{\sum }\limits_{\alpha }\;\langle {{\Psi }_{0}}|{{\tilde{B}}_{Y}}(k)|{{\Phi }_{k^{\prime} ,\alpha }}\rangle \langle {{\Phi }_{k^{\prime} ,\alpha }}|{{\tilde{A}}_{X}}|{{\Psi }_{0}}\rangle \right| \\ \leqslant {{\left( {{\int }_{\parallel k^{\prime} +k\parallel \gt \delta }}{\rm d}k^{\prime} \mathop{\sum }\limits_{\alpha }\;|\langle {{\Psi }_{0}}|{{\tilde{B}}_{Y}}(k)|{{\Phi }_{k^{\prime} ,\alpha }}\rangle {{|}^{2}} \right)}^{1/2}} \\ \times {{\left( {{\int }_{\parallel k^{\prime} +k\parallel \gt \delta }}{\rm d}k^{\prime} \mathop{\sum }\limits_{\alpha }\;|\langle {{\Phi }_{k^{\prime} ,\alpha }}|{{\tilde{A}}_{X}}|{{\Psi }_{0}}\rangle {{|}^{2}} \right)}^{1/2}}. \\ \end{array}\end{eqnarray*}$$

By observing that

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} |\langle {{\Psi }_{0}}|{{\tilde{B}}_{Y}}(k)|{{\Phi }_{k^{\prime} ,\alpha }}\rangle | & = & |\langle {{\Psi }_{0}}|{{B}_{Y}}|{{\Phi }_{k^{\prime} ,\alpha }}\rangle |\left| {{N}_{r}}\mathop{\sum }\limits_{x\in \Lambda }{{{\rm e}}^{{\rm i}(k^{\prime} +k)\cdot x-\frac{\parallel x{{\parallel }^{2}}}{2r}}} \right| \\ {} & {} & \leqslant {{c}_{{\rm gauss}}}{{{\rm e}}^{-\frac{r\parallel k+k^{\prime} {{\parallel }^{2}}}{2}}}|\langle {{\Psi }_{0}}|{{B}_{Y}}|{{\Phi }_{k^{\prime} ,\alpha }}\rangle |, \\ \end{array}\end{eqnarray} \tag{ A.15 }$$

where ${{c}_{{\rm gauss}}}$ is defined in A.2, the first factor can be bounded as

$$\begin{eqnarray*}&&\begin{array}{ccccccccccccccc} {{\int }_{\parallel k^{\prime} +k\parallel \gt \delta }}{\rm d}k^{\prime} \mathop{\sum }\limits_{\alpha }\;|\langle {{\Psi }_{0}}|{{\tilde{B}}_{Y}}(k)|{{\Phi }_{k^{\prime} ,\alpha }}\rangle {{|}^{2}} \\ \leqslant c_{{\rm gauss}}^{2}{{\int }_{\parallel k^{\prime} +k\parallel \gt \delta }}{\rm d}k^{\prime} \mathop{\sum }\limits_{\alpha }\;{{{\rm e}}^{-r\parallel k+k^{\prime} {{\parallel }^{2}}}}|\langle {{\Psi }_{0}}|{{B}_{Y}}|{{\Phi }_{k^{\prime} ,\alpha }}\rangle {{|}^{2}} \\ \leqslant c_{{\rm gauss}}^{2}\;{{{\rm e}}^{-r{{\delta }^{2}}}}{{\int }_{\parallel k^{\prime} +k\parallel \gt \delta }}{\rm d}k^{\prime} \mathop{\sum }\limits_{\alpha }\;|\langle {{\Psi }_{0}}|{{B}_{Y}}|{{\Phi }_{k^{\prime} ,\alpha }}\rangle {{|}^{2}} \\ \leqslant c_{{\rm gauss}}^{2}\;{{{\rm e}}^{-r{{\delta }^{2}}}}\parallel {{B}_{Y}}{{\parallel }^{2}}, \\ \end{array}\end{eqnarray*}$$

whereas for the second factor we use

$$\begin{eqnarray*}&&\begin{array}{ccccccccccccccc} {{\int }_{\parallel k^{\prime} +k\parallel \gt \delta }}{\rm d}k^{\prime} \mathop{\sum }\limits_{\alpha }\;|\langle {{\Phi }_{k^{\prime} ,\alpha }}|{{\tilde{A}}_{X}}|{{\Psi }_{0}}\rangle {{|}^{2}} \\ \end{array}\end{eqnarray*}$$

$$\begin{eqnarray*}&&\begin{array}{ccccccccccccccc} \leqslant c_{{\rm erf}}^{2}{{\int }_{\parallel k^{\prime} +k\parallel \gt \delta }}{\rm d}k^{\prime} \mathop{\sum }\limits_{\alpha }\;{{{\rm e}}^{-qE_{k^{\prime} ,\alpha }^{2}}}|\langle {{\Phi }_{k^{\prime} ,\alpha }}|{{A}_{X}}|{{\Psi }_{0}}\rangle {{|}^{2}} \\ \leqslant c_{{\rm erf}}^{2}\;\parallel {{A}_{X}}{{\parallel }^{2}}. \\ \end{array}\end{eqnarray*}$$

From the above, we thus conclude that

$$\begin{eqnarray}&&\begin{array}{ccccccccccccccc} \left| {{\int }_{\parallel k^{\prime} +k\parallel \gt \delta }}{\rm d}k^{\prime} \mathop{\sum }\limits_{\alpha }\;\langle {{\Psi }_{0}}|{{\tilde{B}}_{Y}}(k)|{{\Phi }_{k^{\prime} ,\alpha }}\rangle \langle {{\Phi }_{k^{\prime} ,\alpha }}|{{\tilde{A}}_{X}}|{{\Psi }_{0}}\rangle \right| \\ \leqslant {{c}_{{\rm gauss}}}\;{{c}_{{\rm erf}}}\;{{{\rm e}}^{-\frac{r{{\delta }^{2}}}{2}}}\parallel {{A}_{X}}\parallel \parallel {{B}_{Y}}\parallel . \\ \end{array}\end{eqnarray} \tag{ A.16 }$$

An identical contribution is obtained from the momentum region, $\parallel k^{\prime} -k\parallel \gt \delta ,$ if the resolution of the identity is inserted in the second term of the right-hand side of (A.14).

For the momentum region, $\parallel k^{\prime} +k\parallel \leqslant \delta ,$ in the third term, we use the same approach with the Cauchy-Schwarz inequality, but we now bound the first factor by

$$\begin{eqnarray*}&&{{\int }_{\parallel k^{\prime} +k\parallel \lt \delta }}{\rm d}k^{\prime} \mathop{\sum }\limits_{\alpha }\;|\langle {{\Psi }_{0}}|{{\tilde{B}}_{Y}}(k)|{{\Phi }_{k^{\prime} ,\alpha }}\rangle {{|}^{2}}\leqslant c_{{\rm gauss}}^{2}\parallel {{B}_{Y}}{{\parallel }^{2}}\end{eqnarray*}$$

and the second factor by

$$\begin{eqnarray*}&&\begin{array}{ccccccccccccccc} {{\int }_{\parallel k^{\prime} +k\parallel \gt \delta }}{\rm d}k^{\prime} \mathop{\sum }\limits_{\alpha }\;|\langle {{\Phi }_{k^{\prime} ,\alpha }}|{{\tilde{A}}_{X}}|{{\Psi }_{0}}\rangle {{|}^{2}} \\ \leqslant c_{{\rm erf}}^{2}{{\int }_{\parallel k^{\prime} +k\parallel \gt \delta }}{\rm d}k^{\prime} \mathop{\sum }\limits_{\alpha }\;{{{\rm e}}^{-qE_{k^{\prime} ,\alpha }^{2}}}|\langle {{\Phi }_{k^{\prime} ,\alpha }}|{{A}_{X}}|{{\Psi }_{0}}\rangle {{|}^{2}} \\ \leqslant c_{{\rm erf}}^{2}{{{\rm e}}^{-q{{E}^{*}}{{(-k,\delta )}^{2}}}}\parallel {{A}_{X}}{{\parallel }^{2}}, \\ \end{array}\end{eqnarray*}$$

where ${{E}^{*}}(k,\delta )={{{\rm min} }_{\parallel k^{\prime} -k\parallel \leqslant \delta }}E(k^{\prime} )$. Again, the contribution of $\parallel k^{\prime} -k\parallel \lt \delta$ in the second term is evaluated completely analogously.

We can thus rewrite (A.14) as

$$\begin{eqnarray}&&|C(k)|\leqslant |\tilde{C}(k)|+{{c}_{{\rm gauss}}}\;{{c}_{{\rm erf}}}\left[ 2{{{\rm e}}^{-\frac{r{{\delta }^{2}}}{2}}}+{{{\rm e}}^{-\frac{q{{E}^{*}}{{(-k,\delta )}^{2}}}{2}}}+{{{\rm e}}^{-\frac{q{{E}^{*}}{{(k,\delta )}^{2}}}{2}}} \right].\end{eqnarray} \tag{ A.17 }$$

To bound the new correlator $\tilde{C}(k)$, we replace the wave packet, ${{\tilde{B}}_{Y}}(k),$ with a completely local version, ${\tilde{\tilde{B}}}_{Y}(k),$ defined as

$$\begin{eqnarray}&&{{\tilde{\tilde{B}}}_{Y}}(k)=\mathop{\sum }\limits_{\frac{x\in \Lambda }{\parallel x\parallel \leqslant \ell }}{{N}_{r}}{{{\rm e}}^{-\frac{\parallel x{{\parallel }^{2}}}{2r}}}{{{\rm e}}^{-{\rm i}k\cdot x}}{{T}^{x}}{{B}_{Y}}{{T}^{-x}},\end{eqnarray} \tag{ A.18 }$$

and thus define yet another correlator , ${\tilde{\tilde{C}}}(k),$ as

$$\begin{eqnarray}&&{\tilde{\tilde{C}}}(k)=\langle {{\Psi }_{0}}|[{{\tilde{A}}_{X}},{{\tilde{\tilde{B}}}_{Y}}(k)]|{{\Psi }_{0}}\rangle .\end{eqnarray} \tag{ A.19 }$$

The error in the operator norm can be bounded by the triangle inequality as

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} \parallel {{\tilde{B}}_{Y}}(k)-{{\tilde{\tilde{B}}}_{Y}}(k)\parallel & = & \parallel \mathop{\sum }\limits_{\frac{x\in \Lambda }{\parallel x\parallel \gt \ell }}{{N}_{r}}{{{\rm e}}^{-\frac{\parallel x{{\parallel }^{2}}}{2r}}}{{{\rm e}}^{-{\rm i}k\cdot x}}{{T}^{x}}{{B}_{Y}}{{T}^{-x}}\parallel , \\ {} & {} & \leqslant {{N}_{r}}\parallel {{B}_{Y}}\parallel \mathop{\sum }\limits_{\frac{x\in \Lambda }{\parallel x\parallel \gt \ell }}{\rm exp} \left[ -\frac{\parallel x{{\parallel }^{2}}}{2r} \right] \\ {} & {} & \leqslant {{c}_{\Lambda }}\parallel {{B}_{Y}}\parallel {\rm exp} \left( -\frac{{{\ell }^{2}}}{2r} \right) \\ \end{array}\end{eqnarray} \tag{ A.20 }$$

where an accurate determination of ${{c}_{\Lambda }}$ requires detailed knowledge about the structure of the lattice, Λ. We can thus write

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} |\tilde{C}(k)| & \leqslant & |{\tilde{\tilde{C}}}(k)|+|\langle {{\Psi }_{0}}|{{\tilde{A}}_{Y}}\left( {{\tilde{B}}_{Y}}(k)-{{\tilde{\tilde{B}}}_{Y}}(k) \right)|{{\Psi }_{0}}\rangle | \\ {} & {} & +\langle {{\Psi }_{0}}|\left( {{\tilde{B}}_{Y}}(k)-{{\tilde{\tilde{B}}}_{Y}}(k) \right){{\tilde{A}}_{Y}}|{{\Psi }_{0}}\rangle | \\ \end{array}\end{eqnarray} \tag{ A.21 }$$

and, using the Cauchy-Schwarz inequality in the different terms and (A.12) and (A.13), we obtain

$$\begin{eqnarray}&&|\tilde{C}(k)|\leqslant |{\tilde{\tilde{C}}}(k)|+{{c}_{\Lambda }}(2{{c}_{{\rm erf}}}+1){{{\rm e}}^{-\frac{{{\ell }^{2}}}{2r}}}\parallel {{A}_{X}}\parallel \parallel {{B}_{Y}}\parallel .\end{eqnarray} \tag{ A.22 }$$

Finally, to bound ${\tilde{\tilde{C}}}(k)$, we separate the time integral in the definition of ${{\tilde{A}}_{X}}$ into two pieces. For $t\lt T$, we obtain

$$\begin{eqnarray}&&\begin{array}{ccccccccccccccc} \left| \frac{1}{2\pi }\int _{-T}^{+T}{\rm d}t\;\frac{\langle {{\Psi }_{0}}|[{{{\rm e}}^{{\rm i}Ht}}{{A}_{X}}{{{\rm e}}^{-{\rm i}Ht}},{{\tilde{\tilde{B}}}_{Y}}(k)]|{{\Psi }_{0}}\rangle }{-{\rm i}t+\epsilon }{{{\rm e}}^{-\frac{{{t}^{2}}}{2q}}} \right| \\ \leqslant \frac{2\parallel {{A}_{X}}\parallel \parallel {{\tilde{\tilde{B}}}_{Y}}(k)\parallel |X|{{{\rm e}}^{-\mu ({\rm dist}(X,Y)-\ell )}}}{2\pi }\int _{-T}^{+T}{\rm d}t\;\frac{{{{\rm e}}^{2\,s|t|}}-1}{-{\rm i}t+\epsilon }{{{\rm e}}^{-\frac{{{t}^{2}}}{2q}}} \\ \leqslant \frac{4\parallel {{A}_{X}}\parallel \parallel {{\tilde{\tilde{B}}}_{Y}}(k)\parallel |X|{{{\rm e}}^{-\mu ({\rm dist}(X,Y)-\ell )}}}{2\pi }\int _{0}^{T}{\rm d}t\;\frac{{{{\rm e}}^{2sT}}-1}{t} \\ \leqslant \frac{4\parallel {{A}_{X}}\parallel \parallel {{\tilde{\tilde{B}}}_{Y}}(k)\parallel |X|{{{\rm e}}^{-\mu ({\rm dist}(X,Y)-\ell )}}}{2\pi }\int _{0}^{T}{\rm d}t\;2\;s{{{\rm e}}^{2sT}} \\ \leqslant \frac{2}{\pi }\parallel {{A}_{X}}\parallel \parallel {{B}_{Y}}\parallel |X|{{{\rm e}}^{2sT-\mu ({\rm dist}(X,Y)-\ell )}}, \\ \end{array}\end{eqnarray} \tag{ A.23 }$$

where we have also used

$$\begin{eqnarray*}&&\parallel {{\tilde{\tilde{B}}}_{Y}}(k)\parallel \leqslant \parallel {{B}_{Y}}\parallel |{{N}_{r}}\mathop{\sum }\limits_{\frac{x\in \Lambda }{\parallel x\parallel \leqslant \ell }}{{{\rm e}}^{-\frac{\parallel x{{\parallel }^{2}}}{2r}}}|\leqslant \parallel {{B}_{Y}}\parallel .\end{eqnarray*}$$

For the contribution of $|t|\gt T$, we use the Gaussian in the integrand to bound

$$\begin{eqnarray*}\begin{array}{ccccccccccccccc} \parallel \tilde{A}_{X}^{\gt }\parallel & = & \left\|\frac{1}{2\pi }{{\int }_{|t|\gt T}}{\rm d}t\;\frac{{{{\rm e}}^{{\rm i}Ht}}{{A}_{X}}{{{\rm e}}^{-{\rm i}Ht}}}{-{\rm i}t+\epsilon }{{{\rm e}}^{-\frac{{{t}^{2}}}{2q}}}\right\| \\ {} & {} & \leqslant \frac{\parallel {{A}_{X}}\parallel }{\pi T}\int _{T}^{+\infty }{\rm d}t\;{{{\rm e}}^{-\frac{{{t}^{2}}}{2q}}}\leqslant \parallel {{A}_{X}}\parallel \sqrt{\frac{2q}{\pi }}\frac{{{c}_{{\rm erf}}}}{T}{{{\rm e}}^{-\frac{{{T}^{2}}}{2q}}}. \\ \end{array}\end{eqnarray*}$$

The bound on ${\tilde{\tilde{C}}}(k)$ is thus given by

$$\begin{eqnarray}&&|{\tilde{\tilde{C}}}(k)|\leqslant \parallel {{A}_{X}}\parallel \parallel {{B}_{Y}}\parallel \left[ \frac{2}{\pi }|X|{{{\rm e}}^{2sT-\mu ({\rm dist}(X,Y)-\ell )}}+\frac{2{{c}_{{\rm erf}}}}{T}\sqrt{\frac{2q}{\pi }}{{{\rm e}}^{-\frac{{{T}^{2}}}{2q}}} \right].\end{eqnarray} \tag{ A.24 }$$

Finally, putting everything together, we obtain

$$\begin{eqnarray}&&\begin{array}{ccccccccccccccc} \frac{|C(k)|}{\parallel {{A}_{X}}\parallel \parallel {{B}_{Y}}\parallel }\leqslant {{c}_{{\rm gauss}}}\;{{c}_{{\rm erf}}}\left[ 2{{{\rm e}}^{-\frac{r{{\delta }^{2}}}{2}}}+{{{\rm e}}^{-\frac{q{{E}^{*}}{{(-k,\delta )}^{2}}}{2}}}+{{{\rm e}}^{-\frac{q{{E}^{*}}{{(k,\delta )}^{2}}}{2}}} \right] \\ +{{c}_{\Lambda }}(2{{c}_{{\rm erf}}}+1){{{\rm e}}^{-\frac{{{\ell }^{2}}}{2r}}}+\frac{2}{\pi }|X|{{{\rm e}}^{2sT-\mu ({\rm dist}(X,Y)-\ell )}}+\frac{2{{c}_{{\rm erf}}}}{T}\sqrt{\frac{2r}{\pi }}{{{\rm e}}^{-\frac{{{T}^{2}}}{2q}}}. \\ \end{array}\end{eqnarray} \tag{ A.25 }$$

In this expression, we can tune the constants r, δ, q, $\ell ,$ and T as a function of ${\rm dist}(X,Y),$ and the characteristics of E(k) around k (i.e., ${{E}^{*}}(k,\delta )$). It is clear that we want to impose the restrictions:

$$\begin{eqnarray}&&\ell \lt {\rm dist}(X,Y),2sT\leqslant \mu ({\rm dist}(X,Y)-\ell ).\end{eqnarray} \tag{ A.26 }$$

Let us start by setting $\ell =\alpha {\rm dist}(X,Y)$ with $\alpha \lt 1$ and $T={\rm dist}(X,Y)/v$ with $v\gt \frac{2s}{\mu (1-\alpha )}$. In addition, we choose $q=\beta \;{\rm dist}(X,Y)$ and $r=\gamma \;{\rm dist}(X,Y)$, so that

$$\begin{eqnarray}&&\begin{array}{ccccccccccccccc} \frac{|C(k)|}{\parallel {{A}_{X}}\parallel \parallel {{B}_{Y}}\parallel } \\ \leqslant {{c}_{{\rm gauss}}}\;{{c}_{{\rm erf}}}\left[ 2{{{\rm e}}^{-\frac{\gamma {{\delta }^{2}}}{2}{\rm dist}(X,Y)}}+{{{\rm e}}^{-\frac{\beta {{E}^{*}}{{(-k,\delta )}^{2}}}{2}{\rm dist}(X,Y)}}+{{{\rm e}}^{-\frac{\beta {{E}^{*}}{{(k,\delta )}^{2}}}{2}{\rm dist}(X,Y)}} \right] \\ +{{c}_{\Lambda }}(2{{c}_{{\rm erf}}}+1){{{\rm e}}^{-\frac{{{\alpha }^{2}}}{2\gamma }{\rm dist}(X,Y)}}+\frac{2}{\pi }|X|{{{\rm e}}^{-[\mu (1-\alpha )-\frac{2\,s}{v}]{\rm dist}(X,Y)}} \\ +\frac{2{{c}_{{\rm erf}}}}{T}\sqrt{\frac{2r}{\pi }}{{{\rm e}}^{-\frac{1}{2\beta {{v}^{2}}}{\rm dist}(X,Y)}} \\ \end{array}\end{eqnarray} \tag{ A.27 }$$

We can now fine-tune the remaining constant to have a similar decay in all exponentials. Assuming the system has reflection invariance (i.e., ${{E}^{*}}(k,\delta )={{E}^{*}}(-k,\delta )$) we can therefore choose the constants such that

$$\begin{eqnarray}&&\mu (1-\alpha )-\frac{2s}{v}=\frac{1}{2\beta {{v}^{2}}}=\frac{{{\alpha }^{2}}}{2\gamma }=\frac{\gamma {{\delta }^{2}}}{2}=\frac{\beta {{E}^{*}}{{(k,\delta )}^{2}}}{2}.\end{eqnarray} \tag{ A.28 }$$

If there is no reflection invariance, the smaller of both ${{E}^{*}}(k,\delta )$ and ${{E}^{*}}(-k,\delta )$ will determine the slowest exponential decay and should be used in the equations above.

Clearly, we should set $\beta ={{[v{{E}^{*}}(k,\delta )]}^{-1}}$. We then obtain $\mu (1-\alpha )-2\;s/v={{E}^{*}}(k,\delta )/(2v),$ from which we determine the optimal velocity, $v=[4s+{{E}^{*}}(k,\delta )]/[\mu (1-\alpha )]$. From the second equality we can fix $\gamma =\alpha /\delta$. The remaining equation is therefore ${{E}^{*}}(k,\delta )/(2v)=\alpha \delta /2$. Inserting the velocity, we obtain

$$\begin{eqnarray*}&&2\mu (1-\alpha ){{E}^{*}}(k,\delta )=\alpha \delta [4s+{{E}^{*}}(k,\delta )]\end{eqnarray*}$$

which determines α as

$$\begin{eqnarray}&&{{\alpha }^{-1}}=1+\frac{\delta }{2\mu }+\frac{2\delta s}{{{E}^{*}}(k,\delta )\mu }\end{eqnarray} \tag{ A.29 }$$

and only leaves δ to be determined. Note that the restriction $\alpha \lt 1$ is satisfied. The fastest exponential decay is obtained by maximizing

$$\begin{eqnarray}&&\alpha \delta =\delta {{\left( 1+\frac{\delta }{2\mu }+\frac{2\delta s}{{{E}^{*}}(k,\delta )\mu } \right)}^{-1}}={{\left( \frac{1}{\delta }+\frac{2s}{\mu {{E}^{*}}(k,\delta )}+\frac{1}{2\mu } \right)}^{-1}}.\end{eqnarray} \tag{ A.30 }$$

If we want to minimize the denominator, there is a clear tradeoff, since increasing δ decreases the first term and increases the second. If k corresponds to a minimum of the dispersion relation E(k), then the function ${{E}^{*}}(k,\delta )$ will be insensitive to δ in some region, and we can choose δ to be as large as possible within this region. However, if k corresponds to a regular point where $\frac{{\rm d}E(k)}{{\rm d}k}(k)\ne 0$, then there is a direct effect from increasing k to decreasing ${{E}^{*}}(k,\delta )$.

A more intuitive result is obtained if we treat the term coming from the Lieb-Robinson bound separately, as the decay properties of this term are specific to the details of the Hamiltonian. Let us assume that we only know about the existence of some maximal velocity of propagation ${{v}_{{\rm LR}}}$, such that for any ${{v}_{{\rm LR}}}|t|\leqslant {\rm dist}(X,Y)$, we can write

$$\begin{eqnarray}&&\parallel [{{A}_{X}}(t),{{B}_{Y}}]\parallel \leqslant {{c}_{{\rm LR}}}\parallel {{A}_{X}}\parallel \parallel {{B}_{Y}}\parallel {\rm exp} \left( -\frac{{\rm dist}(X,Y)}{\xi } \right).\end{eqnarray} \tag{ A.31 }$$

Clearly, choosing, ${{v}_{{\rm LR}}}$ results in a smaller ξ (a quicker exponential decay of the Lieb-Robinson bound), and vice versa. For the Lieb-Robinson bound to be the smallest error in the proof in [3] of the exponential decay of correlations, we need to choose ${{v}_{{\rm LR}}}$ to be at least large enough such that

$$\begin{eqnarray}&&\frac{1}{\xi }\geqslant \frac{\Delta E}{2{{v}_{{\rm LR}}}}.\end{eqnarray} \tag{ A.32 }$$

If we now use the new Lieb-Robinson bound of (A.31) above, we would again choose $\ell =\alpha \;{\rm dist}(X,Y),$ but we would need to fix $T=(1-\alpha ){\rm dist}(X,Y)/{{v}_{{\rm LR}}}$. Hence, we use a fixed velocity of propagation and we do not optimize over it (essentially, the parameter μ in (A.30)). We will again try to have an equal decay in all exponentials, except for the one coming from the Lieb-Robinson bound, which we allow to decay faster. We thus obtain

$$\begin{eqnarray}&&\frac{1-\alpha }{\xi }\geqslant \frac{{{(1-\alpha )}^{2}}}{2\beta v_{{\rm LR}}^{2}}=\frac{{{\alpha }^{2}}}{2\gamma }=\frac{\gamma {{\delta }^{2}}}{2}=\frac{\beta {{E}^{*}}{{(k,\delta )}^{2}}}{2}.\end{eqnarray} \tag{ A.33 }$$

Setting $\gamma =\alpha /\delta$ and $\beta =(1-\alpha )/[{{v}_{{\rm LR}}}{{E}^{*}}(k,\delta )]$ reduces these equations down to

$$\begin{eqnarray}&&\frac{1-\alpha }{\xi }\geqslant \frac{(1-\alpha ){{E}^{*}}(k,\delta )}{2{{v}_{{\rm LR}}}}=\frac{\alpha \delta }{2}.\end{eqnarray} \tag{ A.34 }$$

Clearly, the first inequality is trivially satisfied, since ${{E}^{*}}(k,\delta )\geqslant \Delta E$ and the fixed velocity, ${{v}_{{\rm LR}}},$ satisfies (A.32). From the last equation, we obtain

$$\begin{eqnarray}&&{{\alpha }^{-1}}=1+\frac{{{v}_{{\rm LR}}}\delta }{{{E}^{*}}(k,\delta )}\end{eqnarray} \tag{ A.35 }$$

and the rate of the exponential decay in ${\rm dist}(X,Y)$ is given by

$$\begin{eqnarray}&&\alpha \delta ={{\left( \frac{1}{\delta }+\frac{{{v}_{{\rm LR}}}}{{{E}^{*}}(k,\delta )} \right)}^{-1}}.\end{eqnarray} \tag{ A.36 }$$

We can then optimize over δ to find an optimal decay.

We compute the discrete Fourier transform of a sampled Gaussian by inserting the inverse continuous Fourier transform

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} G(k) & = & {{N}_{r}}\mathop{\sum }\limits_{x\in \Lambda }{{{\rm e}}^{-{\rm i}k\cdot x}}{\rm exp} \left( -\frac{\parallel x{{\parallel }^{2}}}{2r} \right) \\ {} & = & {{N}_{r}}{{\left( \frac{r}{2\pi } \right)}^{d/2}}{{\int }_{{{\mathbb{R}}^{d}}}}{\rm d}k^{\prime} \;{\rm exp} \left( -\frac{r\parallel k^{\prime} {{\parallel }^{2}}}{2} \right)\mathop{\sum }\limits_{x\in \Lambda }{{{\rm e}}^{{\rm i}(k^{\prime} -k)\cdot x}} \\ {} & = & {{N}_{r}}{{\left( \frac{r}{2\pi } \right)}^{d/2}}{{\int }_{{{\mathbb{R}}^{d}}}}{\rm d}k^{\prime} \;{\rm exp} \left( -\frac{r\parallel k^{\prime} {{\parallel }^{2}}}{2} \right){{V}_{\mathcal{B}}}\mathop{\sum }\limits_{K\in \tilde{\Lambda }}\delta (k^{\prime} -k+K) \\ {} & = & {{N}_{r}}{{\left( \frac{r}{2\pi } \right)}^{d/2}}{{V}_{\mathcal{B}}}\mathop{\sum }\limits_{K\in \tilde{\Lambda }}{\rm exp} \left( -\frac{r\parallel k-K{{\parallel }^{2}}}{2} \right) \\ \end{array}\end{eqnarray} \tag{ A.37 }$$

with $\tilde{\Lambda }$ being the reciprocal lattice and ${{V}_{\mathcal{B}}},$ being the volume of the Brillouin zone (which is the unit cell of $\tilde{\Lambda }$). Hence, the Fourier transform of the sampled Gaussian is a sum of Gaussians centered around the different lattice points of the reciprocal lattice, $\tilde{\Lambda }$. We are only interested in the value of G(k) for $k\in \mathcal{B}$, so if r is sufficiently large, the contributions of the Gaussians around the points $K\ne 0$ will be very small. In general, there exists a constant ${{c}_{{\rm gauss}}},$ so that we can bound G(k) by

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} |G(k)| & = & {{N}_{r}}{{\left( \frac{r}{2\pi } \right)}^{d/2}}{{V}_{\mathcal{B}}}{\rm exp} (-r\parallel k{{\parallel }^{2}}/2)\mathop{\sum }\limits_{K\in \tilde{\Lambda }}{\rm exp} \left( rk\cdot K-\frac{r\parallel K{{\parallel }^{2}}}{2} \right) \\ {} & {} & \leqslant {{N}_{r}}{{\left( \frac{r}{2\pi } \right)}^{d/2}}{{V}_{\mathcal{B}}}{\rm exp} (-r\parallel k{{\parallel }^{2}}/2)\mathop{\sum }\limits_{K\in \tilde{\Lambda }}{\rm exp} \left[ -\frac{r\parallel K{{\parallel }^{2}}}{2}+r\parallel K\parallel \mathop{{\rm max} }\limits_{k\in \mathcal{B}}\parallel k\parallel \right] \\ {} & {} & \leqslant {{c}_{{\rm gauss}}}{\rm exp} (-r\parallel k{{\parallel }^{2}}/2). \\ \end{array}\end{eqnarray} \tag{ A.38 }$$