Exact dynamics in dual-unitary quantum circuits

We consider the class of dual-unitary quantum circuits in 1 + 1 dimensions and introduce a notion of “solvable” matrix product states (MPSs), deﬁned by a speciﬁc condition which allows us to tackle their time evolution analytically. We provide a classiﬁcation of the latter, showing that they include certain MPSs of arbitrary bond dimension, and study analytically diﬀerent aspects of their dynamics. For these initial states, we show that while any subsystem of size (cid:96) reaches inﬁnite temperature after a time t ∝ (cid:96) , irrespective of the presence of conserved quantities, the light-cone of two-point correlation functions displays qualitatively diﬀerent features depending on the ergodicity of the quantum circuit, deﬁned by the behavior of inﬁnite-temperature dynamical correlation functions. Furthermore, we study the entanglement spreading from such solvable initial states, providing a closed formula for the time evolution of the entanglement entropy of a connected block. This generalizes recent results obtained in the context of the self-dual kicked Ising model. By comparison, we also consider a family of non-solvable initial mixed states depending on one real parameter β , which, as β is varied from zero to inﬁnity, interpolate between the inﬁnite temperature density matrix and arbitrary initial pure product states. We study analytically their dynamics for small values of β , and highlight the diﬀerences from the case of solvable MPSs.


I. Introduction
The extensive study of isolated quantum matter out of equilibrium carried out in the last two decades reminded us, once again, of how tremendously complex the quantum many-body dynamics can be [1][2][3] . Even though the past decade has witnessed the development of powerful numerical techniques based on matrix product states 4 (MPSs) that are able to determine, quite generally, the dynamics of quantum many-body systems in one spatial dimension [5][6][7][8][9][10] , these methods are usually limited to small or intermediate time scales 10 . This is due to the generic linear growth of the entanglement entropy 11 , which is a major obstacle for the MPS-representation of the time evolving state. For this reason, it is of great interest to find instances where the many-body dynamics can be solved exactly, allowing for an analytic study of interesting physical phenomena, such as the emergence of thermalization 12 .
Integrable models provide a natural arena to search for such solvable examples 13,14 and allowed for great progress in the characterization of large-time properties of manybody systems out of equilibrium [14][15][16][17] . The computation of the full dynamics, however, remains a challenge. In particular, while an impressive number of results have been derived in theories that can be mapped onto free fermionic systems (see Ref. [14] for a comprehensive review), only a few special cases exist where analytic computations could be done in the presence of genuine interactions [18][19][20][21][22][23][24][25][26][27][28][29] . Furthermore, integrable models are by their own nature very special, and can not be representative of generic systems, which are expected to display qualitatively different features.
Recently, increasing attention has been devoted to the class of random unitary quantum circuits, which provide an alternative theoretical laboratory for the study of the many-body dynamics [30][31][32][33][34][35][36][37][38][39][40][41][42][43][44] . The main appeal of these systems is that they represent minimally structured dynamical models where analytic results can be obtained beyond the realm of integrability. One typically considers a set of finite-dimensional Hilbert spaces sequentially updated by unitary gates, which are chosen randomly out of a suitable probability distribution. Analytic predictions for physical quantities are then obtained after averaging over disorder realizations [36][37][38][39][40][41] . This approach allows one to obtain exact results for quantities notoriously hard to compute, such as out-of-time-ordered correlators (OTOCs) 35,39 , operator-space entanglement entropies [45][46][47] and other measures of quantum chaos such as the so-called tripartite mutual information 48,49 . Of particular interest for our work are the settings (sometimes called local random unitary quantum circuits) where the unitary gates couple only neighboring sites, simulating the "local dynamics" of generic many-body quantum systems [34][35][36][37][38][39][40][41] . One can still wonder, however, whether the presence of disorder averages gives rise to qualitative differences when compared to clean, homogeneous systems.
In this respect, an interesting class of quantum circuits without disorder, called "dual-unitary", was recently introduced 50 , for which several dynamical features could be investigated analytically. These systems implement a dynamics in which at each time step the configuration of the system is updated by applying a product of identical unitary gates on neighboring sites. Furthermore, as a defining feature, these gates remain unitary under a reshuffling of their indices. Remarkably, this class was shown to include instances of unitary dynamics both integrable (i.e. with local conservation laws) and chaotic, providing a rare example where the differences between the two can be analyzed to a high degree of control.
Previous works on these circuits focused on the study of infinite-temperature dynamical two-point functions 50 , and on the growth of the operator-space entanglement entropy 51,52 . It is then natural to wonder whether the class of dual-unitary circuits can also provide solvable models for the quantum dynamics arising from given initial states. This question represents the main motivation for the present work.
For one particular dual-unitary circuit, corresponding to the self-dual point of the kicked Ising model [53][54][55][56][57] , this problem has been already partially addressed in Ref. [58], where it was shown that for two special families of initial product states the growth of bipartite entanglement entropy could be computed exactly. In this paper we show that there exists a much broader family of "solvable" initial states, for which the dynamics can be tackled analytically, irrespective of the choice of the dual-unitary gates. This class includes MPSs of arbitrary bond dimension, and allows for the exact computation of several quantities beyond the growth of bipartite entanglement entropy, including the spreading of two-point correlation functions, and the thermalization time of finite subsystems. In this work we will focus on the case of spatially homogeneous systems, with a "Floquet-like" time evolution in which the same product of unitary gates is applied periodically in time. We stress, however, that our constructions also work for circuits with explicit space-time modulations, and can be used to obtain analytic results also in that case.
The rest of this manuscript is organized as follows. In Sec. II we introduce the local quantum circuits considered in this work, while in Sec. III we define the class of solvable initial states with respect to the dual-unitary dynamics. These are classified in Sec. IV and V, while their time evolution is studied in Sec. VI. The dynamics arising from a set of non-solvable states is studied for comparison in Sec. VII, while our conclusions are reported in Sec. VIII. Finally, the most technical aspects of our study are consigned to two appendices.

