Explicit Lieb-Robinson Velocity Bounds

Question

Lieb-Robinson bounds describe how effects are propagated through a system due to a local Hamiltonian. They are often described in the form $$ \left|[A,B(t)]\right|\leq Ce^{vt-l}, $$ where $A$ and $B$ are operators that are separated by a distance $l$ on a lattice where the Hamiltonian has local (e.g. nearest neighbour) interactions on that lattice, bounded by some strength $J$. The proofs of the Lieb Robinson bound typically show the existence of a velocity $v$ (that depends on $J$). This is often really useful for bounding properties in these systems. For example, there were some really nice results here regarding how long it takes to generate a GHZ state using a nearest-neighbour Hamiltonian.

The problem that I've had is that the proofs are sufficiently generic that it is difficult to get a tight value on what the velocity actually is for any given system.

To be specific, imagine a one dimensional chain of qubits coupled by a Hamiltonian $$ H=\sum_{n=1}^N\frac{B_n}{2}Z_n+\sum_{n=1}^{N-1}\frac{J_n}{2}(X_nX_{n+1}+Y_nY_{n+1}), \tag{1} $$ where $J_n\leq J$ for all $n$. Here $X_n$, $Y_n$ and $Z_n$ represent a Pauli operator being applied to a given qubit $n$, and $\mathbb{I}$ everywhere else. Can you give a good (i.e. as tight as possible) upper bound for the Lieb-Robinson velocity $v$ for the system in Eq. (1)?

This question can be asked under two different assumptions:

• The $J_n$ and $B_n$ are all fixed in time
• The $J_n$ and $B_n$ are allowed to vary in time.

The former is a stronger assumption which may make proofs easier, while the latter is usually included in the statement of Lieb-Robinson bounds.

---

Motivation

Quantum computation, and more generally quantum information, comes down to making interesting quantum states. Through works such as this, we see that information takes a certain amount of time to propagate from one place to another in a quantum system undergoing evolution due to a Hamiltonian such as in Eq. (1), and that quantum states, such as GHZ states, or states with a topological order, take a certain amount of time to produce. What the result currently shows is a scaling relation, e.g. the time required is $\Omega(N)$.

So, let's say I come up with a scheme that does information transfer, or produces a GHZ state etc. in a way that scales linearly in $N$. How good is that scheme actually? If I have an explicit velocity, I can see how closely matched the scaling coefficient is in my scheme compared to the lower bound.

If I think that one day what I want to see is a protocol implemented in the lab, then I very much care about optimising these scaling coefficients, not just the broad scaling functionality, because the faster I can implement a protocol, the less chance there is for noise to come along and mess everything up.

---

Further Information

