Structural Stability Hypothesis of Dual Unitary Quantum Chaos

Having spectral correlations that, over small enough energy scales, are described by random matrix theory is regarded as the most general defining feature of quantum chaotic systems as it applies in the many-body setting and away from any semiclassical limit. Although this property is extremely difficult to prove analytically for generic many-body systems, a rigorous proof has been achieved for dual-unitary circuits -- a special class of local quantum circuits that remain unitary upon swapping space and time. Here we consider the fate of this property when moving from dual-unitary to generic quantum circuits focussing on the \emph{spectral form factor}, i.e., the Fourier transform of the two-point correlation. We begin with a numerical survey that, in agreement with previous studies, suggests that there exists a finite region in parameter space where dual-unitary physics is stable and spectral correlations are still described by random matrix theory, although up to a maximal quasienergy scale. To explain these findings, we develop a perturbative expansion: it recovers the random matrix theory predictions, provided the terms occurring in perturbation theory obey a relatively simple set of assumptions. We then provide numerical evidence and a heuristic analytical argument supporting these assumptions.


I. INTRODUCTION
Exact solutions of interacting many body systems are not just monuments to human ingenuity, they are also key instruments in both statistical mechanics and many-body dynamics.In fact, it is often implicitly assumed that each (dynamical) universality class in statistical physics should be endowed with at least one exactly solvable model through which one obtains deeper understanding of the whole universality class.
These are the so-called dual unitary circuits, which are expressed in the form of locally interacting quantum systems in discrete space time, and whose defining feature is that they generate a unitary evolution not only in time but also in space.This fundamental property is most clearly expressed in terms of the so-called space-time duality [36], a space-time swap symmetry of the tensor net-work diagram representing the physical observable of interest.Among other useful features, dual-unitary circuits allow for exact analytical computation of dynamical correlation functions of local operators [18], as well as longrange two-point spectral correlations as expressed in the form of the spectral form factor (SFF) [17,19].
Given this novel class of exactly solvable chaotic systems a fundamental question concerns the stability or robustness of their dynamical features.One may trace the basic motivation for such a question to the notion of structural stability of hyperbolic flows in classical chaotic dynamical systems [37,38].Contrary to integrable systems, which are structurally unstable and where perturbative expansions around them generically diverge (cf.Kolmogorov-Arnold-Moser theory in classical dynamical systems, or divergent Feynman diagram expansions of quantum field theories around their free/non-coupled limits), chaotic systems are expected to be robust against typical perturbations.In this work we formulate the hypothesis of structural stability of Floquet dual-unitary circuits (therefore focusing on discrete time-dynamics with explicit time translation symmetry), and streamline a simple strategy for its verification.Specifically, we conjecture that chaotic/ergodic dual-unitary circuits remain chaotic/ergodic under small, time-translation invariant perturbations.If we assume that chaos is equivalent to a linear-in-t growth of the SFF in the thermodynamic limit, this hypothesis can be reformulated in terms of the socalled spectral Lyapunov exponents (SLE) [39].In that language, the hypothesis states that the leading SLE in each of the t time-translation symmetry sectors decays to 1 exponentially fast in time t.We note that a somewhat simpler approach to structural stability hypothesis has been addressed by two of us and P. Kos in Ref. [40], which studied the robustness of local dynamical correlators.Importantly, however, the decay of local correlators is not sufficient for establishing quantum chaos and ergodicity, hence our current study of SFF, a global (nonlocal) dynamical observable, provides a more stringent characterisation.We also stress that our study goes substantially beyond the one presented in Ref. [41], where the stability of the SFF was investigated for a specific dual unitary circuit at first order in perturbation theory.
Clearly, this question is very challenging, and extremely difficult to address in full mathematical rigour.The purpose of our paper is to establish the minimal set of assumptions -which turn out to be two -needed for a rigorous proof of the stability of the chaotic SFF.We do this by formulating a perturbation theory for the SLE, and identifying sufficient conditions for the expansion to behave in a way that is consistent with ergodicity.We furthermore verify these two assumptions numerically.
We now summarise our assumptions and results in more detail.We study the SFF averaged over an ensemble of locally constructed unitaries E of an extended 1D system with L sites K(t, L) = E U∈E |trU t | 2 .The ensemble is parameterised by ϵ, where the unperturbed point ϵ = 0 corresponds to an ensemble of dual unitary Floquet evolutions.Using the space-time duality noted above, the SFF may be recast as K(t, L) = tr T L for a SFF transfer matrix T acting on 2t sites [17,19].The linear-in-t ramp in the SFF characteristic of ergodic systems (see Eq. ( 10)) is obtained if the leading t eigenvalues of T are are nearly degenerate, i.e., exponentially close to unity, λ j=1,...,t = 1 + e −O (t) .
The unperturbed dual-unitary model provably has this eigenvalue structure; we investigate the circumstances under which it remains true in eigenvalue perturbation theory in ϵ.We show that it follows from two assumptions.The first of these stipulates a non-crossing of the leading eigenvalue of T as a function of ϵ.Recalling the von Neumann-Wigner theorem [42], this corresponds to a genericity assumption on the perturbation.The second, more substantial assumption can be expressed in terms of an exponential bound on certain multi-point correlation functions involving the unperturbed SFF transfer matrix.We numerically demonstrate the validity of these assumptions in a particular family of perturbed Floquet dual unitary circuits, and corroborate those results with a heuristic analytical argument involving the spectral decomposition of the resolvent of the unperturbed SFF transfer matrix.
The rest of the paper is structured as follows.In Sec.II we recall the precise setting considered, introduce a simple minimal model that we use for the numerical tests, and briefly describe the numerical methods used in our computations.In Sec.III we present a numerical survey suggesting that the dual unitary behaviour is indeed structurally stable.In Sec.IV we present our perturbative argument.In Sec.V we discuss the validity of our second assumption, presenting numerical tests and an heuristic analytical argument.Finally, in Sec.VI we report our conclusions and discuss the outlook of our research.Some technical details and proofs are reported in the three appendices.

