The Value of the Early-Time Lieb-Robinson Correlation Function for Qubit Arrays

: The Lieb-Robinson correlation function is one way to capture the propagation of quantum entanglement and correlations in many-body systems. We consider arrays of qubits described by the tranverse-ﬁeld Ising model and examine correlations as the expanding front of entanglement ﬁrst reaches a particular qubit. Rather than a new bound for the correlation function, we calculate its value, both numerically and analytically. A general analytical result is obtained that enables us to analyze very large arrays of qubits. The velocity of the entanglement front saturates to a constant value, for which an analytic expression is derived. At the leading edge of entanglement, the correlation function is well-described by an exponential reduced by the square root of the distance. This analysis is extended to arbitrary arrays with general coupling and topologies. For regular two and three dimensional qubit arrays with near-neighbor coupling we ﬁnd the saturation values for the direction-dependent Lieb-Robinson velocity. The symmetry of the underlying 2D or 3D lattice is evident in the shape of surfaces of constant entanglement, even as the correlations front expands over hundreds of qubits.


I. INTRODUCTION
Entanglement in many-body quantum systems and the spread of quantum correlations has been of considerable interest for both fundamental reasons and for the possible applications in quantum computing.Focus on quantum information naturally employs entropic measures derived from the von Neumann entropy of a state.Multiple measures have been proposed to capture the quantum entanglement of spatially separated systems and this remains an area of great activity [1][2][3].Information so conceived is a property of the state of the system, characterized by either a state vector for pure states or a density operator for pure or mixed states.
Another avenue, not tied to the state of the system, considers the commutation relations between local operators that act on spatially separated subsystems.Preeminent here is investigation of the Lieb-Robinson operator [4], which quantifies the quantum correlation between a Heisenberg operator Âj (t) which acts on subsystem j of the composite many-body system, and an operator Bk which acts on a distinct subsystem k.The operator is given by ĈAj,B k (t) = [ Âj (t), Bk ] ( The Lieb-Robinson correlation function is given by the norm of the operator. At t = 0, the operators Â and B, having different support, commute and the correlation function is zero.We can think of the operator Âj (t) as spreading out into the system and as its support comes to include that of B the correlation increases.This captures the quantum correlation, between the two systems which is the source of entanglement.The Lieb-Robinson correlation function depends only on the choice of the operators and the Hamiltonian of the system, independent of any initial state.
Lieb and Robinson established, for quite general conditions, that the correlation so defined propagates outward with a finite speed and with an exponentially bounded leading edge in space: where d(j, k) is an appropriate distance measure between the two subsystems and v LR is the Lieb-Robinson velocity.
Because of the importance of this result, there has been a substantial effort to refine the bound in particular circumstances [5][6][7][8][9][10] and to understand the connection between entropic measures and correlation measures [11,12].Spin models and equivalent qubit arrays have played a prominent role in studying bounds for specific systems [13] and can also be explored experimentally [14].
Following the work of Luitz et al. [15][16][17][18], we focus here on a specific class of Hamiltonians and compute the correlation function itself, rather than a bound on the correlation function.This comes at the cost of generality, of course, but yields insight into the way correlations propagate.In Section II we consider a linear array of near-neighbor coupled qubits-the 1D tranverse field Ising model [13], and calculate the Lieb-Robinson correlation function numerically.This is only practical for small systems.We then derive an analytical expression for the correlation function which we can compare with the numerical results.The analysis is extended to an arbitrary array of qubits in Section III, and then applied to regular square lattices in two and three dimensions in Section IV.
For the purposes of stating a bound, the norm in Eq. (3) need not be specified, but to calculate the value of the correlation function a particular norm must be chosen.Here we use the normalized Frobenius norm where N is the dimension of the Hilbert space.The normalization ensures that the value is independent of the size of the space.
We consider first a linear array of N q qubits with an energy ∆ coupling adjacent z-components, as described by the following Hamiltonian: Here σ(k) is the Pauli spin operator with α = (x, y, z) operating on the k th qubit.Each operator is understood to be embedded in the full direct-product Hilbert space of the array with dimension N = 2 Nq .This system is the 1D transverse field Ising model.The first term represent the internal dynamics which allow each qubit to flip its state.The energy γ sets the characteristic time for the dynamics We have the usual commutation relations for Pauli operators.
We consider the Lieb-Robinson operator between the z components of the first qubit and qubit k.
The Lieb-Robinson correlation function is given by the norm of the operator.
There is no loss of generality here by choosing the first qubit as the reference because the coupling is the same between any adjacent qubits.The correlation function C k captures the quantum correlation between two qubits that are (k − 1) apart.
The time-development of the Heisenberg operator σ(z) The correlation function can then be written Figure 1 shows the Lieb Robinson correlation function for a line of nine qubits by a direct numerical evaluation of Eq. ( 11) for the Hamiltonian in Eq. (5).
We are interested in the early-time behavior as correlations initially spread down the line.By "early-time" behavior we mean times such that the correlation with the reference qubit 1 is small.Our focus is on the leading order term in time as the correlations between the first site and the k th site initially grow.Figure 2 shows a view of the correlations zoomed in to the lower left-hand corner of Fig. 1.
The early time behavior of the Lieb-Robinson correlation function has a power-law dependence, as is evident in Fig. 2 and in the log-log plot shown in Fig. 3.The dots in the figure represent the numerical solution of equation (11) for ∆/γ = 1.The lines are the results of the analytical expression derived below.Figure 4 similarly shows the numerical result for the case when ∆/γ = 5.The larger interaction strengths makes the correlations  11) for the Hamiltonian in Eq. ( 5) with ∆/γ = 1.The lines represent the analytic result of Eq. (30).spread more rapidly down the line.

