Arbitrary relaxation rate under non-Hermitian matrix iterations

We study the exponential relaxation of observables, propagated with a non-Hermitian transfer matrix, an example being out-of-time-ordered correlations (OTOC) in brickwall (BW) random quantum circuits. Until a time that scales as the system size, the exponential decay of observables is not usually determined by the second largest eigenvalue of the transfer matrix, as one can naively expect, but it is in general slower -- this slower decay rate was dubbed"phantom eigenvalue". Generally, this slower decay is given by the largest value in the pseudospecturm of the transfer matrix, however we show that the decay rate can be an arbitrary value between the second largest eigenvalue and the largest value in the pseudospectrum. This arbitrary decay can be observed for example in the propagation of OTOC in periodic boundary conditions BW circuits. To explore this phenomenon, we study a 1D biased random walk coupled to two reservoirs at the edges, and prove that this simple system also exhibits phantom eigenvalues.

Observables O(t) which are propagated using a transfer matrix tend to decay towards their asymptotic value O(∞) exponentially In recent work [26,[31][32][33][34] it was shown that purity and out-of-time-ordered correlations (OTOC) [35][36][37] in random circuits can be propagated using a transfer matrix approach.It was noted [25,26] that these quantities exhibit a two stage exponential relaxation towards their asymptotic values, namely where λ 2 is the second largest eigenvalues of the transfer matrix (the largest is equal to 1 and corresponds to the stationary state).Note that the initial decay λ ph persists until extensive times in the system size, meaning that the effective decay λ ph ̸ = λ 2 will be present for all times in the thermodynamic limit.This behavior was dubbed phantom eigenvalue in previous papers [25].
Typically, the decay in processes involving a non-Hermitian transfer matrix is determined by the largest value λ ps in the pseudospectrum [27,28].In this paper we instead observe that the decay can be an arbitrary value between λ 2 and λ ps .This happens for example in the OTOC relaxation in brickwall (BW) periodic boundary conditions (PBC) random quantum circuit.We work out the details behind this arbitrary relaxation and find that, surprisingly, the behaviour is not determined exclusively by the properties of the transfer matrix, but depends on the specifics of initial vectors used in the iteration.
Such a decay was observed by examining the behavior of OTOC in a chain of 2n qudits of dimension q where X i and Y j are two local Pauli-like matrices located at qudit i and j, respectively.Note that O(t) depends on both i and j, but we will leave the dependence implicit.
The matrix X i is propagated in time with the BW random quantum circuit U of depth t (2t rows of a BW geometry), namely X i (t) = U † X i U .The random matrix U is composed of two sites independent Haar unitaries.Because the unitaries in the circuit are random, O(t) is independent on the choice of the matrices X and Y .First, we will look into PBC random quantum circuits.Specifically, we will derive a transfer matrix propagation approach to study the OTOC.Consequently, we observe an intriguing behavior in its decay towards its asymptotic value.Remarkably, the decay in the thermodynamic limit (TDL) is neither given by the second largest eigenvalue of a transfer matrix, nor by the pseudospectrum, but See Fig. 1 for a graphic representation.To understand better this phenomenon, we will delve in the simpler case of OTOC in BW open boundary condition (OBC) circuits .Surprisingly, we will find that the transfer matrix used to describe OTOC dynamics can be also used to describe other processes, such as a biased 1D random walk with dissipation at the edges, see App.C 4. In App.C 4, we leverage the equivalence between OTOC dynamics and a biased random walk to obtain a closed form solution for O(t) when j = 1.Focusing on general iterations of this transfer matrix, we will explain why λ ph can be a value between λ 2 and λ ps .
Periodic boundary conditions -The average time evolution of OTOC in BW PBC circuits is equivalent to a partition function of an Ising-like model [38].Shortly, the partition function is obtained by summing certain domains of up-spins that evolve on a 1D Ising chain of n sites [39], one site for each random gate, see App.A for a brief overview.In App.B 1 we derive a Markovian evolution of OTOC by composing a transfer matrix that propagates all domains in the previously mentioned 1D Ising chain [40].To calculate the average OTOC at time t, we iterate the transfer matrix A on an initial vector |v⟩ t times, t ∈ N. The domains relevant for the partition function from [38] are extracted by projecting A t |v⟩ on the vector ⟨p| Note that ⟨p| depends on the qudits j.
Once we obtain a transfer matrix approach for computing OTOC, we can compare the actual exponential decay with the second largest eigenvalue λ 2 of A. Fig. 2 illustrates that the actual decay λ ph of OTOC between qubits n and j = 1, q = 2, lies between λ 2 = (4/5) 4 (computed in App.B 2) and λ ps = 1, namely λ ph = (4/5) 2 [41].Exact results were computed also for general q, namely λ 2 = 16q 4 /(1 + q 2 ) 4 (see App. B 2), λ ps = 1 and λ ph = λ 2 ph .This intriguing behavior is retained for all possible choices of q and j.The underlying reason for this behavior can be attributed to our specific selections of ⟨p| and |v⟩.For general or randomly chosen initial vectors, the expression ⟨p|A t |v⟩ decays as λ t ps until t ∼ n, see Fig. 2b for an example.shows the iteration of A on two different sets of initial vectors: blue symbols denote the choice of vectors that give OTOC dynamics (black dashed line from (a)), black symbols denote the choice of the physical ⟨p|, but p n(n−1)+1 = 0, and random v k , k v k = 1, that decays as λ t ps .The green dotted line denotes the decay given by λ2, the red dotted line denotes λ ph and the orange one the decay for λps.(b) also compares data for n = 80 (light symbols) with data for n = 40 (dark symbols), which shows that the initial decays persist until t ∼ n.
Open boundary conditions -Let us continue by examining the propagation of OTOC for BW OBC random quantum circuits.Similar to the case of PBC, OTOC can be obtained by evolving all domains on a 1D Ising chain of n sites and summing only the relevant domains that are described in [26,38] or in App. A. To analyze the correlation between the first and the j-th qudit, we focus just on the right edge of our domain since the left edge remains fixed at the left boundary of the system.Consequently, at every time there are only n possible domains, one for each domain width, in contrast to the ≈ n 2 domains in the PBC scenario.The initial vector |v⟩, v k = q 4 q 4 −1 δ 1,k , is propagated with the transfer matrix (1+q 2 ) 2 .See App.C 1 for more details.Note that the transfer matrix A exhibits a distinct structure.Namely, the first n − 1 components are propagated using a tridiagonal Toeplitz matrix T [42].If we are interested in the OTOC between the first qudit and the j-th one, the summation over relevant domains is represented with the inner product with the vector ⟨p|, p k = 1 if k ≥ j and 0 otherwise, leading to O(t) = ⟨p|A t |v⟩.
In the physical systems that we analyze, the vector |v⟩ is localized at the leftmost position.Thus, for times t < n, the iteration with T is equivalent to the iteration with A. Therefore, when we want to analytically study the iteration until extensive times t ∼ n (the interesting domain) we can focus on the properties of the iteration O(t) = ⟨p|T t |v⟩.However, note that O(∞) is still defined as lim t→∞ ⟨p|A t |v⟩.
Before moving to the properties for general δ,τ and σ, we will first examine the classical limit q → ∞.In this case, only σ ̸ = 0 and T corresponds to a transposed Jordan kernel.We will generalize this transposed Jordan kernel to a transposed Jordan block with diagonal elements equal to δ (which is zero in for q → ∞) and lower diagonal elements σ.Note that the n times degenerate eigenvalue of T is δ and the largest value in the pseudospectrum is δ + σ.The physical case q → ∞ is recovered by setting δ = 0 and σ = 1.When |v⟩ is leftlocalized, the iteration for times t < n can be expressed as follows The second row is obtained by replacing the t-th power of the Jordan block with the t-th power of a diagonal matrix of entries δ plus the nilpotent matrix with lower diagonal elements equal to σ.The vectors e j with k-th component δ j,k are generalized eigenvectors of the nilpotent matrix.The quantity C(r) from the last row is the convolution between the initial vectors ⟨p| and |v⟩ .If C(r) is independent of r, for example for p k = 1 and v k = δ k,1 , then the sum over r can be evaluated, and O(t) ∝ (δ + σ) t for t < n.It is worth noting that δ + σ corresponds to the largest value in the pseudospectrum of T .Interestingly, In the TDL, if we wish to normalize ⟨p|, we must choose µ > 1.In all the examples O(∞) = 0, consequently, the decay of O(t) to its asymptotic value becomes an arbitrary number between δ, the largest eigenvalue, and δ + σ, the pseudospectrum of T .This brief derivation illustrates that we cannot determine the decay of O(t) solely by examining the properties of the transfer matrix.In the case C(r) = µ −r , the decay depends on the specific form of the initial vectors.Now, let us move back to the case of general q.The largest value λ ps = 1 in the pseudospectrum of T is larger than its largest eigenvalue λ 2 = δ + 2 √ στ , see App.C 2 for the properties of T .We would expect that quantities like O(t) − O(∞) = ⟨p|T t |v⟩ decay as λ t ps .However, it is known that OTOC decays as λ t 2 [26].It turns out that the initial vectors for OTOC are special, see App.C 4 for an argument made by considering OTOC dynamics as a biased random walk or App.C 3 for an exact computation of O(t).When considering general (or random) ⟨p| and |v⟩, the decay of O(t) would indeed be given by the largest value in the pseudospectrum, see Fig. 3.
ps .The green dotted line denotes λ2, the red dotted line the value λps.The plot also compares data for n = 80 (small symbols) with data for n = 40 (big symbols), which shows that the initial decays persist until t ∼ n.
Arbitrary decay -We have observed that the decay λ t 2 is achieved through the use of special vectors.Is it possible to get an arbitrary decay rate with an appropriate choice of ⟨p| and |v⟩, as we saw for Jordan blocks?Fig. 4a) shows how Although the largest value in the pseudospectrum is 1, this quantity decays as ≈ 0.85 t , which is between λ 2 and λ ps .
The arbitrary decay can be intuitively understood in terms of the pseudospectrum.Instead of propagating O(t) = ⟨p|T t |v⟩ for p k = µ −k and v k = δ k,1 , we instead express it like ⟨p|D −1 T t D|ṽ⟩, where pk = 1, so D The quantity ⟨p|A t |v⟩ decays towards its asymptotic values as λ t ph ≈ 0.85 t , and λ2 < λ ph < λps.The red line denotes λ ph , the orange line the value λps > λ ph and the green line the asymptotic decay of finite system sizes λ2.The plot shows data for two system sizes, namely n = 80 for light symbols and n = 40 for dark symbols, which shows that λ ph persists until times extensive in the system size.(b) Cartoon representation of the pseudospectrum of the matrix from Eq. 8.For different exponential localization µ of the initial vector we get different pseudospectra.As we increase µ, we get from the pseudospectrum of A for µ = 1 to the spectrum of A for µ = 4.The solid black line on the real axis corresponds to the spectrum of the matrix.is a diagonal matrix with diagonal elements D k,k = µ k .The action of D −1 retains the right vector's form but rescales it.Recall that non-unitary similarity transformations can alter the matrix's pseudospectrum [43].In fact, the largest value in the pseudospectrum of is λ ps = δ + σ µ +τ µ which coincides with the decay of O(t) from Fig. 4a for µ = 1.35.A schematic representation of the pseudospectrum of D −1 T D for different values µ is shown in Fig. 4b.It turns out that λ ph can only be between λ 2 and λ ps .To see this one must solve exactly all the sums in the expression O(t) = ⟨p|T t |v⟩, see App.C 3 for the computation of λ ph .We conclude that the decay of O(t) towards its asymptotic value is not determined solely by the properties of the transfer matrix, but it is highly dependent on the initial vectors used in the iteration.
Note that with Eq. 8, we could get a non-Hermitian matrix from an initially Hermitian matrix (σ = τ ).Does this mean that a decay slower than λ t 2 is possible in Hermitian systems?It turns out that achieving λ ph > λ 2 is only possible if we select p k = µ −k with µ < 1.In the TDL, this implies that the vector ⟨p| is not normalizable, causing O(t) to be exactly equal to zero already at t = 0, see App.D for the proof.Since the OTOC is exactly zero in the TDL, λ ph > λ 2 in Hermitian systems becomes a finite size effect.
Discussion -We studied the decay of observables O(t) = ⟨p|A t |v⟩ propagated with a non-Hermitian transfer matrix A. Such systems can be found when studying purity or OTOC propagation in random quantum circuits or other Markovian systems.The peculiarity of these systems is that they can exhibit phantom eigenvalues, that is, the convergence of O(t) towards its asymptotic value is not determined by the second largest eigenvalue λ 2 of A but rather by the largest value λ ph in the pseudospectrum of A. In this paper we have shown that this is not always the case.The exponential rate at which O(t) relaxes can be an arbitrary value between λ 2 and λ ph .Such "arbitrary" decay happens in physical systems, for example in OTOC relaxation in PBC BW random quantum circuits.To compute the actual decay rate of O(t) one must not just look at properties of the transfer matrix A, but rather at the whole system.Namely, we found that when ⟨p| was exponentially localized, this changed the convergence rate of O(t) to an arbitrary value, which depends on the localization length of ⟨p|.
Although the decay is not solely determined by the transfer matrix properties, the pseudospectrum still plays a crucial role in the computation of the phantom decay λ ph .However, it is not the pseudospectrum of the transfer matrix A which we have to explore, but rather the pseudospectrum of the transformed matrix D −1 AD, where D is chosen so that ⟨p|D is exponentially localized.Currently, it is understood that the pseudospectrum of A (or D −1 AD) determines the exponential relaxation of observables for general initial vectors.However, for special choices of initial vectors the observable can decay as λ 2 .Ultimately, it would be useful to develop a general technique to determine whether an observable decays as λ ph > λ 2 solely from looking at properties of A and initial vectors.This will be ground for future studies.Each domain has its own weight, which is determined by the number of + spins in the top edge of the grid and by length of the domain walls.For each + in the top edge we multiply the weight of the domain by q 2 .Each opposite neighboring spins (horizontally) contribute with q/(q 2 + 1).For the example in Fig. 5, the weight of the domain is (q 2 ) 3 • ( q q 2 +1 ) 16 , because the domain has 3 qudits in the top edge and 16 differently oriented neighboring spins.In an infinite system, the number of differently oriented neighboring spins is equal to 4t, but in a finite system it could be smaller if the domain hits the boundary of the system or if at some point in time the + spin span over all sites.
OTOC in BW circuits can be expressed as a weighted sum of domains on a 2D grid that start at qudit i and contain the qudit j in the top edge.Instead of looking at the partition function of a 2D grid, as the authors of [38] did, in the following Sections we will compute OTOC by evolving all possible domains of the 1D Ising chain, obtained by considering the vertical axis of the 2D grid as time.To obtain O(t) we then just sum the weight of relevant domains at time t.These weights can be now though as probabilities in this Markovian propagation.These Markovian propagation of domain wall to compute OTOC is preferable when dealing with finite systems, because it is easy to take in account the boundaries of the system.Computing partition functions of the 2D grid of Ising spins works well for infinite systems, but gives complicated solutions expressed with recursion for finite systems [26].Using the Markovian propagation of domain weights described in this paper, in App.C 4 we will see that we can obtain a simpler close-form solution of OTOC in OBC circuits.Here we derive a Markovian propagation of all domains for PBC circuit.As described in App.A, OTOC is then recovered by summing only certain domains that contain the site j at time t.
To begin, we shall encode all domains in the vector |v⟩ such that the domain beginning at the i z -th spin and having width w ∈ {1, . . ., n − 1} is written in the n • i z + w component of |v⟩.Note that n counts the number of spins, which is half of the number of qudits.The only domain of width n will be put in the last position of |v⟩, so |v⟩ will be a vector of n • (n − 1) + 1 components.
To propagate |v⟩ in time, we can construct a transfer matrix that propagates all domains of the 1D chain of Ising spins (see App. A. Following App.A we get where T is a block circulant matrix where δ = 2q 2 (1+q 2 ) 2 , τ = 1 (1+q 2 ) 2 and σ = q 4 (1+q 2 ) 2 .A row of T is thus composed of n blocks with size (n − 1) × (n − 1).As last, the last row b of A is a vector of n(n − 1) components.The non-zero components of b are b (n−1)i−1 = σ 2 and b (n−1)i = δσ + q 2 σ for i ∈ {1, . . ., n}.
We shall always begin with the initial vector |v⟩ located on the last spin.This choice is irrelevant because of the periodic boundary conditions.However, we must be careful when we choose ⟨p|.Depending on the position of the qudit j, the non-zero components are p (n−1)(i−1)+k = 1 for i ∈ {1, . . ., n} and k ∈ {j − i + 1 mod n − 1, . . ., n − 1} and p n(n−1)+1 = 1.The vector |v⟩ is obtained by propagating the domain on the last site for a half-time step.The only nonzero components of |v⟩ are v 1 = q 2 , v (n−1) 2 +1 = q 2 and v (n−1) 2 = q 4 .

Diagonalization
To diagonalize the transfer matrix A, one can begin by diagonalizing the block circulant matrix T .If λ k and |v k ⟩ are an eigenpair of T , then λ k is also an eigenvalue of A whose eigenvector is obtained from |v k ⟩ by adding one last component ⟨b|v k ⟩ λ k −1 .There is also an eigenvector of A which is not derived from the eigenvectors of T , namely (0, . . ., 0, 1) with eigenvalue 1.To diagonalize A thus one needs first to diagonalize T .
The block circulant matrix can be diagonalized by applying a block Fourier transform F † T F , where F is a matrix of n × n blocks of size n − 1 × n − 1.The block at the i-th row and j-th column of F is a diagonal matrix with constant elements exp(2πijk/n)/ √ n.After applying this similarity transformation, we will end up with a block diagonal matrix with n blocks, where the k-th block is where d −2 = q 4 exp(2πik/n), d −1 = 2q 2 (1 + exp(2πik/n)), d 0 = 2 cos(2πk/n), d 1 = 2/q(1 + exp(−2πik/n)) and d 2 = 1/q 4 exp(−2πik/n), k ∈ {1, . . ., n}.Following [44] we can rewrite the matrix above as a product of two commuting matrices where where The tridiagonal matrices Ãk and Bk can be diagonalized simultaneously.The eigenvalues of T are, after simplifications, where k ∈ {1, . . ., n} is the momentum from the Fourier transformation and j = {1, . . ., n − 1} runs through the components in the k-th Fourier mode.The corresponding left and right eigenvectors (l /n e πikm/n sin(jmπ/n), (B11) for l ∈ {1, . . ., n} and m ∈ {1, . . ., n − 1}.

Pseudospectrum
In this Subsection, we shall compute the pseudospectrum of T .In Fig. 6 we plot the ϵ-pseudospectrum for system sizes 80,120 and 180 and ϵ = 10 −5 .The largest value in the pseudospectrum approaches 1 as we increase the system size.We conjecture that in the thermodynamic limit the pseudospectrum T is the union over all values k of the product of the pseudospectra of Ãk and Bk from Eqs.B9.This coincides with c 1 e iϕ + a 1 + b 1 e −iϕ c 2 e iϕ + a 2 + b 2 e −iϕ , for ϕ ∈ [0, 2π] and k/n ∈ [0, 1].The conjectured region is shown in Fig. 6 in black and it seems to match with the plots from Fig. 6 for n → ∞.As we increase n, the ϵ-pseudospectrum starts filling the black region, which is our conjecture for the pseudospectrum of T .The black region is obtained as the union of all curves from the conjecture for every possible value of k/n.
the system, so the only part of the domain that can move is its right edge.This means that we have to keep track of just n different domain widths.We begin the derivation of the transfer matrix by encoding all the domains by increasing width in a n dimensional vector |v⟩.The transfer matrix propagating the domains is 2 ) 2 and σ = q 4 (1+q 2 ) 2 .The initial vector must contain the domain on just the first spin, so |v⟩ = ( q 4 q 4 −1 , 0, . . ., 0).The vector used to extract O(t) depends on the qudit j and it is p k = 1 for k ≥ j/2 and 0 otherwise.

Transfer matrix properties
Similarly as for the matrix A for PBC, for OBC we can also compute the spectrum of A by first computing the spectrum of the Toeplitz tridiagonal matrix T .Namely, if λ k and |v k ⟩ are an eigenpair of T , then λ k is also an eigenvalue of A whose eigenvector is obtained from |v k ⟩ by adding one last component There is also an eigenvector of A which is not derived from the eigenvectors of T , namely (0, . . ., 0, 1) with eigenvalue 1.The eigenvalues of T are [42] where k ∈ {1, . . ., n − 1}.The right eigenvectors |r k ⟩ and left eigenvectors The pseudospectrum of T is [43] δ Note that the spectrum of T lies on the real values, meanwhile the pseudospectrum is an ellipse in the complex plane with largest value δ + σ + τ , see Fig. 4b.

Evaluation of O(t) through spectral decomposition of T
In this Appendix we will evaluate the expression for OTOC As we shall see, the OTOC behaves as ∼ λ t ps = 1, but if we subtract the steady state of the whole transfer matrix A the leading term cancel out making O(t) decay as λ t 2 , so we get O(t) − O(∞) ≍ λ t 2 as expected.Even though O(t) does not exhibit phantom behaviour, the evaluation of O(t) is interesting because the phantom appears in the solution, but it gets cancelled with the steady state that comes from the eigenvalue λ 1 = 1.This shows that the physical case of OTOC is special, because the choice of vectors ⟨p| and |v⟩ makes the phantom cancel exactly with O(∞).
We begin by writing Eq.C6 with the help of the spectral decomposition of T where To simplify the expression we shall first evaluate the sum over the index k.We get Plugging Eq.C8 into Eq.C7 and rearranging some terms we get where we labeled the largest value in the pseudospectrum δ + σ + τ with λ ps .To simplify Eq.C9 we will assume τ < σ and n ≫ 1, so that we can neglect the term τ σ n 2 .We are interested in the time dependence of O(t), so we can forget about the factors outside the sum, because they do not contribute to the time dependence of O(t).After these simplifications we get We can replace the fraction , because λ ps > λ h .Omitting constant terms we get To prove that O(t) decays as λ t ps we will show that the terms in the sum over the index k are zero for k ≈ n, meaning that the sum over k runs from k ≈ n to ∞ and it is independent on t.This in turn implies that the only time dependence in the expression Eq.C11 is λ t ps .We continue by summing over the index h where is the Dirichlet kernel [45].The Dirichlet kernel is composed of a term (−1) x+1 , which cancels with the other functions D n (x) and an infinite series of Kronecker delta functions δ x/2,πp , p ∈ Z.In our case x = r−2s+n+2 n in the first Dirichlet kernel.The only possible values of p are 0 and 1, because for other values there are no s and r to satisfy x/2 = πp.Moreover, for p = 0 and p = 1 not all values of r give non-zero contribution.For example x = r − 2s + n + 3 and p = 1 gives s = (r − n + 1)/2.The variable s runs from 0 to n, so (r − n + 1)/2 ≥ 0 and (r − n + 1)/2 ≤ r.The former bound gives r ≥ n − 1, the latter r ≥ −n.Computing these bounds for all three Dirichlet kernels and for both p = 0 and p = 1 we see that r ≳ n for each term.This in turn implies that the sum over r runs from ≈ n to k, so k ≳ n to have non-zero contributions.For k < n all terms in the sum over k will be zero, so we can substitute ∞ k=t with ∞ k≈n , which makes the only time dependent part of the leading term in O(t) λ t ps .
We saw that O(t) behaves as λ t ps , but OTOC decay to O(∞) as λ t 2 .This comes from the fact that the leading term in Eq.C9 (the one with (−1) h ) exactly sums to 1 for t < n (not shown)).This means that by subtracting 1 from O(t) we cancel out the leading term of O(t).The subleading term (the neglected term in Eq.C9) decays as λ t 2 , hence OTOC decays as given by the largest eigenvalue of A.
The exact evaluation of ⟨p|T t |v⟩ above can be also used for the case p k = µ −k .The choice of ⟨p| reflects in the evaluation of ⟨p|r h ⟩, where |r h ⟩ is an eigenvector of T .We get The Eq. C15 is exactly Eq.C8 if we substitute τ with τ µ 2 .We can repeat all calculations from before until we get with λ(µ) = δ + σ/µ + τ µ.Eq.C16 is analogous to Eq. C10 if we replace λ(µ) with λ ps .We can repeat all calculations below Eq.C10 to show that λ(µ) is indeed the true decay of O(t).The value of λ(µ) is greater than λ 2 , because otherwise we cannot repeat the steps in Eq.C11.Moreover, λ(µ) ≤ λ ps = 1 otherwise ⟨p| is not normalizable.We conclude that λ 2 ≤ λ(µ) ≤ λ ps .

Biased random walk
In this Section we prove that OTOC dynamics in a OBC BW circuit is equivalent to a biased 1D random walk coupled to reservoirs at the edges.Using this equivalence, we will be able to compute the decay of OTOC to , where λ 2 is the second largest eigenvalue of the transfer matrix that propagates OTOC (see App. C 2).
The transfer matrix from Eq. 6, used to propagate OTOC, can be also used to describe a 1D biased random walk.To get a proper Markov chain transfer matrix, the elements in each column should sum to 1 to conserve probabilities.To achieve this and keep the tridiagonal transfer matrix to propagate probabilities one can make the random walk dissipate on the left and right boundary, as shown in Fig. 7. Doing so, the Markov chain transfer matrix is where τ + σ + δ = 1, and where the first and last site correspond to the left and right bath, respectively.Note that we could propagate OTOC using this transfer matrix, with p k = 1 − δ k,1 and v k = δ 2,k .Because p 1 = 0, the OTOC can be equivalently propagated by forgetting about the first column and row in Eq.C17.Using the random walk interpretation, the OTOC O(t) = ⟨p|A t |v⟩ is equivalent to the probability of not being in the left bath by beginning at the second position (left edge of the bulk).
Let us label the probability of being at site j in the bulk with r j+1 (so the leftmost site has probability r 2 ), and the probability of being in the left bath with r 1 , then OTOC is O(t) = 1 − r 1 .To prove that OTOC decays to its asymptotic value O(∞) = 1 as λ t 2 , we will compute the probability r 1 .Because we can enter the left bath just once, the probability r 1 of being in the left bath at time t can be found by summing the probabilities of entering the left bath for all times T + 1 ≤ t.The probability of entering the left bath at time T + 1 is obtained as τ times the probability of being at the leftmost site in the bulk 1 at time T without ever being in the left reservoir before.The number of paths of length T beginning at 1 and ending at 1 without going to the left reservoir corresponds to the number of Dyck words [46] with σ and τ .The number of these Dyck words can be expressed with the Catalan number C T /2 [47], which is T T /2 /(T /2 + 1) if T is even and 0 otherwise.The corresponding probability for such a Dyck word is (τ σ) T /2 .To get the probability of being at 1 at time T without ever being in the left bath, we should also include the moves where we stay at the same place; this moves have probability δ and can be places anywhere in the Dyck word.To sum up, the probability of being at site 1 at time T is obtained as the sum over all possible combinations of moves to left/right and staying at the same position.At the end we get OTOC is obtained as 1 − r 1 .For simplicity, we will now compute O(t) for qubits, q = 2.We get where 2 F 1 (a, b, c; z) is the hypergeometric function.With the help of the biased random walk picture, in Eq.C19 we obtained a closed form solution of OTOC in OBC BW circuits with Haar random 2-site gates.A closed form solution was obtained also in [26], but it is much more complex, because it is expressed with recursion.In [38] there is a simple result for OTOC, but it is for infinite systems, whereas our solution holds for any system size n.
Because of the simple form of Eq.C19, we can compute the rate at which OTOC decays to O(∞).Assuming OTOC decays to the asymptotic value exponentially, we can get the exponent as λ eff (t) = O(t+1)− .

