Correlations in Perturbed Dual-Unitary Circuits: Efficient Path-Integral Formula

Interacting many-body systems with explicitly accessible spatio-temporal correlation functions are extremely rare, especially in the absence of integrability. Recently, we identified a remarkable class of such systems and termed them dual-unitary quantum circuits. These are brick-wall type local quantum circuits whose dynamics are unitary in both time and space. For these systems the spatio-temporal correlation functions are non-trivial only at the edge of the causal light cone and can be computed in terms of one-dimensional transfer matrices. Dual-unitarity, however, requires fine-tuning and the degree of generality of the observed dynamical features remained unclear. Here we address this question by introducing perturbations. First we show that if the deviation from dual-unitarity is random and independently distributed at each space-time point, dynamical correlations maintain the dual-unitary form. Then, considering fixed perturbations, we prove that for a particular class of unperturbed elementary dual-unitary gates the correlation functions are still expressed in terms of one-dimensional transfer matrices. These matrices, however, are now contracted over generic paths connecting the origin to a fixed end point inside the causal light cone. The correlation function is given as a sum over all such paths. Our statement is rigours in the"dilute limit", where only a small fraction of the gates is perturbed, and in the presence of random longitudinal fields, but we provide theoretical arguments and stringent numerical checks supporting its validity even in the clean case and when all gates are perturbed. As a byproduct, in the case of random longitudinal fields -- which turns out to be equivalent to classical Markov circuits -- we find four types of non-dual-unitary interacting many-body systems where the correlation functions are exactly given by the path-sum formula.

Understanding the dynamics of extended quantum many-body systems with local interactions is the core problem of nonequilibrium statistical mechanics, with a wide range of applications ranging from condensed matter physics to high energy theory and quantum gravity. In particular, the set of two-point spatio-temporal correlation functions of local observables can be considered as the prime quantifier of the dynamics. For example, they can be used in the framework of linear response theory [1,2] to express coefficients, such as conductivities and kinematic viscosities, that describe macroscopic transport properties.
A major obstacle is that computing dynamical correlation functions in interacting systems is notoriously hard. This is true for numerical simulations of correlations in real time, which are typically exponentially hard (in physical simulation time) [3,4] (see also [5]), and even more so for analytical computations. Even in Bethe-ansatz integrable systems the task turns out to be daunting: while recent breakthroughs allowed for the calculation of dynamical correlations in the late-time, hydrodynamic, regime [6][7][8], the computation of two-point correlations at arbitrary intermediate times remained out of reach. This leaves us with non-interacting, quasi-free, or Gaussian theories, as the only general class of systems where dynamical correlations can be analytically computed in general. Dynamical correlations in coupled (interacting) theories are then usually formulated in terms of (Keldysh) diagrammatic many-body perturbation theory [1]. In these approaches correlations are written as power series in the coupling constant around the underlaying free theory. Such perturbative series, however, generically have a vanishing radius of convergence -especially in the thermodynamic limit -and their relevance for determining the actual physical behaviour of correlation functions in extended systems is questionable. One typically finds an even qualitative change in the behaviour of correlations and transport coefficients at finite temperatures when going from free to strongly coupled theories. While the former usually behave as ballistic conductors, the latter are expected to display diffusive transport. The latter is thus a non-perturbative effect and establishing its microscopic origin remains one of the main challenges of statistical physics [9,10].
Very recently, we proposed a new class of locally interacting many-body systems, which allows for the exact computation of spatio-temporal correlation functions of local observables [11] and of some other indicators of quantum chaos and scrambling of quantum information [12][13][14][15][16][17][18][19][20][21]. These systems are special brick-wall type quantum circuits (see, e.g., Ref. [22][23][24]) called "dualunitary" quantum circuits and are characterised by unitary evolution, not only in the direction of time propagation, but also in the orthogonal direction of space propagation. Owing to this property, dual-unitary systems allow for a conceptual exchange of space and time axes and, in a vague sense, are discrete-space-time analogues of conformal field theories. In particular, the correlations in these systems can propagate only at the maximal speed and along one-dimensional straight paths in the space-time. One may think of these models as being 'statistically exactly solvable ', and yet, their dynamics is possibly ergodic [11] and quantum chaotic, both from the point of view of spectral statistics [12] and dynamical complexity [13,15]. This is similar to the situation in classical dynamical systems [25], where correlation func-tions of certain strongly chaotic single-particle models, such as e.g. Arnold cat maps or Baker maps, are exactly computable [26,27] despite their deterministic trajectories being unpredictable and uncomputable in the algorithmic complexity sense. Similarly, for generic quantum chaotic dual-unitary dynamics, the Heisenberg evolution of local operators is exponentially hard to simulate (classically) [13,15] even though the two-point functions at infinite temperature are exactly computable.
Dual unitarity, however, requires fine tuning of coupling parameters in the elementary two-body gates defining the systems and it is thus clearly not stable under generic perturbations. Nonetheless, one might wonder whether the properties of ergodicity and quantum chaos of dual-unitary circuits are structurally stable in a similar way as in classical ergodic theory. There, one can prove that linear automorphisms on the tori (Arnold cat maps) are topologically (structurally) stable under perturbations [28]. A quantitative measure of such stability may be a uniform continuity of two-point time correlation functions with respect to the perturbation parameter.
While one may imagine integrable and free extended quantum many-body systems as isolated points or low dimensional sub-manifolds (clearly structurally unstable in the above sense) in some high dimensional manifold of all appropriate (say, translationally invariant, locally interacting, etc.) systems, the chaotic and ergodic models may represent a finite domain with positive measure. As a matter of fact, the prevailing opinion of the scientific community is that almost all models are ergodic, depending on the precise definition of the ergodic class.
In this paper, we make the first step in addressing the question of structural stability by studying correlation functions in perturbed dual-unitary circuits. We consider three distinct classes of perturbations, which can be thought of as increasingly more realistic models of the generic situation: (i) Each of the elementary two-body gates composing the circuit is perturbed by an independent identically distributed random multiplicative gate (which can be considered the exponential of a hermitian random matrix); (ii) The perturbation is fixed (nonrandom, systematic dual-unitarity breaking) but the system is noisy: there are random local longitudinal fields at each space-time point over which we average; (iii) The perturbation is fixed and the system is clean. In the case (i), we show that the light cone structure of dualunitary dynamical correlations [11] is preserved (stable), the only effect of the perturbations is an additional exponential damping along light rays. In the case (ii), we rigorously show that there exists a class of generic (not fine-tuned) dual-unitary circuits that is "stable" under perturbations. More precisely, when these systems are perturbed, correlations continue to propagate along onedimensional paths in the space-time, though, these paths do not need to be straight anymore. As a result, all points of the causal light cone acquire non trivial correlations. To find the correlations between any two points one has to sum the contribution of all "skeleton diagrams", i.e., all allowed space-time paths connecting them. Although our rigorous result holds only in the limit of vanishing density of perturbed gates (or, at fixed, but arbitrary large order in perturbation theory), we provide further theoretical arguments and numerical evidence that it holds even when all gates are perturbed by a small enough perturbation. Moreover, in this setting we discover four new types of circuits where the two-point correlation functions are exactly (without approximations) given by the sum of skeleton diagrams. Finally, we formulate two simple analytical conditions for identifying the stable subclass of dual-unitary circuits, and we show numerically that the stable subclass exists also in the generic case (iii). Considering qubit gates we explicitly show that, for sufficiently small perturbation's strength, the sum of skeleton diagrams approximates exact two-point correlation functions with arbitrary precision.
Alongside the dynamics of unitary quantum circuits we also discuss those of a set of classical bistochastic many-body Markov chains, which can be written as local "Markov circuits". Indeed -as we show in the paper -the formal treatment of these two classes of systems is completely analogous. In particular, in Markov circuits the property of "dual-unitarity" is replaced by that of "dual-bistochasticity", meaning that the circuit is bistochastic when propagating in both space and time directions. Remarkably, the average over random longitudinal fields discussed at point (ii) maps quantum unitary circuits into Markov circuits. This means that our aforementioned rigorous results find a direct application also in the case of interacting classical stochastic dynamics.
The rest of the paper is structured as follows. In Sec. II we describe the technical setting of the problem and introduce the necessary diagrammatic notation, while in Sec. III we present the basic strategy and the main results of the paper. In Sec. IV we detail our findings in the case of random perturbations while in Sec. V we consider fixed perturbations in noisy or classical Markov systems. The validity of our conclusions for clean (non-random) perturbed dual-unitary circuits is demonstrated Sec. VI by combining analytical arguments with numerical results. Finally, Sec. VII contains our conclusions. Some technical details and proofs relegated to the appendices.

II. SETTING OF THE PROBLEM
We consider one-dimensional many-body systems composed of 2L qudits, each of them with d internal states. In other words, we examine a chain where at each site there is a quantum system with Hilbert space H C d and denote an orthonormal basis of H by The time evolution in this system is generated by discrete applications of "Floquet operators" U of the form Here we considered periodic boundary conditions, we assumed translational invariance, and we represented diagrammatically the unitary operator U x,x+1/2 = U ("local gate"), x ∈ {−L/2, . . . , −1/2, 0, 1/2, 1, . . . , Finally, we adopted the convention of time running upwards, namely in a product AB the symbol for the operator B is depicted below A. Introducing the diagrammatic representation we can express the unitarity of the gates in terms of the following diagrammatic rules A. Dynamical Correlations in the Folded Picture In this work we will adopt the so called "folded" representation of the circuit, a standard trick in tensornetwork theory [29,30] that has recently found many applications in studies on local quantum circuits [5, 13-15, 22, 23, 31-33]. Here we consider the simplest possible folding mapping that consists of turning operators into states of a quantum circuit with larger local Hilbert space by "folding" them in such a way that their upper and lower legs lie on top of each other. More specifically, the folding mapping turns operators in End((C d ) ⊗2L ) into states in (C d 2 ) ⊗2L . For example, for local operators we have Note that in this paper we will always consider operators that are Hilbert-Schmidt normalised, which lead to normalised states. Moreover, introducing we can rewrite the r.h.s. of (6) as The folding mapping can also be applied to time-evolved operators in Heisenberg picture. Considering again the example of local operators we have where we introduced the "doubled gate" (or the Heisenberg operator gate) and denoted by (·) T the matrix transposition. Note that, on the level of the doubled gate, the relations (5) read as namely, they impose a unitality condition on W ensuring that it maps the "identity state" | into itself. In this setting we consider dynamical correlation functions in the infinite-temperature state, which in the folded picture are represented as Whenever t < 2L − |x + 2y|, namely when we are effectively in the thermodynamic limit, non-trivial correlations are contained in a causal light cone. This is readily seen by using the rules (11), which allow us to rewrite (12) as follows (13) Here we conveniently rotated the picture by 45 • clockwise and, for concreteness, we depicted correlation functions for integer coordinates y ∈ Z and x ∈ Z. The cases with half-integer coordinates can be treated in an analogous fashion: if y ∈ Z + 1/2 the states at the bottom ( and ) are exchanged, while when x + y ∈ Z + 1/2 those at the top ( and ) are exchanged.
We see that in this representation correlation functions correspond to partition functions of a statistical mechanical model (with complex weights determined by the tensor W ) defined on a rectangular lattice of dimensions x where we introduced the ceiling function · , such that x ∈ Z and x ≤ x < x + 1 for any x ∈ R (this definition applies also for negative x). Moreover, we note that for values of x such that x ± ≤ 0 the correlations vanish identically. For definiteness, from now on we will always consider y = 0 unless otherwise stated.
In the representation (13), it is natural to think of the correlations in terms of horizontal and vertical transfer matrices Note that the above transfer matrices fulfil the following two properties i) A ab x and C ab x are contracting, i.e. their eigenvalues lie within the unit circle in the complex plane, for all a, b. This is a consequence of the unitarity of W and can be established following the derivation in Appendix A of Ref. [13].
ii) The state | ⊗x is an eigenvector of A x and C x with eigenvalue one. This is a direct consequence of the unitality relations (11).
The folding mapping described in this subsection turns the evolution of operators in the quantum circuit defined by the elementary gate U into that of states in a larger (super) quantum circuit defined by the elementary gate W . In this language the correlation functions are nothing but matrix elements of powers of the evolution operator between two specific states. In particular those of local operators are matrix elements of W t between two "one particle states" composed by the tensor product of 2L−1 copies of | and one state orthogonal to it (|a and |b in (12)). Matrix elements of this kind can be brought to the form (13). Even though the gate W is unitary by construction (cf. (10)), the unitarity of W is not needed for the simplification: one just needs the unitality property (11). Therefore, this setting can be used to study more general problems than that of computing correlations in unitary quantum circuits. An example of it is given in Sec. V 2 where we use it to study correlations in a class of classical Markov chains.