A. Physical System
We consider a unitary quantum circuit acting on a chain of 2L qubits (the local Hilbert space has dimension d = 2) at half-integer positions that are evolved by discrete applications of the Floquet operator Here {U x } x=0,1/2,...,L−1/2 are the local gates, i.e., unitary matrices acting on two adjacent qubits, at positions x and x + 1/2.Matrices acting at different positions are generically different and we denote by the subscript x the leftmost site where the matrix acts non-trivially.The local gates can be parameterised as with where σ = (σ (1) , σ (2) , σ (3) ) is a vector of Pauli matrices, and σ x is the corresponding local embedding in (C 2 ) ⊗2L , while (u x ⊗ v x ) is a tensor product of two one-site unitaries u, v positioned at sites x and x + 1/2 respectively.
Since the operator U is a special kind of matrix-productoperator, it can be depicted using the standard diagrammatic representation of tensor networks [43].In this rendition tensors with n indices are represented as two dimensional shapes with n protruding legs.For instance, the local gate is represented as The sum over indices is represented by connecting the corresponding legs and matrix multiplication acts bottom to top.For example, see Fig. 1 for a portray of tr[U t ] where each local gate is represented by and different shades denote in principle different matrices.
Note that, choosing J 1,x = J 2,x = π/4 the quantum circuit becomes dual-unitary [18], and also the left-toright contraction of the diagram can be thought of as the trace of a power of unitary operator [17,19].
For simplicity from now on we focus on the case where the interaction term is the same at each half step.Namely, we consider and set while the one-site gates u x and v x are position dependent, i.e., θ x , ϕ x in (3) are explicitly x−dependent.