(C20)
The term 3/2+t 3+t behaves as 1 − 3 2t + O(1/t 2 ), and the ratio of the hypergeometric functions decays to 1 faster than exponentially (Fig. 8), so we conclude that O(t) decays to 1 as λ eff = 16/25 = λ 2 .We computed the time evolution of a random walk, where we start at the leftmost position in the bulk and we are interested in the probability of staying in the bulk.This random walk coincide with the OTOC evolution in OBC random quantum circuits.When computing the relaxation of OTOC to their asymptotic value O(∞), we subtract the leading term from O(t), which results in the relaxation given by λ 2 .For different initial conditions subtracting O(∞) from O(t) does not cancel exactly the leading term in O(t), meaning that O(t) initially does not relax, but stays constant, as we expect by looking at the largest value in the pseudospectrum of the transfer matrix, λ ps = 1.OTOC in this sense can be considered a special case of initial conditions.

Appendix D: Hermitian transfer matrix
In this Appendix, we shall prove that the phantom decay λ ph > λ 2 of quantities O(t) = ⟨p|T t |v⟩ to their asymptotic value O(∞) is a finite size effect if T is a symmetric tridiagonal Toeplitz matrix with upper and lower diagonal elements β.Increasing the system size, O(t) approaches 0 for every non-zero time as seen if Fig. 9.In the TDL, O(t) will be exactly equal to zero for every non-zero time.The quantity ⟨p|T t |v⟩ decays towards its asymptotic values as λ t ph ≈ 0.85 t , and λ2 ≈ 0.72.The red line denotes λ ph ≈ 0.85 and the green line the asymptotic decay of finite system sizes λ2.Even though the initial decay λ ph is larger than λ2, the plot for different system sizes clearly shows that O(t = 0) approaches 0 by increasing the system size n.
Let us begin with the initial vectors p k ∝ µ −k and v k = δ k,1 .The expression O(t) can be computed using p k = 1 if we appropriately transform T with a similarity transformation, similarly as we did in the main text in Eq. 8.In this case we get an effective iteration with a non-Hermitian matrix, where a decay λ ph > λ 2 should not be surprising.We have where D is a diagonal matrix with diagonal elements equal to D k,k = µ k .The choice µ > 1 is the only possible choice if we wish to normalize ⟨p|.In this case the upper diagonal of D −1 T D is larger than its lower diagonal, resulting in a transfer matrix of the same form as the one used to propagate OTOC.Even the initial vectors p k = 1 and v k = δ k,1 are the same as those used for OTOC, so we know that O(t) will decay to its asymptotic value O(∞) = 1 as λ 2 .
Another possible choice for ⟨p| is when µ < 1.In this case λ ph > λ 2 , however if we wish to normalize ⟨p|, then O(t = 1) will scale as ⟨p|v⟩ = The normalization and right-localization of ⟨p| implies that O(t) is exactly zero for every non-zero time in the TDL.
To conclude, we explored the case of Hermitian T and saw that the rate at which O(t) decays to O(∞) is λ ph > λ 2 only ⟨p| localized at the right edge of the system, i.e. µ < 1.In this case, however, O(t) decays to 0 for every non-zero time as we increase the system size n, which makes the phantom eigenvalue a finite size effect in Hermitian systems.In contrast, when T is non-Hermitian, λ ph > λ 2 until extensive times also when O(t) is not zero in the TDL, which makes λ ph the only true decay in the TDL.