B. Dual-unitary Gates
As observed in Ref. [11], the correlations drastically simplify whenever the gate, together with (11), also fulfils We use a different (orange) color to denote doubled gates fulfilling these two additional conditions. The conditions (20) originate from the following requirements on the single (unfolded) gates and essentially imply that the evolution of the system remains unitary when one exchanges the roles of space and time. Gates with these properties, called "dual-unitary" gates, have recently been used to obtain exact results in a number of different problems concerning non-equilibrium dynamics of quantum many-body systems and quantum many-body chaos [11][12][13][14][15][16][17][18][19]. The simplification of dynamical correlations can be seen, for instance, by applying the first of (20) to the rightmost corner of (13) We see that the relation can be applied further and ulti-mately allows to bring the diagram in the following form b x |a 0 (t) du = b a . (23) By applying the second of (20) we see that the correlation factorises, and reduces to These two terms are zero for a and b orthogonal to the identity operator, i.e. if they are traceless. Following this derivation it is easy to see that the only non trivial correlations are obtained in diagrams with no corners with two neighbouring . Namely, when the rectangle in (13) reduces to a line with two bullets at the ends. For integer starting points, y ∈ Z, this happens when In this case the diagram reads as where in the last step we used the transfer matrix (16) for a simple one-dimensional system (x = 1). Moreover, we added the subscript "du" to stress that A du,x is made with dual-unitary gates. Analogously, when the starting point is a half oddinteger, say y = 1/2, the correlation is non-zero only for and reads as In summary, in dual-unitary circuits the correlations are entirely determined by 1d transfer matrices (or equivalently 1-qudit maps) and take the following simple form b x+y |a y (t) du =δ x+t b|(C du,1 ) 2t |a where δ x denotes the Kronecker delta function. As discussed in Ref. [11], depending on the spectrum of the transfer matrices A du,1 and C du,1 , these correlations can show four increasing degrees of ergodicity ranging from non-interacting behaviour -where correlations are all constant -to the ergodic and mixing one -where all correlations decay exponentially. In particular, by providing a complete parametrisation of dual-unitary gates for d = 2, Ref. [11] showed that the ergodic and mixing case is typical (i.e. it has measure one in the parametrisation space).
Finally, we stress that if the double gate is defined as in (10) and U fulfils (21), then also W and W † fulfil an analogous diagrammatic relation. This is however not needed to obtain the results of this subsection. We only need the conditions (20), which we dub "dual-unitality". Even though when the gate W comes from a folded quantum circuit the two conditions are equivalent, (20) is less restrictive and can hold in a more general setting (cf. Sec. V 2).

III. STRATEGY AND RESULTS
The goal of this paper is to develop a perturbative expansion of correlation functions around the dual-unitary point. The idea is to consider circuits with a number of non-dual-unitary gates U ε composed of a dual-unitary term U du and a non-dual-unitary correction. To have gates that are manifestly unitary we consider perturbations of the form where V is a generic hermitian 2-qudit operator and nonnegative real parameter ε sets the strength of the perturbation. The folded gate W ε can then be written as When representing correlation functions as partition functions in a lattice of doubled gates (cf. (13)) one can think of the gates W ε as defects. To better control the perturbation theory it is useful to also modulate the number of defects in the lattice. To this aim we introduce an additional parameter: the density δ of defects. We fixed the density because the actual arrangement of the defects does not affect appreciably the physics: one can imagine to randomly place δx + x − defects among the x + x − gates in the lattice (13). For simplicity, however, sometimes it will be useful to imagine the defects covering a regular sublattice of (13) with lattice spacings ν + and ν − , It follows directly from the above definitions that in both limiting cases, ε = 0 and δ = 0, we recover a dual-unitary circuit. The two perturbations are highly inequivalent.
In particular, the case of small density δ 1 is substantially easier to treat than that of small strength ε 1 and allows for rigorous results. This difference can be appreciated through a simple combinatoric argument: While for small δ one can work with a single partition function with a small number of defects and large dual-unitary islands, the expansion in ε of (31) generates a complicated sum of terms. In particular, the number of contributions at a given order ε n corresponds to the number of ways to dispose n identical objects in x + x − identical drawers and becomes exponentially large in the volume for large enough n.
Remarkably, in this paper we find that -under certain conditions on the unperturbed dual-unitary gate W duthe leading order contribution to the correlations can be directly computed in both cases and, surprisingly, takes the same form. Specifically, we observe that -at the leading order in time -correlations are still determined by the 1d transfer matrices A du,1 and C du,1 (cf. (26)- (28)). The difference is that, instead of being contracted along straight lines as in (26) and (28), now the maps can also be contracted along zig-zag lines like b a , (33) which we dub "skeleton diagrams". In particular, the correlation between two arbitrary points in the causal light cone is obtained by summing the contributions of all skeleton diagrams connecting the two points. The turns in the diagrams are generated by the defects, this means that for δ < 1 all the possible positions of the turns are restricted to a sub-lattice, while for δ = 1 the turns can be anywhere in the lattice. Note that all skeleton diagrams with down-or left-turns are forbidden. Indeed, these diagrams are cut by the rules (11) and do not contribute to the correlation. Such decomposition of the correlation function can be interpreted as a discrete path-integral on the 2d lattice (32).
More specifically, here we study this problem for three increasing levels of difficulty. We start by considering random defects, then we move on to fixed defects on "reduced gates" which are effectively 4 × 4 matrices (see below), and at the end we consider the general case.

A. Results on Random Defects
We begin by considering the simplest possible scenario: uncorrelated random defects. We prove that GUErandom defects can produce no turns, but, at the same time, maintain non-trivial correlations. Averaging over random GUE defects, i.e. where V in (30) is a GUE matrix with unit variance of matrix elements, produces correlations that are still in dual-unitary form. The only effect of the defects is an additional damping factor that causes or enhances -depending on the degree of ergodicity of the unperturbed dual-unitary circuit (see the discussion at the end of Sec. II B) -the exponential decay of the correlations along the light cone edge (light ray). In Sec. IV we find (36)). The function returns values in where the transfer matrices A du,1 and C du,1 are defined as in (16) and (17) but in terms of the unperturbed dualunitary gate W du , while we introduced the rescaled variablesx and the function Finally, the symbol · in Eq. (35) denotes the floor function, such that Fig. 1), the dual-unitary form of the correlations is stable against GUE-random defects for any strength ε and density δ. Note that, instead, Haar random defects, i.e. replacing U ε in (30) by a Haar random unitary matrix, produce f (x) = −1 for all x, i.e., they destroy the correlations (i.e. correlation function is a Kronecker delta in spacetime) for any positive density δ > 0 (see Sec IV).

B. Results on Reduced Gates
Secondly, we study folded gates that are effectively 4-dimensional (rather than 16-dimensional which is the minimal achievable dimension for W defined as in (10)). More precisely, we consider the case where each wire in (32) is effectively a two-state system with a basis where | is orthogonal to | . In the tensor product of two bases (37) the reduced two-body gate reads as where we used thin lines to highlight that the wires are now two-dimensional. As explained in Sec. V, this situation can originate from averaging the doubled gates (10) with respect to a single-site Haar U (1) measure (e.g. using local random magnetic field in the z-direction), or  Table I (gate  2). We see a reminiscence of the dual-unitary behaviour with a strong peak along the light-cone edge. As opposed to the pure dual-unitary result, however, the correlations are nonzero within the whole light cone.
from the interpretation of (12) in terms of a classical Markov chain. These two physically very different interpretations generate two distinct parameterisations for the elements of w, which are explicitly reported in Appendix B. Interestingly, we find that in both cases the gate becomes a bistochastic matrix (see the definition (66)) when conjugated with H ⊗ H, where is the Hadamard transformation. Most importantly, we see that (38) -even though not unitary -fulfils the conditions (11). Moreover, we see that the gate also fulfils the conditions (20) if we set ε 1 = ε 2 = 0, namely For reduced gates all non-trivial correlations are proportional to x | 0 (t) and we find the following three main results:

Exactly solvable cases
Apart from the dual-unitary point, ε 1 = ε 2 = 0, the gate (38) has four additional non-trivial exactly solvable points: where the elements that are not set to zero can take arbitrary values. As shown in Sec. V B in all these cases the correlation functions are exactly given by the sum of skeleton diagrams. In particular, considering a regular sub-lattice of defects as in Eq. (32) we find and the rescaled light cone coordinatesx ± are defined in Eq. (35) [34]. An example of the correlation pattern described by (41) is depicted in Fig. 2.
We stress that, generically, the gates fulfilling either of (i)-(iv) generate highly complex dynamics, e.g., they efficiently scramble quantum information. For example, we numerically computed a standard dynamical chaos indicator for locally interacting systems -the so called local-operator entanglement [13,14,31,[35][36][37][38][39][40] observing a linear growth. The key property leading to the simple form (41) is that the dynamics generated by the gates (i)-(iv) are not time-reversal symmetric. In particular, it is true that evolving forward (backward) in time the support of local operators can grow, forming larger and larger strings of local operator products. Yet these large strings cannot shrink back and do not contribute to the overlap with ultra-local operators. The fact that very large strings do not contribute much to the correlations of local operators is expected to hold quite generally. For example, a similar idea has recently been invoked in Ref. [5] to devise a numerical method able to access the late-time regime. The key point is that this property becomes exact in the cases (i)-(iv).

Low density limit
Two simple conditions on the spectrum of horizontal and vertical transfer matrices constructed with the reduced dual-unitary gate (40) (see Sec. V C 1) allow us to prove that the expression (41) is the dominant contribution to the correlation function at low density δ and for any ε 1 , ε 2 . More precisely, if the horizontal (vertical) transfer matrix fulfils the aforementioned conditions, (41) is dominant in the limit whereν +(−) is the minimal separation among the defects in the horizontal (vertical) direction, and the relative error decays exponentially inν +(−) . In Sec. V C 1 we prove that these conditions hold if the parameters of the gate w 0 (cf. (40)) fulfil where α and β denote the largest singular values of the sub-matrices c e d g and a f b g .
3. Small strength at density one Based on the rigorous results described at Sec. III B 2 we argue that skeleton diagrams give the dominant contribution also for δ = 1 and ε 1 ∝ ε 2 1 if both conditions (43) are satisfied. As explained in Sec. V C 2, the main idea is that, even though (41) does not provide a complete perturbative expansion in ε 1 (or ε 2 ), at each fixed order in perturbation theory the skeleton diagrams dominate for x + and x − both large. This argument is tested in Sec. V C 2 by comparing (41) with the exact correlations evaluated numerically: the observed agreement is excellent.

C. Results on Generic Gates
The contribution of skeleton diagrams in the generic case (considering again a regular sub-lattice of defects) reads as the "defect-maps" and the primed sums are subject to the constraints The expression (45) is a direct generalisation of (41), where one replaces numbers a, c, ε 1 , The conditions on the spectrum of horizontal and vertical transfer matrices mentioned in Sec. III B 2 are sufficient to rigorously prove the dominance of skeleton contributions at low densities also in the generic case (where the matrices are constructed with non-reduced dual-unitary gates). Although we cannot analytically determine the family of gates W du for which the conditions are fulfilled, we numerically identify such a family for quantum circuits of qubits (d = 2) in Sec. VI. For gates in this family the argument discussed in Sec. III B 3 remains valid as well. Accordingly, we find that if x + and x − are both large, (45) gives the most relevant contribution to correlations also for δ = 1 and ε 1.

IV. RANDOM DEFECTS
To analytically treat the problem outlined in Sec. III we first make a very drastic simplification and consider correlations averaged over random defects. Specifically, we take V in Eq. (30) to be random matrices independently distributed at each space-time point according to the Gaussian Unitary Ensemble, i.e. parameter ε and all physical quantities depend on ε and σ only through ε/σ.
Under these assumptions the average of operator gate W ε reads as with f (x) defined in (36). Eq. (50) is determined in Appendix A employing standard group-theoretic arguments.
Since both terms on the r.h.s. of (50) fulfil the "dualunitality" conditions (20) and their coefficients sum to one, the averaged gate fulfils (20)  Here we distributed the defects along the regular sublattice of Eq. (32) for a more convenient representation. The treatment of more general (e.g. random) dispositions of defects is completely analogous.
Expanding each averaged gate as in (50) we produce a sum of 2x ± terms (cf. (35)). As can be directly verified using (11), only the term with no | | gives a nontrivial contribution. Considering that term we directly find (34).
We conclude by observing that the GUE average considered here is very different from the flat average over the Haar measure on U (d 2 ). Taking a Haar random perturbation R in (30) instead of e iεV immediately gives where in the second step we used the right multiplication invariance of the Haar measure, and the last step is a trivial instance of the known result on integrals over unitary groups [41]. A direct consequence of this result is that all correlations are trivialised by the Haar average for any density δ of defects. Interestingly, (50) does not reproduce (53) for any value of ε (and hence for any value of the variance of the GUE distribution).

V. FIXED DEFECTS IN REDUCED GATES
In the previous section we proved that dual-unitary physics is stable under random perturbations distributed according to the GUE ensemble (for any density). In essence, this stability is due to the trivial form of the "mean defect" -which is in the GUE case proportional to the identity -and the fact that the probability measure is unbiased. In this case the mean defect does not effect the dynamics of two-point correlations and all deviations from the dual-unitary dynamics average out, resulting in a simple rescaling of the unperturbed gate.
Here we want to move on to defects with a nontrivial "mean". To study this very difficult problem, and yet make some analytical progress, we consider a minimal setting, where the wires in (32) are effectively twodimensional. This can be achieved by considering two seemingly very different physical problems. Even though, as we will see, these two problems lead to the same mathematical formulation, they come from very different physical contexts and it is interesting to consider them both.

U (1)-Random Unitaries
Let us begin considering the quantum circuit defined in Sec. II in the following special case.
where U is a generic U (4) matrix, σ z is a Pauli matrix, and the phases φ 1 , φ 2 , φ 3 , φ 4 are independent random variables uniformly distributed in [0, 2π]. Once again, the random variables {φ j } at each gate/space-time point are considered independent.
(ii) Consider two-point correlations averaged with respect to all {φ j }. This corresponds to averaging over the singlesite Haar measure of the U (1) group. On the level of the folded gate W (cf. (10)) the average corresponds to a simple projection, i.e. where In other words, the average over the Haar U (1) measure projects each wire onto the subspace spanned by the diagonal matrices {| ≡ |1 , | ≡ |σ z } (cf. (7)), and, from now on, we conveniently use this reduced basis notation. This means that the averaged 2-qubit folded gates are effectively 4 × 4 matrices (while they are 16 × 16 without the average).
We see that, even though the U (1) average simplifies the problem by reducing the size of the relevant local Hilbert space, it does not completely trivialise the twofolded gate and in turn the two-point correlations. This is in contrast to most averages over random unitary gates that have been considered over the last few years [22-24, 31, 32, 42-61]: in the latter cases the average of the doubled gate W is trivial (as in (53)) and so are twopoint functions. The minimal folded gate with a nontrivial average is the one composed of four copies of the time evolution operator (two copies of U and two of its hermitian conjugate U † ).
To be more specific, let us count the number of free parameters in the gate (55). This can be conveniently done using the following parametrisation for the generic unitary U (cf. (54)) [62,63] Here φ ∈ R, and u 1 , u 2 , u 3 , u 4 are elements of SU (2) in the fundamental representation. They can be expressed in terms of Euler angles From this parametrisation it is easy to see that the averaged folded gate (55) depends on 11 angles: two Euler angles (β j and γ j ) for each single-wire unitary u j , and the parameters J 1 , J 2 , J 3 . The explicit form of the averaged gate in the basis (37) is given by Eq. (38), where the explicit parametrisation of the coefficients ε 1 , . . . , g in terms of the above angles is given in Appendix B.
It is important to stress that the averaging procedure described above can be carried on for any fixed value of the angles {β j , γ j , J j }. In particular, setting J 1 = J 2 = π/4 one averages over gates in the dual-unitary subclass [11]. In this specific case, the averaged gate, now depending on 7 parameters (see Appendix B), takes the form (40). Comparing with (38) we see that the effect of requiring dual unitarity is to set ε 1 and ε 2 to zero. Note that after the averaging the gate continues to fulfil the dual-unitality conditions (20).
Before turning to the analysis of correlation functions, it is interesting to consider an alternative classical stochastic interpretation for the gates (38) and (40).

Classical local Markov chains
The averaged gates introduced in the previous subsection admit an alternative interpretation in terms of classical stochastic processes. To see this, let us consider (12) as the fundamental object, forgetting its origin in terms of a folded quantum circuit and focus on the following setting: (i) Each wire in the tensor network diagram (12) has generic dimension N ∈ N, not restricted to squares of positive integers. In particular, we choose a certain basis and interpret each state as a possible state of a classical spin (or any abstract configuration of a classical system). Thus, we now view our physical system as a chain of 2L classical spins. The probability distributions over configurations of such chain can be formally expanded in the product basis Namely, we can write where the coefficient p(t; {α j }) gives the probability that the system is in the configuration α 1 , . . . , α 2L at time t, and hence (ii) The state | has the following expansion in the basis (60) Namely, | is the flat sum of all possible configurations of a single spin and, apart from a normalisation factor, it represents the flat probability distribution. In the context of classical stochastic processes, such state is known as the maximal entropy state and is typically denoted by |ω . Therefore, restoring the correct normalisation we have (iii) The local two-body gate W is not unitary but bistochastic in the tensor product of two bases (60). Specifically, Therefore, W in (19) is now a Markov chain propagator evolving probability distributions of the configurations of a chain of 2L classical spins (or any other discrete variable degrees of freedom). More precisely, one can define as the time-evolved probability distribution.
The bistochastic property (66) of the elementary local gate W implies that the maximal entropy state is stationary: W |ω ⊗ |ω = |ω ⊗ |ω . This means that the unitality relations (11) are satisfied even if the gate is not unitary and we recover the natural light-cone causal structure leading to (13) in the thermodynamic limit.
In this language, (12) is the correlation function in the maximal entropy state of two (diagonal) local observables such that where |a = α A α |α and |b = β B β |β are the states appearing in (12). In other words, we can rewrite the r.h.s. of (12) as b x+y |a y (t) = N 2L ω . . . ω|B x+y W t A y |ω . . . ω , (70) which is the expression for correlation functions in classical Markov chains and classical cellular automata [64][65][66][67].
Moreover, defining the "dual" local Markov gateW by and requiring it to be bistochastic in the basis {|αβ ≡ |α ⊗ |β }, the dual-unitality relations (20) are also satisfied. This immediately implies that -if both W andW are bistochastic -the correlations take the simple form (29). Indeed, the simplification of the correlation functions is only based on the diagrammatic relations (11) and (20) without utilising the unitarity of the gates. We will refer to bistochastic gates W du with this special property as dual bistochastic. Hereby, we can now establish a direct connection with the previous subsection. Consider the minimal case where each classical spin takes only N = 2 values and define | as the state orthogonal to | , i.e.
Using the second requirement in (66) it is immediate to verify that To sum up, one can think of the gates (38) and (40) in two different ways: (i) as (dual-)unitary double gates averaged according to the single site Haar U(1) measure, or (ii) as (dual-)bistochastic gates conjugated with H ⊗ H, i.e., written in Hadamard transformed basis. Even though both procedures yield gates of the form (38) and (40), the matrix elements ε 1 , . . . , g are parametrised differently. They can in principle generate two different families of gates with unclear relation. Using the explicit parameterisations reported in Appendix B we numerically established that the families produced by both strategies, (i) and (ii), are equivalent.

A. Correlation functions as sums over paths
Let us now examine the r.h.s. of Eq. (12) in such minimal setting: two-dimensional wires and gates of the form (38). Inserting a resolution of the identity | at each wire, we can explicitly decompose (12) into the sum of 2 x+x− terms. Each term admits a simple interpretation as a (weighted) path on the two-dimensional space, or, analogously, as the spatial configuration of a certain polymer. To see this, let us introduce a different diagrammatic representation of the weights in terms of "tiles" where we connect s with solid lines and ignore s. For example, The complete set of tiles (corresponding to non-zero coefficients of the gate (38)) is given by Then, the correlation function is expressed as We used that all non-trivial states (i.e. those orthogonal to the "identity" | ) are proportional to | . Namely, for all a, b such that tr[a] = tr[b] = 0. Moreover, in writing (75) we considered the case of x, y ∈ Z. The other three possibilities correspond to different positions of the initial and final lines. Specifically, for y ∈ Z + 1/2 the initial line enters the bottom left tile from below and for x + y ∈ Z + 1/2 the final line exits the top right tile from above. As suggested by the diagrammatic representation in Eq. (75), the correlation function can be thought of as the sum of paths with fixed end points, where different configurations have different (possibly negative) associated weights. The paths are allowed to split and merge, with weights e, f and b, d, respectively, but they cannot "jump": dangling ends of lines inside the diagram are forbidden. Finally, since downward and leftward turns are forbidden (cf. (74)), the paths can form loops only if they split (see, e.g., the last two diagrams on the r.h.s. of (75)).

B. Exactly solvable cases
By inspecting the set of allowed tiles (74), we identify four non-trivial solvable cases, which we will analyse in the following two subsections.

Almost dual-unitary cases
As observed before, the case ε 1 = ε 2 = 0 corresponds to the dual-unitary point and all correlations propagate along straight paths, with no turns allowed. Depending on the initial conditions (i.e. on whether y is integer or half-odd integer) the straight lines are going either upward or rightward. This is simply a path-sum reformulation of the general dual-unitary result (29) (cf. Ref. [11]).
Remarkably, however, the correlations remain exactly solvable also when only one of ε 1 or ε 2 vanishes. Indeed, in this case we only allow for paths with a single turn and the rules (74) do not permit any "dressing", i.e. any thickening of the lines due to loops (see Sec. V C 1 for a more precise definition). For example, choosing ε 1 = 0 and ε 2 = 0, the correlations take the form given in (41), where only the first terms in each of the two lines can be non-zero (depicted in (99)). We see that, apart from the straight paths described above, (41) establishes nontrivial correlations between integer and half-odd integer points.

Non-dual-unitary solvable cases
Surprisingly, the tiles (74) contain another non-trivial solvable limit that is not dual-unitary. Indeed, there can be no loops in the paths whenever the "split weights" f and e or the "merge weights" b and d are zero. Consequently, when (e, f ) = (0, 0) or (b, d) = (0, 0) the only allowed paths are the "skeleton diagrams" described in Sec. III (see, e.g., the first two diagrams on the r.h.s. of Eq. (75)). In this case we can directly evaluate the correlation function (75) by summing all skeleton diagrams. In particular, if the defects cover a regular sub-lattice as in (32), one obtains the expression (41).
The formula (41) can be derived by straightforward combinatorics. Let us detail its derivation considering the case x ∈ Z as an example. We begin by noting that, since the path goes from the bottom left corner to the upper right one, it must involve the same number n of right-and up-turns, with weights respectively given by ε 1 and ε 2 . The minimum number of such turns is 0 (it contributes only for x − = 1), while the maximum, min(x + ,x − − 1), is set by the size of the defect's sub-lattice (the rescaled light cone coordinatesx ± are defined in Eq. (35)). For each fixed n, one has n turns, x + − n horizontal segments, contributing with a factor a x+−n , and x − − 1 − n vertical ones, contributing with a factor c x−−1−n . To count all the possible ways in which the elementary pieces can be combined, we can consider the horizontal and the vertical directions separately. In the horizontal direction we need to distribute n indistinguishable pairs of turns (first up and then right) inx + positions, leading to the combinatorial factor In the vertical direction we are instead more constrained. Indeed, the first and last turns must be in the first and last row, respectively. The other n − 1 pairs can be distributed freely in the remainingx − − 2 positions, leading to Putting it all together, we obtain the desired result. For defects on irregular sub-lattices the reasoning is similar but one has different combinatorial coefficients depending on the actual shape of the sub-lattice.

C. Perturbation theory around the dual-unitary point
Let us now consider the case of circuits that are perturbed away from the dual-unitary point. As discussed in Sec. III we use two parameters to control the perturbations: strength ε and density δ. We begin by considering the case of perturbations in the density of defects, which allows for a more rigorous analysis. Later we will see that, surprisingly, most of the rigorous conclusions drawn in that case apply also for small ε and arbitrary δ ≤ 1.

Low density, unit strength
We begin our analysis by focusing on defects with strength ε = 1 placed on a regular sublattice as in Eq. (32). In this case there are regular strips (vertical and horizontal) composed only of dual-unitary gates. Whenever the widths ν ± (cf. (32)) of these strips become large enough, we can simplify the contribution by considering only the leading eigenvectors of the strips' transfer matrices a du,x and c du,x . These are defined as in Eqs. (16)-(17) but using the 4 × 4 dual-unitary (or rather, dual bistochastic) gate w du (cf. (40)), namely To treat these matrices we make use of the following rigorous result, proved in Appendix C.
Property 2. The matrices a du,x , c du,x take the following block diagonal form where we defined and the "reminders" r 1,x , r 2,x are non-zero only in the subspaces defined by the projectors 1 − x j=0 p x,j . Moreover, and analogously where α, β ∈ [0, 1] are respectively the operator norms (largest singular values) of the matrices Consider now a vertical strip x − × ν + with horizontal transfer matrix a du,x− . Property 2 guarantees that for we can make the replacement with an error bounded in operator norm r 1,x− by a 2 + |bf |/(1 − α). The replacement (88) makes the calculation of correlations extremely easy. Let us illustrate this considering the diagram on the r.h.s. of (32), which becomes where in the second step we used the unitality relations (11) to contract the vertical lines. Using the explicit form of the one-dimensional transfer matrices and of the ones with orange gates that are obtained by setting ε 1 = ε 2 = 0 in the above, we see that (89) is nothing but the sum (41) of skeleton diagrams. Taking into account the bound on the norm of r 1,x , we find that the error associated with the replacement (88) is bounded by so that (89) becomes exact in the limit x + , ν + → ∞ with x + fixed. This bound is obtained by replacing just one of the vertical dual-unitary strips with the remainder r 1,x . Specifically, considering the relative error where x | 0 (t) is the exact result and x | 0 (t) sk is the skeleton expression (Eq. (41)), we get Thus, (41) gives a good approximation to the full result when ν + is large enough. A completely analogous reasoning holds for horizontal strips ν − × x + , whenever and ν − is large enough. In particular, for x − , ν − → ∞ withx − fixed we again find an exact statement. We stress that in the above argument we never used the fact that the defects are disposed along a regular sublattice: we just used that their minimal separation,ν ± , in one of the two directions becomes large. Provided that the latter condition applies, the correlations are sums of skeleton diagrams. In the case of irregular disposition of defects, however, one has to sum only over the skeleton diagrams connecting the (irregular) subset of lattice sites containing defects. Hence, the combinatorial factors differ from those in (41).
Finally, we note that all this admits a simple interpretation in terms of the paths introduced in the previous subsection. The different eigenvectors of the transfer matrix can be seen as increasingly thicker horizontal lines of length 1. Their contributions -the eigenvalues of a du,x -are obtained by complicated combinations of the tiles (74), which can be interpreted as a complicated "dressing". In this language, Property 2 implies that if (87) holds, the dominant weight is carried by the "bare" propagator . (96)

Unit density, small strength
Now let us consider a different perturbative limit: small strengths ε 1 and fixed density. In particular, for definiteness, we focus on the case of unit density δ = 1, which is the most interesting from the physical point of view. Namely, we take all the gates in the diagram (32) to be green and set with the matrix elements such that, when the gate is conjugated with H ⊗H (cf. (39)), it becomes bistochastic. The advantage of the parametrisation (97) is that in this way -in the language of Sec. V A -parameter ε "counts the turns". Namely, each fixed order in perturbation theory is determined by the sum of all allowed paths with a fixed number of turns. Even though this is useful for classifying possible terms of different perturbative orders, the explicit evaluation of each order appears far from trivial.
Firstly, consider the simplest case: the leading order that is proportional to ε, when the x ∈ Z + 1 2 , y ∈ Z. It is given by a single skeleton diagram: , and there are no other diagrams of the same order. Therefore, the relative error for ε → 0 goes as O(ε), even for small x − or x + . In contrast, the first non-trivial order in the case x, y ∈ Z and x + x − > 1 is of order ε 2 and has much more complex structure. A simple set of contributions to this order is given by skeleton diagrams of the form . (100) Still, many more diagrams contribute at the same order in perturbation theory. For example, we have , , .
The crucial observation at this point is that there exists a regime for the parameters of the gate w du , where all these "complicated" diagrams give negligible contribution when both x − and x + become large. Specifically, this happens when both of the conditions discussed in the previous section -(87) and (95) -hold. This observation can be understood as follows. First, we note that the diagrams in (101) can be thought of as skeleton diagrams with "dressed" horizontal one-particle propagator. Then, we observe that, as we are working at a fixed order in perturbation theory, these dressed propagators are composed of dual-unitary tiles only. We can then make use of Property 2 to bound their contribution by where x 1 is the length of the corresponding segment of the dressed propagator. We see that for large enough x this contribution is exponentially suppressed with respect to the "bare" line (∝ a x ). If both conditions, (87) and (95), hold, this reasoning can be repeated at any fixed order in perturbation theory to show that the skeleton diagrams are leading when x + and x − are both large.
To check the above reasoning, we performed a direct numerical evaluation of the correlation function (13) and compared it with the prediction (41), obtained by summing all skeleton diagrams. The comparison is extremely encouraging: the deviations are typically undetectable on the scale of the plot, see Fig. 3 for a representative example.
Turning to a more quantitative analysis of the agreement, we considered the relative error (93) where x | 0 (t) is calculated numerically with no approximations while x | 0 (t) sk is calculated using (41). The results are reported in Figs. 4 -6. Specifically, Fig. 4 concerns the case of large x + = x − . We see that when the "unperturbed gate" w du fulfils the conditions (87) and (95), the relative error is always very small and appears to vanish with ε. Interestingly, even if the predictions (41) for x ∈ Z and x ∈ Z + 1/2 are of different orders in ε, we observe almost the same relative errors in the two cases. Lastly, an important point highlighted by Fig. 4 is that the relative error is of order one for any ε when the conditions (87) and (95) are violated.
Our argument above relied on the fact that x − and x + are both large. When one of the two, say x − , is fixed, we expect the relative error to be O(ε 0 ) in the case x, y ∈ Z (O(ε) for the case of integer and half-integer endpoints), and bounded by (102) with y = x − . For small x − this can be directly verified by computing the exact correlations through Eq. (18). For example, in the case of x − = 2 and x ∈ Z we find This result shows that the leading correction is a "dressing" of the skeleton diagram in the "short direction" and gives This point is confirmed by our numerical results, as illustrated in Fig. 5. For increasing values of x − the order zero contribution becomes increasingly small and is eventually dominated by the O(ε) contributions. Finally, an interesting question concerns how many orders in ε one has to keep in the skeleton-diagram expansion (41) in order to get an accurate result. This question is considered in Fig. 6, where we compare the relative errors obtained by (i) keeping all orders in the expansion (41), (ii) keeping only the leading order in ε.  93)) for x− x+ as a function of ε. The gate is fixed to gate 2 from Table I   We see that the prediction (i) is much more accurate for larger ε, even though the two predictions coincide for small enough ε.