B. Analytic solution for early-time behavior
We adopt the following notation for an n-nested commutator.
It will also be useful to define a notation for a nested commutator with a sequence of operators.
(13) Using the notation of Eq. ( 12), we can write the time dependence of any Heisenberg operator as Applying this to σ(1) z (t) we have The early time Lieb-Robinson operator for qubit k is the lowest-order Ĉ(n) k which is not zero.All other terms in Eq. ( 17) will be vanishingly small for early enough times.
For convenience, we define so that, To calculate the self-correlation (bit 1 with itself), we first evaluate ( Ĥ) (1) , σ(1) z = (−γ)[σ (1)  x , σ(1) z ] and The early-time Lieb-Robinson self-correlation is therefore where in this case the n = 1 term is the first nonzero term which dominates at early times.The linear increase in C 1 (t) is evident in Figs.
For the near-neighbor correlation C 2 (t), the first nonzero term in Eq. ( 18) occurs for n = 3.

Ĝ(3)
2 = ( Ĥ) (3) , σ(1) The term ( Ĥ) (3) , σ z contains many cross terms from the three iterated commutators with Ĥ in Eq. ( 5).But only one term survives the commutator with σ(2) z in Eq. ( 24), with the result that Ĝ x (25) and the early-time result for the correlation function is For C 3 , the n = 5 term is the leading term.The structure of the only surviving term for Ĝ(5) 3 (Eq.( 18)) in this case is illustrative.We write it using the notation of Eq. (13).Ĝ( 5) x , σ(1) ) traversing the shortest path (directly) between qubit 3 and qubit 1.The alternation between terms requires 5 nestings of commutators with the the Hamiltonian.A nesting of this form is necessary to avoid the whole expression becoming zero.The general feature is that to connect each pair of qubits along the path requires two nested commutators: a σ(k) x term followed by a σ(k) z term.Therefore, only odd-order terms contribute.
Each commutator in Eq. ( 27) with σx generates a factor of (−2i)γ and each commutator with the product σz σz generates a factor of (−2i)(−∆/2) so Ĝ(5) x σ(1) x . (28) The norm of the Pauli string is 1.Using (19), we then obtain In general, the lowest-order nonzero Ĝ(n) k is proportional to a Pauli string of k σx terms, one for each qubit between 1 and k, and n = 2k−1.The result for the earlytime Lieb-Robinson correlation function between the reference qubit 1 and qubit k is given by: Correlations calculated with Eq. (30) yield the solid lines in Figs. 3 and 4 and the dots show the values from numerical solution of Eq. (11).The numerical calculations include all orders in time, and thus are correct at all times, not just early times, as illustrated in Fig. 1.Equation (30) shows the leading order term, which is arbitrarily close to the exact result for early enough times.The excellent agreement between the two seen in Figs. 3 and 4 shows that the higher order terms are indeed negligible.

