Asymmetric butterfly velocities in 2-local Hamiltonians

The speed of information propagation is finite in quantum systems with local interactions. In many such systems, local operators spread ballistically in time and can be characterized by a ``butterfly velocity", which can be measured via out-of-time-ordered correlation functions. In general, the butterfly velocity can depend asymmetrically on the direction of information propagation. In this work, we construct a family of simple 2-local Hamiltonians for understanding the asymmetric hydrodynamics of operator spreading. Our models live on a one dimensional lattice and exhibit asymmetric butterfly velocities between the left and right spatial directions. This asymmetry is transparently understood in a free (non-interacting) limit of our model Hamiltonians, where the butterfly speed can be understood in terms of quasiparticle velocities.


I. INTRODUCTION
Understanding the quantum dynamics of thermalization in isolated many-body systems is a topic of central interest. While memory of a system's initial conditions is always preserved under unitary dynamics, this information can get "scrambled" and become inaccessible to local measurements, thereby enabling local subsystems to reach thermal equilibrium [1][2][3][4]. This scrambling can be quantified by studying the spatial spreading of initially local operators under Heisenberg time evolution. Under dynamics governed by a local time-independent Hamiltonian H, an initially local operator near the origin, A 0 , evolves into A 0 (t) = e iHt A 0 e −iHt . As A 0 (t) spreads in space, it starts to overlap with local operators B x at spatially separated locations x. The effect of scrambling is thus manifested in the non-commutation between A 0 (t) and B x , which can be quantified via an out-of-time-ordered correlator (OTOC): where A 0 , B x are local unitary operators, ℜ represents the real part, and the expectation value is with respect to the infinite temperature thermal ensemble. The OTOC is expected to exhibit the following features in systems with scrambling dynamics [6,7,11,17,19,24,[53][54][55][56][57][58][59][60][61][62][63][64][65][66][67]: At early times, A 0 (t) approximately commutes with B x and the OTOC is nearly equal to one. At late times, A 0 (t) becomes highly non-local and spreads across the entire system, and the OTOC decays to zero [11,14,21,25]. At intermediate times, the operator has most of its support within a region around the origin defined by left and right operator "fronts" that propagate outwards, and generically also broaden in time [58,59]. As the operator front approaches and passes x, the OTOC C(x, t) decays from nearly one to zero. We will restrict ourselves to translationally invari-ant systems where operators spreads ballistically with a butterfly speed v B , which is similar in spirit to the Lieb-Robinson speed [68] characterizing the speed of information propagation. In these cases, the operator fronts define a "light-cone" within which the OTOC is nearly zero.
A set of recent papers illustrated that the butterfly velocity can depend on the direction of information spreading [69,70]. In one dimension, the asymmetry between the different directions can be quantified by the butterfly speeds v r B and v l B , where the superscript r (l) represent propagation directions to the right (left). While Ref. [69] showed how this asymmetry could be induced by anyonic particle statistics, Ref. [70] constructed models of asymmetric unitary circuits, and Hamiltonians inspired from such circuits.
In this work, we present a complementary and physically transparent way for constructing a family of Hamiltonians with asymmetric information propagation. Our construction does not rely on particle statistics, nor is it inspired by unitary circuits. Instead, we start with noninteracting integrable spin 1/2 models where the butterfly speed is related to quasiparticle propagation velocities and can be analytically calculated [66,67,71,72]. We show how the butterfly speed can be made asymmetric in such models, before generalizing to non-integrable Hamiltonians by adding interactions.