II. The dual-unitary dynamics
We consider a chain of 2L sites, each one associated with a d-level system. The corresponding Hilbert space In the following, we will denote the corresponding basis vectors by |j , j = 1 , . . . d. We are interested in local quantum circuits, which implement a particular periodically driven unitary evolution. Specifically, given the initial state |Ψ 0 we study the quantum dynamics defined by where t ∈ N, while Here U j,k are two-site unitary operators acting on the product of local spaces h j ⊗ h k , and periodic boundary conditions are implemented in Eq. (3). We will mainly consider the case of infinite systems, namely we will study the above dynamics in the limit L → ∞. We focus on the special class of dual-unitary quantum circuits recently introduced in Ref. [50]. These circuits are defined as follows. Let U be a unitary gate, and definẽ U as the operator given by the following reshuffling of indices k| ⊗ |Ũ |i ⊗ |j = j| ⊗ |U |i ⊗ |k .
We say thatŨ is the dual gate of U . Then, dual-unitary circuits are defined by local unitary gates U such thatŨ is also unitary, namely This definition is reminiscent of those of perfect tensors and 2-unitary matrices, introduced in Refs. [59] and [60] respectively. The former are defined as tensors which are isometric operators for any bipartition of their indices into two subsets (not necessarily of the same size). The latter are matrices, acting on the tensor product of two qudits, that remain unitary for any bipartition of indices into two pairs. We note, however, that dual-unitary gates are generally neither 2-unitary matrices nor perfect tensors, since they only correspond to unitary operators for two special bipartitions of the indices. Finally, we mention that the partitioning of indices in the definition of dual-unitarity is the same appearing in the one of the cross-norm and realignment criteria for separability [61][62][63] . The dual-unitary condition can be naturally expressed using the customary tensor-network graphical representation of quantum circuits. In this language, matrix elements of local operators are depicted as boxes with a number of incoming and outgoing legs. To each leg corresponds an index associated with one of the local sites on which the local operator acts on. In particular, for the operators U and U † we have the representation , .
When legs of different operators are joined together it is understood that one should sum over the index of the corresponding local space. Finally, an explicit label for the legs can be omitted when it does not generate confusion. Using this notation, Eq. (6) can be rewritten as while Eq. (7) reads = , = , (10) where continuous solid lines represent the identity operator.
In the case of qubits, i.e. circuits with local dimension d = 2, it was shown in Ref. [50] that the most general dual-unitary gate reads where φ, J ∈ R, u ± , v ± ∈ SU(2) and Furthermore, it was shown that this family includes both integrable [64][65][66] and non-integrable cases. In particular, it contains a full parameter line of the integrable trotterized XXZ chain 65,66 and a quantum circuit representation of the self-dual kicked Ising (SDKI) model with We recall that the dynamics defined by the gate (14) is integrable for h = 0, while it is chaotic otherwise 50 , and, accordingly, in the latter case its spectral form factor is described by Random Matrix Theory 57 .

III. The solvable initial states
Even for dual-unitary quantum circuits, the computation of the time-evolution from arbitrary initial states appears to be extremely hard. Still, in the special case corresponding to the self-dual kicked Ising Floquet dynamics, it was found that the evolution of the bipartite Rényi entropies could be computed exactly for a particular family of product states 58 . This result relied on a specific mathematical property of the latter, which were called "separating". In this section we see that a similar logic can be followed for arbitrary dual-unitary circuits; in particular, we show how one can introduce a natural notion of "solvability" for a given initial state, and how this relates to the "separating" property in the special case of Ref. [58].
We consider a generic initial state |Ψ L 0 in the form of an MPS  where A ij j are matrices of dimensions χ j × χ j+1 . Since we are interested in the limit of infinite system sizes, it is natural to restrict to initial states that are invariant under translation of p sites, where p is some integer number. For reasons that will become clear later, we choose in the following p = 2, which also contains the class of translationally invariant states. In this case, we can rewrite so that now |Ψ L 0 only depends on two sets of matrices of dimensions χ × χ and χ × χ respectively. Finally, we consider MPSs that are normalized in the thermodynamic limit, namely In analogy with the previous section, we can make use of a standard graphical notation and represent the individual tensors as and so that MPSs are represented by a sequence of circles connected by lines, with additional out-coming legs corresponding to the physical local spaces.
To isolate the special property that the MPSs (17) should have in order to generate an exactly-solvable dynamics, it is instructive to consider the computation of the time evolution of one-point functions. Specifically, we consider Ψ L t |O 1 |Ψ L t , where O 1 is an arbitrary local operator acting on site j = 1. A pictorial representation of this expectation value is reported in Fig. 1, where we employed the graphical notations introduced above. Borrowing standard ideas from the literature on tensor networks 4,9,67 , one can write where E(t) and E O1 (t) are appropriate transfer matrices acting on the tensor product of 2t + 2 local sites along the "transverse direction". There are several (in general inequivalent) ways to define these operators. For the purpose of the present discussion, it is useful to define E(t) and E O1 (t) directly in terms of the elementary tensors in Eqs. (8), (19) and (20). In particular, using a graphical notation and focusing for concreteness on the case where t is even, we can define  (22) Here the right out-coming 2t + 2 legs correspond to the input space on which E(t) and E O1 (t) act on. The validity of Eq. (21) is straightforwardly established. Indeed, one can simply note that the graphical representation for tr E k (t)E O1 (t)E k (t) , which is obtained by placing 2k+1 transfer matrices side by side, is the same as for the expectation value of Ψ t |O 1 |Ψ t in a chain of 4k + 2 sites (where periodic boundary conditions are implemented).
Next, suppose that the largest eigenvalue λ 0 of E(t) is non-degenerate; more precisely, suppose that its algebraic multiplicity (namely, the number of diagonal elements in the Jordan form of E(t) that are equal to λ 0 ) is 1 and that there are no other eigenvalues with the same absolute value. Then, for large L where we used that the length of the system is 2L. Since Ψ L t |Ψ L t = Ψ L 0 |Ψ L 0 , Eq. (18) implies λ 0 = 1. In turn, this yields where we denoted by |L and |R the left and right eigenstates of E(t) associated with λ 0 , with the normalization L|R = 1. In general, the evaluation of Eq. (24) can only be done numerically for small times. However, for dual-unitary circuits there exist a class of states for which |L and |R can be determined exactly. Consider in particular an initial two-site shift invariant MPS, as defined in Eq. (17), and suppose that there exists a χ-dimensional matrix S such that Then one can show that is a right eigenstate of E(t) with eigenvalue λ 0 = 1. Here S α,β = α|S|β are the matrix elements of S in the basis {|α } χ α=1 of the auxiliary space associated with the initial MPS. The proof can be carried out graphically by noting that Eq. (25) can be represented as where we introduced the following notation Indeed, as shown in Fig. 2, this allows one to compute the action of E(t) on |R by making use of the diagrams in Eqs. (10) and (27). The above discussion motivates us to introduce the notion of solvable initial states for the quantum dynamics, and take Eq. (27) as a defining property of solvability. A priori, however, this condition alone is not sufficient to guarantee the uniqueness of the leading eigenvalue, which is necessary, for instance, to obtain Eq. (24). Accordingly, we say that a two-site shift invariant MPS [as defined in Eq. (17)] is solvable with respect to the class of dual-unitary quantum circuits if the following two conditions are satisfied: C1. The transfer matrix E(t) has a unique eigenvalue λ 0 with largest absolute value, λ 0 = 1 and λ 0 has algebraic multiplicity 1 ∀t ∈ N; C2. There exists a non-zero χ-dimensional matrix S satisfying Eq. (25) .
As we will see in the following, condition C1 could be removed. Indeed, as it will be clear from our derivations, the sum of k MPSs satisfying C1 and C2 gives us another MPS whose transfer matrix E(t) has k known eigenvectors with maximum absolute value, and for which analytic results could be derived. However, condition C1 allows us to classify the most "elementary" solvable states, from which all the others can be built simply out of linear superposition.
It is worth to discuss a connection between the present work and the findings of Refs. [9,67], where a folding technique to contract infinite tensor networks was introduced. Indeed, it can be seen that the solvable initial states are such that the leading right eigenvector |R of E(t) is a product state in the folded tensor network introduced in [9], in analogy with what happens from different initial states in the toy model studied in Ref. [67]. Note, however, that the latter was manifestly non-interacting, while the dual-unitary circuits generally implement a chaotic time evolution 52 , making the dynamics of solvable initial states non-trivial.
We also note that the logic of the present paper is similar to that of Refs. [68], [69], where it was shown that for any Bethe-Ansatz integrable Hamiltonian it is always possible to find "integrable" MPSs for which the quantum dynamics can be tackled analytically. In that case, these MPSs were defined as the initial states for which the transverse transfer matrix E(t) (obtained after a discretization of time through an appropriate Trotterization procedure 65 ) becomes Bethe-Ansatz integrable, so that its eigenvectors can be determined exactly. In our case, we stress that the transfer matrix E(t) corresponding to solvable MPSs is in general not Bethe-Ansatz integrable and, accordingly, we only have access to the eigenstate with largest absolute value.
Finally, before leaving this section, we recall that there is a close relation between the above definition of solvable and that of "separating" initial states, as defined in Ref. [58]. This is reported in Appendix A, to which we refer the interested reader.

IV. Classification of the solvable initial states
In this section, we proceed by providing a complete classification of the states that can be described by solvable MPSs in the thermodynamic limit. This is achieved by Theorem 1, which is stated below. The proof of the latter, which is based on well established techniques in quantum information and tensor network theory 4,70,71 , is rather technical and the interested reader can find it in Appendix B. In the rest of this section we only pro- vide the definitions that are needed in order to present its statement.
We begin by introducing the following parametrization for two-site shift invariant MPSs is a single set of χ-dimensional matrices, which encode all of the information stored in the sets {A j } d j=1 and {B j } d j=1 . We note that the solvability condition (25) can be written in terms of the tensors Next, we provide two additional definitions that are needed in order to state our main result. First, let {|Φ L 0 } L be a class of states defined on systems of increasing (even) sizes. We say that where O R is any observable acting non-trivially only on a finite product of local Hilbert spaces h j1 ⊗ · · · ⊗ h j R . In this case, we also say that |Φ L 0 and |Ψ L 0 (M) are equivalent in the thermodynamic limit.
Second, we recall the well-known notion of injectivity 4 : we say that a two-site shift invariant MPS is injective. It can be proven 4 that if |Ψ L 0 (M) is injective for L > 0, so is for L > L. This means that we can define the class Using the above definitions, we are now in a position to state one of our main results, i.e. the following theorem.
As we already mentioned, the proof of this theorem is rather technical, and is therefore reported in Appendix B.
At this point, it is important to note that Eq. (33) also allows us to write down the left eigenvector of the from which it follows that right and left eigenvectors are equal, namely Altogether, Theorem 1 provides us with a useful criterion to construct solvable MPSs, as we now explain. Given the tensors N (i,j) , we begin by defining the matrix W (N ) acting on the tensor product where i, j = 1, . . . d, α, β = 1, . . . χ, so that Eq. (33) can be straightforwardly rewritten as This equation is particularly useful: when combined with Theorem 1 it tells us that solvable states can be completely parametrized by matrices W ∈ End(C d ⊗C χ ) that are unitary. Next, note that if |Ψ L 0 (N ) is an MPS satisfying Eq. (33), then the same is true for the MPS where u j , v j ∈ U (d) are arbitrary unitary operators acting on the local Hilbert space h j . This means that solvable MPSs can be classified up to products of local unitaries. On the level of the matrix W (N ), this transformation reads as where 1 χ is the identity on C χ .
In the next section we see how the above consideration can be used to construct explicitly solvable MPSs.

V. Qubit systems: explicit solutions
In this section we address the explicit construction of solvable MPSs for the simplest case of a qubit system, corresponding to d = 2. In particular, we provide formulae for solvable MPSs with bond dimensions χ = 1 and χ = 2. We recall that here χ denotes the bond dimension of the MPS |Ψ L 0 (N ) obtained by grouping together two sites, namely it is the dimension of the matrices N (i,j) .

A. Bond dimension χ=1
The simplest example of solvable MPSs is given by product states, corresponding to bond dimension χ = 1. In this case, up to products of local unitaries, it is always possible to choose W (N ) ∝ 1 2 , which leads to the following explicit form On the other hand, it is trivial to verify that the state |Ψ L 0 defined above is injective. Putting all together, we obtain a single solvable MPS with bond dimension χ = 1.

B. Bond dimension χ=2
The case of bond dimension χ = 2 is more interesting. Now the unitary matrix W defined in Eq. (36) acts on the tensor product C 2 ⊗ C 2 . Then, we can use the following known parametrization 72 for W ∈ U (4), which is complete up to a global irrelevant phase: where t, u, v, z ∈ SU (2), while with K j ∈ R. By performing the matrix exponential V {K j } 3 j=1 , and rearranging the indices as in Eq. (36), we arrive at the following general parametrization (up to products of local unitaries) where u, v ∈ SU (2) and K ± = K 1 ± K 2 ∈ R. Finally, the MPS |Ψ L 0 (N ) is injective, except for a set of measure zero in the space of parameters {K j } 3 j=1 . This simply follows from the fact that for generic choices of {K j } 3 j=1 the operators N (i,j) span the whole set of 2 × 2 matrices.

VI. The exact quantum dynamics
In this section we finally explore the dynamics arising from the solvable MPSs. We focus in particular on three aspects: the thermalization time of local observables, the two-point correlation functions, and the entanglement growth.
A. Local thermalization and the quasi-particle picture As we have already discussed, the knowledge of the left and right eigenvectors of the transfer matrix E(t) allows us to compute the expectation values of observables that are supported over finite regions of space. However, it follows from Eq. (24) that the dynamics of observables localized at one site is quite trivial. Precisely, for a solvable MPS (evolved at time t) |Ψ L t (N ) , we have lim where we used Eq. (35). Namely, one-point functions remain constant at the infinite-temperature value.
In general, we can extend this result to local observables that are supported over finite regions of space containing more than one site. Consider the generic operator where O α j ∈ End (h j ) and where tr O α j = 0. First note that due to even-odd effects, its expectation value on time-evolved solvable MPSs will be always zero unless j R − j 1 ≡ 1 (mod 2). Second, one can prove where The proof of this can be easily performed graphically. Indeed, for a given observable O {α1:α R } {j1:j R } , one can write down the tensor network corresponding to its expectation value. By a repeated use of the graphical identities in Eqs. (9), (10) and (27) it is then straightforward to see that if t > t * the tensor network simplifies to a new one, where all the "upper legs" of O {α1:α R } {j1:j R } are connected to its "lower legs". In turn, this is exactly the tensor network corresponding to the infinite-temperature expectation value in the r.h.s. of Eq. (49). From the graphical representation it is also straightforward to see that the correlation function in Eq. (49) displays an even-odd effect in time, i.e., it is zero for all even (odd) times if j 1 is even (odd) (where we are again considering traceless operators).
Eq. (49) implies that, while the time evolution of solvable MPSs is non-trivial, finite regions reach infinite temperature in a time which is proportional to their sizes. It is interesting to observe how this behavior is exactly predicted by the standard conformal quasi-particle picture 11 . Indeed, t * in Eq. (50) is exactly the time needed for a pair of quasi-particles produced at the center of the region R and moving at maximal speed (which is 1 for our choice of units) to exit from it. Finally, we see that local thermalization to infinite temperature takes place irrespective of the presence of local conserved quantities (namely, of the integrability of the unitary dynamics) and that t * does not depend on the unitary gates chosen. In the next section, we will see that qualitative differences emerge in the study of the light-cone spreading of two-point correlation functions for integrable and nonintegrable circuits.

B. The two-point correlation functions
We now move on to examine two-point correlation functions. For concreteness we focus on the case of qubits (although our treatment is valid for general physical local dimension d), and consider where σ α j , α = x, y, z are Pauli matrices acting at site j, while |Ψ L t (N ) is a solvable (injective) initial MPS satisfying Eq. (33), evolved at time t. From the results of the previous subsection, we have that for fixed j, r the function C α,β (j, r, t) will become zero after a time t * = (r + 1)/2. Here, however, we are interested in its full time evolution, which is non-trivial for generic unitaries U .
It turns out that it is possible to compute exactly C α,β (j, r, t) for arbitrary values of j, r and t, following an approach similar to the one developed in Ref. [50] for infinite-temperature dynamical two-point functions. The first step consists in simplifying the graphical representation for C α,β (j, r, t) by means of the dual unitarity conditions (9), (10) and the solvability relation (27). By doing this, we obtain the formula C α,β (j, r, t) = δ r(mod 2),1 δ j−t(mod 2),1 C α,β (r, t) , (52) where Here D α,β 1 (t) and D α,β 2 (r, t) are functions which admit a simple graphical representation. As an example, D α,β 2 (r, t) is depicted in Fig. 3 for r = 7 and t = 2. The simplified tensor networks associated with D α,β 1 (t) and D α,β 2 (r, t) can be contracted efficiently, by slightly generalizing the method employed in Ref. [50]. In particular, based on its graphical representation, one can derive the following formula for the function D α,β 1 (t) where χ is the bond dimension of the initial solvable MPS |Ψ L 0 (N ) and (·) T denotes matrix transposition. Here we introduced the following definitions. First, the functions F, F are maps acting on the space of linear operators End(C 2 ), reading where d = 2 for qubits. Second, the function P is a map on the same space which is defined in coordinate space as where repeated indices are summed over. From the above definitions, it is clear that Eq. (54) can be easily evaluated numerically, with a little computational effort. An expression with a similar structure can be derived for the function D α,β 2 (r, t). In particular, we have Here F, F are given in Eq. (55) and (56), while E N reads Next, the functions G and G are maps acting on the space of linear operators End(C 2 ), and are defined in coordinate space as where repeated indices as summed over. Note that also the function E N admits a convenient graphical representation, which reads (E N [a]) m,n = a n m .
Once again, the proof of Eq. (58) can be easily obtained by graphical inspection, and since it presents no difficulty we omit it here. Putting all together, Eq. (53) provides an efficient formula for the computation of two-point correlation functions, which allowed us to produce numerical results for different solvable initial states |Ψ L 0 (N ) and dual-unitary operators U . Examples of our findings are displayed in Fig. 4 and 5, where we report data for different choices of the latter, and the same solvable MPS. Note that in order to remove even-odd effects, we have plotted correlations averaged over neighboring sites, namely C α,β 0 (r, t) = 1 4 x,y=±1 which only depend on the relative distance r of the two Pauli matrices. In Fig. 4(a), we report the two-point correlation function C x,x 0 (r, t) for the integrable point of the self-dual kicked Ising chain, cf. Eq. (14): we see that there is no decay along the main light-cone. Conversely, in Fig. 4(b),(c) we report data for unitary gates where integrability is broken by increasing the value of the magnetic field h and in this case the plots clearly display an exponential decay in time. In fact, from the explicit formula (54) and the results of Ref. [50] this behavior is expected. Indeed, in Ref. [50] unitary gates were classified in four classes of increasing level of ergodicity, depending on the behavior of dynamical infinite-temperature correlation functions. Even though the formula for C α,β 0 (r, t) along the main light-cone is different from the one found in Ref. [50], both of them involve a repeated application of the function (1-qudit channel) F defined in Eq. (55), and hence depend essentially on the spectrum of F. Accordingly, the classification of unitary gates reported in Ref. [50] not only distinguishes the behavior of dynamical correlations at infinite temperature, but also applies to predict the qualitative features of two-point functions during the quantum dynamics arising from solvable initial MPSs .
In Fig. 5 similar data are reported for C z,z 0 (r, t) and different dual-unitary quantum circuits. In particular, Fig. 5(a) displays the case of the integrable evolution corresponding to the XXZ gate in Eq. (13), while in Fig. 5(b) and Fig. 5(c) integrability is broken by choosing nontrivial matrices u ± , v ± in Eq. (11). As for Fig. 4, we see that a breaking of integrability corresponds to an exponential decay of the correlations along the light-cone.
Finally, it is important to stress that the analytic formula in Eq. (53) does not predict a broadening of the light-cone during the time evolution, as we also observe from Figs. 4 and 5. Once again, this is consistent with an exact CFT picture of non-interacting quasi-particles that are created in pairs at time t = 0 and spread ballistically for t > 0 [11].

C. The entanglement growth
As a final aspect of the dynamics arising from solvable MPSs, we discuss the spreading of entanglement. In particular, we consider the setting already studied in Ref. [58] for the special case of the self-dual kicked Ising chain. We focus on a system of size 2L with periodic boundaries, and compute the time evolution of the entanglement of the subsystem A -composed by a connected block of sites -and the rest, in the limit L → ∞ (see also Ref. [73], where a similar analysis was carried out for the entanglement of half chain). The initial state is given by a solvable MPS |Ψ L 0 (N ) and the system is evolved using a dual-unitary quantum circuit. As it is customary, in order to measure the spreading of entanglement we study the growth of the Rényi entropies S (α) where ρ A (t) is the density matrix reduced to the subsystem A = [j 1 , j 1 + −1] containing neighboring sites. At this point we note that the entanglement growth displays an even/odd effect in space and time. In order to simplify our discussion, in the rest of this section we will consider the case of t even, and choose j 1 to be odd, although a very similar treatment can be carried out, with minor modifications, when either t is odd or j 1 is even. In the case of t even and j 1 odd, we can further assume, without loss of generality, to be even. Indeed, if = 2k + 1, it is easy to see that where A k is the connected region containing sites from j 1 to j 1 + k − 1, while 1 k is the identity acting on site j 1 + k − 1. Using (27) and the unitarity of the gates we can write the thermodynamic limit of the reduced density matrix as where we used that, thanks to Theorem 1, the matrix S in (27) can be chosen equal to the identity. This circuit can be further contracted by repeated use of the initial state's solvability condition (27) combined with the dual-unitarity property (10). Two different results are obtained depending on whether or not 2t is larger than . Specifically, for 2t ≤ we find while for 2t > the circuit simplifies to The explicit form (66)-(67) of the reduced density matrix directly gives where θ H [x] is the step function (with θ H [0] = 0) and we introduced the d x -dimensional matrix The result (68) follows directly from the observation that the reduced density matrix (66) is unitary equivalent to the tensor product of the matrices O −2t and 1 2t /d 2t .
To analyze (68), it is instructive to consider its prediction for S A (t)/ in the scaling limit , t → ∞ where the ratio /t = ζ is kept fixed. Let us start considering the second term on the r.h.s. of (68). This term can be written in terms of the state transfer matrix [cf. (B1)] (71) In the scaling limit, for any ζ < 2 (the term is trivially zero in the opposite case), the number of matrices τ (N ) in (71) becomes infinite and we can make the following replacement where we used that, thanks to Theorem 1, the unique eigenvector of τ (N ) corresponding to eigenvalue 1 is a maximally entangled pair of qudits. The replacement (72) leads to In particular, note that the r.h.s. of (73) is finite: this means that only the first term on the r.h.s. of (68) contributes to the scaling limit when dividing by . As a consequence, S In the last few years this form has been observed in conformal invariant models 11,74 , in generic isolated systems [75][76][77] , and in local random unitary circuits [34][35][36]39 . The ubiquitousness of Eq. (74) has been recently explained by introducing the so called "minimal membrane picture" 34 . In essence, the idea is to estimate the entanglement between a subsystem A and the rest by measuring the length of the minimal membrane in space-time that separates the subsystem from the rest. Ref. [58] showed that (74) is exact at any time in the self-dual kicked Ising model evolving from "separating states". Considering the special case χ = 1 of (68) we see that this statement carries over for generic dual-unitary circuits. Indeed, in this case we have tr [O α x ] = 1 for all α. From (68), however, we also see that for more general initial "solvable" MPSs there appear some corrections that can make the entanglement spectrum non-trivial. These corrections are encoded by the matrix O x and depend solely on the initial state.
We remark that the corrections contained in (68) can be generically calculated numerically with high efficiency because the rank of the matrix that encodes them is constant and equal to χ. A simple analytically tractable limit is that of infinite length of the block A. Indeed, proceeding as we did to obtain (73), we readily find lim →∞ S (α) We see that in this case the entanglement spectrum is again completely flat. Finally, we note that the form (66)-(67) of the reduced density matrix can also be used to compute the evolution of the entanglement entropy of disjoint blocks. In this case one has to connect some of the central uppermost and lowermost lines of the circuits (66)-(67) to divide A into two disconnected parts. The resulting quantum circuit is again very simple in the case (67) and implies that for 2t > the reduced density matrix is again proportional to the identity. The case (66), however, becomes more complicated: the density matrix is not unitary equivalent to O −2t ⊗1 2t /d 2t anymore, and the result depends on the specific dual-unitary gate considered.

VII. Non-solvable initial states
It is natural to wonder how the features of the quantum dynamics studied in the previous section depend on the solvability of the MPSs chosen as initial states. In general, for instance, one may expect that while the system will still locally approach an infinite temperature density matrix for generic initial states and unitary gates, local expectation values will display an exponential decay to zero, rather than approaching zero in a finite number of time steps.
The dynamics arising from generic states can be studied using numerical MPS techniques 10 . Here, in order to compare the solvable and non-solvable cases, we follow a different approach, and consider a family of initial states which depend on one real parameter β, interpolating between the infinite temperature density matrix and arbitrary initial pure product states, as β is varied from zero to infinity. This allows us to study analytically the dynamics for small values of β, highlighting some qualitative differences that arise with respect to the solvable dynamics.
For concreteness, we focus once again on the case of a qubit system and consider the initial mixed state where we introduced the single-site density matrix Here a ∈ End(C 2 ), with tr(a) = 0 and a 2 = 1, while Z(β, a) = tr [exp (−βa)], so that ρ (j) (β, a) admits the expansion In this section, we will focus on the computation of timedependent one-point correlation functions. First, note that according to Eqs. (1) and (2), the initial density matrix (76) is time-evolved as Next, making use of Eq. (78), one can formally write In the following, we show how dual unitarity allows us to compute the coefficients c a n (t) exactly up to n = 2, providing a perturbative knowledge of the one-point correlation functions.
As for the solvable initial states, one-point functions for the density matrix (76) display an even/odd effect in time, due to the discrete nature of the dynamics. In order to simplify the following discussion, in analogy to Sec. VI C, also in this section we restrict for concreteness to even values of times and choose an operator placed at an odd site j = 2k + 1.
We begin by noting that Eq. (81) can be pictorially represented as in Fig. 6, where a small black rectangle placed at site j denotes the single-site density matrix ρ (j) (β, a) defined in Eq. (77). Next, from Eq. (78) it is clear that each operator a bears a factor tanh(β), so that at the Pictorial representation for the first coefficient c a 1 (t), for τ = 2 (namely, t = 4). From the picture, it is clear that it can be computed using the analytic formula for the infinite-temperature dynamical two-point functions derived in Ref. [50]. The small green rectangle corresponds to the operator −1/2 tanh(β)a first order in x = tanh(β) Fig. 6 simplifies in a sum of tensor networks, each one corresponding to setting all operators ρ (j) (β, a) equal to 1/2, except for one. By means of the usual graphical identities, it is straightforward to see that only one term in this sum is non-zero. In the case of Fig. 6, for instance, this corresponds to the tensor network displayed in Fig. 7, where a small green rectangle now corresponds to the operator −1/2 tanh(β)a. We recognize that the latter is proportional to an infinitetemperature two-point correlation function, and can thus be computed using the results of Ref. [50], thus obtaining Here F[a] is the map defined in Eq. (56), (·) T denotes, as before, matrix transposition, while τ was introduced in Eq. (82). The same logic can be followed to compute the secondorder coefficient c a 2 (τ ), which can be expressed as a sum of some particular three-point correlation functions. The computation of the latter is slightly more involved, but one can once again exploit dual unitarity to evaluate them in an efficient way, in complete analogy with what we did in Sec. [introduced in Eqs. (55) and (56)], but with U and U † exchanged, namely with d = 2 for qubits. Finally, the function S is a map acting on the space of linear operators End(C 2 ⊗ C 2 ) which is defined in coordinate space as Dual unitarity does not seem to provide a substantial simplification in the computation of higher-order coefficients c a k>2 (t) when compared to a generic evolution. Still, the cases k = 1, 2 already allow us to highlight some of the expected differences with respect to the dynamics arising from solvable MPSs. In Fig. 8 we present results for the evolution of the expectation value of O j = σ x j in a chaotic quantum circuit corresponding to the self-dual kicked Ising chain (we chose U SDKI as defined in Eq. (14), and set the magnetic field to h = 0.15). We initialized the system in the state (76) with a = σ x . We see that, as expected, both coefficients c x 1 (τ ), c x 2 (τ ) approach zero, although in an exponential fashion rather than in a finite number of time steps. In general, based on this result, we also expect that our analytic formulae for two-point correlation functions, valid for solvable MPSs, will acquire some "exponential corrections" in the case of generic initial states and dual-unitary gates.