VI. FIXED DEFECTS IN GENERIC GATES
In this section we use a combination of numerical and analytical arguments to show that the main mechanisms observed in the minimal setting of the previous section carry over to the case of generic clean quantum gates. Under certain conditions on the unperturbed gate W du the correlations computed by summing all skeleton diagrams still agree strikingly well with the exact numerical results. To simplify our numerical studies we consider the case of qubits (d = 2), where the double gate W is a 16 × 16 matrix.
We begin by briefly describing how to compute the contribution (45) of all skeleton diagrams for the case of defects on regular sub-lattices. Firstly, we construct a single skeleton diagram by drawing a zig-zag line connecting the two operators on the lattice (32), and multiply the appropriate one-dimensional transfer matrices along the line (cf. (33)). The total contribution is obtained by summing the contributions of all possible paths on the lattice connecting the two endpoints and containing no left and down turns. In particular, l + j (l − j ) in (45) are of the lengths of each horizontal (vertical) segment and the constraints (48) come from the simple requirement that the sum of all horizontal (vertical) segments is equal to x + (x − ) minus the number of turns. Note that, because of the non-commutativity of the matrix product, Eq. (45) is considerably harder to evaluate than Eq. (41). The problem of its evaluation is addressed in Appendix D.
Following the analysis of the previous section we now consider small densities of defects and investigate the spectra of horizontal and vertical transfer matrices, A du,x and C du,x , composed only of dual-unitary gates W du,ε . We remind the reader that, as discussed in Sec. II A, these transfer matrices have the spectrum contained in the unit circle and a "trivial" eigenvector -| ⊗x -corresponding to the eigenvalue 1. This state, however, does not affect the correlation functions. One can see this by considering (18) and observing that | ⊗x has zero overlap with the states |a . . . , . . , C a x | . . . . The relevant quantity for correlations is then the next-to-leading eigenvalue. In analogy with the discussion in Sec. V C 1, we conclude that if (i) The next-to-leading eigenvalue λ 1 of A du,x has eigenvectors with support one, i.e. it is of the form (ii) There is a finite gap between λ 1 and the rest of the spectrum; the skeleton diagrams give a good approximation to the correlation functions for large enoughν + (minimal separation of the defects in the horizontal direction). Relative errors are bounded by where λ 2 is the third leading eigenvalue of A du,x . Note that, in the language of the previous section, this still corresponds to saying that the bare propagator carries the dominant weight. Analogous statements hold for the vertical direction, if one replaces A du,x by C du,x andν + byν − . In this more complicated setting we are not able to prove a rigorous statement, such as Property 2, to determine the gates for which (i) and (ii) hold. Nevertheless, investigating the spectrum numerically, we find that a class of gates, for which the leading non-trivial eigenvectors are of the form (105), exists, see the left panel of Fig. 7. Moreover, focussing on this class, we isolated a subclass that shows a finite gap: an example of the probability distribution for the gap is reported in the right panel of Fig. 7 for three increasing values of x (i.e. the length of the transfer matrix). Finally, we remark that the left panel of Fig. 7 also shows a rapid drop in the probability that a certain eigenvector is leading with the size of its support. This suggests that if (i) or (ii) does not hold, the dominant contribution to the correlations is still given by skeleton diagrams. In general, however, these diagrams will have "thickened" lines corresponding to "dressed" propagators: the width of the line corresponds to the support of the leading eigenvector.
Finally, we consider δ = 1 and small ε. In complete analogy with the discussion in Sec. V C 2 we decompose the gate W ε in two parts as follows where W du,ε is dual-unitary and the only non-zero elements of δW ε (now 18) are the "turns" Note that both, the elements of W du,ε and those of δW ε , will in general depend on ε. The difference is that for small ε the former are O(ε 0 ), while the latter are O(ε).
Equipped with the definition (107), we are now in a position to repeat the argument of Sec. V C 2 by formally considering a perturbative expansion in the number of turns (which at the first order is equivalent to that in ε). We then conclude that, if the points (i) and (ii) above hold for both A du,x and C du,x , Eq. (45) gives the leading contribution at each fixed order in ε when x + and x − are both large. This is in agreement with our numerical findings: representative examples are reported in Figs. 8 and 9. In particular, we again see that for large x + and x − the relative error decreases with ε (cf. right panel of Fig. 8) and for fixed ε and large x + it decreases with x − (cf. Fig. 9). Specifically, in the latter case we expect the error to be bounded by (cf. (102) and (106)) This is confirmed by the numerical results in the figure 9, which show an even faster decay, suggesting that the bound is not tight. Finally, from the right panel of Fig. 8 we see that if the conditions (i) and (ii) are violated, the agreement immediately becomes very poor.

VII. CONCLUSIONS AND PERSPECTIVES
In this paper we studied correlation functions in perturbed dual-unitary circuits, or, equivalently, in perturbed dual-bistochastic Markov chains. We considered the problem for three increasingly more realistic (increasingly complex) settings -clean systems with random defects, noisy systems with fixed defects, and clean systems with fixed defects -and identified a class of dual-unitary circuits that is stable under perturbations. More precisely, when these systems are perturbed, their correlations continue to be given in terms of one-dimensional transfer matrices (single qudit maps) and hence to (generically) decay exponentially. These exponentially decaying correlations are continuous with respect to dual-bistochasticity-breaking perturbations demonstrating structural stability. The main qualitative change induced by the perturbations is that they allow the correlations to spread through the whole causal light cone and not just along its edge, as it occurs in the pure dual-unitary case. One can intuitively think of the perturbations as defects that permit the correlations "to turn" in spacetime. Remarkably, the quantitative features of the correlations are extremely well captured by a simple "path integral" formula, corresponding to the sum of correlations over all one-dimensional paths within the causal light cone that connect the endpoints. In conclusion, this paper presents, to our knowledge, the first set of theoretical results on the stability of the ergodic regime of quantum many-body systems under generic perturbations and opens the research arena of quantum manybody ergodic theory.
The case of noisy systems with fixed defects turned out to be an intriguing minimal model. It shows all the physical features of the generic setting but, simultaneously allows for rigorous derivations. This is not an isolated case (see, e.g., [15]): introducing a small degree of randomness to isolate minimal settings appears to be a fruitful route for accessing non-trivial information about interacting many-body systems. In the minimal setting we identified four additional classes of circuits -beside the dual-unitary ones -where correlation functions are exactly solvable -i.e. the aforementioned "path integral" formula applies exactly. Systems in these classes are generically strongly interacting and generate highly complex dynamics. An obvious direction for future research is to study them further, understanding, for example, their spectral statistics and their non-equilibrium dynamics. Moreover, it would be very interesting to understand whether these classes can be defined in generic quantum circuits.
Another outstanding question raised by our work con-cerns what happens when one perturbs dual-unitary circuits that are not in the stable class. One possibility, which seems to be hinted by our numerical results, is that the correlation functions continue to have a "path integral" form but the paths are "thickened". In other words they are computed in terms of n-dimensional transfer matrices rather than one-dimensional ones. If n does not scale with time -as we seem to observe in general -the physical picture remains very similar to the one studied here and, in particular, all correlations continue to decay exponentially. However, it would be interesting to understand whether there exists a class of circuits for which n grows with time. This could lead to a phase transition in the behaviour of correlations and therefore to richer physics. Importantly, this would allow for power-law decaying correlations signalling non-trivial conservation laws.
We can immediately propose two generalisations of our setup. First, instead of unitary quantum circuits, one can consider circuits of completely positive maps. These include the class of circuits with projective measurements that is currently attracting substantial attention as it displays measurement driven phase transitions [68][69][70][71]. Specifically, it is straightforward to generalise the concept of dual bistochasticity to dual quantum bistochasitity. Second, our treatment can be directly extended to perturbed dual-unitary (or dual-bistochastic) circuits in higher spatial dimensions where generic circuits display non-trivial complexity transitions [72]. In this case we again expect the correlations to be written as sums over one-dimensional paths.
The coefficients are readily computed by multiplying M by P sym and P anti−sym , respectively, and taking the trace. This leads to Expanding α in powers of ε we find The above expectation values can be calculated exactly and give Plugging back into (A5) we find Appendix B: Parametrisation of the gates (38) and (40) In this appendix we present an explicit parametrisation of the gates (38) and (40) depending on whether they are obtained as U (1)-averaged (dual-) unitary double gates or as conjugated (dual-) bistochastic ones.

