The Lieb-Robinson correlation function for the quantum transverse field Ising model

,


I. INTRODUCTION
How does quantum information or influence spread?A significant advance in understanding this important question was famously provided in 1972 by Lieb and Robinson [1].They focused on manybody systems composed of localized but interacting sub-systems.They considered the norm of the commutator where Âk operates on subsystem k and Bm (t) = e i Ĥt/ℏ Bm (0)e −i Ĥt/ℏ (2) is a Heisenberg operator which acts on subsystem m at time t.They showed that under very general conditions, the norm of the commutator in Eq. ( 1) is bounded, where d k,m is an appropriate distance between the two subsystems.The Lieb-Robinson velocity v LR is understood to be an effective speed limit on the propagation of quantum information or influence.This bound holds even for non-relativistic quantum mechanics and is often envisioned as a sort of effective "light-cone" of possible interaction outside of which influence is exponentially suppressed.
There is by now a large literature calculating and improving the values of v LR in the Lieb-Robinson bound for many different systems and including the effects of finite temperatures and disorder [2][3][4][5][6].The Lieb-Robinson result has been shown to imply finitelength correlations in the ground state of any gapped many-body system [7,8].Considerable work has been done on specific lattice spin systems, with nearneighbor coupling and those with power-law interactions [9,10].Work on propagation in such spin models is particularly relevant to current developments in quantum computing and quantum information processing.Recent reviews of the Lieb-Robinson bound and its connection to fundamental questions of operator locality and quantum information are given in references [11] and [12].The bounds have been explored experimentally [13]; a review is in [14].
The Lieb-Robinson analysis is related to issues of quantum entanglement [15].Equations ( 1) and (2) involve only Âk , Bm , and the system Hamiltonian Ĥ. Entanglement by contrast is a feature of the quantum state of a composite system, which develops precisely because of the interaction between subsystems.The Lieb-Robinson commutator captures an aspect of the effective, possibly mediated, interaction between the two subsystems which is the precondition for them becoming entangled with one another.As emphasized by Chen, Lucas, and Yin [11], the quantity C A,B (t) represents a fundamental constraint on purely quantum non-locality and entanglement.
Here we focus not on the bound in Eq. ( 3) in the general case but on one specific model-the quantum transverse field Ising (QTFIM) model for a onedimensional chain of qubits (or equivalently, spins).Following the work of Colmenarez and Luitz [10], we define the Lieb-Robinson correlation function by where the operators are the usual Pauli operators on qubits 1 and k.
The QTFIM is the simplest model that includes both intrinsic local dynamics and coupling between subsystems, and it has other advantages as well.There is a well-known quantum phase transition that occurs when the strength of the local dynamics is equal to the inter-qubit coupling.A disordered system for small coupling becomes ferromagnetically ordered at a critical value of the coupling.The Hamiltonian can be diagonalized using now-standard methods which map the manyspin system onto free fermion quasiparticles [16,17].Though simple in form, the QTFIM is also experimentally and practically relevant; it is used in D-Wave quantum-annealing based computers.Recent results have been reported for one-dimensional Ising chains consisting of 2000 qubits [18] and higherdimensional configurations with 5000 qubits [19].Quantum-dot cellular automata dynamics have been mapped to the QTFIM [20].
We study the behavior of C k (t) as a function of space (qubit index k) and time.Our goals are to understand in detail the propagation of these quantum correlations down the chain and in particular (a) the effects of varying the Hamiltonian parameters, (b) the fingerprints of the quantum phase transition, and (c) the relation of the propagation speed to the Lieb-Robinson velocity.
Directly evaluating C k (t) has been practically challenging for all but relatively short chains even for the QTFIM.The reasons are familiar: the size of the state space grows exponentially with the length of the chain, and the matrix exponentials in Eq. ( 2) are costly to compute.Heroic methods may be required for chains of 22 qubits [10].For the QTFIM, we devise an operator Pauli walk method which scales only linearly with the system size, and so allows evaluation of Eq. ( 4) for chains of hundreds of qubits and long times.This is sufficient to establish key features of the propagation for chains of any length.For the special case of the QTFIM at the quantum critical point, the method is extended to yield a closed form expression for the semi-infinite chain.
The term correlation function often denotes an expectation value, for example of the form and so depends on the system state [21].By contrast, the Lieb-Robinson correlation function defined in Eq. ( 4) is a correlation between operators independent of state.Moreover, because of the commutator, it captures only explicitly quantum mechanical effects not classical correlations.It is also fairly associated with bounding the spread of information in a quantum channel [22], though it is not defined in Shannon information theory terms.The quantity C k (t) could also be conceived as quantifying influence, or even potential influence of one system on another (and reciprocally).We will refer to it as a correlation function; its precise meaning is clear from the definition.We consider the dynamics of the QTFIM in the sense of calculating C k (t) using the time-dependent Heisenberg operators.The Hamiltonian itself is time-independent.We are not addressing the situation of the time-dependent state of the system when the Hamiltonian changes abruptly-a so-called "quench." The organization of the paper is as follows.Section II describes the QTFIM and the corresponding operator space.Section III presents the operator Pauli walk method which is at the heart of our calculations.The special case of evaluating C k (t) at the quantum critical point is described in Section IV.At the leading edge of the correlation front which propagates down the chain, a simpler analysis is possible.This was described in [23] in one, two, and three dimensions.We review the main results for one-dimensional chains in Section V in order to connect them with the present calculations.In Section VI, we examine the effect of the inter-qubit coupling strength on C k (t), its saturation value, and propagation speed.A discussion of the main results follows.