B. Spectral Form Factor and Space Transfer Matrix
Our aim is to characterise the spectral statistics of the Floquet operator (1) for generic choices of the local gates.Namely, we want to understand the general features of the distribution of the eigenvalues of U, i.e., where the quasienergies φ j can be taken in [0, 2π).To this end we compute the spectral form factor (SFF) which measures spectral correlations over arbitrary distance and, over the last few years, has emerged as the standard spectral-correlation measure in extended systems, see, e.g., Refs.[16,17,19,[44][45][46][47][48][49][50][51][52][53][54].Here E[•] denotes an expectation value over an ensemble of similar systems, which we conveniently generate by taking the local gates v x , u x in Eq. ( 2) to be independent and identicallydistributed random matrices (equivalently one can take the angles {θ x , η x } to be i.i.d.random).The specific distribution of {v x , u x } is irrelevant for the discussion below, and we will make a concrete choice when discussing the parameterisation used in our numerical studies (cf.Sec.II D).Introducing the average allows us to filter out the system-specific details and obtain a universal result which is expected to only depend on gross properties of the system.Note that the SFF is not self averaging [55], therefore without disorder averaging one should not expect universal results (see, e.g., Refs.[56]).We expect our results to apply for any finite -no matter how small -magnitude of disorder in the thermodynamic limit.Specifically, in ergodic systems the SFF is expected to take a universal form that coincides with that observed in random matrices of the same size.The specific randommatrix ensemble to compare with depends on the antiunitary symmetries of the Floquet operator (e.g.time reversal symmetry).In the generic case of no anti-unitary symmetries, which is the one of interest here, the relevant prediction is that of Dyson's Circular Unitary Ensemble (CUE), which reads as [6,57] showing a characteristic ramp-like shape.On the other hand, whenever the system is strongly non-ergodic (e.g.integrable, or localised) the energy levels are expected to be statistically independent.This means that the SFF should reproduce the Poissonian-distribution result As discussed in Refs.[17,19] (see also Refs.[49,50,52] where this approach is applied to non-dual-unitary systems), the spectral form factor of a quantum circuit can be rewritten in terms of the trace of the L-th power of a transfer matrix acting along the space direction.The idea is to exploit the symmetry of the diagram in Fig. 1 under a 90°degree rotation (space-time swap), and the fact that the disorder is uncorrelated in space.The main steps go as follows.First we observe that, for t > 0, the SFF in Eq. ( 9) can be written as Using now the parameterisation (6, 7) we can depict the quantity inside the average on the r.h.s. as in Fig. 2 where we introduced the symbols 2. Graphical representation of K(t, L) (cf.Eq. ( 12)).
The symbols represent the local gates as per Eqs.( 13)-( 14).Random one-site gates along the same column coincide, while those on different columns are uncorrelated.
and random one-site gates of the same shade of colour are the same.As it is clear from the figure, this quantity can be equivalently represented as the trace of a matrix acting on the vertical lattice (of 2t qubits).
More precisely, we denote by Ṽ and W the many-body operators composed by vertical columns of gates acting on the time lattice, i.e., where Ṽ and W are obtained by reshuffling the indices of the local gates to propagate from left to right.Namely, the are obtained from V and W through the mapping ( Using these objects we define the transfer matrix which is highlighted in the shaded box of Fig. 2.Here (•) T denotes transposition, Π n is the one-site shift operator in a chain of n qubits, and the tensor product ⊗ r combines operators acting in the forward evolving time-sheet with those acting on the backward evolving one (cf.Eq. ( 9)): from now on we will refer to the lattices over which these operators act respectively as the forward and backward time lattice.Using the transfer matrix in Eq. ( 17) we can express the SFF as Recalling now that random gates at different positions (i.e.those along different columns in the figure) are uncorrelated, we can bring the average E [•] inside the trace to obtain Explicitly, we have where O is a non-expanding map implementing the average over the local disorder [58], i.e., Note that, although we cannot generically prove that T is diagonalisable, Eq. (19) implies that the SFF is solely determined by its spectrum and the size of its Jordan blocks.To see this we write the Jordan decomposition of T as follows where R is invertible, D is diagonal, and K is strictly upper triangular with zeros on the diagonal.The eigenvalues of D coincide with those of T and the degeneracy of a given eigenvalue λ is given by where j labels all the Jordan blocks corresponding to λ (their total number is N λ ), designated as J j,λ , while dim(A) denotes the dimension of the matrix A. Plugging the decomposition ( 22) into (19) we find where we used that products of D and at least one K are traceless (using fact K strictly upper diagonal).
A notable property of T is that it has a global Z t × Z t symmetry under independent two-site translations in the forward and backward lattices, i.e.
As a result, we can block-diagonalise it by considering a fixed double-momentum sector labelled by (ν, ν ′ ), with ν, ν ′ = 0, . . ., t − 1.Therefore, the eigenvalues of T (or D) can be labelled as where a = 0, . . ., N (ν,ν ′ ) − 1 with N (ν,ν ′ ) the size of the (ν, ν ′ ) sector.Noting that the projector onto the sector (ν, ν ′ ) can be written as Y (ν) ⊗Y (ν ′ ) , where we introduced such that , we find Ref. [19] proved that in the dual unitary limit of the models considered here and away from the trivial noninteracting point (specifically for J ′ 1,2 = J 1,2 = π/4 and J ′ 3 , J 3 ̸ = π/4 in Eq. ( 7)), the transfer matrix T has exactly t eigenvalues λ = 1, corresponding to onedimensional Jordan blocks, while all other eigenvalues are contained in a disk of radius r < 1.In fact, there is exactly one maximal-magnitude eigenvalue in each diagonal sector (ν, ν), and their corresponding eigenvectors read as In the r.h.s. of Eq. ( 29) we represented the operator in Eq. ( 27) as a state of a Hilbert space with doubled dimension using to the operator-to-state mapping The presence of t dominant eigenvalues with unit magnitude allows one to simply show that K(t, L) is indeed described by the CUE form in Eq. ( 10) for large enough L. As pointed out in Ref. [50], this t-fold degeneracy of T , and the fact the eigenvectors preserve only the diagonal part of its Z t × Z t symmetry, indicates that the ramp in the spectral form factor is a manifestation of a spontaneous breaking of symmetry Z t × Z t to Z t .The goal of this paper is to understand if and why this spontaneous symmetry breaking is stable as one moves away from the dual unitary point, from the perspective of perturbation theory.

C. Minimal Example
Even though our theoretical analysis can be carried out for the full family of circuits, Eqs. ( 13)-( 14), for our numerical investigations it is useful to fix some of the parameters and consider a minimal toy model example.Specifically, we consider x = 0, and average over θ (1) x , θ x , ϕ x , ϕ x with a Gaussian measure with zero mean and infinite variance, i.e., maximal disorder strength.This gives where we set s = 0, 1/2.These choices simplify drastically our numerical analysis as they reduce the number of parameters to only two, however, they still contain a rich phenomenology.The two extremal cases are found for ϵ 1 = ϵ 2 = 0, when the model corresponds to an ensemble of ergodic dual unitary circuits, and for ϵ 1 = ϵ 2 = π/4 when the model is trivially localised as there is no coupling between different sites.We verified that including the term Jσ (3) ⊗ σ (3) in the local gate (J has been set to 0 in Eq. ( 31)), or a finite disorder variance, does not qualitatively modify the numerical results.For instance, in Fig. 3 we show the gap in the 0,0) momentum sector of T at the dual unitary point as a function of J and t.We note that the gap is largest when J = 0, and closes (∆ 0 → 0) when J → π/4, i.e., when the local unitary U equals the swap matrix and the circuit becomes non-interacting.Moreover, in Eq. ( 33) we took σ → ∞.To confirm that this choice does not introduce non-generic behaviour we investigated the gap of the space transfer matrix for various values of finite variance.A representative example of this is reported in the bottom panel of Fig. 3, where we took J = 0. We see that σ 2 /2 = 1 we are already close to the σ = ∞ limit.For example at t = 8, J = 0, σ 2 /2 = 1 the gap is found to be ∆ 0 (8) ≈ 0.660709.Instead taking σ → ∞ one gets ∆ 0 (8) ≈ 0.6607102.That is, the correction to the gap is observed in the fifth decimal place.We also emphasise that the choice of an infinite variance σ 2 → ∞, i.e., a maximal disorder strength, should be the most challenging for the survival of dual unitary behaviour away from the dual unitary point.Indeed, this is the optimal regime in which one might expect Floquet many-body localisation (MBL).
The number of parameters can be reduced further by fixing the ratio between ϵ 1 and ϵ 2 .The value of the ratio, however, affects rather drastically how dual-unitarity is broken and has to be chosen with care.Here we choose the two values of ϵ 2 /ϵ 1 corresponding to weakest and strongest breaking of dual-unitarity.To find them we note that ρ Ũ = Ũ Ũ † /4 is a positive matrix with unit trace [59] and can be interpreted as a quantum state.Therefore, the dual-unitarity breaking can be estimated by computing the fidelity between ρ Ũ and the maximally FIG. 3. The gap ∆0(t) = 1 − |λ1| at the dual unitary point for various t and J, all data is collected from (ν, ν ′ ) = (0, 0) symmetry sector.The top plot takes α → ∞ as in Eq. ( 33) while the bottom plot sets J = 0 and takes σ finite.The bottom plot is the only data set in this article where σ is taken to be finite, all other datasets and figures take σ → ∞.
The data showcased in the top plot with J = 0 is featured in Fig. 9.
mixed state ρ ∞ = 1/4.In particular, an explicit calculation gives where |A| = √ AA † .This means that, considering without loss of generality the two cases corresponding respectively to the weakest and strongest breaking of dual unitarity are ϵ 2 = 0 and ϵ 2 = ϵ.In the following we consider both these cases referring to them, respectively, as Case I and Case II.

D. Numerical Approaches
Besides standard exact diagonalisation of (1) for small systems, in this paper we provide three numerical tests to characterise the space transfer matrix T .The first is an Arnoldi exact diagonalisation method to converge several leading eigenvalues simultaneously.We use this to study the full spectrum of T without resolving the translation symmetry in Eq. (25).Due to the number of eigenvalues we want to compare and the behaviour of the spectrum this method is limited to t ≤ 6 (our full forward-backward time lattice comprises of 4t qubits).
Next, we isolate a specific double-momentum sector (ν, ν) and characterise the action of the transfer matrix within the sector (a detailed discussion of how this is achieved is presented in Appendix C).This removes degeneracy near the spectral edges of T and in most cases allows us to study the spectrum with the power iteration method, allowing us to access times t ≤ 9 (for some data sets we had to use the Arnoldi method also within a given sector as the gap was too small).Moreover, we use this representation to evaluate the coefficients in our pertubative analysis of Sec.IV.
Our third method is a Monte Carlo-based approach that approximates the maximal eigenvalue in a given double-momentum sector by stochastically unraveling the average in Eq. ( 18) and allowing us to approach times t ≤ 12.The idea is to observe that, before the average, the transfer matrices in Eq. ( 18) are written as the tensor product of two operators (related by complex conjugation) acting only on the forward and backward lattices.Therefore, if we do not perform the average explicitly we can consider only one of the lattices, halving the number of qubits we need to simulate.More concretely we proceed as follows.Instead of considering directly Eq. ( 19) we construct local gates on the time lattice as where u = exp iθ (1) σ (1) exp iθ (3) v = exp iϕ (1) σ (1) exp iϕ (3) σ (3) .
Therefore, the transfer matrix is characterised by four angles ⃗ θ = (θ x , θ x , ϕ x , ϕ x ), which are uniformly generated random numbers, and can be written as We then observe that for large enough N where |ψ⟩ is a random state on the forward lattice, c 0 = ⟨λ 0 |(|ψ⟩ ⊗ r |ψ⟩ * ), and the average is performed over all choices of ⃗ θ n , n = 1 . . .N .We now estimate the the left hand side of Eq. ( 41) by sampling over the choices of ⃗ θ n .Letting we have where m labels the samples of ⃗ θ 1 , . . ., ⃗ θ n and Λ denotes the sample size.Therefore, plugging back into (41) we obtain The left hand side of this equation can be easily computed for a large number of samples and moderately large t because it is defined only on the forward lattice (the objects involved live in a vector space of half the size).A simple way to extract λ 0 from this relation is to take the logarithm of both sides and find the slope of the data as a function of N .In what follows we take Λ = 10 7 for small t and reduce our sample size to Λ = 10 5 for t = 11, 12.This method works best when λ 0 is significantly larger than one.Finally we mention that if one is only interested in characterising the evolution in time of the leading eigenvalues of T one can use the time-evolution-based method recently introduced in Ref. [49].This method assumes a certain structure for the spectrum of T and evaluates the leading eigenvalues by performing an evolution in time in a quantum circuit with finite length and judiciously selected twisted boundary conditions.Here, however, we are interested in an assumption-free characterisation of T which goes beyond its leading eigenvalues (cf.Sec.IV).Therefore, we do not use this approach.

III. STABILITY OF THE ERGODIC PHASE: NUMERICAL SURVEY
As recalled in Sec.II B, on the dual-unitary manifold the quantum circuit is chaotic in the sense that its spectral form factor exhibits the linear ramp characteristic of random matrix theory (cf.Eq. ( 10)).This property is found to correspond to the space transfer matrix T having t eigenvalues equal to one.In this section we investigate how the spectrum of T behaves when we move away from the dual unitary point using the model in Eq. (31).
Our expectation is that the ϵ-dependence of the spectrum of T should be smooth for small enough ϵ.This means that there should exist a phase where the spectrum of T is qualitatively similar to the that of the dualunitary point: The leading eigenvalue, or SLE, should have an approximate t-fold degeneracy near unity (which becomes increasingly exact as t → ∞) and the corresponding eigenvectors should lie in the diagonal momentum sectors (i.e.(ν, ν)).This phenomenology would be consistent with the spontaneous symmetry breaking scenario proposed in Ref. [50] (cf.Sec.II B).More concretely, we expect that at late enough times Eq. ( 24) is dominated by the leading eigenvalues so that For ϵ = 0 this approximation is exact, the eigenvalues are all equal to 1, so that K(t, L) = t for large enough L. When ϵ ̸ = 0, the requirement is that a linear ramp ensues on time scales in excess of the so called Thouless time (expected to be sub-polynomial in L in the absence of conservation laws).Combined with Eq. ( 45), this suggests that the leading eigenvalues in each sector are split from 1 by an amount that decays exponentially in time, i.e., On the other hand, an MBL phase should be characterised by a single leading eigenvalue going to λ 0,(0,0) = 4 for large t, while all remaining eigenvalues having magnitude less than unity.In the language of Ref. [50] this corresponds to a symmetry-unbroken phase.
Let us now proceed to substantiate these expectations using our exact numerical approaches.Even though maximal-magnitude eigenvalues exist in all the (ν, ν) symmetry sectors, we begin by focussing our attention on the (0, 0) sector (and temporarily drop the sector label) as this gives us the ability to investigate longer times (larger time-lattice sizes).At the end of this section we will discuss all symmetry sectors simultaneously.We observe that in the (0, 0) sector, ϵ ̸ = 0 implies λ 0 (t) > 1 for all the observable times t, that is, perturbing away from the dual unitary point increases the leading eigenvalue (Fig. 6).Conversely, in the other sectors we detect a decreasing in the size of the maximal eigenvalue, i.e., |λ 0,(ν,ν) | < 1 for ν ̸ = 0 and ϵ ̸ = 0.
In Fig. 4 we report ln(λ 0 (t) − 1) versus t for both Cases I and II and variety of values of ϵ, all relatively small.We find good agreement with the following scaling form where γ(ϵ) > 0 and c(ϵ) are constants depending solely on ϵ.In general, we observe that γ(ϵ) decreases monotonically with increasing ϵ.In particular, we find that γ(ϵ) is substantially larger in Case I than in Case II (cf.Sec.II C).More specifically, we observe that taking ϵ twice as big in Case I compared to Case II produces similar γ(ϵ).For example, γ(0.05)| CaseII ≈ 0.6191 while γ(0.1)|CaseI ≈ 0.5962.This persists for all ϵ tested.We similarly observe γ(0.3)CaseII ≈ 0.0687 and γ(0.6)CaseI ≈ 0.0805.For larger values of ϵ the exponent γ(ϵ) becomes too small to be reliably determined in the time window accessible by our exact methods.However, our numerics indicate that for Case I γ ϵ → π 4 is non-zero.This suggests that for Case I ergodicity is stable for the entire parameter regime.In principle one can explore larger times by using the time-evolution based approach introduced in Ref. [49] (cf.Sec.II D), however, here our focus  is chiefly on small ϵ.We further note the apparent zigzag pattern in the data of Fig. 4 differentiating even and odd t.We associate this behaviour with the fact that for even t there is an additional parity symmetric sector of T besides (0, 0), i.e., (π, π).This changes the structure of the eigenspaces (and their dimensions) also affecting the value of λ 0,(0,0) .Having discussed the behaviour of the leading eigenvalue in the (0, 0) sector we now move on and study the gap between the latter and the rest of the spectrum.In particular, in Fig. 5 we report as a function of time and for different choices of ϵ (again for both Cases I and II).In all the cases explored we find that the gap satisfies ∆(t) > 0 in (0, 0) sector and, therefore, λ 0 (t) is sufficient to characterise the large L behaviour of the SFF.This means that, for small ϵ, the ν = 0 component of Eq. ( 45) is confirmed.
To conclude our numerical test of Eqs.(45) and (46) it remains to check that this behaviour is mirrored in other symmetry sectors.We achieve this by using the Arnoldi method to check all sectors simultaneously.Our results are reported in Fig. 6.Interestingly, the leading eigenvalue in the (0, 0) sector is always observed to be the dominant one.Along with this property, it is observed to converge to 1 more slowly than the other leading eigenvalues.For example the t = 6, ϵ = 0.6 data set for Case I has λ 0,(0,0) = 1.45898 while the second furthest from unity is λ 0,(1,1) = λ 0,(5,5) = 0.838665.Leading eigenvalues are always found to be real numbers, while sub-leading values in general are real or complex with magnitude smaller than unity.The leading eigenvalue in the (0, 0) sector is the only one to be consistently greater than unity, other leading eigenvalues can oscillate around 1 as functions of time.
In summary, our numerical analysis is consistent with the expectation that for small ϵ the spectrum of T is a smooth deformation of that at the dual unitary point: We have t real eigenvalues close to unity while the rest FIG. 5.The gap ∆ = |λ0| − |λ1| versus time t for various choices of perturbation.This data was retrieved using the Arnoldi method in the (ν, ν ′ ) = (0, 0) symmetry sector.
have significantly smaller magnitudes.We stress that this is the case despite the fact that our minimal model (31,32) involves maximal disorder strength.Interestingly, in all our numerical experiments we also observed that the (0, 0) sector has the largest eigenvalue.This is consistent with the symmetry breaking picture of Ref. [50]: as the system is perturbed away from the maximally ergodic point, i.e., the dual unitary point, a single symmetric eigenvalue becomes the dominant one.Further numerical surveys suggest that this phenomenology continues also when moving away from our minimal model (31,32).

IV. PERTURBATION THEORY
In this section we propose an analytical explanation for the main observation of Sec.III.Namely, that the dualunitary eigenvalue structure is maintained for finite ϵ suggesting structural stability of the dual-unitary phase.The idea is to fix ϵ and write the maximal eigenvalue of T in each diagonal double-momentum sector (ν, ν) -we denote it by λ (ν,ν) ≡ λ 0,(ν,ν) -as a perturbative series in an auxiliary parameter.Studying this series we then show that, if two assumptions are fulfilled, then |λ (ν,ν) −1| is bounded by a term that is exponentially small in t, implying that our circuit models are ergodic at the considered value of ϵ.Remarkably, this happens even for the minimal model with maximal disorder strength discussed in Sec.II C. Specifically, our two assumptions are (i) No "maximal level" crossing occurs in the perturbative expansion, i.e., the evolution of the maximal eigenvalue in each sector can be followed by tracing the smooth deformation of the maximal eigenvalue at the dual-unitary point.
(ii) The (ϵ-dependent) coefficients of our perturbative series are bounded by an exponentially decaying function of t and grow at most exponentially in n, where n is the perturbative order.
A more precise formulation of these assumptions is given in the upcoming derivation.The first assumption is safe for generic enough perturbations: if the perturbation couples all the eigenvectors in a given sector, all adjacent level encounters are avoided crossings.Therefore, Assumption (i) can be rephrased by saying that we assume our perturbation to be sufficiently generic.This assumption is consistent with our numerical survey of Sec.III (cf.Fig. 5).The second assumption is our main one and, as we discuss in Sec.V, can be partly justified by an analytical argument.In Sec.V we also show that Assumption (ii) fails at the trivially localised point (i.e. for Case II and ϵ = π/4, cf.Sec.II C) while we give numerical evidence of it holding for small ϵ.
We now proceed to show that when (i) and (ii) hold the quantum circuit (1) is ergodic.We begin considering the eigenvalue equation for the maximal eigenvalue λ of the transfer matrix resolved to the double momentum sector (ν, ν), namely Here and in the following, we drop the dependence on the subscript (ν, ν) whenever it is not ambiguous to do so.Since for ϵ = 0 only the diagonal sectors contain the maximal eigenvalues, and the latter are unique, Eq. ( 49) contains all the necessary information to characterise the spectral form factor for small enough ϵ.It is not obvious, however, that this will continue to hold also for finite ϵ.
Here we assume this to be the case, namely we make the following postulate Assumption 1 (No leading eigenvalue crossing).The leading eigenvalues of T for small enough ϵ are obtained by smooth deformation of those of T | ϵ=0 .Namely, λ is a smooth function of ϵ.FIG. 6. Full spectrum λ n,(ν,ν) analysis (all diagonal double momentum sectors) of T for various ϵ and t.Points are plotted in polar coordinates, with the radius being the magnitude of the eigenvalue (this plot therefore covers up some degeneracy in the sub-leading eigenvalues).The polar angle is 2πν/t, where ν labels the symmetry sector (ν, ν).Results were obtained by an Arnoldi method converging n = 12 eigenvalues at the edge of the spectrum.
By virtue of Assumption 1 we can limit our treatment to diagonal double-momentum sectors and only consider Eq. ( 49), which we treat using a "dressed" perturbative approach.Specifically, we proceed as follows.First we set and rewrite Eq. ( 49) as Next we solve in perturbation theory in x for fixed ϵ.Finally we recover a perturbative solution of Eq. ( 49) by setting x = ϵ in the end.The advantage of this approach is that it involves only a first order correction to the transfer matrix as it happens in standard time-independent perturbation theory in quantum mechanics, at the cost of making the perturbation manifestly ϵ-dependent.
Explicitly, we expand both |λ⟩ and λ in x and impose (52) order by order in x.As shown explicitly in Appendix A, this yields (55) and where we set Note that G is the resolvent R(z) = (z1 − T 0 ) −1 of T 0 projected away from the leading-eigenvalue subspace and evaluated at z = 1.The projection makes this operator well defined.Using these expressions one can write λ (n) in terms of the following expectation values of products of the perturbation Tϵ and the projected resolvent Explicitly, the first few orders read as Continuing to arbitrary order we have where δ[x] and θ[x] are equal to one when x = 0 and x ≥ 0 respectively and to 0 otherwise, while C[{k ij }] counts the combinatorial multiplicity of a given term (we do not need its explicit expression).Here q denotes the number of [k 1 , . . ., k m ] symbols appearing in a given term, p m is the length of the m-th symbol.Each symbol contains p m + 1 factors of the perturbation Tϵ , while each [ ] contains one factor.Therefore, the final line of Eq. ( 62) contains in total n factors of the perturbation, consistent with it being an n-th order perturbative term.
To find the constraints we used that p i +1 is the number of Tϵ in [k i1 , . . ., k ipi ] and pq j=1 k ij is the number of Gs.Since the total number of Tϵ in each term at order n equals n we have where m is the number of [ ]'s in the term.On the other hand from the last of ( 56) and the last of ( 59) we have that each term contributing to the n-th order in perturbation theory contains n − 1 occurrences of G. Therefore we find Our goal is to bound from above the magnitude of the perturbative correction λ (n) to the maximal eigenvalue.
To this end we make the following assumption on the scaling of the symbols [k 1 , . . ., k m ] Assumption 2. The symbols [k 1 , . . ., k m ] and [ ]can be bounded as follows where β > 0, α, γ, t 0 are independent of m, k j and t.
Using Assumption 2 in Eq. ( 62) we obtain the following bound where N n is defined as This number can be computed with three simple observations.First we note that (56) implies Next, we observe that N n>2 can be alternatively written as where # k1,...,k ℓ denotes the number of terms in the ex- we immediately find N 3 = 3 and the following recursive relation where we used N 1 = N 2 = 1.This means that N n+1 fulfils the recursive relation of the Catalan numbers with the same initial condition.Therefore Plugging back into (54) we find x n e (α+γ)n N n . (72) To conclude we observe that the sum on the r.h.s. is always convergent for small enough x.Namely we have convergence whenever x ≤ ϵ ≤ e −(α+γ) /4.
For all the values of ϵ fulfilling the above bound we then have This expression recovers Eq. ( 46) and shows that whenever Assumptions 1 and 2 hold the ergodic phase is stable.
V. DISCUSSION OF ASSUMPTION 2 Our Assumption 2 on the behaviour of the perturbative coefficients in Eq. ( 60) can be justified by an analytical argument assisted by numerical observations.For definiteness we again focus on the sector (ν, ν) = (0, 0), although other double momentum sectors show similar behaviour.We begin by considering the simplest of the coefficients in Eq. ( 60), i.e., and compute it numerically for n = 1, . . ., 30 and t = 3, . . ., 8. Some representative examples of our results, for both Cases I and II, are reported in Figs.7 and 8. Overall we see that, in agreement with Eq. ( 65), the term increases exponentially as a function of n and is exponentially suppressed as a function of t.
A more refined analysis is provided by fixing t, varying n, and performing a linear fit.Namely we set and find α(t, ϵ) and δ(t, ϵ) providing the best fit [60].Studying these coefficients (cf.Fig. 9) we find that, very interestingly, the slope α(t, ϵ) is roughly independent of ϵ.Moreover -and this is a key observation -it matches remarkably well the logarithm of the inverse of the spectral gap calculated at the dual-unitary point, i.e., In fact, since in this case the largest sub-leading eigenvalue is unique and real (cf.Tab.I) we have ∆ 0 (t) = 1 − λ 1 (t).In fact, our Assumption 2 requires δ(t, ϵ) to decrease linearly in t.Here our data are less convincing and sensitive to the parity of t, but have the correct overall trend.This result can be reproduced by making two assumptions on the structure of the dual-unitary transfer matrix T 0 .First, we assume that T 0 is diagonalisable: this seems a reasonable assumption given that T 0 is an average over matrices, and that defective (i.e.non-diagonalisable) matrices are non-generic.In fact, our upcoming reasoning continues to hold also when there are non-trivial Jordan blocks but only for the eigenvalue 0. The latter requirement is easier to check numerically, and it is fulfilled in all our numerical observations [61].Therefore, we write G as [62] where |λ j , l⟩ and ⟨λ j , r| are respectively the (orthogonal) right and left eigenvectors corresponding to the j-th subleading eigenvalue λ j .Next, assuming 1 − λ 1 ≪ |1 − λ j>1 | we truncate the spectral decomposition (78) to the leading eigenvalue  9. We in general observe the value to be real for this choice of parameters.We also include λ2(t) which was extracted from the symmetry resolved Arnoldi method.
[n] terms in perturbation theory as a function of t and n.Solid dots represent the natural logarithm of data retrieved through exact evaluation of [n] in the (ν, ν ′ ) = (0, 0) sector.Dotted lines are numerical fits indicating exponential dependence on the independent variables t, n.
This gives which is consistent with the numerical observation (77).
In fact, the approximate form (79) of the resolvent can also be used to explain the exponential decay in time observed in the upper panels of Figs. 7 and 8.We begin noting that for a large enough ℓ ∈ N one can make the following approximation where Q is the projector on the leading eigenspace of T 0 (cf.Eq. ( 59)) and we neglected terms that are prima facile O((λ j /λ 1 ) ℓ ).This approximation is in fact more subtle than it might appear because the terms |λ j , r⟩⟨λ j , l| /(⟨λ j , r|λ j , l⟩) can have large operator norm (possibly even exponentially large in t).This means that one might need to consider ℓ = O(t) to safely neglect higher order terms.Here we assume that this is not the case and take ℓ to be O(t 0 ).
Using (81) we can rewrite Eq. (80) as Next, we rewrite the last line in terms of the original time evolving gates, undoing the space-time duality transformation discussed in Sec.II B. The idea is to represent the expectation value of a product of n transfer matrices on a state in terms of the time-evolution operator of a chain of n qubits.The state translates into the boundary conditions imposed on the time-evolution operator.In our case the boundary conditions will pair forward and backward evolution, giving a non-unitary boundary term to the evolution operator.More concretely, considering for FIG. 8.
[n] terms in perturbation theory as a function of t and n.Solid dots represent the natural logarithm of data retrieved through exact evaluation of [n] in the (ν, ν ′ ) = (0, 0) sector.Dotted lines are numerical fits indicating exponential dependence on the independent variables t, n.
instance the first term we have where |O⟩ is the state corresponding to the operator O under the mapping in Eq. ( 30).In the second step we dropped terms that are at most O(2 −t ) by using and noting that the sum on the r.h.s. is bounded by 2 t .The terms on the r.h.s. of Eq. ( 83) can be easily translated in the time evolving picture because the states |Π 2j 2t ⟩ implement simple pairings between backward and forward evolution.For instance, the term with τ = 0 is written as where we introduced the 4 2n−1 × 4 2n−1 matrix The latter is written in terms of the time evolution operators for 2n − 2 qubits (∈ End(C 2 2n−2 )) and the boundary matrices (∈ End(C 4 )) Introducing a convenient diagrammatic representation for objects acting on both the forward and backward time we can depict B n,0 as in Fig. 10a.Analogously, a generic term with τ ̸ = 0 is written as where we introduced the matrices Note that B n,τ ∈ End(C 4 2n+τ −1 ) and b τ ∈ End(C 4 2 ).Introducing the following diagrammatic representation for b τ , we can depict B n,τ as in Fig. 10b.The traces of B n,τ can be treated following Ref.[63].In particular, using Theorem 1 of the aforementioned reference we have that if there are no x ≤ y such that for some local operators a, a ′ , b, b ′ , then tr B t n,0 = 1 + O(e −βt ) , β > 0. (99) Note that, although β > 0 for all values of n, Ref. [63] gives no information on its n dependence.With a similar reasoning we prove in Appendix B that if (97) does not hold for any x, then ρ(B n,τ ) < 1, where we used ρ(•) to denote the spectral radius.This implies Since for (ϵ 1 , ϵ 2 ) ̸ = (π/4, π/4) the gates fulfilling (97) and (98) have measure zero in the disorder average, we conclude that where the constant A(ϵ) vanishes for ϵ = 0 because the l.h.s. is trivially equal to one for ϵ = 0 while Fig. 8 suggests that β(0) is finite.Note that, since we do not control the β dependence on ℓ, we cannot exclude that it approaches 0 in the limit of infinite ℓ.This is why we had to assume ℓ = O(t 0 ) in Eq. ( 81).Proceeding analogously we find with B(ϵ) ≃ A(0) ′ ϵ/2 for small ϵ.Putting all together in Eq. ( 82) we then have where C 2 (ϵ) is O(1) for small ϵ.If we compare Eq. ( 103) with our two parameter fit for [n] (Eq.( 76)), we predict that −β is the slope of δ with respect to t.The right panel of Fig. 9 then suggests that β(ϵ) depends weakly on ϵ for small enough ϵ.
Proceeding along similar lines we can estimate all the coefficients in Eq. ( 60).In particular, using the approximations (79) and (81) we have that the generic coefficient While applying (99) and (100) we have with C m (ϵ) = O(1) for small ϵ.This provides an analytical justification to Assumption 2. In Figs.11 and 12 we provide an independent numerical test of this assumption by considering the behaviour of the flat coefficients These contributions are those for which the approximation (79) is the least justified as the latter becomes exact only when the resolvent is taken to infinite power.Nevertheless, from Fig. 11 we clearly see an exponential decay in time in agreement with Eq. ( 105) and hence with Assumption 2. This despite the relatively short times accessible in our numerical simulations.Instead, Fig. 12 reports the behaviour of [1, 1 . . .1] as a function of m, i.e. the number of ones in the coefficient.We see that the coefficient decays exponentially in m.This suggests that C m+1 (ϵ) in Eq. ( 105) is bounded by an exponentially decaying constant.In fact we note that, at least for Case I and small enough ϵ, Fig. 11 shows that the exponential decay in time becomes stronger when m increases.This can be explained by our asymptotic form Eq. (105) if we admit that the factor λ 1 (1 − λ 1 ) appearing in the denominator increases as a function of time.Note that this growth cannot be unbounded as λ 1 (1 − λ 1 ) ≤ 1/4.Finally we stress that Eq. (105) do not hold at the trivially localised point ϵ 1 = ϵ 2 = π/4.Indeed, in that case Eqs.(99) and (100) do not apply as Eqs.(97) and (98) are clearly satisfied for all x and y.As a result, in this case there is no exponential decay in time of the coefficients (60).For instance, using the simple form of T at the localised point (see, e.g., Ref. [51]) it is easy to show

VI. CONCLUSIONS
In this work we laid down a general framework to investigate the structural stability of dual-unitary spectral correlations, which can be loosely thought of as the quantum many-body analogue of the theory of structural stability of hyperbolic flows in classical chaotic dynamical systems [37,38].
Our guiding principle has been that, contrary to integrable systems, dual-unitary systems should be robust under typical perturbations as they are quantum chaotic.Therefore, the spectral correlations of a perturbed dualunitary system, or at least their universal part, should be accessible by devising an appropriate perturbation theory.Here we formulated such a perturbation theory and identified two key assumptions needed for a rigorous proof of its convergence.We then provided a compelling numerical evidence for the validity of these assumptions in a particular family of perturbed Floquet dual-unitary circuits, and corroborated them with a heuristic analytical argument (supported by numerical evidence) involving the spectral decomposition of the resolvent of the unperturbed transfer matrix.
Besides their consequences on the structural stability of dual-unitary correlations, our findings have an important consequence on the interplay between ergodicity and disorder.Indeed, since spatial disorder does not affect dual-unitarity breaking (only two-body couplings can break dual-unitarity), our results imply that whenever a quantum many-body system is close to an interacting (non-SWAP) dual unitary point its spectral correlations are always random-matrix-like, irrespective of the disorder strength -for instance, the particular family used in our numerical analysis has maximal disorder strength.This rules out the possibility of Floquet MBL in the thermodynamic limit for systems close enough to the dual-unitary point.
Although these findings are remarkable, they are in many ways only a stepping stone to the development of a comprehensive theory of structural stability for Floquet dual-unitary circuits, and many key questions remain open.In particular, it would be important to explore the degree of generality of the structural stability we found in dual unitary circuits by investigating more general perturbations.Moreover, further research is required to corroborate and expand our perturbative approach.To this end, we identify two main directions.
The first is to find quantitative estimates or bounds for the radius of convergence of our perturbative expansion.Indeed, at the moment we have merely shown that, under our two assumptions, the radius is finite.However, we gave no information on its value.A quantitative estimate of the radius of convergence could potentially lead to the identification of the point of transition to the nonergodic (i.e.localised) regime, which might be occur for a finite value of two-body coupling or only at the trivially localised point where the qubits are disconnected.A related question is whether one can identify the transition point expanding around the trivially localised point.The second direction is, of course, to provide rigorous mathematical proofs of our assumptions, in particular to the second one that appears the more substantial.We believe that the heuristic analytical argument we provided in support of that assumption can be used as blueprint for such a proof.In this section we briefly describe how to reduce the numerical analysis to a given double-momentum sector.As mentioned in the main text the transfer matrix T is invariant under two-site translations in the forward and backward lattices [Π 2τ1  2t ⊗ Π 2τ2 2t , T ] = 0, τ 1 , τ 2 = 0, . . ., t − 1. (C1) For the forward and backward lattice the most natural basis to work in is therefore the eigenbasis of the two-site shift operator.Considering the forward time lattice, we have 2t qubits, and we know we have Π 2t 2t = 1, giving eigenvalues e 2πiν/t with ν = 0, . . .t − 1.To generate the eigenbasis we select a set of reference states |f ⟩ (taken to be product states in the computational basis) and write where L f is the period of the reference state Π 2L f 2t |f ⟩ = |f ⟩.Typically L f = t however some special states will have periods that are integer multiples of t.Computationally this representation reduces our overall storage from 2 2t → 2 2t /t approximately, with the largest sector being the ν = 0 sector (see, e.g., Ref. [64] for more details).For the purposes of discussing complexity later in this section we will call D = 2 2t and D ν = 2 2t /t.
An important observation about T is that it is made up of a product of operators which independently are translation invariant under two shifts.This is evident from Eq. ( 20) and implies that we can update vectors in distinct steps.First let us treat the case where we do not couple the forward and backward lattice, i.e., we consider an operator of the form where Ũ is a two-local operator (it acts non-trivially only on a pair of nearest neighbours).We will use the above as an example to illustrate working in the momenta basis for operators with this general structure.Since Where the sum over f is taken over the set of representation basis states.One memory inefficient way to evaluate this expression is to simply evaluate Ũo |f ⟩ in the full Hilbert space, and then compress the state back into the translation invariant representation.This approach is computationally costly and likely memory bound, severely slowing down the code.It also removes the advantage of working in the symmetry resolved basis by increasing storage requirements to D, which, once we couple the forward and backward lattice, will eliminate our advantage with this approach.In fact, an open question we were not able to answer is how to evaluate Ũo |ψ ν ⟩ faster than O(D 2 ν ).The operations similarity to a discrete Fourier transform indicates this may be possible.
Instead of updating the vector directly from the expression (C5) we focus our efforts on computing matrix elements of the operator This allows us to store it in a D ν ×D ν dimensional matrix, the same size as we will see, as the many body vectors once we combine the forward and backward lattice.A vector on the full space can be represented by where the sum is taken over the set of representation basis states.Updating the full forward and backward lattice state with the symmetry resolved operator is now given by The final piece of the puzzle is to understand the operation We focus here on O (3) 0 as it is diagonal in the computational basis.Other non-diagonal terms can be evaluated with a simple basis rotation, and then following the steps we will outline.The action of this operator on the translation invariant basis is trivial, we have,

FIG. 1 .
FIG. 1. Diagrammatic representation of tr U t .The boxes represent local gates and different legs act on different spatial sites.Matrix product is represented by joining legs and goes from bottom to top.The lines at the left and right edges are joined because we consider periodic boundary conditions (L ≡ 0), while those at the top and bottom are joined because of the trace.The background grid specifies the space-time lattice.The gates in the same vertical column are identical.

30 FIG. 4 .
FIG.4.ln(λ0(t) − 1) as a function of time for various ϵ1, ϵ2.Dotted lines indicate numerical fits.Circle data points are retrieved with the power method isolated in the (ν, ν ′ ) = (0, 0) symmetry sector.Diamond data points are calculated using the Monte Carlo.Monte Carlo data consists of 10 7 samples for t ≤ 10 and t > 10 use 10 5 data points.In the top two panels we supply exact values and Monte Carlo estimates for all times t, while in the bottom two panels we simply plot exact values for t ≤ 9 and the remaining points are Monte Carlo estimations.
⟨m ν | Ũo |f ν ⟩ = The above equation is simpler to evaluate.To see this without loss of generality take r = p = 0.Because Ũo is made up of a product of commuting terms the expression factorises⟨m| Ũo |f ⟩ = t j=0 ⟨m 2j m 2j+1 | Ũo |f 2j f 2j+1 ⟩, (C7)where we have broken the representative state into its computational basis form for individual time lattice qubits.Note that Eq. (C6) can be computed in O(t) steps due to repeated computations.For convenience we will call this new symmetry resolved operator Ũ(ν) 0,(m,f ) ≡ ⟨m ν | Ũo |f ν ⟩.(C8)