VIII. Conclusions
In this work we have considered the dynamics of the recently introduced class of dual-unitary quantum circuits, and exhibited a family of initial states for which different physical quantities can be computed exactly. We have characterized these states in terms of a particular "solvability" condition, and provided a complete classification of them. In particular, we have shown that this family includes MPSs of arbitrary bond dimension χ, and provided explicit examples for χ = 1 , 2. (cf. Sec. V).
For the dynamics arising from these initial states, we have derived three main results. First, we were able to compute explicitly the time it takes for an observable supported on a finite region to approach its infinitetemperature value, and shown that t ∝ (Sec. VI A). Second, we have provided an exact formula, efficient to evaluate, to compute two-point equal-time correlation Here we chose the local operator Oj = σ x j , while for the initial Gibbs state we set a = σ x . The dynamics is driven by the quantum circuit corresponding to the self-dual kicked Ising chain, namely we chose USDKI as defined in Eq. (14), and set the magnetic field to h = 0.15.
functions of local observables, and we have shown that they display different qualitative features depending on the ergodicity of the quantum circuit (Sec. VI B). Third, we have derived a closed formula for the time evolution of the entanglement entropy of a connected block of length , computing explicitly the limit → ∞ (Sec. VI C). This generalizes the exact result of Ref. [58] to all dualunitary quantum circuits, extending it for more general initial states. Remarkably, we showed that solvable MPSs with bond dimension larger than one produce non-trivial finite-time corrections. Finally, we have also considered a family of nonsolvable initial mixed states depending on one real parameter β, which interpolate between the infinite temperature density matrix and arbitrary initial pure product states, cf. Sec. VII. We have studied analytically their dynamics for small values of β, and highlighted the main differences from the case of solvable MPSs.
In the light of our results, there are several natural directions to be explored. For example, the class of solvable MPSs introduced in this paper was defined by the property that the leading eigenstate of the transverse transfer matrix is a product state in the "folded picture" of Refs. [9,67]. One can wonder whether it is possible to find different types of initial states or unitary gates for which such leading eigenstates are written instead in the form of non-trivial MPSs (with fixed bond dimension). In this case, the dynamics would still be solvable, but a richer phenomenology would emerge for the evolution of local observables.
Finally, it is certainly interesting to wonder whether some aspects of the present paper can be generalized to the case of higher spatial dimensions, where the application of numerical techniques is known to be much harder with respect to the one-dimensional case. A successful description of the dual-unitary dynamics in higher spa-tial dimensions would provide a rare benchmark, for instance, for the development of numerical computational methods.
Acknowledgments LP acknowledges support from the Alexander von Humboldt foundation. BB and TP acknowledge support by the EU Horizon 2020 program through the ERC Advanced Grant OMNES No. 694544, and by the Slovenian Research Agency (ARRS) under the Programme P1-0402. JIC acknowledges support by the EU Horizon 2020 program through the ERC Advanced Grant QENO-COBA No. 742102, and from the DFG (German Research Foundation) under Germany's Excellence Strategy -EXC-2111 -390814868. Part of this work has been carried out during a visit of LP to the University of Ljubljana, whose hospitality is kindly acknowledged.

A. Comparison between the solvability and separating conditions
In this section we discuss the relation between the "separating" initial states of Ref. [58] and the solvable MPSs introduced in this paper. First, we observe that the work [58] focused on one particular realization of a dualunitary quantum circuit, corresponding to the self-dual kicked Ising model 57 . In this specific case, the separating property was also introduced as a technical condition on initial product states, allowing for an analytic determination of the right eigenstate of an appropriate transverse transfer matrix.
There are two main differences between the present work and Ref. [58]: first, here we allow for more general initial states in the form of MPSs, and do not restrict to product states; second, we look for initial conditions that are analytically tractable for any dual-unitary circuit, and not just for one specific case. In fact, we show in the following that the separating states of Ref. [58] are not solvable MPSs for a generic dual-unitary evolution, but become so after applying one layer of unitary operators, which corresponds to the first time step of the self-dual kicked Ising Floquet dynamics.
We start by recalling that the Floquet evolution associated with the self-dual kicked Ising chain is defined by 50,73 where t ∈ N. Here U e SDKI (U o SDKI ) is the transfer matrix defined by applying the local gate (14) to each pair of neighboring qubits, where the first one is chosen at an even (odd) position, while U SDKI = U e SDKI U o SDKI . The operators U e/o I and U I are defined analogously, with two-site gates given by Next, the separating states were defined in Ref. [58] as the two subclasses of product states defined by where 1 denotes a vector of length 2L with all entries equal to 1. Note that we do not need to specify the value of φ in Eq. (A7). In order to compare with the present article, where we restricted to two-site shift invariant initial states, we can assume We see that the above class of states does not belong to the family of solvable MPSs of dimension χ = 1, which was characterized in Sec. V A. However, from Eq. (A1), (A2), we have that the Floquet dynamics of the self-dual kicked Ising model is made up of different steps. In particular, given an initial state |ψ , the first step consists in multiplying it by either U e I or U e SDKI U e I : after that the evolution is dictated by the dual-unitary quantum circuit encoded in U t SDKI . Accordingly, one should check that if |ψ ∈ T ∪ L, then U e I |ψ and U e SDKI U e I |ψ are solvable. This is indeed the case, as one can verify by direct computation.
In conclusion, the Floquet dynamics of the separating product states introduced in Ref. [58] can be described in terms of a dual-unitary circuit where the initial configuration is indeed a solvable MPS with bond dimension χ = 1.