A. System Hamiltonian
We consider a linear array of N q qubits (or spins) with near-neighbor interactions.The system is described by the transverse-field Ising model with Hamiltonian where the operators σ{x,y,z} k are the usual Pauli operators on site k.The first term in Eq. ( 6) represents the internal dynamics of each qubit and produces a characteristic time for these dynamics given by The second term represents the interaction between adjacent qubits and results in an energetic cost for neighboring opposite spins.We can re-write Eq. ( 6) in dimensionless form as with J ′ = J/γ now characterizing the system.The quantum phase transition between the ordered and disordered phase occurs at J ′ = 1.

B. Operator space
The set of operators on the qubit array form a Hilbert space with the inner product with N = 2 Nq being the dimension of the operator state space.The norm induced by this inner product is known as the normalized Frobenius norm.The operator norm, which is equal to the maximum modulus of the singular values, is frequently used in defining the Lieb-Robinson bound.Different norms can in general produce somewhat different results [3], but for the problem at hand the two norms yield identical values (see Section III F below and the Supplementary Material [24]).Pauli operators on the same site k obey the commutation relations and Pauli operators on different sites commute.It is sometimes convenient to use the common alternative notation We consider Pauli strings of the form For example, using the simpler notation, or The set of all 4 Nq Pauli strings we denote P.These Pauli strings form an orthonormal basis for Hermitian operators on the system such that

III. CALCULATING THE LIEB-ROBINSON CORRELATION FUNCTION
A. Heisenberg time dependence In the Heisenberg picture, the time dependence of any operator can be written where Q = Q(0), and we use the usual notation for n iterated commutators Using Eq. ( 7), the time dependence of σz 1 can then be written The Lieb Robinson correlation function is defined by and we take this as the most direct method of computing it.Making explicit use of the norm defined in Eq. (10) we have From the definition, we see C k (0) = 0 because σz 1 commutes with itself and with the Pauli operators on all other sites.The time dependence of σz 1 (t) can be envisioned as the operator spreading out from the first qubit progressively down the array, incorporating components of Pauli operators on other qubits.As the expanding operator encounters site k, the quantum correlation between the qubit 1 and qubit k, as quantified by C k (t), increases.It is this spread and growth that we want to understand in detail.
We focus first on the operator formed by the commutator [σ z k , σz 1 (t)] and using Eq. ( 19) obtain Because Pauli strings span the operator space, the commutator can be expanded as a weighted sum of Pauli strings as To solve for the coefficients C n,σs , we multiply both sides of (23) by σs ′ Ĥ′ n , σz and take the trace Now using Eq. ( 16), we obtain or By substituting Eq. ( 23) into Eq.( 22), we have By virtue of Eq. ( 16), one can show that the operators σ † k,s also have an orthogonality relationship, From Eq. ( 28) can write (31) It is necessary to calculate the trace of the absolute value square of this commutator so that we can evaluate the Lieb Robinson correlation function using Eq. ( 21).We proceed by defining two quantities that will be temporarily useful, and Note that both the quantities R n (t) and B σs (t) are numbers rather than operators.Using Eqs.(32) and (33), we can rewrite Eq. ( 31) in the compact form The absolute value squared of the commutator is then where we have used Eq. ( 30) and define Re-inserting the expressions for B σs (t) and R n (t) from Eqs. (32) and (33) yields and so from the definition for the correlation function in Eq. ( 21), we have To calculate this correlation function we need to determine the quantity C n,σs from Eq. ( 27) and D k,σs from Eq. (39).