U (1)-averaged (dual-) unitary double gates
To find a convenient parametrisation we note that the elements of the gate (38) can be computed by evaluating Here we definedŪ where u(α, β, γ) is defined in Eq. (59) and V [{J i }] in Eq. (58). Eq. (B1) follows directly from the parametrisation (57) and the form (55) of the averaged gate.
Thus, applying triangle inequality to (C3), we find where we used that a du,x = 1 (this can be proven by reasoning as in Appendix A of Ref. [13]). Furthermore, we denoted by α the operator norm of the 2 by 2 matrix in the first term on the r.h.s. of (C3). Explicitly, α can be computed by considering the largest singular value of the matrix and reads as α = c 2 + e 2 + g 2 + d 2 + (c 2 + e 2 − g 2 − d 2 ) 2 + 4(cd + eg) 2 2 . (C12) By means of the parameterisations (B19)-(B23) and (B31)-(B33), one can explicitly verify that α ∈ [0, 1]. Using (C11), we see that λ • x is bounded by y x , which is defined as the solution of y x = αy x−1 + |f |, and reads explicitly Hence, where in the last step we assumed α < 1. Combining everything we get Considering now the range of parameters for which we have that for λ 1,x−1 ≤ |a|, it follows λ 1,x ≤ |a|. In particular, since λ 1,1 = |a|, we conclude Moreover, from (C15) we conclude aa du,x−1 − ap x−1,0 + ba du,x−1 ≤ a 2 + |bf | 1 − α , contributions, which becomes eventually impractical for x ± , n 0 1. This problem can be overcome by (i) truncating the expansion (45) (this leads to accurate results for small enough ε) or (ii) by approximating the higher orders by substituting A 0,1 and C 0,1 with the projectors on their largest non-trivial eigenvalues, which we respectively denote by λ 1,a and λ 1,c . Thus we can rewrite the contribution of the higher orders in (45) exactly in the form Eq. (41), with the sums running from a certain n = n c (n = n c − 1) instead of 1 (0). The parameters are expressed as ε 1 = λ 1,c | E 1 |λ 1,a , ε 2 = λ 1,a | E 2 |λ 1,c , a = λ 1,a , c = λ 1,c , and there is an additional factor of b||λ 1,a λ 1,a |a ( b|λ 1,c |λ 1,a |a ) in front of it. If all l ± j 1 this replacement would give a negligible error, however, in the sum there are also terms with l ± j ≈ 0 and the projection is strictly speaking unjustified. Nevertheless, this gives a useful approximation, especially in the limit x ± → ∞.