B. Proof of Theorem 1
In this appendix we provide a full proof of Theorem 1. Let us start by proving the first part of the statement, namely that each solvable MPS is equivalent in the thermodynamic limit to an injective MPS satisfying Eq. (33).
Consider the following transfer matrix associated with the MPS |Ψ L 0 (M) which acts on the tensor product C χ ⊗C χ of two auxiliary spaces. It is easy to show that condition C1 is equivalent to requiring that the transfer matrix τ (M) has a unique eigenvalue λ 0 with largest absolute value, with λ 0 = 1 and algebraic multiplicity equal to 1. This follows from Ψ L t |Ψ L t = Ψ L 0 |Ψ L 0 , and the identity where we used that the system size is 2L. Defining now the state where V is the χ-dimensional matrix with entries V i,j . Hence, the transfer matrix τ (M) has a unique largest eigenvalue λ 0 = 1 if and only if the linear map has a unique largest eigenvalue λ 0 = 1. On the other hand, setting i = j in Eq. (30), and summing over j, we obtain d j,k=1 Namely, combining conditions C1 and C2, we obtain that if |Ψ L 0 (M) is a solvable MPS, then S is the only fixed point for the linear map E M (X) in Eq. (B5).
It is important to observe now that E M (X) is a completely positive, linear map on the space of matrices End(C χ ) and that its spectral radius (defined as the largest among the absolute values of its eigenvalues) is equal to 1. This allows us to exploit known results in the theory of positive maps on C * -algebras. In particular, using Theorem 2.5 in Ref. [78], it follows that E M (X) has a positive fixed point 4 , namely there exists a Hermitian matrix T with a non-negative spectrum such that E M (T ) = T . However, due to condition C1, it follows that T ∝ S, since E M (X) has a unique fixed point. Hence, up to a proportionality factor, we just showed that for a solvable MPS, the matrix S in condition C2 must be positive.
Suppose now that the matrix S is invertible. Then S must be strictly positive and we can make use of the following lemma to conclude that the MPS |Ψ L 0 (M) is injective. 2. The corresponding eigenvector Λ is a strictly positive operator, namely it is Hermitian with strictly positive spectrum.
The proof of this Lemma is non-trivial, as it requires some technical tools in the theory of quantum channels. The basic idea underlying the Lemma is to exploit the fact that injectivity for sufficiently large L implies that the channel E M (X) is primitive 70 . In turn, one can show that this property is equivalent to the above conditions 1 and 2. We omit further details of the proof here, and refer the reader to Ref. [70] where it is explicitly worked out, cf. Prop. 3 therein. For a pedagogical introduction of the relevant quantum information tools, see instead Ref. [79].
Conversely, suppose that the matrix S is not invertible. Then, we can assume that there exists a basis of C χ in which the matrices M (i,j) can be written in block diagonal form, so that without loss of generality we can assume where {M (i,j) α } d i,j=1 are χ α -dimensional matrices with χ α < χ, α = 1, 2. In order to see this, we proceed as follows.
Since S is positive, we can write S = χ α=1 µ α |α α|, where χ < χ and µ α > 0. Following [4], and defining P R to be the projector onto the space R spanned by |α 's, we can prove that M (i,j) P R = P R M (i,j) P R , namely M (i,j) |α ∈ R , ∀(i, j), α. Indeed, suppose that this is not true. Then, there exists (m, n), β such that which is a contradiction. Thus, the matrices M (i,j) must be block-triangular. However, it is now immediate to show that we can write the same state as an MPS |Ψ L 0 ( M) with tensors M (i,j) that are block-diagonal and obtained from M (i,j) by setting to zero off-diagonal terms. Putting all together, we conclude that we can indeed assume Eq. (B7) without loss of generality, with χ 1 being the rank of the matrix S.
We note now that Eq. (B7) implies namely |Ψ L 0 (M) can be written as the sum of two distinct MPSs. Next, we observe that the norm of |Ψ L 0 (M 2 ) must be vanishing in the thermodynamic limit, since the map E M2 (X) has spectral radius r < 1. Indeed, suppose that this is not the case and there exists S such that E M2 (S) =λ 0S with |λ 0 | ≥ 1. Then, are both eigenstates for E M (X) with eigenvalue |λ 0 |, |λ 0 | ≥ 1, which contradicts the hypothesis. In conclusion, one can find an injective MPS |Ψ L 0 (M 1 ) , with bond dimension χ 1 < χ, which is equivalent to |Ψ L 0 (M) in the thermodynamic limit.
Putting all together, up to equivalence in the thermodynamic limit, we can always assume that the fixed point S of E M (X) is strictly positive, and define N (i,j) = S −1/2 M (i,j) S 1/2 , where we also have S 1/2 † = S 1/2 . It is now straightforward to verify that |Ψ L 0 (N ) is equivalent to |Ψ 0 (M) in the thermodynamic limit, is injective and satisfies Eq. (33). This proves the first part of the statement of Theorem 1.
Finally, let us prove the second part of Theorem 1, namely that if a state is equivalent to an injective MPS satisfying (33), then it is also equivalent to a solvable MPS satisfying conditions C1 and C2. Clearly, using Lemma 1, we only need to prove that the algebraic multiplicity of the leading eigenvalue is equal to 1. This can be done once again by invoking a general result in the theory of quantum channels 80 : if E M (X) is a positive linear map with spectral radius r = 1 satisfying Eq. (33) (i.e. it is unital ) and λ 0 is an eigenvalue with |λ 0 | = 1, then its algebraic and geometric multiplicities coincide, namely all blocks in the Jordan form corresponding to λ 0 are onedimensional. Here we omit the proof of this statement, for which we refer the reader to Ref. [80] (see Proposition 6.2 therein). This concludes the proof of Theorem 1.