B. Calculating C n,σs using operator Pauli walks
In this section, we develop an algorithm for calculating C n,σs using algebraic methods, graph theory, and finally matrix methods.Specifically, we aim to determine the subset of the 4 Nq Pauli strings that actually contribute to the sum in (41), and to calculate the inner products in Eq. (27).
We now write the Hamiltonian in the more compact notation We need to calculate the projection of Pauli string σs on the iterated commutator We first consider a chain of N q = 4 qubits and then will generalize later.Our notation for Pauli strings will suppress identity operators so, for example, X 1 Y 2 will be understood to mean X 1 Y 2 I 3 I 4 .We focus first on the iterated commutator We begin by calculating [H ′ , Z 1 ] and generate the other commutators that arise as a result of the repeated iteration of commutation with H ′ in Eq. ( 44).
The Pauli operator commutation relations yield In each of these equations, the term underlined on the right gives us the next commutator to evaluate.
For each equation, the commutator between H ′ and the non-underlined term can be evaluated using the previous results.The sequence terminates in his case because there is no fifth qubit.One result of this process is that the sequence allows us to order the relevant Pauli strings as FIG. 1. Operator node graph for a 4 qubit line.
The Pauli strings σ m enumerated in Eqs. ( 54)-(61) form a closed set in the sense that the iterated commutator Ĥ′ n , Z 1 will be a linear combination of these Pauli strings and no others for all values of n.All other Pauli strings σs will yield an inner product C n,σs = 0. Rather than considering all 4 4 = 256 Pauli strings, we will only need to consider C n,σm for m = [1,8].
Next we can use Eqs.( 46)-( 53) to calculate the iterated commutators in Eq. ( 45) and the inner products in (43).We evaluate the first few using operator algebra.For the n = 1 iterated commutator we have therefore from Eq. ( 43) The next level of iteration can be evaluated using Eq.(47) to obtain so evaluating the inner product in Eq. (43) yields Similarly To evaluate the sum in Eq. ( 41), we need to evaluate C n,σm for arbitrarily large values of n.
We now construct a finite directed graph which represents this algebraic process.Figure 1 shows an array of operator nodes (in blue) corresponding to the 8 Pauli strings listed in (54)-( 61).The connectivity of the nodes is determined by Eqs. ( 46)-(53).Two nodes are connected by an edge if the two corresponding operators appear on the left and right side of one of these equations.The direction of the edge (indicated in the figure by an arrow) is from the operator that appears on the left side of the equation to the operator that appears on the right side.The edge weight is given by the coefficient of the operator on the right side of the equation.Edges are color-coded in Figure 1 according to their weight and direction.
For example, consider the node in Figure 1 corresponding to the operator X 1 Z 2 .There are four edges connecting this node to other nodes: • The right-directed edge going from node Y 1 with weight 2iJ ′ given by Eq. ( 47) • The right-directed edge going to node X 1 Y 2 with weight 2i given by Eq. ( 48) • The left-directed edge going to node Y 1 with weight −2iJ ′ given by Eq. ( 48) • The left-directed edge going from node X 1 Y 2 with weight −2i given by Eq. ( 49).
We define an operator Pauli walk on this graph as an ordered list of nodes (Pauli strings) such that each consecutive pair of nodes is connected by a directed edge.Each node may be visited more than once and each edge traversed more than once.The length of such a Pauli walk is the number of edge traversals, or equivalently, one less than the number of (not necessarily unique) nodes visited along the walk.The weight product of a Pauli walk is the product of the edge weights along the walk [25].
The sequence of iterated commutators generated by Eqs. ( 46 Some examples: From Figure 1 we see that the weight product for this walk is in agreement with Eq. (72).
• There are 2 Pauli walks of length 3 connecting node Z 1 to node Y 1 : From the figure we see that the sum of the corresponding weight products is in agreement with Eq. ( 70).
• There are 3 Pauli walks of length 4 connecting node From the figure we see that the sum of the corresponding weight products is In general, C n,σm is a polynomial in J ′ .As the walk length n becomes large, the number of walks becomes very large, even for a modest N q , because each walk can go back and forth between nodes many times.
The process of computing the sum of the weight products for all Pauli walks of a given length can be automated by defining the adjacency matrix A whose (ℓ, p) element is the edge weight from the ℓ th node to the p th node.In the 4-qubit example illustrated in Figure 1 the 8 × 8 adjacency matrix is (81) It is convenient to also define A ′ ≡ A/(2i).The sum of the weight products of all walks of length n that begin at node j and end at node k is the (j, k) th element of the n th power of A. In our case, this means we can evaluate the sum of all weight products for Pauli walks starting at the first node σ1 = Z 1 and ending at node σm by evaluating Using (82), we can now rewrite Eq. ( 41) with the sum over the set of Pauli strings σ s now expressed as a sum over the ordered set {σ m } as where We have seen that for the 4-qubit line, only the 8 Pauli strings enumerated in Eqs. ( 54)-(61) contribute to the correlation function in Eq. ( 41).Examining the structure of this set, we see that the Pauli strings in {σ m } appear in pairs associated with each qubit k.The first element of each pair ends with Z k , and the second ends with with Y k , after which all the succeeding σ m 's have the operator X k .The quantity D k,m is zero for any Pauli string which includes Z k or I k because both commute with Z k .Only if σ m includes Y k or X k is D k,m ̸ = 0, and in both of those cases it has a value of 1.The result is that D k,m is 0 for m < 2k and 1 otherwise.Therefore the effect of D k,m is to limit the sum over m in Eq. ( 83) to start at m = 2k, so we can write E. Generalizing the method of operator Pauli walks Generalizing from the case of N q = 4 to any number of qubits is now straightforward.In place of the 8 equations relations in Eqs. ( 46)-(53), we have 2N q commutation relations.For each qubit with index j = 1 to N q , there are two relevant equations: one is the commutator between the Hamiltonian and a Pauli string that ends in Z j , and the second is a commutator between the Hamiltonian and a Pauli string that ends in Y j .Specifically, we have and These connect the Pauli strings that contribute to Eq. ( 41).These 2N q operators σ m also occur in pairs associated with qubit index j ∈ [1, N q ].
The method of enumerating Pauli walks on the nodes labeled by σ m as described in the previous section is naturally extended to any length chain.The adjacency matrix A, determined by the coefficients in Eqs. ( 87) and (87), can be written in the general case as and we again define A ′ ≡ A/(2i).
The value of the Lieb-Robinson correlation function can be calculated by calculating powers of A ′ and evaluating ..
(90) This can be rewritten in terms of the matrix exponential of the adjacency matrix: Importantly, whereas the size of the Hamiltonian matrix for a system of N q is 2 Nq × 2 Nq , the size of the connectivity matrix A is only 2N q × 2N q .This makes the method tractable for much larger systems.To evaluate C q (t) with this method only requires evaluating Eq. (91).

F. Choosing a norm
We can now revisit the question of which norm to use in computing the Lieb-Robinson correlation function.
The operator norm and Frobenius norm of an operator Q can be expressed as The operator norm is the square root of the maximum eigenvalue of QQ † and the Frobenius norm is the square root of the average of the eigenvalues of QQ † .In general these are of course not equal to each other.
For the QTFIM, we have shown above that the operator Q = [σ z k , σz 1 (t)] can be written as a linear combination of just the 2N q Pauli strings in Eq. (88).One can show that therefore in this case, for all times QQ † is a (time-dependent) multiple of the identity operator.For the identity operator I, the eigenvalues are all 1, so the largest eigenvalue and the average eigenvalue are identical.The two norms then yield the same result.The detailed proof is included as Supplemental Material [24].As a consequence of this identity, the Lieb-Robinson correlation function defined by Eq. ( 4) can be calculated with either definition of the norm; we have chosen to use the Frobenius norm.20).The points show the results calculated from the method of Pauli walks using Eq. ( 91).The Pauli walk method in this case is more than 100 times faster and can be extended to much larger chains.