II. INTEGRABLE HAMILTONIANS
In this section, we construct time-independent integrable Hamiltonians for spin 1/2 degrees of freedom living on an infinite one dimensional lattice. The Hamiltonians only have local terms acting on 2 spins at a time. These models are exactly solvable, so the butterfly velocities can be analytically calculated, and demonstrated to be asymmetric. This family of Hamiltonians parameter-ized by λ takes the form: where σ x j , σ y j , σ z j are the Pauli spin 1/2 operators located at site j, J > 0, h zz , h x , h yz , h zy are constants, and the parameter λ lies in the range [0, 1]. This model can be mapped to a system of free fermions via a Jordan Wigner representation. When λ = 1, the Hamiltonian is the well known transverse Ising model with inversion symmetry about center of the chain. On the other hand, for λ < 1, the Hamiltonian does not have inversion symmetry when h yz = h zy .
In order to detect the ballistic light cone and asymmetric butterfly velocities, we consider the OTOCs where µ, ν ∈ {x, y, z} and β = 0 represents the infinite temperature thermal state. We note that the mapping to free fermions allows Pauli operators to be written in The asymmetry in butterfly speeds for 0 < λ < 1 can be directly understood using the quasiparticle description of the free model. It is known that the butterfly speed in an integrable model is the maximum quasiparticle group velocity [66,67,72], and the operator fronts generically broaden either diffusively or sub-diffusively depending on whether the integrable system is interacting or not [72].

III. NON-INTEGRABLE HAMILTONIANS
In this section, we construct non-integrable Hamiltonian by adding longitudinal fields to the free Hamiltonian H λ [64,65]. The asymmetric butterfly velocities are estimated from a variety of measures including out-of-timeordered correlations, right/left weight of time-evolved operators, and operator entanglement.
The interacting Hamiltonian on a one dimensional lattice with open boundary conditions is where L is the system size, and h z is a longitudinal field strength. We select the particular parameters λ 0.5, although none of our results are fine tuned to this choice. The longitudinal field breaks integrability and is expected to thermalize the system. For non-integrable Hamiltonians with thermalizing dynamics, the level statistics is consistent with the distribution of level spacings in random matrix ensembles [73]. Let E 0 < · · · < E n < E n+1 < · · · be the sequence of ordered energy eigenvalues and s n = (E n+1 − E n ) be the level spacings. One defines the ratio of consecutive level spacings r n = s n /s n−1 , and the distribution of r n can be described by the Wigner-like surmises for non-integrable systems [74,75] where W = 1, Z 1 = 8/27 for Gaussian Orthogonal Ensemble (GOE), and W = 2, Z 2 = 4π/(81 √ 3) for Gaussian Unitary Ensemble (GUE), while they are Poisso-nian for integrable systems. As shown in FIG. (3), the ratio distribution provides evidence supporting the nonintegrability of the Hamiltonian. When λ = 0.5, the Hamiltonian is complex Hermitian, and its ratio distribution agrees with the Gaussian Unitary Ensemble (GUE). When λ = 1, the Hamiltonian is real, symmetric and has the inversion symmetry with respect to its center, and its level statistics in the sector with even parity agrees with the Gaussian Orthogonal Ensemble (GOE) [64,65]. We now characterize the asymmetric spreading of quantum information in this model using a variety of different diagnostics that were discussed in [70].

A. Asymmetric butterfly velocities from OTOCs
In this subsection, we estimate the asymmetric butterfly velocities from OTOCs.
As discussed earlier, as the time-evolved operator spreads ballistically, OTOCs can detect the light cone and butterfly velocities. The saturated value of OTOCs equals approximately one outside the ballistic light cone and zero inside it. Near the boundary of the light cone, the OTOCs decay in a universal form C(j, t) = 1 − f e −c(j−vB t) α /t α−1 [66,67], where c, f are constants, v B describes the speed of operator spreading, and α controls the broadening of the operator fronts. In a generic "strongly quantum" system (i.e. away from large N / semiclassical/weak coupling limits) the operator front shows broadening which corresponds to α > 1 so that the OTOC is not a simple exponential in t [67].
Nevertheless, the decay can still look exponential along rays j = vt in spacetime, C(j = vt, t) = 1 − f e −λ(v)t , defining velocity-dependent Lyapunov exponents (VDLEs) which look like [67]. The VDLEs provide more information about the operator spreading than the butterfly velocities alone.
First, as shown in FIG. (4), we observe asymmetric butterfly velocities in relatively large systems with L = 41 spins. In our numerical calculations, we use the time-evolving block decimation (TEBD) algorithm after mapping matrix product operators to matrix product states [76][77][78], which is able to efficiently simulate the evolution of operators in the Heisenberg picture. The OTOCs shown in FIG.(4) clearly demonstrate asymmetric butterfly velocities. Second, we estimate the asymmetric butterfly velocities from the extracted VDLEs λ Because of the limited computational resources for exact diagonalization, the right and left butterfly velocities are measured by setting the initial local operator at the boundary j = 1 and j = L respectively. In FIG. (5), the OTOCs exponentially decay along the rays with different speed . For a given velocity v, λ(v)/v is the slope of logarithm of the left and right propagating OTOCs versus the distance x. After extracting the VDLEs λ(v) from the OTOCs, here we give a rough estimation of the butterfly velocities via fitting the curve λ(v) ∼ −c(v−v B ) α . In FIG. (6), we obtain the results v r B ∼ 0.29J and v l B ∼ 0.66J.

B. Asymmetric butterfly velocities from right/left weights
Now we turn to the analysis of asymmetric butterfly velocities directly measured from right and left weights of the spacial spreading operators.
To define the right/left weight, note that every operator in a spin 1/2 system with length L can be written in the complete orthogonal basis of 4 L Pauli strings S = ⊗ L j=1 S j , i.e. O(t) = S a S (t)S, where S j = I, σ x , σ y or σ z . Unitary evolution preserves the norm of operators, so S |a S (t)| 2 = 1 holds for a normalized operator. The information of operator spreading is contained in the coefficients a S (t). In order to describe the spatial spreading, the right weight is defined by where the left weight is defined analogously. Because of the conservation of operator norm j ρ r(l) (j, t) = 1, the weight can be interpreted as an emergent local conserved density for the right/left fronts of the spreading operator. Recent studies [58][59][60][61] showed that the hydrodynamics for the right/left weight can be characterized by a biased diffusion equation in non-integrable systems, which means that the front is ballistically propagating with diffusively broadening width. Thus, when the time-evolved operator spreads, ρ r moves to the right with velocity v r B , and ρ l moves to the left with velocity v l B . Here in the numerical calculations of exact diagnolization, the right and left weights are obtained by setting the initial local operator at the boundary j = 1 and j = L respectively. The right weight ρ r (1 + x r , t) of σ x 1 (t) is calculated in order to compare the left weight ρ l (L − x l , t) of σ x L (t), where x r (x l ) is the distance between the right (left) end and the location of initial operator. As shown in FIG. (7), the estimated velocities are v r B ∼ 0.34J and v l B ∼ 0.77J by fitting the times when the weights reach half of the maximum peak for given distances. The symbols ×/+ mark the times when the right/left weights reach half of the maximum peak for given distances. Lower panel: time of the half-peak versus the distance. The solid and dashed lines are the results of linear fitting.

C. Asymmetric butterfly velocities from the operator entanglement growth
In this subsection, the growth of operator entanglement is investigated to estimate the asymmetric butterfly velocities.
The entanglement of time-evolved operators encodes information about its spatial spreading. Refs. [70,79] discussed a coarse-grained hydrodynamic description for the entanglement dynamics of a spreading operator under Heisenberg time evolution. LetŜ(t, x) be the operator entanglement across bond x at time t. To leading order, its growth depends on the local entanglement gradient where s eq is the equilibrium entropy density at infinite temperature, and Γ(s) is a non-negative growth function. Let us discuss the properties of the growth function Γ(s). First, the growth function equals zero in equilibrium. In thermal equilibrium at infinite temperature, the operator entanglementŜ(t, x) has a pyramid shapê S(t, x) = 2s eq min{x, L − x} in a one-dimensional system of length L. Thus, Γ(s eq ) = Γ(−s eq ) = 0 is satisfied. Second, the butterfly velocities are given by the derivatives v r B = −s eq Γ ′ (s eq ) and v l B = s eq Γ ′ (−s eq ). Thus, this function Γ(s) provides yet another way of obtaining butterfly velocities, and it follows that any Γ(s) with asymmetries at its endpoints will have asymmetric butterfly speeds. FIG. (8) shows the growth of operator entanglement across different spatial cuts in time, showing the system approaching the late time equilibrium pyramidal shape.
The statements above are obtained from a 'minimal curve picture of Ref. [79], i.e. the entanglementŜ(t, x) is the minimum of the sum of the initial entanglement S(0, y) and the integral of a velocity-dependent line tension E(v). The minimal curve is a straight line with constant velocity v = (x − y)/t, so thatŜ(t, x) = min y ts eq E( x−y t ) +Ŝ(0, y) . Then one can get the growth function is Γ(s) = min v E(v) − v s/s eq , where s is between −s eq and s eq . Applying the inverse Legendre transformation, the line tension can be expressed as E(v) = max v Γ(s) + v s/s eq . The line tension captures the leading-order hydrodynamics of operator spreading, and the left/right butterfly speeds can be obtained from the intersections of the line tension and the velocity Thus, the butterfly velocities can be estimated if one knows the growth function or the line tension. According to Reference [79], the line tension can be estimated via an effective one E ef f (v) = 1/(s eq t) * S U (vt/2, −vt/2, t) obtained from the entanglement S U (x, y, t) of the time evolution operator U (t) = exp(−iHt), where U (t) acting on L spins is treated as a state on 2L spins. Looking at the intersections in FIG. (9) indicates v r B ∼ 0.5J, v l B ∼ 0.9J. This method gives a very rough estimation. We expect that when the system size is large enough, there is a single speed characterizing information propagation in every direction which agrees across different methods of estimation.

IV. CONCLUSION AND DISCUSSION
In summary, we have constructed 2-local integrable and non-integrable Hamiltonians with asymmetric ballistic information spreading. Exact solutions of the butterfly velocities are obtained in the integrable models. The asymmetric butterfly velocities are estimated from different quantities characterizing operator spreading including out-of-time-ordered correlations, right/left weight of time-evolved operators, and operator entanglement.
Given the constructions and studies in this paper, several open questions would be interesting to explore in the future work. Here we have focused on the information spreading at infinite temperature. How does the asymmetric spreading change at finite temperature? Additionally, it is worth studying how asymmetries encoded in various quantities are intertwined with each other. For example, does the transport of conserved quantities (like energy) inherit the same signatures of asymmetry as the spreading of local operators? Is it possible to disentangle them? Finally, probing the asymmetry of information propagation may also interesting to explore in manybody localized systems or disordered systems with Griffiths effects, where the butterfly velocities are zero and the light cones are logarithmic or subballistic.  The Jordan Wigner mapping allows spin operators to be written in terms of free Majorana fermions as follows: : σ x j = iγ 2j γ 2j+1 , σ z j = ( j−1 k=−∞ iγ 2k γ 2k+1 )γ 2j and σ y j = ( j−1 k=−∞ iγ 2k γ 2k+1 )γ 2j+1 . Below, we obtain analytic solutions for time-evolved operator for this Hamiltonian [Eq. (2)] within the Heisenberg picture. Denoting γ 0 (t) = n f n (t)γ n and γ 1 (t) = m h m (t)γ m , the time-evolved operator is σ x 0 (t) = iγ 0 (t)γ 1 (t) = i n<m F n,m (t)γ n γ m , where F n,m (t) = f n (t)h m (t) − f m (t)h n (t), and the out-of-time-ordered correlations are C xz (j, t) = 1 − 2 n≤2j,m≥2j+1 |F n,m (t)| 2 , C xx (j, t) = 1 − 2 n<2j |F n,2j (t)| 2 + |F n,2j+1 (t)| 2 + m>2j+1 |F 2j,m (t)| 2 + |F 2j+1,m (t)| 2 .
Next, we get the analytic solution of time-evolved operators γ 0 (t) = n f n (t)γ n and γ 1 (t) = m h m (t)γ m in the integrable Hamiltonian [Eq. (2)]. Plugging the candidate solution into the Heisenberg equation, it is straightforward to get the differential equations for the coefficients f n (t) where the initial condition is f n (0) = δ n,0 . After applying the Fourier transformation f 2n (t) = 1 2π π −π dq e −inq A(q, t), f 2n+1 (t) = 1 2π π −π dq e −inq B(q, t),