C. Asymptotic behavior
Figures 1 to 4 show the Lieb Robinson correlation C k (t) for individual qubits with index k in the linear array as a function of time.Figure 5 shows snapshots of the correlation function, computed from Eq. ( 30), as a function of k at particular times.This lets us see the "correlation wave" propagating from the reference bit down the line.Values of C k larger than 10 −2 are clipped and not shown here because our focus is on early-time behavior for each qubit.At the times shown in the figure, the correlation front decays rapidly in space-faster than an exponential.
We examine the asymptotic shape of the correlation front propagating down the line quantitatively by taking the large k limit of Eq. (30).The factorial can then be replaced by Stirlings approximation to obtain where v LR is the Lieb-Robinson velocity given by This velocity is expressed in terms of the natural dimensionless time t/τ .In dimensional form it would be v LR = (e/ √ 2) √ ∆γ/ qubits/s.We can examine how the propagation of the correlation front converges to a constant speed.Let C thresh be a threshold value of the correlation function and let t k /τ be the time at which qubit k crosses that threshold of correlation, so C k (t k /τ ) = C thresh .We define the backwards finite-difference velocity of the propagating front by  .By this point the correlation front, calculated from Eq. ( 30), is well-approximated by the exponential dependence of Eq. ( 34) and propagates at the saturated Lieb-Robinson velocity given by Eq. ( 32).The dotted line shows the correlation front at one time evaluated in the limit of very large k using Eq.(34).
Figure 6 shows this velocity v k as a function of qubit index down the line.The velocity saturates to the value v LR given by Eq. (32).The line shows values calculated from Eq. (30) for ∆/γ = 1 and C thresh = 10 −25 .The saturated Lieb-Robinson velocity is then eπ/ √ 2 ≈ 6.04 in terms of the dimensionless time t/τ .Figure 6 also shows the values of v k calculated numerically directly from the defining Eq. (11).
Figure 7 shows the correlation front much farther down the line (qubits with index more than 10, 250) and at considerably longer times (t/τ more than 1360).By this point, the shape of the leading edge of the front is approximately exponential and moving at a steady speed.Taken together, Figs. 5 and 7 also illustrate that "earlytime" does not refer to small values of t/τ , but rather to Center for Nanoscience and Technology the span of time before the quantum correlations between qubit k and the reference qubit 1 become large.
In the region of the advancing correlation front for very large values of k, Eq. ( 31) can be written: This matches the slope of the lines and the spacing in Fig. 7.
We note that the magnitude of early-time correlations become smaller as the front propagates down the chain because of the k −1/2 dependence in Eq.( 34).This also means the shape of the front differs slightly from a simple e −2kx in each snapshot shown in Fig. 7.The attenuation of the Lieb-Robinson correlation function was also noted in Ref. [18] for the Heisenberg XXX model with shortrange interactions.

III. CORRELATION FUNCTION FOR AN ARBITRARY QUBIT ARRAY
The result in Eq. (30) can be generalize for any network with arbitrary interaction strengths ∆ j,k between qubits j and k.The Hamiltonian for the network is Figure 8 shows a network of nine qubits with the strength of the coupling ∆ j,k between each indicated.The dots in Figure 9 show the result of numerical calculation of the Lieb-Robinson correlation function between the first qubit and the k th qubit in the array using The solid lines are the result of the analytic expression derived below.
For the linear array studied in the previous section, there was one shortest path between qubit 1 and qubit k, the direct path along the line, that led to the dominant term in Ĝ (Eq.18).This was seen, for example, in the calculation for G 3 in Eq. ( 27), for which the path was simply {1, 2, 3}.We must now generalize this analysis to include more than one path through the network.A path is then a list of indices that begin at qubit j and end at qubit k.The length of a path is the number of edges that connect successive qubits.
For an arbitrary array of interacting qubits, there can be several paths with the same minimum length L that contribute to the leading order early-time correlation function.We denote each minimum-length path by P α (j, k) where α is the index of the particular path.The index of the m th qubit along path P α (j, k) is denoted q α m .The lowest order nonzero operator Ĝ(2L+1) then becomes a weighted sum of Pauli strings of σx operators for each qubit traversed in the path.The weights are given by the products of the couplings ∆ q α m ,q α m+1 along the path.We then obtain for the early-time Lieb-Robinson correlation function between qubits j and k where the sum is over all paths with the same minimum length L. Equation (38) makes no assumption about the locality of the coupling ∆ between qubits, which could extend far beyond near-neighbors.Equation (38) reduces to Eq. ( 30) when all the interaction strengths ∆ are the same and connect only near neighbors.In that case, there is just the single minimum path between qubit 1 and qubit k which has length L = k − 1.The results shown in Fig. 9 are in excellent agreement with the numerical calculations.