G. Comparison with the direct time-exponential method
Figure 2 shows the Lieb-Robinson correlation function C k (t) for a chain of N q = 10 qubits.The timescale is in units of τ as defined by Eq. (7).The solid curves show the results of a direct calculation using the operator exponential time dependence in Eq. (20).The points show the result of the operator Pauli walk method using Eq.(91).The agreement between the two methods is essentially exact.The results are shown for qubit coupling J ′ = 1/2.Because large-matrix exponentials are costly to calculate, the Pauli walk method here provides a speed-up of more than a factor of 100 for the calculation.To assure accuracy, we use an extended-precison arithmetic package [26].
Several general features of C k (t) are visible even in this relatively short line.The value of the correlation function is zero at t = 0 for all k because operators on the first qubit and the k th qubit initially commute-they cannot know about each other yet.As the influence of the first qubit propagates down the chain, the correlation function for each site starts to rise.Quantum oscillations on the characteristic time scale of τ are due to the internal dynamics of each qubit.The maximum value of C k (t) is C sat = 2.This has its origin in the Pauli commutation relations of Eq. ( 11)-the maximum amount of "noncommutation" is 2. The correlation of qubit 1 with itself, C 1 (t) is initially zero because [σ z 1 (0), σz 1 ] = 0, and rises because σz 1 (t) starts to mix in components of σx 1 and σy 1 .The nesting of C k (t) for successive values of qubit index k is apparent in Figure 2. Note also that for each pair of successive qubits k and k + 1, there are particular times for which the values of C k (t) and C k+1 (t) are nearly identical.For example C 2 (t) and C 3 (t) are nearly the same for t/τ ≈ 0.7.Calculations on a fine grid of times around such points reveal that the values do not in fact become identical, merely quite close.
The computational cost of using direct operator exponentiation as in Eq. ( 20) becomes prohibitive for chains of even 15 qubits.Moreover quantum reflections from the end of the chain propagate back and so relatively short chains can give a misleading impression of the way in which correlations would propagate in a longer chain.In the case shown in Figure 2, for example, after about t/τ = 2 the behaviour observed is strongly influenced by interference due to the finite chain length and reflections from the end.The method of operator Pauli walks gives us the ability to examine much larger systems and observe behavior far from the ends.We will explore this behavior more in Section VI.

IV. THE SEMI-INFINITE CHAIN AT THE PHASE TRANSITION
The results for C k (t) given by Eq. ( 91) can be simplified even further if we consider the case of a semi-infinite chain at the critical inter-qubit coupling J ′ = 1.The adjacency matrix in Eq. ( 89) becomes It will prove helpful in this context to renumber the operator nodes in Eq. (88) (not the qubits) so the first node has the index 0. We then define and write the correlation function The quantity G n,m is the sum of the weight products for all walks of length n from node 0 to node m.The weight product for each of the walks will be ±1.Each negative step (a step to the left) in the walk generates a negative sign in the weight product.The overall sign of the weight product for a walk is determined by the number of negative steps.All of the walks of length n from node 0 to node m must have the same number of negative steps, and thus the same overall sign.One corollary is that in computing G there is no interference between walks to m that have the same length.Consider a walk from node 0 to node m of length n.Let n p be the number of positive steps and n n be the number of negative steps.The weight product for the walk will be (+1) np (−1) nn .We have n = n p + n n and m = n p − n n , so Since n n is an integer, we must have n − m ≡ 0 (mod 2), i.e., n and m must have the same parity.A walk with an even (odd) number of steps cannot end on an odd (even) node.
Let N w (n, m) be the number of walks from node 0 to node m of length n.A well established result in combinatorics (Bertand's ballot problem) is that in the case of a line of nodes unbounded to the right the number of such walks is (See, for example, the expression in reference [27] with a = n p + 1 and b = n n .)Note that for m > n, N w (n, m) = 0.
Combining the information about the number of walks and the weight product for each walk, we can write So, if we consider a semi-infinite line of qubits (N q → ∞) we can write for the correlation function Because only values of n with the same parity as m contribute to the sum in Eq. (99), we make the change of variables n ′ ≡ (n − m)/2, n = 2n ′ + m , with the result where we have then substituted the symbol n for n ′ .The sum over n can be connected to a Bessel function of the first kind through the relation (8.440 in [28]) Substituting Eq. (101) into Eq.(100) yields where we let We divide the sum S into two parts The first term can be written Now we can use the relations (11.4.7 in [29]) and (8.513.4 in [28]) to write Eq. ( 106) as Substituting Eq. ( 109) and Eq. ( 105) into Eq.( 103), we obtain Substituting in Eq. ( 104), we arrive at a closed-form solution for the Lieb-Robinson correlation function for the J ′ = 1 case,

V. THE LEADING EDGE OF CORRELATIONS
The complicated behavior of the Lieb-Robinson correleation function shown in Figure 2 is much simpler if we focus on the the early turn-on of correlation in each successive bit-the leading edge of the correlation.Figure 3 shows C k (t) for J ′ = 1/2 as correlations first start to grow and propagate along the qubit chain.The solid lines show the results computed directly from Eq. ( 20) and the dots indicate the results from the Pauli walk method of Eq. (91).For each qubit k, as the correlation with qubit 1 just begins to grow, and C k (t) is still small, the growth  20) and black dots are the results of the Pauli walk method given by Eq. ( 91).
follows a power-law behavior in time ∼ t 2k−1 .At this leading edge, we have shown in reference [23] that the Lieb-Robinson correlation function is given by: (112) Only odd orders contribute to the initial growth so the next most important term is proportional to t 2k+1 .
Figure 4 shows the early growth of the Lieb-Robinson correlation function for the ten-qubit chain on a log-log scale.Three methods of calculating the correlation are shown: the direct exponential time dependence of Eq. ( 20) (solid lines), the operator Pauli walk method of Eq. (91) (solid dots), and the leading-edge expression of Eq. ( 112) (open circles).The agreement between the three methods is excellent.
The leading edge correlation response quantifies the way in which the correlation between qubit k and qubit 1 initially turns on.How fast does this leading edge of correlation propagate down the chain?For qubits far enough down the chain that we can replace the factorial in the denominator of Eq. ( 112) with Stirling's approximation, we obtain where the value of v LR is given by which we identify as the Lieb-Robinson velocity.A comparison of the leading edge expression of Eq. ( 114) with the results of the operator Pauli walk method is shown in Figure 5 for the J ′ = 2 case.Snapshots of C k (t) are shown for various times as the correlations propagate down a chain of 200 qubits.The solid lines are from the Pauli walk method, which includes all higher orders in time.The calculation using the leading edge expression of Eq. (114) matches well for early times, but not as well later as the correlations moves down the chain and higher order terms become more important.As we will see in the next section, the correlation front moves down the chain and the leading edge description is accurate only well out ahead of that front.Note that the magnitude of C here is very small.
To see something more like a simple exponential decay in front of the correlation front, one must go much further down the chain.In that limit, where k is quite large and in the region of the leading edge, where k ≈ v LR t, Eq. ( 113) can be further simplified to The details of the connection between Eq. (113) and Eq.(115) are in Appendix A. We recognize that this nearly-exponential leading edge of correlation corresponds to the classic results of Lieb and Robinson.For the TFIM, we have here the additional information of the specific prefactor and its dependence on J ′ , and the modification of the wavefront by the factor 1/ √ k (also noted in [30]).Figure 6 shows a comparison between the results   of Eq. (112) (lines), Eq. (114) (red dots), and the exponential form given by Eq. (115) (black dots).For C k (t) to assume this simple form-an exponential front moving down the chain at velocity v LRrequires looking very far down the line, here more than 10,000 qubits.

VI. PROPAGATION OF CORRELATIONS
We now have three techniques that let us characterize the way the correlation quantified by C k (t) propagates in larger systems.The operator Pauli walk method described in Section III is our primary tool, augmented by the critical point (J ′ = 1) results from Section IV, and by the expression for the leading edge of correlations from Section V.
A. Overall behavior and effect of coupling strength Several features of the way in which the correlation spreads down the chain now become clearer.As one might expect, away from the ends, the turning on of correlations in each qubit becomes more regular, with a similar pattern for successive values of qubit index k down the line.In the time-frame shown, reflections from the end of the chain play no role.The characteristic quantum oscillations remain, but now we can see that in the absence of end effects, C k (t) saturates to the maximum value of 2 eventually for all values of k.This feature was not as obvious in the shorter chain of Figure 2.While Figure 7 plots the time dependence of C k (t) for specific qubits down the chain, Figure 8 shows the same process in a different way.The figure shows C k (t) for all the qubits in the chain at particular snapshots in time: t/τ = [1, 3, 5, . . ., 39].This displays the correlation front moving down the chain.Behind the correlation front, C k saturates at the maximum value of 2 (for this case when J ′ = 1/2).The regular spacing of the curves in both Figures 7 and 8 suggests that the front moves with a constant velocity down the chain, creating the so-called "light cone" of propagating influence.(It will turn out that this is not the Lieb-Robinson velocity.)At a time far in advance of the front's arrival, the state of qubit k is necessarily independent of whatever is happening, or has happened, with qubit 1 so C k (t) is small.The structure of the Hamiltonian has not yet allowed qubits 1 and k the possibility of becoming quantum entangled with each other.
We note that Figure 8 shows an initially perplexing feature-what appears to be small plateaus in C k as a function of k at particular times.As the inset shows, this is simply a manifestation of the abovedescribed near-equalities that occasionally occur for neighboring qubits.The plateau in the small box is simply due to the fact that C 39 (t) and C 40 (t) are very nearly equal at t/τ ≈ 13.
Results for the propagation of correlations for the critical coupling J ′ = 1 are shown in Figures 9 and  10.In this case, we can use the analytic result for  a semi-infinite chain of qubits given by Eq. ( 111).
The behavior is qualitatively similar to the J ′ = 1/2 case discussed above, although the curves are smoother.The values of C k (t) still saturate at 2. There are no near-equalities for neighboring qubits and so no plateaus.One noteworthy difference that is immediately visible is that the correlation front moves down the chain much faster than the case with J ′ = 1/2 (compare Figures 8 and 10).We will return to analysis of the velocity of propagation below.The correlation function in the case of stronger coupling, here J ′ = 2, is shown in Figures 11 and  12, for N q = 200 qubits.This stronger coupling produces a ferromagnetically ordered ground state with a double degeneracy.These results are calculated with the method of Pauli walks.Unlike either the weak-coupling or critical coupling case, the correlation function now saturates at less than the maximal value of 2-in this case C k (t) saturates at the value 1.The velocity of the correlation front is the same as that in the J ′ = 1 case.

B. Saturation value
For a long chain, where reflection from the end can be neglected, we see that C k (t) for any qubit eventually reaches a saturation value which we denote C sat .When the correlation front that is moving down the line is well passed a particular qubit k, the correlation function C k for that qubit will approach arbitrarily closely to C sat .
The calculated values of C sat for different values of J ′ are shown in Figure 13.For J ′ ≤ 1, corresponding to the weakly-coupled disordered phase, C sat = 2 independent of J ′ .We observe that for J ′ ≥ 1, corresponding to the ordered phase, the value of C sat decreases inversely with increasing J ′ : The agreement between the numerical calculations using the Pauli walk method and this conjectured, though unproven, analytical expression is excellent.

C. The velocity of correlation propagation
We have seen in Figures 7-12 that the Lieb-Robinson correlation front propagates down the chain of qubits.How fast does this front travel?We quantify this in the obvious way by picking an arbitrary threshold level C threshold and finding the time t k when C k (t) = C threshold .Using finite differences, the velocity of the front can then be calculated.Near the beginning of the chain, the velocity varies, but subsequently saturates to a stable value, as is clear in the figures.We find that this value v front is This front velocity can be connected to a wellestablished result for the QTFIM.For the infinite qubit chain, or one with periodic boundary conditions, the Hamilonian of Eq. ( 6) can be mapped onto an equivalent fermion problem using the Jordan-Wigner transformation [16,17,31].A subsequent Bogoliubov transformation rotates the fermion creation and annihilation operators to yield a free (noninteracting) fermion Hamiltonian with quasiparticle excitation spectrum Here q is the wavenumber of the harmonic excitation and g ≡ 1/J ′ .The corresponding group velocity can be calculated via the usual expression The maximum value of this group velocity v max g is given by precisely the expression for v front in Eq. ( 117).The derivation is included for completeness in Appendix B. Thus, the correlation front quantified by C k (t) and depicted in Figures 7-12 travels down the qubit chain at the speed of the maximum quasiparticle group velocity.
There is another velocity in the problem, as we have already described in Section V, the Lieb-Robinson velocity in Eqs. ( 113) and (114).The front velocity v front and the Lieb-Robinson velocity v LR are plotted as a function of coupling strength J ′ in Figure 14.The points are obtained from numerical calculations using the operator Pauli walk method.The solid lines are from the analytic expressions in Eqs. ( 114) and (117).The phase transition at J ′ = 1 is clearly evident in v front , but is not manifested in v LR .For any value of J ′ > 0, The relationship of the two velocities in the problem can also be seen in the light-cone plot of C k (t) in Figure 15.The color in the figure represents the logarithm of the correlation function for a 200-qubit chain with J ′ = 2.The calculation is done using the Pauli walk method of Eq. (91).Isocontours for the correlation function are shown in white (slight oscillations are an artifact of the contouring algorithm and the discreteness of the grid).The black line segment corresponds to the velocity v front , and the red line segment corresponds to v LR .Because time is on the vertical axis, a lower slope on the graph corresponds to a higher velocity.The main correlation front expands with a linear light cone given by v front .The velocity of the leading edge is initially faster than that given by Eq. (114), and saturates at the Lieb-Robinson velocity from above [23].As the primary front expands down the chain, the slope of the isocontours at a particular qubit shift to correspond to v front .The region that is sufficiently far enough ahead of the primary front that it is well-described by C k ∼ e −2(k−vt) is more than 100 qubits down the chain and has an extremely small value of the correlation function.
A note about terminology is perhaps in order here.The original Lieb-Robinson paper [1] used the term "group velocity" to describe the velocity that appears in Eq. (3).Their result was very general and did not depend on a specific energy dispersion relationship, or indeed on the precise Hamiltonian.The meaning of this velocity is provided by the context in which it appeared-a bound on correlations that depends on the quantity (d−vt).Here we use the term "group velocity" in the standard sense provided by Eq. (119).In this case, the group velocity v g is understood as the speed of a wavepacket composed of a superposition of single-quasiparticle Hamiltonian eigenstates with some distribution of wavenumbers centered about q.We identify our observed (calculated) front velocity v front with the maximum group velocity so defined.