Figure 2 .
Figure 2. OTOC time dependence for a BW PBC circuit, q = 2. (a) shows the propagation of O(t) − O(∞) for different final positions j for i = n, n = 50.The gray curve denotes the critical time tc when λ eff = 0.55.(b) shows the iteration of A on two different sets of initial vectors: blue symbols denote the choice of vectors that give OTOC dynamics (black dashed line from (a)), black symbols denote the choice of the physical ⟨p|, but p n(n−1)+1 = 0, and random v k , k v k = 1, that decays as λ tps .The green dotted line denotes the decay given by λ2, the red dotted line denotes λ ph and the orange one the decay for λps.(b) also compares data for n = 80 (light symbols) with data for n = 40 (dark symbols), which shows that the initial decays persist until t ∼ n.

Figure 3 .
Figure3.Iteration of A (OBC, q = 2) on two different sets of initial vectors: blue symbols denote the choice of vectors that give OTOC dynamics, black symbols denote the choicep k = 1 − δ n,k and random v k , k v k = 1, that decays as λ tps .The green dotted line denotes λ2, the red dotted line the value λps.The plot also compares data for n = 80 (small symbols) with data for n = 40 (big symbols), which shows that the initial decays persist until t ∼ n.

Figure 4 .
Figure 4. (a) Iteration of A (OBC, q = 2) for p k = µ −k and v k = δ k,1 , µ = 1.35.The quantity ⟨p|A t |v⟩ decays towards its asymptotic values as λ tph ≈ 0.85 t , and λ2 < λ ph < λps.The red line denotes λ ph , the orange line the value λps > λ ph and the green line the asymptotic decay of finite system sizes λ2.The plot shows data for two system sizes, namely n = 80 for light symbols and n = 40 for dark symbols, which shows that λ ph persists until times extensive in the system size.(b) Cartoon representation of the pseudospectrum of the matrix from Eq. 8.For different exponential localization µ of the initial vector we get different pseudospectra.As we increase µ, we get from the pseudospectrum of A for µ = 1 to the spectrum of A for µ = 4.The solid black line on the real axis corresponds to the spectrum of the matrix.