IV. CORRELATION FUNCTION FOR 2D AND 3D REGULAR LATTICES
We consider a two-dimensional square lattice of qubits with uniform near-neighbor coupling.
Figure 10 shown snapshots of isocontours of C evaluated using Eq. ( 40) for several values of t/τ .The Lieb-Robinson correlation front expands outward from the origin.All the minimum-length paths are Manhattan paths, with the most direct route for expanding correlations along the coordinate axes.The correlation function along coordinate axes reduces to the one-dimensional result of Eq. ( 30).(Note that the reference qubit is here indexed (0, 0) rather than 1.) As with the one-dimensional case, the speed of the correlation front motion saturates to the Lieb-Robinson velocity, but for for the square 2D lattice the velocity depends on direction.We can express this in terms of the angle θ with the x-axis and find: Calculating the value of the Lieb-Robinson correlation function itself, rather than a bound, for a specific Hamiltonian illuminates some aspects of the way quantum information propagates.We are able to see how the Lieb-Robinson velocity emerges as the saturation velocity of the correlation fronts propagation.A specific value for this velocity, and its angular dependence, is also obtained for regular lattices.We also see that the leading of the correlation front is super-exponential as the propagation starts, and acquires the exponential dependence in the bound of Eq. (3) only later.
Equation ( 38) is the main result of this work and applies to a large class of ZZ-coupled Hamiltionians given by Eq. (35).We have specialized this result to address the case of regular square qubit arrays in 1D, 2D and 3D, with near-neighbor couplings.
There is some competition between the strength of coupling between qubits and the number of paths through intermediate qubits.But for early times the cost of adding path length is very high because of the 2L+1 exponent in the time exponent of Eq. ( 38), and because of the factorial in the denominator.This latter factor also results in the decay of correlations in space as the front propagates outward.This is reflected the k −1/2 in Eq. (31) for the Ising-coupled qubit line.

Figure 1 .Figure 2 .
Figure1.The Lieb-Robinson correlation function between the reference qubit 1 and qubit k for a linear array of 9 qubits with near-neighbor interactions.The calculation is done by numerical evaluation of Cq using Eq.(11) and the Hamiltonian of Eq. (5) with ∆/γ = 1.

Figure 3 .
Figure 3.The Lieb-Robinson early-time correlation function for the nine-qubit line demonstrates the power-law dependence on time.The qubit index is k and the dots represent numerical evaluation of Eq. (11) for the Hamiltonian in Eq. (5) with ∆/γ = 1.The lines represent the analytic result of Eq. (30).

Figure 4 .
Figure 4.The Lieb-Robinson early-time correlation function for the nine-qubit line with qubit index k.The dots represent the direct numerical evaluation of Eq. (11) with ∆/γ = 5 and the lines represent the analytic result of Eq. (30).The stronger coupling between neighboring qubits yields a faster propagation of correlations down the chain.

Figure 5 .
Figure 5. Snapshots of the early-time Lieb-Robinson correlation at different times for a line of qubits as a function of qubit index k.The correlation front propagates down the line and, for the times shown, the shape of the front is attenuated faster than exponentially.The solid lines are the results of Eq.(30) with coupling ∆/γ = 1.The dotted line shows the correlation front at one time evaluated using Eq.(31).

Figure 6 .
Figure 6.Velocity saturation to the Lieb-Robinson velocity for the one-dimensional qubit line.The finite-difference velocity defined by Eq. (33) is calculated from Eq. (30) for for ∆/γ = 1 and C thresh = 10 −25 .The values marked by an × are the results of numerical solution of Eq.(11)

Figure 7 .
Figure 7. Snapshots of the Lieb-Robinson correlation as in Fig.5, but now for much later times and much farther down the qubit line.Snapshots are shown for even values of t/τ between 1360 and 1400.Note that the qubit index k is offset by 10 4 .By this point the correlation front, calculated from Eq. (30), is well-approximated by the exponential dependence of Eq. (34) and propagates at the saturated Lieb-Robinson velocity given by Eq. (32).The dotted line shows the correlation front at one time evaluated in the limit of very large k using Eq.(34).

Figure 8 .
Figure 8.An arbitrary network of qubits with different coupling parameters ∆.

Figure 9 .
Figure 9.The Lieb-Robinson correlation function for the qubit network shown in Fig. 8. C 1,k is the correlation between qubit 1 and qubit k.The dots are the result of numerical evaluation of Eq. (36) and the lines represent the general analytic result of Eq. (38).Results for some values of k nearly overlap: k =2, 3; k = 4, 5, 6; k =7, 8.

Figure 10 .
Figure 10.Isocontours of the Lieb-Robinson correlation function for a square two-dimensional array of qubits with nearneighbor coupling.Contours are shown for correlation values of 10 −60 at several values of time.

Ĥ
qubit indices k x and k y go from from −N to +N .The reference qubit is at the origin (k x , k y ) = (0, 0).Clearly there is no distinction between positive and negative directions, so the value of the Lieb-Robinson correlation function at a particular qubit depends only on (n, m) = (|k x |, |k y |).The minimum path length from the origin to (n, m) is simply m+n, and the number of paths with that length is (n+m)!/(n!m!).Applying the general result of Eq. (38), we obtain for the correlation between the reference qubit at the origin and the (n, m) qubit (in any quadrant):

Figure 11 PFigure 11 .
Figure 11 shows snapshots of the isosurfaces of the Lieb-Robinson correlation function C for three successive times.The influence of the Manhattan paths is clear in the rounded octahedral shape.The direction-dependent Lieb-Robinson velocity for the 3D square lattice is expressed in terms of the polar angle φ and azimuthal angle θ: