Two-step phantom relaxation of out-of-time-ordered correlations in random circuits

We study out-of-time-ordered correlation (OTOC) functions in various random quantum circuits and show that the average dynamics is governed by a Markovian propagator. This is then used to study relaxation of OTOC to its long-time average value in circuits with random single-qubit unitaries, finding that relaxation in general proceeds in two steps: in the first phase that lasts upto an extensively long time the relaxation rate is given by a phantom eigenvalue of a non-symmetric propagator, whereas in the second phase the rate is determined by the true 2nd largest propagator eigenvalue. We also obtain exact OTOC dynamics on the light-cone and an expression for the average OTOC in finite random circuits with random two-qubit gates.


I. INTRODUCTION
Complexity of quantum evolution is of wide theoretical and practical interest. It can be captured in different ways, one common idea is to quantify what is colloquially called scrambling of quantum information [1,2]. Details of what scrambling means depend on a particular situation, but broadly speaking one can quantify it in two ways: either by properties of the evolved state, for instance its entanglement under random circuit evolution [3], or by complexity of time-evolved operators as measured for instance by the out-of-time-ordered correlation (OTOC) functions [4][5][6]. In the present paper we shall study dynamics of OTOC O β (i, j, t), being equal to the Hilbert-Schmidt norm of the commutator, O β (i, j, t) = 1 2 AA † = 1 2 n+1 tr(AA † ), where A = [σ α i (t), σ β j ] is a commutator between two local traceless operators, one of them being evolved in time. It therefore measures how fast correlations spread from the spatial position i to j, and, more importantly for our discussion, also how fast σ α i (t) becomes "random". Due to its relative simplicity and relevance for quantum information OTOCs have been studied in very many different contexts. Limiting just to homogeneous systems, these include field theory [7,8], Luttinger liquids [9], and (chaotic) many-body systems [10][11][12][13][14][15]. In studies of quantum many-body systems one has to usually resort to numerics and that is why any exact results are greatly appreciated. Simplification that allows for analytical results often comes due to symmetries, either exact ones (e.g., integrable systems) or for instance effectively increasing symmetry by some averaging procedure. One such case when averaging brings simplification is that of random circuits where it has been shown that the average dynamics can be described by a Markov chain [16], leading to exact entanglement generation speeds for specific circuits [17]. Today many other theoretical approaches to random circuits are known. A prominent example is generic hydrodynamic description of operator spreading and OTOC dynamics [18,19], explicitly verified by exact results for random U(4) circuits. For a slightly different Brownian Hamiltonian evolution see Ref. [20]. Another powerful example of exactly solvable dynamics are so called dual-unitary circuits [21], among them also random dual-unitary circuits [22]. One of the distinguished features of dual-unitary circuits is that 2-point spatio-temporal correlations are nonzero only on the light-cone boundary [21], and one can use a powerful finite transfer matrix formalism. Similar is the case in integrable circuits [23] in which the gates satisfy the Yang-Baxter equation. For dual-unitary circuits OTOC decay exponentially with time close to the lightcone boundary [24]. Some simplification is possible also for certain small perturbations to dual-unitarity [25].
Many studies try to find some indication of chaoticity. Remembering that the notion of chaos is in classical systems defined as a property of the long-time limit, one might ask if such long-time complexity is somehow reflected also in the OTOC dynamics. The answer is not clear, what one can say however is that for lattice systems with finite local Hilbert space dimension it is not clear how to distinguish chaoticity from integrability via OTOCs. A possible approach is to get some measure of instability, like "quantum" Lyapunov exponents, from OTOCs dynamics. However, this is bound to fail for several reasons. One is that one might get an exponential behavior that is unrelated to chaos, for instance simply due to unstable fixed points [26][27][28]. Hydrodynamic behavior of the operator front might also look the same in chaotic and integrable systems [29] (for free systems see Ref. [30]). On top of it, in lattice models with finite local dimension, like chains of qubits, there is no obvious small parameter and so any possible exponential Lyapunov-like growth of OTOCs can hold only upto finite (short) times [31,32].
We are going to study OTOC dynamics in random quantum circuits, mostly in one-dimensional geometry and for qubits. In random circuits there is no dichotomy between integrability and chaos -random circuits can be thought of as being models of chaotic systems -and so we are not aiming at coming up with some chaoticity criterion. What we shall focus on is the long-time dynamics of OTOCs, specifically on how fast OTOC relaxes to its asymptotic value reached at long times that corresponds to a completely scrambled evolution. As we shall see, this will reveal interesting mathematical and physical properties.
Deriving a Markovian description of the average OTOC dynamics in random circuits we shall show that the relaxation rate typically exhibits a discontinuity at a specific time linear in the number of qubits. What is more, the relaxation time in this first phase, which is dominant in the thermodynamic limit, is not given by the gap of the Markovian matrix. Instead, it is given by a so-called phantom eigenvalue -a fake "eigenvalue" that is not in the spectrum. Illustration of such phantom relaxation is in Fig. 1. Looking at a particular OTOC O(1, 4, t) (see Eq. (2) and (11) for definitions), whose dynamics is given by a particular Markovian matrix M , we study its relaxation towards O(1, 4, t → ∞) ≈ 1. One can see that the relaxation proceeds in two steps: asymptotically at large t > t c ≈ n/2 one has the expected exponential decay ∼ |λ 2 | t , where λ 2 is the second largest eigenvalue of M ; however, for t < t c the exponential relaxation goes as ∼ λ t ph , where λ ph is a "phantom" eigenvalue that is larger than any true eigenvalue of M . Because in the thermodynamic limit t c diverges, the correct relaxation rate that one will observe at any finite time is not given by the spectrum of M , but instead by the phantom λ ph . Curiously, it turns out that in some cases λ ph is equal to the 2nd largest eigenvalue of a different circuit not related in any obvious way to M . Similar phenomenon of phantom relaxation has been recently observed also in purity dynamics [33]. Perhaps also related is an observation that in nonequilibrium dynamics described by the Lindblad equation the gap does not necessarily give the correct relaxation time [34,35], and of non-Hermiticity of transfer matrix describing integrable circuits [23].

II. RANDOM QUANTUM CIRCUITS
In this paper we deal with random quantum circuits defined on a system of n qubits. The unitary propagator U is a product of local elementary gates U i,j acting on qubit pairs (i, j), that is U = i,j U i,j . Every elementary step is, in turn, defined as a product of two independent one-site random unitaries V i and V j and a two-site uni- Two examples of random quantum circuits, where the product of elementary gates is ordered in a brick-wall (BW) pattern and in a staircase (S) pattern, can be seen in Fig. 2. As can be deduced from the name, the BW protocol is defined as a configuration where in each unit of time we first couple nearest-neighbor qubits (i, i + 1) with an odd i, then all pairs with even i. Apart from being widely studied for its simplicity, we mainly focus on this protocol because it turned out to be the fastest possible local scrambler of entanglement [33]. Another configuration that we will encounter in this paper is the S configuration. The S configuration consists of operators U i,i+1 , where at each step we increase i by 1. In the main part FIG. 1. Phantom relaxation of OTOCs. (a) In the thermodynamic limit OTOC relaxes to its long-time value as λ t ph , instead of |λ2| t , where |λ2| is the 2nd largest eigenvalue of the OTOC transfer matrix M , and λ ph is a phantom eigenvalue (which is not an eigenvalue of M ). (b) Spatio-temporal plot of OTOCs O(1, k, t) for n = 34 where one can see regions of relaxation with λ ph , and the asymptotic region with λ2. White vertical line marks the cross-section shown in (a). All is for a random circuit with the XXZ gate with az = 0.2 and the brick-wall protocol with periodic boundary conditions. we shall mostly focus on random quantum circuits acting on 1-dimensional (1D) chains of qubits with either open boundary conditions (OBC) or periodic boundary conditions (PBC); that is, qubits are distributed on a line (OBC) or on a circle (PBC).
One obtains various random circuits by different choices of a fixed two-site gate W i,j and the ordering of elementary steps. To distinguish various choices of W i,j we shall parametrize it in the following canonical form [36][37][38] where V α k are one-site unitary operators, σ x,y,z are Pauli matrices and a = (a x , a y , a z ) are three real parameters, which can be constrained to 0 ≤ a z ≤ a y ≤ a x ≤ 1. In this paper, we will be interested in the average dynamics of OTOCs generated by random quantum circuits. Due to randomness on single qubits (at every elementary step Illustration of a brick-wall (BW) protocol (a) and a staircase (S) protocol (b) on a qubit chain of size n = 8 with periodic boundary conditions. Blue boxes represent elementary steps Ui,j. Red dotted lines represent integer times, which are measured so that one unit corresponds to the action of one period of the random quantum circuit. One period of a BW protocol consists of local operators Ui,i+1 where we first act on qubits with odd i, then on qubits with even i. In a S protocol in one period we subsequently act with elementary steps Ui,i+1 starting from i = 1 and increasing i by one for each local operator.
we act with random unitaries V i and V j ) the choice of local operators V α k does not affect our averaged dynamics, so only the choice of the three real parameters (a x , a y , a z ) is what matters. To conclude, without loss of generality we can take our fixed two-site unitary to be w i,j , which is in turn parametrized by only three constrained real parameters 0 ≤ a z ≤ a y ≤ a x ≤ 1.

III. OTOC MARKOV CHAIN
We shall study out-of-time-order correlations defined as with σ γ k denoting the Pauli matrix at position k, and γ ∈ {x, y, z}. The time-evolved Pauli matrix is obtained as σ α i (t) = U † σ α i U . OTOC thus measures how correlations between two initially localized operators spread in the system. Its minimal value is 0 for operators that commute, i.e., until σ α i (t) begins to overlap with σ β j , whereas its maximal value is 2 reached for e.g., σ β j = σ x j and σ α i (t) = σ y j . If σ α i (t) at large times randomly spreads over all available operator space the average OTOC will converge towards its thermal value O ∞ ≈ 1 (see Appendix A 1). We are going to study how OTOCs converge to this long-time stationary value. Note that often the name OTOC is used just for the 2nd term in Eq. (2), whose asymptotic value goes towards 0. Because we will be interested in relaxation we will in fact study this 2nd term.
It has been shown that averaging over one-site random unitaries leads to a Markov chain description of the evolution of the average purity [39,40]. Because OTOCs are, similarly as purity, also quadratic in the time-evolved operator, their average evolution can also be written in terms of a Markov chain. This has been done for the special case of a random U(4) elementary step U i,j in Ref. [18], whereas we derive the Markovian matrix description for a protocol consisting of an arbitrary twoqubit W i,j conjugated by independent single-qubit unitaries.
The derivation relies on the fact that it is possible to express OTOCs as a linear combination of all possible purities of a system of n qubits. Writing the operator σ α i (t) in the basis of Pauli strings with coefficients a σ (t) where we use the label σ = (σ 1 , σ 2 , . . . , σ n ) while the product of Pauli matrices on all sites is denoted by σ = σ 1 σ 2 · · · σ n , and where for brevity we defined two sets that will be useful in specifying various summations. For instance the sum in Eq. (4) runs over all Pauli strings σ = (σ 1 , σ 2 , . . . , σ n ) except for those having σ β j or 1 at the site j. We wish to relate a vector containing all possible OTOCs O(i, j, t) for every position j to a vector of purities through a linear transformation. To obtain purity we write the density operator in terms of Pauli strings coefficients c σ , ρ(t) = 1 √ 2 n σ c σ σ. Purity I A , which measures pure-state entanglement between two complementary subsets of qubits denoted by A and B (consisting of n A and n B qubits, respectively), is then Expression (6) is invariant with respect to an arbitrary permutation of the three Pauli matrices at any site. In other words, it is only the totally symmetric sum of c 2 σ for all three Pauli matrices that matters for purity. For instance, for a system of two qubits with subsystem A being the 1st qubit, we have I A = 2c 2 (1,1) + 2(c 2 (σ x ,1) + c 2 (σ y ,1) + c 2 (σ z ,1) ). So instead of bookkeeping all 4 2 coefficients c 2 (σ1.σ2) it is enough to keep track of only 2 2 combinations of them, which we can neatly pack into a two-site vector (for definition of S 1 see Eq.(5)) We can obtain purities for all possible bipartitions of two qubits from components of Φ, specifically, if the 1st qubit is in A we have I A = 2Φ 0 + 2Φ 1 , whereas if the 2nd qubit is in A one has I A = 2Φ 0 + 2Φ 2 , where we labeled the 4 components in Eq. (7) by Φ 0,1,2,3 . Generalizing Φ to n qubits it will have 2 n components that we label by bit strings s = (s 1 , . . . , s n ), where s j ∈ {0, 1}, with the components being To shorten the notation we shall occasionally also use the integer value of the bit string s instead of specifying the full s = (s 1 , . . . , s n ), as s ≡ n j=1 2 j−1 s j . Purity for an arbitrary bipartition is now given by a particular component of vector Φ I obtained as Specifically, the component [Φ I ] s is equal to the purity for a bipartition in which the subsystem A consists of qubits for which s j = 0, i.e., the bit s j encodes the subsystem in which the j-th qubit is.
Ref. [39] showed that it is possible to write the evolution of purities Φ I averaged over one-site Haar random unitaries as a Markov chain. Abusing notation and from now on using Φ I (t) to denote the average purity after t steps of our random circuits (2), one has The transfer matrix M describing one period of our circuit is a product of matrices M i,j , one for each elementary step U i,j [39]. For example, a transfer matrix describing t 2 periods of a BW PBC circuit on n = 4 qubits would be Note that because the two-site gates W i,j are the same for all steps all transfer matrices are independent of time.
Looking at the expressions for OTOC in Eq. (4) and Φ in Eq. (8) we can see that they look rather similar. Because we know how average purities are evolved (10), we also know how to evolve Φ(t), namely defining Φ(t) = . This will in turn lead us to the evolution of OTOC.
To achieve that let us rather look at the OTOC averaged over three possible σ β j , Note that the dependence on site index i is implicitly hidden in the expansion coefficients a σ (t) (3) of the initial σ α i . Using Φ for a vector defined as in Eq. (8) but for coefficients a σ , and formally defining a vector Φ O by where s j = 1 and s k =j = 0. Therefore, n components of Φ O are equal to OTOCs while the other 2 n − n components are some other combinations of a 2 s not related to OTOCs. Note that the choice of A O is not unique; the two 1 in the top row take care of summing over both sets S 0 and S 1 for sites k = j in Eq. (11), while the 4 3 in the 2nd row accounts for an overall prefactor accounted by a single bit s j being 1, ie., summation only over where the vector e k has components [e k ] s = δ s,k . The vector Φ containing coefficients of σ α i (t), instead of ρ(t), is propagated in exactly the same way as for average purity [41], that is, because Φ = A −1 where M is the transfer matrix propagating purities. Using M calculated for random circuits and an arbitrary W i,j parameterized by (a x , a y , a z ), Eq.(13) from Ref. [33], we immediately get the transfer matrix M describing the evolution of average OTOC under one elementary step, Here c ± = 1 12 (9 ± 2u − v) and d = 1 6 (v − 3), with u = cos (πa x ) + cos (πa y ) + cos (πa z ) and v = cos (πa x ) cos (πa y ) + cos (πa x ) cos (πa z ) + cos (πa y ) cos (πa z ). Note that e 0 is a trivial eigenvector of M corresponding to the eigenvalue 1, i.e., is a stationary state. However, M has another nontrivial eigenvector with λ = 1 containing the asymptotic stationary values of OTOC O ∞ .
Summarizing, the matrix M i,j (15) acts nontrivially only on 2 sites i and j and is written in the basis of bit strings ordered as {0 i 0 j , 1 i 0 j , 0 i 1 j , 1 i 1 j }. To get the matrix propagating OTOCs for the complete circuit for one unit of time one must multiply appropriate 2-site M i,j in the same order as the gates are applied in the protocol, for instance, for a n = 4 site BW PBC circuit one has M = M 4,1 M 2,3 M 3,4 M 1,2 . The state Φ O on which full M acts has 2 n components. Fixing the site i in the initial OTOCs O(i, j, 0), the initial vector (13) The transfer matrix description of the average OTOC dynamics (14) that we obtained offers several advantages. First, it gives a neat analytical description on which one can use standard tools of analyzing Markov chains, like for instance trying to connect the spectral properties of M to the asymptotic relaxation of OTOC to its infinitetime values. Second, it also greatly simplifies numerical simulations of OTOCs -instead of, e.g., explicitly simulating the dynamics of operators, averaging over different realizations, one can directly simulate the average OTOCs dynamics.

IV. EXACT DYNAMICS ON THE LIGHT-CONE
In this section we will obtain the exact dynamics of OTOCs on the light-cone, from which we will be able to determine in a very simple way the set of two-site W that result in maximum velocity circuits. Maximum velocity quantum circuits were defined in Ref. [24] as circuits where the butterfly velocity v B [7,[42][43][44] equals the Lieb-Robinson velocity v LR [45]. The Lieb-Robinson velocity determines the causality light-cone of which boundaries are at positions k = i ± v LR t (i is the location of σ α i (t = 0)), so that OTOC O(i, j, t) with parameters (j, t) outside the light-cone vanish. In a random circuit its value is determined solely by the circuit geometry and is for instance v LR = 2 for the BW configuration. The butterfly velocity on the other hand is determined as v B = |j − i|/t min , where t min is the minimal time when O(i, j, t) ∼ 1 at fixed large |j − i|. Contrary to v LR the butterfly velocity depends both on the geometry and on the choice of the gate W . For most random quantum circuits, v B = v LR . An illustration of these two velocities can be found in Fig. 3.
Let us focus on a random quantum circuit with a brick-wall configuration of gates acting on an infinite system, n → ∞. We would like to compute quantities O(i, i ± v LR t, t). Due to symmetry it is enough to consider OTOC only at the right light-cone. We also limit ourselves to odd i such that the light-cone edge is at k(t) = i + 2t (for even i one would have It is important to note that due to causality all val- Iterating this by half-steps to smaller times until we reach O(i, i, 0) = 4/3 one obtains the OTOC on the right light-cone. Similar procedure works also on the left light-cone and even i, resulting in Looking at the left light-cone at odd i, or the right lightcone and even i, one instead gets The additional term c + in Eq. (18) comes from the interaction OTOC on the light-cone therefore decay exponentially with the rate 2 ln c − , hence one gets v B = v LR = 2 iff c − = 1. Solving c − = 1 for a x , a y , a z we obtain a x = a y = 1 and an arbitrary a z , which corresponds to dual-unitary circuits [21] and for which one can explicitly calculate all 2-point correlations which are nonzero only on the light-cone boundary [21]. This means that taking W from the dual-unitary set of gates (i.e., so-called XXZ gates) is the only choice leading to the maximum velocity 0 10 20 Average OTOC for the BW random circuit with PBC and the XY gate (a = (1, 1, 0)), n = 28. Operator relaxation towards long-time thermal asymptotics will be studied by observing how a fixed-position O(i = 1, j = 6, t) converges towards O∞ with time (inset).
random circuits (of the type studied in this paper), i.e., circuits for which OTOC do not decay along the lightcone. The same set of maximum velocity gates was also obtained in Ref. [24] for circuits without one-site random unitaries. Besides identifying maximum velocity gates our simple derivation also gets us the exact dynamics of OTOC on the light-cone for arbitrary gates W . Note that for circuits with a dual-unitary two-qubit gate W one can also get a closed expression for the OTOC decay in the vicinity of the light-cone, for non-random circuits see [24], for random [22]. We also observe that the same set of gates, except at a z = 1, results in the maximal possible entanglement scrambling speed [33]. The gate with a z = 1 is the SWAP gate and is special. The OTOC dynamics for a random circuit with the SWAP gate is trivial because the transfer matrix M i,j itself (15) is equal to a SWAP gate resulting in O(i, j, t) that is non-zero only on the light-cone, while at the same time such W produces no entanglement.

V. CONVERGENCE RATE
Under the application of a random quantum circuit the initially localized operator will spread in space, causing OTOC to increase from being zero outside of a lightcone to a nonzero value inside it. Often one is interested in this ramp-up of OTOC as for instance measured by the butterfly velocity. We shall instead investigate the late-time convergence rate of OTOC O(i, j, t). That is, we are interested in how fast O(i, j, t) at some fixed i and j relaxes towards its final value, see Fig. (4) for an illustration.
The asymptotic value O ∞ is reached at long times when the time evolved operator σ α i (t) becomes a uniform mixture of all possible Pauli strings (identity excluded) on n qubits, i.e., when the propagator U resembles a random unitary. A derivation of O ∞ can be found [46] in A 1 and gives us If the eigenvalues of the transfer matrix M are gapped away from 1, which indeed is the case, we expect that OTOC exponentially relax to their asymptotic value O ∞ as Our main object of study is the convergence rate r. Considering that OTOC are propagated by M t one might think that the convergence rate will be determined by the 2nd largest eigenvalue λ 2 of M . However, as we shall see, this is not always the case.
In section V A we shall first discuss protocols in which at each step we randomly pick a pair of qubits on which we act, that is protocols with a random ordering of gates. For those we will see that indeed OTOC decay exponentially with the convergence rate r given by the second largest eigenvalue of the transfer matrix r = − ln |λ 2 |.
In sections V B and V C we shall on the other hand study protocols with a nearest-neighbor deterministic order of gates, mostly the BW or the S configuration (Fig. 2), and find, similarly as for purity [33], that the convergence rate can be either equal or smaller than − ln |λ 2 |. In Appendix C we also numerically demonstrate that in the thermodynamic limit one does not need explicit averaging over single-qubit unitaries, i.e., dynamics is self-averaging and therefore one will obtain the same results also for a single circuit realization.
Relying on a map from OTOC to a partition function of an Ising-like model found in [18] we analytically compute OTOC for the BW PBC and BW OBC in the case where every elementary step of the circuit is independently drawn from the Haar measure on U(4). Contrary to previous literature, where the analytic expression for OTOC was obtained in the TDL [18,22,24], in Appendix B we present new results in finite systems with either OBC or PBC.

A. Random protocols
Random protocols studied here are defined as random quantum circuit where at every elementary step we couple two qubits chosen randomly. We have two different possibilities: a) at each step we uniformly choose one of the all possible n qubits, for example the i-th one, and we act with gates M i,i+1 , and b) at each step we randomly choose two qubits i and j and we act with M i,j . We call the former case the random nearest neighbor protocol (r.n.n.) and the latter scenario the all-to-all coupling.
In the r.n.n. case, the average elementary step can be written as the average over all possible choices of i, with L = n − 1 or L = n, depending on the boundary conditions. Similarly, for the all-to-all case we obtain with L = n(n − 1)/2. The transfer matrix propagating OTOC for one unit of time is M =M L . Because each M i,j is just a similarity transform of the purity M i,j , Eq. (14), the spectrum of M is identical to the spectrum of purities transfer matrix M . Furthermore, as was shown in [33], the average elementary steps M i,j propagating purity can be linearly transformed to a real symmetric matrix. Therefore the spectrum of M is equal to the spectrum of a symmetric purity matrix. This is important because the spectrum is real with orthogonal eigenvectors. The spectral decomposition of a Hermitian M takes the form M = k λ k |v k v k |, and we can expand |Φ as |Φ = k c k |v k with c k = v k |Φ being bounded by |c k | 2 ≤ Φ|Φ . The time iteration is For a Hermitian M there are therefore no surprises; if the 2nd largest eigenvalue λ 2 is gapped away from other eigenvalues the asymptotic decay rate will be given by λ 2 and will kick-in at a system size independent time. For the r.n.n. protocol the 2nd largest eigenvalue has been computed numerically for arbitrary gates [47] and analytically for a few Clifford gates [17]. For λ 2 in the all-to-all case and Clifford W see Ref. [17], for arbitrary W Ref. [33].

B. Brick-wall protocol with PBC
For protocols with a deterministic order of gates things can and will be completely different. The decay rate will not necessarily be given by λ 2 . Remember that for a deterministic order of gates the transfer matrix is just a product of corresponding two-site M i,j , for instance, for a 4 qubit BW protocol it is M = M 4,1 M 2,3 M 3,4 M 1,2 . The difference compared to random protocols is that a product of symmetric matrices needs not to be symmetric. As a consequence, the eigenvectors of such M are not orthogonal, c k are not upper bounded, and, as has been seen in purity evolution [33], the relevant decay can differ from |λ 2 | t .
In the following we shall plot how the value O(1, j, t) behaves for a fixed position j. We will always fix i = 1, because OTOC in PBC circuits depend only on |j−i|. We will plot values of |O(1, j, t)−O ∞ | and the time derivative There is a phantom eigenvalue: initially, the rate is given by λ2 for a BW OBC circuit (red dashed line for n = 30). At late times the rate is instead equal to − ln |λ2| for a BW PBC circuit (green dashed line for n = 20). The inline plot shows a transition in the exponential decay of the same data, including red and green dashed exponential functions corresponding to red and green rates in the main plot.
in order to investigate OTOC convergence rate (note that r(t → ∞) = r from Eq. 20).
Let us start with a generic two-qubit gate Data for O(1, j = 7, t) is shown in Fig. 5 and demonstrates that OTOC converge to their final value with a rate different than − ln |λ 2 | (for W g one has |λ 2 | ≈ 0.72 for n = 20). The rate is (initially) smaller, as if there would be an eigenvalue larger than λ 2 -a phantom eigenvalue. Such slower decay persists up to times that are proportional to the system size. The value of the phantom eigenvalue is equal to the second largest eigenvalue of the transfer matrix for the BW OBC circuit. Remember that we are looking at a circuit with PBC, not OBC, nevertheless, it is perhaps expected that for initially localized quantities and until the boundary conditions (PBC) influence OTOC dynamics, the convergence rate is given by λ 2 of BW OBC (see also next section). Namely, choosing the initial vector localized roughly equally far from the left and right boundary (i ≈ n/2), the dynamics generated by BW PBC or OBC circuit is identical up to times t ≈ n/4. Therefore, what might be surprising is that λ 2 of M for BW with PBC and OBC are different. Looking at OTOC on a different site, j = 7, one might observe a slightly different graph from Fig. 5, however the behavior remains qualitatively the same: at early times that scale as ∼ n the dynamics is always determined by a phantom eigenvalue, which is the same for every j, whereas at late time the dynamics is given by the second largest eigenvalue of M . Of special interest are gates with canonical parameters a = (1, 1, a z ), a z < 1 (dual unitary W ). We purposely skip a z = 1 because of its trivial dynamics. Contrary to the generic gates W g , for dual unitary gates we will see that early-time dynamics is always determined by λ 2 of an S PBC circuit (shown in Fig. 2(b)), which though is always larger than λ 2 for BW PBC. One will therefore again have a situation where the relevant relaxation rate is not given by λ 2 of the BW PBC circuit.
We will first take a look at a circuit with the dualunitary gate with a z = 0.2. Because there are some differences between even and odd j at later times, essentially due to even/odd effects of the light-cone boundary position (see Sec.IV), we show in Fig. 6 how O(1, j, t) converge for j = 7 in (a) and (b), as well as for j = 8 in (c) and (d). Looking at Fig. 6(c) that focuses on short times we can see that r is zero until the right light-cone boundary hits the site j = 8. We assume that j − i is odd and j − i < n − (j − i), i.e., the first information that hits the site j comes from the right light-cone boundary, not from the wrapped-around (PBC) left light-cone boundary. OTOC and the rate are therefore zero until t ≈ (j − i)/2. After that r stays at a value that is not given by |λ 2 | of the BW PBC transfer matrix, but rather by λ 2 of the transfer matrix for the S PBC (red dashed line) and for which we have a conjectured analytical form, see Ref. [48]). At the time t c = (n + 1)/2 − (j − i)/2, determined by the time when the left light-cone boundary hits the site j, the rate suddenly transitions to its ultimate asymptotic value given by λ 2 of the BW PBC (green line). In Fig. 6(d) we can see that this rate stays roughly constant upto small modulations at times larger than t c , e.g. at t ≈ 20 for n = 34. They happen at times of successive light-cone boundary wrappings (for more details see next paragraph). There are some interesting differences for odd j (Fig.6(a,b)). Specifically, because for odd i = 1 the left light-cone boundary is at even sites and therefore never overlaps with an odd j, the rate has a transition to its asymptotic form only when the right light-cone boundary hits the odd site j = 7 for the 2nd time (due to PBC). This happens at t c = (j − i)/2 + n/2, e.g., t c = 20 for shown n = 34 and j = 7. As one can see from O in Fig.6(b), the rate itself does not change; rather the OTOC exhibits a jump. As we shall see in the next paragraph, the ultimate asymptotic decay is nevertheless still determined by λ 2 of the BW PBC circuit.
In order to be able to better explore those spikes we shall next look at a dual unitary gate with a z = 0.6, because OTOC decay slower and we are able to simulate longer times (rounding errors of double precision floating point numbers ultimately limit the smallest O we can calculate). Results are shown in Fig. 7. From the figure we learn that the convergence of the initial rate to that given by the eigenvalue of the S PBC protocol is rather slow with n; smaller system sizes have rates that do not yet converge to − ln |λ 2 | S−PBC . There are also small kinks in the decay of |O(1, j, t) − O ∞ | that are due to light-cone wrapping boundaries. Overall though the rate changes only once from the initial one to the asymptotic − ln |λ 2 | BW−PBC at the already dis- cussed time that is proportional to n (see frames (b) and (d)). We can now also clearly see several spikes at times when the light-cone boundary wraps around the system multiple times. Specifically, starting with an odd i = 1 the right light-cone boundary will hit a site at odd j at times t = (j − i)/2 + kn/2, where k in an integer (blue vertical lines in the Figure), whereas the left light-cone boundary will hit it at times t = kn/2 − (j − i − 1)/2 (black vertical lines). There is a slight asymmetry between the effects of left and right light-cone boundary: spikes due to the left one are prominent only for even j which comes due to an asymmetry in the behavior of OTOC on the light-cone boundary, Eqs.(1718)for even j the left light-cone has an additional factor c + = 1 3 (1 + cos (πa z )) = 1 − |λ 2 | S−PBC . We also observe that r decreases with increasing a z ( [48]). This means that the fastest relaxation of OTOC among dual-unitary gates is obtained for the circuit with the XY gate, i.e., a = (1, 1, 0), when one has r = log 3 at t n .
We have seen that in all BW circuits with PBC, for generic gates as well as dual unitary gates, the relevant relaxation rate of OTOC that holds upto times of order ∼ n, i.e., until OTOC become exponentially small in system size, is not given by the 2nd largest eigenvalue of the BW PBC transfer matrix. For the generic gate W g the rate was given by |λ 2 | BW−OBC which is larger than |λ 2 | BW−PBC -a phantom eigenvalue phenomenon. One might be inclined to justify this result based on a trivial fact that the choice of boundary conditions of course does not matter up to times that are proportional to ∼ n. Until boundary effects kick in OTOC evolve as they would in a BW OBC system (if the Pauli matrix at time t = 0 is positioned "far enough" from the boundaries). This however is not really a full explanation; remember also Operators in the same period of the S OBC circuit are represented with the same color, meanwhile operators in the same period of BW OBC are labeled by the same parameter t (see Fig. 2). Due to causality, all gates outside the future lightcone starting from i = 6, and the past light-cone originating from j = 8, do not matter (crossed out gates). By stacking together S OBC protocols one obtains the same set of gates as for BW OBC.
that for dual unitary gates the rate (phantom eigenvalue) was given by |λ 2 | S−PBC and not |λ 2 | BW−OBC , despite the evolution still being the same as it would be in the BW OBC circuit.
In the next section we shall study circuits with OBC. Based on results presented so far we can predict that for BW OBC with generic gates one will have no phantoms, whereas we expect to see a phantom rate given by |λ 2 | PBC−S for BW OBC circuits with dual unitary gates.

C. OBC protocols
Let us first stress one important property of a family of OBC protocols that comes about due to the locality of the initial vector Φ O (t = 0) (Eq. 13). We conjecture that OTOC dynamics is not influenced by permutations of elementary gates in one period of the BW OBC circuit. For example, looking at a BW OBC protocol one could permute the order of elementary steps in one period and obtain an S OBC circuit without affecting the average OTOC dynamics, e.g., its decay rate.
To support this claim, we will show that for a fixed W the BW OBC protocol generates the same OTOC dynamics as the S OBC protocol upto a constant time-shift.
We will rely on Fig. 8 to explain the equivalence. One can easily see that stacking together S protocols one obtains a circuit of the form shown in Fig. 8, i.e., a brick wall protocol in the middle (between t = 2 and t = 3 in the Figure), and two "triangles", one at the top right (gates after time t = 3 in Fig. 8) and one at bottom left (gates before time t = 2 in Fig. 8). Let us focus on the calculation of O(i, j, t). Due to causality the evolved local operator vanishes outside the light-cone starting from the i-th qubit. We are interested in the component of Φ O (t) representing O(i, j, t), i.e. the 2 j−1 -th one, this means that also operators in the past light-cone starting from the qubit j vanish. The relevant gates are therefore those inside the two light-cones, i.e., in Fig. 8 the gates that are not crossed. The same set of relevant gates would be obtained acting with a BW OBC protocol. The only difference between S OBC and BW OBC circuits is a time-shift that comes from the difference between v LR of the two circuits. This is reflected in the fact that O(i, j, 1) = 0 for arbitrary j in S OBC circuits, whereas using a BW OBC protocol we have O(i, j, t) = 0 at all times smaller than ∆t ≈ |j − i|/2, which is equal to the time-shift between the two protocols.  (6, 8, 6). The time-shift is constant and depends only on the value j − i (and can be a half-integer). This can be seen also in explicit numerical data in Fig. 9 where O(1, 8, t) for BW OBC circuit (triangles) is the same as O(1, 8, t − 3) (squares) obtained for the S OBC.
Using similar arguments one can see that if one iterates an arbitrary OBC configuration, that is a protocol in which each nearest-neighbor gate is applied exactly once per unit of time, one always gets a brickwall pattern of gates. Therefore one can show that the OTOC of local operators and any OBC protocol is upto a time-shift equal to the one in say BW OBC circuit. We have also checked numerically on a few examples of random gate permutations that this is indeed the case. We remark that in Ref. [33] it has been shown that the spectra of transfer matrices M for a single iteration are the same for all OBC protocols.
This equivalence though holds only for OBC. For instance, for the XXZ gate S PBC and BW PBC protocols can behave rather differently, see circles and stars in Fig. 9, what is more, the BW PBC circuit exhibits a phantom relaxation. On the other hand, the S PBC with the XXZ gate does not exhibit a phantom relaxation, while the S PBC with the generic gate W g does (data not shown).
Regarding possible phantoms in the OBC setting we can see in Fig. 9 that for dual unitary gates BW OBC does exhibit a phantom (the initial rate is given by |λ 2 | of the S PBC), while for generic gates it does not (data not shown), which is expected (we have seen in Fig. 5 that the rate for BW PBC was given by |λ 2 | of BW OBC, and the two O(i, j, t) should agree until t ∼ n). Let us have a closer look at the dual unitary gate a z = 0.5 and BW  OBC protocol. From data in Fig. 10 we indeed see that there is a phantom -the initial rate is smaller -and that there are, similar as in the PBC case (Fig. 7), again spikes in the rate. Those spikes are associated with jumps in the relaxation of OTOC (frame (b)) that happen every time the reflected right light-cone returns to site j, i.e., at times kn − (j − i − 1)/2. All different cases of random circuits and their relaxation are summarized in Table I [33], i.e., λ2 is the same for BW and S). The 2nd largest eigenvalue for XXZ gates and PBC protocols are conjectured to be equal to |λ2| Spbc = (2 − cos(πaz))/3 and |λ2| BWpbc = (2 − cos(πaz)) 2 /9 [48].
for the OBC in the case of generic gates, or to |λ 2 | S−PBC in the case of XXZ gates. Predicting when does one have a phantom relaxation and what is this λ ph equal to is not simple. In some cases λ ph is equal to the second largest eigenvalue for a circuit with OBC, like for the generic gate W g with PBC, while for the XXZ gate with S OBC protocol it is instead equal to |λ 2 | for the PBC.
There is also the case of XXZ gates in the BW configuration where the phantom eigenvalue is equal to |λ 2 | for a whole different circuit, namely the S configuration. Understanding in detail the physics of phantom relaxation therefore remains an open problem.

VI. CONCLUSION
We have derived a Markovian propagator for the average out-of-time-ordered correlations of local operators in random quantum circuits in which each two-qubit transformation is composed of a fixed two-qubit gate W and two random single-qubit unitaries. This allows us to get an exact expression for OTOC on the light-cone and any W .
We then focus on the asymptotic relaxation rate at long times with which OTOC relaxes to its long-time average corresponding to a completely scrambled evolution. Similarly as in the case of purity relaxation [33] we find that this OTOC relaxation rate is in many cases not given by the second largest eigenvalue of the Markovian matrix M that governs dynamics of OTOCs. One has a so-called phantom relaxation -a relaxation where the approach to the steady state asymptotically goes as λ t ph , with λ ph being some number that is, opposite to "expectations", not equal to any of the eigenvalues of M . Be-cause λ ph is in fact larger than any nontrivial eigenvalue |λ j | of M we call it a phantom eigenvalue. In short, M encodes all the information about OTOCs evolution but its eigenvalues do not give the correct relaxation rate in the thermodynamic limit, despite |λ 2 | being gapped away from λ 1 = 1.
Such phantom relaxation proceeds in two steps, where in the first step that lasts upto times that are linear in system size the rate is given by the phantom eigenvalue, while in the second it is eventually given by the second largest eigenvalue λ 2 . Because the transition time between the two regimes diverges in the thermodynamic limit one has a situation where at a fixed system size and t → ∞ one gets the naively expected (but thermodynamically incorrect) relaxation as |λ 2 | t , while in the correct thermodynamic limit of first taking the system size to infinity and only then time to infinity one will observe the relaxation rate given by the phantom eigenvalue. The phenomenon occurs because the limits t → ∞ and n → ∞ do no commute, while mathematically it comes about because the transfer matrix M is not symmetric, resulting in spectral expansion coefficients that blow up with system size [33], see also Refs. [23,34,35,49,50] for other situations where that occurs.
We find such two-step phantom relaxation for brickwall circuits with dual unitary (i.e., XXZ type) as well as with generic two-qubit gates, and for periodic or open boundary conditions. Phantoms are also found for the staircases configuration with open boundary conditions and dual unitary type gates; see Table I for an overview. We numerically observe that the phantom eigenvalue λ ph is equal to the 2nd largest eigenvalue of M of a different circuit that can have different boundary conditions as well as different gates configuration (theoretical reasons for that are at present not understood).
For circuits with open boundary conditions we demonstrate that up to a time-shift all different circuit geometries, i.e., brick-wall, staircases, etc., have the same OTOC dynamics. We also numerically verify that the dynamics is self-averaging, that is, one will get a phantom relaxation even for a single random circuit realization, and even without spatial or time independence of single-qubit random unitaries. An explicit randomness therefore seems not to be essential. This leaves an interesting possibility that a similar phenomenon could be observed also in in other systems, for instance in Floquet models.
The important message therefore is that: (i) when one deals with finite non-Hermitian matrices the leading eigenvalue might not give the correct asymptotic dynamics, and (ii) that this leads to a two-step relaxation process with a sudden discontinuous transition in the relaxation rate at a time when the light-cone hits the site in question for the second time (either due to a reflection from a boundary for open boundaries, or due to a wrapping around for periodic boundary conditions). On the mathematical level it is therefore due to the fact that boundary conditions apparently can affect the leading relevant eigenvalue in a nontrivial way. While we do obtain some exact properties of the Markovian matrix, like a conjectured exact expression for λ 2 in the case of periodic boundary conditions, much remains to be under-stood, in particular under which physical conditions one gets such a two-step relaxation. Support from Grants No. J1-1698 and No. P1-0402 from the Slovenian Research Agency is acknowledged.

Appendix B: U(4) exact results
This appendix is dedicated to analytic solutions of OTOC time-dependence in a random quantum circuit with either BW OBC or BW PBC configurations and the choice of random two-site gate. Choosing a random two-site gate means that every elementary step is composed by a random gate uniformly drawn from the unitary group U(4) according to the Haar measure. All shall be calculated independently on the local Hilbert space dimension q -instead of qubits we shall now work with general qudits. The analytic solution when q = 2 can be used as a non-trivial check of the exactness of our newly derived Markov chain. Namely, average dynamics for random two-site gate can be obtained by setting u = 0, v = −3/5 in the transfer matrix from Eq. 15.
In order to analytically solve the time evolution of OTOC, we will heavily rely on a reduction of OTOC dynamics to an Ising-like partition function found in [18]. Before we continue with solutions for BW OBC and BW PBC, we will give a brief overview of the reduction obtained in the previously mentioned paper. We won't give a detailed description of the reduction, but rather explain what one has to do in order to obtain the final result.
The Ising-like model of which we will calculate the partition function is obtained by replacing all two-site elementary steps U i,i+1 in the random quantum circuit with two-level spins (s ∈ {↑, ↓}). When dealing with a BW protocol, one obtains a grid of spins (see Fig. 11). Contrary to the main text, here we measure time τ such that one unit corresponds to a row of the BW circuit. The choice of OTOC O(i, j, τ ) reflects itself in the upper and lower boundary of the spin grid: at the lower boundary one must place an up-spin at position (i, i + 1) (for simplicity we shall assume that i is odd), at the upper boundary one must place an up-spin at position (j, j + 1) or (j, j − 1), depending on time τ and parity of j. Moreover, boundary conditions dictate that all other spins at the lower boundary must be ↓.
Reference [18] gives a detailed description on how to obtain such constraints, together with a derivation of the interaction laws between spins at different positions. Namely, spins interact with a three-body interaction. Take three spins, (i, i + 1) and (i + 2, i + 3) at time τ , and (i + 1, i + 2) at time τ + 1. The weight of the interaction is 1 if all spins are equally oriented, q/(q 2 + 1) if the spins at time τ are different and 0 otherwise. It is now clear that in order to obtain an up-spin at position (j, j + 1) at the upper boundary the spin grid must contain two domain walls, determining the boundary between up-spins and down-spins (see Fig. 11). The authors of [18] showed that OTOC are now recovered with the following equation Interaction between spins are colored according to their weight (1 for same spins and q/(q 2 + 1) for different spins). In the case shown here is represented a system with 12 qudits up to time τ = 8. At the lower boundary, we start with the up-spin located at the 5th (or equivalently 6th) qudit. The domain walls, describing the boundary between up and down-spins, must contain the qudit j at time τ = 8 (for the example shown j must be taken from the interval [3,8] in order to the picture to represent a possible configuration). For the configuration shown in the appropriate term of the sum in Eq. B1 one should take "int. = 14" and "width = 3".
The first term in Eq. B1 comes from the lower boundary condition and will be present in all configurations (PBC and OBC) that we study here. The sum runs over all possible ends of the domain walls on the upper boundary with width "width". With "int." we denoted the number of interactions between differently oriented spins. In the sum we must count all possible domain walls with given "ends" (#walls). For a more graphical explanation of Eq. B1 see Fig. 11.
In [18], OTOC were calculated for an infinite size system (infinite number of qudits). The authors obtained  corresponding to the convergence rate of OTOC in an infinite system (see Eq. B7).
• All configuration with width always smaller than n/2. In this case ( q q 2 +1 ) int. = ( q q 2 +1 ) τ −1 . The number of different spin configurations will be computed as the sum over all domain wall endpoints (so that up-spins contain the qudit j) that do not reach 0 minus all domain walls that at arbitrary times hit the right boundary n, that is H(τ 0 , n). The factor (q 2 ) width depends on the endpoint of the domain wall. Such configurations contribute to the partition function as A detailed explanation on how the numbers of paths were obtained will not be given here, because it deviates too much from the topic. OTOC O(1, j, τ ) for BW OBC with random 2-site gates is equal to the sum of Eq. B12 and Eq. B13. The obtained result was used to plot the time dependence of OTOC O(j = 10, t) for different values of q for a system with 20 qudits, see Fig. 12.

OTOC in a finite system with periodic boundary conditions
The case of BW PBC is more complicated. In this case our analytical result does not give any computational advantage over the Markov chain iteration method. We will see that the final results will be given by a recursion, which is time consuming during numerical evaluations. However, it is still useful to obtain an analytical result for completeness. As in the OBC case, here we also differentiate configuration that are wide n/2 and configurations that are never wider than n/2.
• The number of domain walls with width n/2 at time τ c will be computed recursively. LetN τc (u 0 → u 1 , v 0 → v 1 ) denote the number of domain walls pairs, starting at (u 0 , v 0 ) and ending at (u 1 , v 1 = u 1 + n/2). All domain walls inN τc (u 0 → u 1 , v 0 → v 1 ) must never be wider than n/2 at times τ < τ c . To shorten the notation, we shall always take u 0 = 0 and v 1 = 1, i.e., i = 1, and writeN τc (u 0 → u 1 , v 0 → v 1 ) asN where binomial coefficients n k = 0 for n < k or k < 0. The recursion in Eq. B14 ends at τ c = n/2 by takingN (B16) • Domain walls with width always smaller than n/2 will be computed with the help ofN (τc) u1 . The number of desired paths will be calculated by subtracting domain walls that reach width n/2 from the total number of all possible paths with v b = u + mod (D − mod (−τ + 2u + n/2, n), n)/2 + 1 and D = n/2 − 1 + j (note that we always take i = 1).  (26)). One can see that even though n = 24 is not yet in the thermodynamic limit, at tc (vertical dashed line) the relaxation rate does change also for individual circuit realizations.
The final result is obtained by summing Eq. B16 and Eq. B17. Note that due to recursive terms, this analytical result does not give any substantial computational advantage over the Markov chain iteration method, however we do not claim that there is no way to count the number of domain walls in a simpler way.