4 Figure 5 .
Figure 5. Schematic representation of a domain that must be summed in the partition function for Oi,j(t = 4).The domain starts around the site of the initial Pauli-like matrix i and contains the site j at time t = 4.

Appendix B: PBC case 1 .
Transfer matrix derivation vector b describes the transition probabilities to the steady state, i.e. a 1D chain of Ising spins with all +.C dictates how the domain width changes for fixed i z .U and D describe how to obtain domains with i z − 1 or i z + 1 from i z , respectively.The matrices D, U and C are tridiagonal matrices because of locality in the random walk of domain width and have dimension (n − 1) × (n − 1).T has a block circulant form because of the locality in the changes in i z .These matrices are . . .0 δσ τ σ 0 0 . . .0 σ 2 δσ τ σ 0 . . .0 . . . . . . . . . . . . . . .0 . . .σ 2 δσ τ σ 0 0 . . .0 σ 2 δσ τ σ

Figure 6 .
Figure 6.The colored dots represent the ϵ-pseudospectrum for n = 80 (green), n = 120 (orange) and n = 180 (blue) and ϵ = 10 −5 .As we increase n, the ϵ-pseudospectrum starts filling the black region, which is our conjecture for the pseudospectrum of T .The black region is obtained as the union of all curves from the conjecture for every possible value of k/n.

Figure 7 .
Figure 7. Cartoon picture of a biased random walk on 4 sites with dissipation on the boundary.Above the sites (boxes) are shown possible moves (arrows) with the corresponding probability.One the random walk enters the reservoirs it cannot return back to the bulk of the chain.

25 Figure 8 .
Figure 8. Ratio of the two hypergeometric functions from Eq. C20.The ratio decays fast to 1.

Figure 9 .
Figure 9. Iteration of the tridiagonal matrix T with δ = 8/25, τ = σ = 5/25 for p k ∝ µ −k and v k = δ k,1 , µ = 0.4.The vector ⟨p| is normalized such that k p k = 1.The quantity ⟨p|T t |v⟩ decays towards its asymptotic values as λ tph ≈ 0.85 t , and λ2 ≈ 0.72.The red line denotes λ ph ≈ 0.85 and the green line the asymptotic decay of finite system sizes λ2.Even though the initial decay λ ph is larger than λ2, the plot for different system sizes clearly shows that O(t = 0) approaches 0 by increasing the system size n.