$\newcommand{\bra}[1]{\left<#1\right|}\newcommand{\ket}[1]{\left|#1\right>}\newcommand{\bk}[2]{\left<#1\middle|#2\right>}\newcommand{\bke}[3]{\left<#1\middle|#2\middle|#3\right>}\newcommand{\proj}[1]{\left|#1\right>\left<#1\right|}$There are some nice features of this Hamiltonian which I assume make calculation easier. In particular, the Hamiltonian has a subspace structure based on the number of 1s in the standard basis (it is said to be excitation preserving) and, even better, the Jordan-Wigner transformation shows that all the properties of higher excitation subspaces can be derived from the 1-excitation subspace. This essentially means we only have to do the maths on an $N\times N$ matrix $h$ instead of the full $2^N\times 2^N$ matrix $H$, where $$ h=\sum_{n=1}^NB_n\proj{n}+\sum_{n=1}^{N-1}J_n(\ket{n}\bra{n+1}+\ket{n+1}\bra{n}). $$ There is some evidence that the Lieb-Robinson velocity is $v=2J$, such as here and here, but these all use a close to uniformly coupled chain, which has a group velocity $2J$ (and I assume the group velocity is closely connected to the Lieb-Robinson velocity). It doesn't prove that all possible choices of coupling strength have a velocity that is so bounded.

I can add a bit further to the motivation. Consider the time evolution of a single excitation starting at one end of the chain, $\ket{1}$, and what its amplitude is for arriving at the other end of the chain $\ket{N}$, a short time $\delta t$ later. To first order in $\delta t$, this is $$ \bra{N}e^{-ih\delta t}\ket{1}=\frac{\delta t^{N-1}}{(N-1)!}\prod_{n=1}^{N-1}J_n+O(\delta t^{N}). $$ You can see the exponential functionality that you would expect being outside the 'light cone' defined by a Lieb-Robinson system, but more importantly, if you wanted to maximise that amplitude, you'd set all the $J_n=J$. So, at short times, the uniformly coupled system leads to the most rapid transfer. Trying to push this further, you can ask, as a bit of a fudge, when can $$ \frac{t^{N-1}}{(N-1)!}\prod_{n=1}^{N-1}J_n\sim 1 $$ Taking the large $N$ limit, and using Stirling's formula on the factorial leads to $$ \frac{etJ}{N-1}\sim 1, $$ which suggests a maximum velocity of about $eJ$. Close, but hardly rigorous (as the higher order terms are non-negligible)!

$$ \newcommand{\bra}[1]{\left<#1\right|}\newcommand{\ket}[1]{\left|#1\right>}\newcommand{\bk}[2]{\left<#1\middle|#2\right>}\newcommand{\bke}[3]{\left<#1\middle|#2\middle|#3\right>}\newcommand{\proj}[1]{\left|#1\right>\left<#1\right|} $$

Answer 1

Let me first answer the general question how to get a reasonably tight Lieb-Robinson (LR) speed when you are facing a generic locally interacting lattice model, and then I'll come back to the 1D XY model in your question, which is very special to be exactly solvable.

---

General Method

The method to obtain the tightest bound to date (for a generic short-range interacting model) is introduced in Ref1=arXiv:1908.03997. The basic idea is that the norm of the unequal time commutator $\|[A_X(t),B_Y(0)]\|$ between arbitrary local operators can be upper bounded by the solution to a set of first order linear differential equations living on the commutativity graph of the model. The commutativity graph, as introduced in Sec.II A of Ref1, can be easily drawn from the model Hamiltonian $\hat{H}$, and is designed to reflect the commutation relations between different local operators presented in $\hat{H}$. In translation invariant systems, this set of differential equations can be easily solved by a Fourier transform, and an upper bound of the LR speed can be calculated from the largest eigenfrequency $\omega_{\max}(i\vec{\kappa})$ using Eq.(31) of Ref1. In the following I'll apply this method to the 1D XY model as a pedagogical example. For simplicity, I'll focus on the time-independent and translation invariant case $|B_n|=B>0$, $|J_n|=J>0$ (the resulting bound doesn't depend on signs of $B_n,J_n$). For the translation non-invariant, time-dependent case, you can either solve the differential equation numerically (which is an easy computational task for systems of thousands of sites), or you can use an overall upper bound $|J_n(t)|\leq J,~|B_n(t)|\leq B$ and proceed to use the method below (but this slightly compromises tightness compared to the numerical method).

1. First we draw the commutativity graph, as below. Each operator in the Hamiltonian~($X_{n}X_{n+1}$, $Y_nY_{n+1}$, $Z_n$) is represented by a vertex, and we link two vertices if and only if the corresponding operators don't commute (or, in the current case, anti-commute).

2. Then write down the differential equations Eq.(10) of Ref1: \begin{eqnarray} \dot{\bar{\gamma}}_{\alpha,n}&=&J[\bar{\gamma}_{\alpha,n-1}(t)+\bar{\gamma}_{\alpha,n+1}(t)]+B[\bar{\gamma}_{3,n}(t)+\bar{\gamma}_{3,n+1}(t)],~~\alpha=1,2,\nonumber\\ \dot{\bar{\gamma}}_{3,n}&=&J\sum_{\alpha=1,2}[\bar{\gamma}_{\alpha,n-1}(t)+\bar{\gamma}_{\alpha,n}(t)]. \end{eqnarray}

3. Fourier transforming the above equation, we have \begin{equation} \frac{d}{dt}\begin{pmatrix} \bar{\gamma}_{1,k}\\ \bar{\gamma}_{2,k}\\ \bar{\gamma}_{3,k} \end{pmatrix}=\begin{pmatrix} 2J\cos k & 0 & B(1+e^{ik})\\ 0 & 2J\cos k & B(1+e^{ik})\\ J(1+e^{-ik}) & J(1+e^{-ik}) & 0\\ \end{pmatrix}\begin{pmatrix} \bar{\gamma}_{1,k}\\ \bar{\gamma}_{2,k}\\ \bar{\gamma}_{3,k} \end{pmatrix}. \end{equation} The eigenfrequencies are $2J\cos k,J\cos k\pm \sqrt{(J\cos k)^2+2BJ(1+\cos k)}$. The LR speed is given by Eq.(31) of Ref1: $$v_{\text{LR}}\leq \min_{\kappa>0}\frac{\omega_{\max}(i\kappa)}{\kappa} =\mathcal{Z}_{\frac{B}{2J}}J,$$ where \begin{equation}\label{eq:vLRFHWy} \mathcal{Z}_y\equiv \min_{\kappa>0}\frac{\cosh\kappa+\sqrt{\cosh^2\kappa+4y(1+\cosh\kappa)}}{\kappa}. \end{equation}

Note: This bound diverges when $B/J\to \infty$, while the physical information propagation speed stays finite. We can get rid of this problem by using the method in Sec. VI of Ref1. The result is $v_{\text{LR}}\leq 4\mathcal{X}_0 J$ in this limit, where $\mathcal{X}_y$ is defined as the solution to the equation $x\mathrm{arcsinh}(x)=\sqrt{x^2+1}+y$.

---

Velocity bounds for some classic models

The above method is completely general. In case you are interested in more, I listed the velocity bounds for some classic models in the following table, obtained in a similar way. Notice that the LR velocity $v_{\text{LR}}$ is upper bounded by the smallest of the all the expressions listed (so in different parameter regions different expressions should be used). The function $F(J_x,J_y,J_z)$ is defined as the largest root of $x^3-(J_xJ_y+J_xJ_z+J_yJ_z)x-2J_xJ_yJ_z=0.$ All parameters are assumed positive (just take absolute value for the negative cases).

$$\begin{array}{|c|l|} \hline \text{Model} & v_{\text{LR}} \\ \hline d\text{-dimensional TFIM} & 2\mathcal{X}_0\sqrt{dJh}= 3.02\sqrt{dJh} \\ \hat{H}= J\sum_{\langle mn\rangle} X_m X_{n}+h\sum_n Z_n & 4\mathcal{X}_{\frac{d-1}{d}}dJ \leq 8.93 d J \\ & 4\mathcal{X}_0dh = 6.04 d h \\ \hline d\text{-dimensional Fermi-Hubbard} & 2 \mathcal{X}_{\frac{3U}{4dJ}}dJ \\ \hat{H}=J\sum_{\langle mn\rangle, s=\uparrow,\downarrow} (a^\dagger_{m,s}a_{n,s}+\mathrm{H.c.}) & 8\mathcal{X}_{\frac{d-1}{d}}dJ\leq 17.9 d J \\ ~~+U \sum_n a^\dagger_{n\uparrow}a_{n\uparrow}a^\dagger_{n\downarrow}a_{n\downarrow} & \mathcal{Z}_{U/J}J~(d=1) \\ \hline \text{1D Heisenberg XYZ} & 4\mathcal{X}_0 F(J_x,J_y,J_z) \\ \hat{H}=\sum_n (J_x X_n X_{n+1}+J_y Y_n Y_{n+1}+J_z Z_n Z_{n+1}) & 34.6 \max\{J_x,J_y\}\\ \hline \end{array}$$

As for how good these bounds are, I haven't investigated in general, but for the 1D TFIM at critical point $J=h$, exact solution gives $v_{\text{LR}}=2J$, while the above bound gives $2\mathcal{X}_0J\approx 3.02 J$. Similarly, at the $U=0$ point of FH and the $J_x=J_y,J_z=0$ point of Heisenberg XYZ, the above bound are all larger than exact solution by a factor of $\mathcal{X}_0\approx 1.50888$. [Actually at these special points the latter two are equivalent to decoupled chains of TFIM, as can be directly judged from their commutativity graph.]

---

Tighter bound for 1D XY by mapping to free fermions

Now let's talk more about the 1D XY model. As you noticed, it's exactly solvable by mapping to free fermions: $$\hat{H}=\sum_nB_n(a_n^\dagger a_n-1/2)+\sum_n J_n(a_n^\dagger a_{n+1}+\mathrm{H.c.}). $$ For general $B_n(t),J_n(t)$ you need to solve the free-fermion problem numerically, but let me mention two special cases that are analytically tractable.

1. $B_n(t)=B,J_n(t)=J$ are fixed and translation invariant. Then the exact solution is \begin{eqnarray}\label{eq:aisigmat} a_{n}(t)&=&\frac{1}{\sqrt{2\pi}}\int_{-\pi}^\pi \tilde{a}_{k}e^{-i2Jt\cos k }e^{ikx}d k\nonumber\\ &=&\sum_{m}\mathcal{J}_{|n-m|}(2Jt) a_{m}(0), \end{eqnarray} where $\mathcal{J}_{|n-m|}(2Jt)$ is the Bessel function of order $|n-m|$. So the LR speed is $\color{red}{v^{\text{XY}}_{\text{LR}}=2J}$.

2. $B_n,J_n$ are fixed in time but are completely random (quenched disorder). Then due to many-body localization (or Anderson localization in the fermion picture), information don't propagate in this system, so $v_{\text{LR}}=0$. More rigorously, in arXiv:quant-ph/0703209, the following bound is proved for disordered case: $$[A_X(t),B_Y(0)]\leq \mathrm{const.}~ t ~e^{-d_{XY}/\xi},$$ with a decelerating, logarithmic light cone $d_{XY}=\xi \ln t$.