VII. DISCUSSION
The primary result of this investigation is to highlight the fact that for the 1D QTFIM two different velocities characterize the propagation of the Lieb-Robinson correlation function.A correlation front moves down the qubit chain with velocity v front , ahead of which C k (t) is exponentially smaller.This velocity increases linearly with coupling strength in the disordered Ising phase (J ′ < 1) and is independent of coupling strength in the ordered phase.The velocity is given by the maximum group velocity of the single quasiparticle excitation band.Behind the propagating front, C k (t) saturates to a value of 2 in the disordered phase, and 2/J ′ in the ordered phase.
Well ahead of the correlation front, the exponential leading edge of the correlation front moves faster.This velocity saturates from above to the Lieb-Robinson velocity v LR = e √ Jγ/ℏ (equivalent to Eq. ( 114)).This velocity is unaffected by the phase transition at J ′ = 1.The Lieb-Robinson velocity describes a region well out ahead of the main correlation front in which the value of C k (t) is extremely small.
It is not clear that the operator Pauli walk method developed here for the QTFIM can be extended to other models.The key feature that enabled this method to work was that the number of Pauli strings that were required to support the spreading Heisenberg operator was reduced from 4 Nq to 2N q .For longer-range interactions, that may well not be the case.Further investigation will be required to reveal whether the general features of the propagation of C k (t) as described here hold for such systems.
)-(53) correspond precisely to a Pauli walk among the nodes associated with the Pauli strings σ m in Eqs.(54)-(61).The quantity C n,σm is equal to the sum of the weight products for all Pauli walks of length n which start at the first node (σ 1 = Z 1 ) and ends at the node labeled σm .

FIG. 2 .
FIG.2.The Lieb-Robinson correlation function for a ten-qubit chain with J ′ = 0.5.The solid curve shows the results for C k (t) calculated using the exponential time evolution operator in Eq.(20).The points show the results calculated from the method of Pauli walks using Eq.(91).The Pauli walk method in this case is more than 100 times faster and can be extended to much larger chains.

FIG. 3 .
FIG.3.Leading edge of Lieb-Robinson correlations for a ten-qubit line with J ′ = 1/2.Solid lines calculated directly from Eq. (20) and black dots are the results of the Pauli walk method given by Eq. (91).

FIG. 4 .
FIG.4.Leading edge of Lieb-Robinson correlations for a ten-qubit line with J ′ = 1/2.The solid curves shows the results for C k (t) calculated using the exponential time evolution operator in Eq. (20).The black dots show the results calculated from the method of Pauli walks using Eq.(91), and the open circles are calculated from the leading edge expression of Eq. (112).

FIG. 5 .
FIG.5.Leading edge of the Lieb-Robinson correlation function for J ′ = 2.The solid lines show snapshots at particular times calculated with the full operator Pauli walk method of Eq. (91).The black dots are calculated using the simplified leading edge expression of Eq. (113).

FIG. 6 .
FIG.6.The leading edge of the Lieb-Robinson correlation function far down a qubit chain for J ′ = 2. Snapshots of C k (t) are shown for even values of t/τ between 828 and 862.The solid lines are calculated using Eq.(112).Red dots are the results of Eq. (113), and black dots are calculated using Eq.(115).

Figure 7
Figure 7 shows the Lieb-Robinson correlation function calculated using the method of operator Pauli walks for a line with N q = 200 qubits and J ′ = 1/2.The figure shows C k (t) as a function of time for qubit indices k = [10, 15, 20, . . ., 70] and times up to 20τ .Several features of the way in which the correlation spreads down the chain now become clearer.As one might expect, away from the ends, the turning on of correlations in each qubit becomes more regular, with a similar pattern for successive values of qubit index k down the line.In the time-frame

FIG. 9 .
FIG.9.Time dependence of the Lieb-Robinson correlation function for a chain with coupling J ′ = 1.This is the critical point for the quantum phase transition between the disordered (J ′ < 1) and ferromagnetic (J ′ > 1) ground states.

FIG. 10 .
FIG.10.Snapshots of the Lieb-Robinson correlation function at different times for a chain of qubits with critical coupling J ′ = 1.The correlation front is moving down the chain much faster than in the J ′ = 1/2 case shown in Figure8.

FIG. 11 .
FIG. 11.Time dependence of the Lieb-Robinson correlation function for a chain of Nq = 200 qubits with coupling J ′ = 2.

FIG. 12 .
FIG. 12. Snapshots of the Lieb-Robinson correlation function at different times for a chain of Nq = 200 qubits with coupling J ′ = 2.

FIG. 14 .
FIG. 14. Correlation front velocity and the Lieb-Robinson velocity for the leading edge.Points are from numerical solutions using the operator Pauli walk method.Solid lines are from Eqs. (114) and (117).
so the exponential leading edge of correlation is traveling down the chain considerably ahead of, and faster than, the main correlation front seen in Figures7-12.