Influence matrix approach to many-body Floquet dynamics

Recent experimental and theoretical works made much progress towards understanding nonequilibrium phenomena in thermalizing systems, which act as thermal baths for their small subsystems, and many-body localized ones, which fail to do so. The description of time evolution in many-body systems is generally challenging due to the dynamical generation of quantum entanglement. In this work, we introduce an approach to study quantum many-body dynamics, inspired by the Feynman-Vernon influence functional. Focusing on a family of interacting, Floquet spin chains, we consider a Keldysh path-integral description of the dynamics. The central object in our approach is the influence matrix (IM), which describes the effect of the system on the dynamics of a local subsystem. For translationally invariant models, we formulate a self-consistency equation for the influence matrix. For certain special values of the model parameters, we obtain an exact solution which represents a perfect dephaser (PD). Physically, a PD corresponds to a many-body system that acts as a perfectly Markovian bath on itself: at each period, it measures every spin. For the models considered here, we establish that PD points include dual-unitary circuits investigated in recent works. In the vicinity of PD points, the system is not perfectly Markovian, but rather acts as a bath with a short memory time. In this case, we demonstrate that the self-consistency equation can be solved using matrix-product states (MPS) methods, as the IM temporal entanglement is low. A combination of analytical insights and MPS computations allows us to characterize the structure of the influence matrix in terms of an effective “statistical-mechanics” description. We finally illustrate the predictive power of this description by analytically computing how quickly an embedded impurity spin thermalizes. The influence matrix approach formulated here provides an intuitive view of the quantum many-body dynamics problem, opening a path to constructing models of thermalizing dynamics that are solvable or can be efficiently treated by MPS-based methods, and to further characterizing quantum ergodicity or lack thereof.

Recent experimental and theoretical works made much progress towards understanding nonequilibrium phenomena in thermalizing systems, which act as thermal baths for their small subsystems, and many-body localized ones, which fail to do so. The description of time evolution in many-body systems is generally challenging due to the dynamical generation of quantum entanglement. In this work, we introduce an approach to study quantum many-body dynamics, inspired by the Feynman-Vernon influence functional. Focusing on a family of interacting, Floquet spin chains, we consider a Keldysh path-integral description of the dynamics. The central object in our approach is the influence matrix (IM), which describes the effect of the system on the dynamics of a local subsystem. For translationally invariant models, we formulate a self-consistency equation for the influence matrix. For certain special values of the model parameters, we obtain an exact solution which represents a perfect dephaser (PD). Physically, a PD corresponds to a many-body system that acts as a perfectly Markovian bath on itself: at each period, it measures every spin. For the models considered here, we establish that PD points include dual-unitary circuits investigated in recent works. In the vicinity of PD points, the system is not perfectly Markovian, but rather acts as a bath with a short memory time. In this case, we demonstrate that the self-consistency equation can be solved using matrix-product states (MPS) methods, as the IM temporal entanglement is low. A combination of analytical insights and MPS computations allows us to characterize the structure of the influence matrix in terms of an effective "statistical-mechanics" description. We finally illustrate the predictive power of this description by analytically computing how quickly an embedded impurity spin thermalizes. The influence matrix approach formulated here provides an intuitive view of the quantum many-body dynamics problem, opening a path to constructing models of thermalizing dynamics that are solvable or can be efficiently treated by MPS-based methods, and to further characterizing quantum ergodicity or lack thereof.

I. INTRODUCTION
Describing non-equilibrium quantum matter and harnessing it for quantum technology is one of the central challenges in modern physics. The problem of highly non-equilibrium dynamics of many-body systems, both isolated and open, has been attracting intense experimental and theoretical interest over the past years [1][2][3]. Ergodic isolated systems are believed to thermalize as a result of their quantum evolution; qualitatively, such a system can act as an efficient thermal bath for its sufficiently small subsystems. Recent breakthroughs identified classes of systems that do not reach thermal equilibrium [4][5][6][7] and therefore may exhibit new phenomena not envisioned within the framework of statistical mechanics.
Floquet systems, where the Hamiltonian is periodically varied in time, play a special role in the family of non-equilibrium systems, thanks to their natural experimental realizations. Although periodic driving sequences have been utilized in nuclear magnetic resonance for decades [8], recent works revealed a range of new surprising phenomena in Floquet systems. In particular, it was shown that many-body Floquet systems may exhibit new topological properties [9]. Many-body localization in Floquet systems can protect them from heating [10- * These two authors contributed equally to this work 12], enabling new non-equilibrium states of matter with remarkable properties not attainable in thermal equilibrium [13][14][15][16][17].
The central difficulty in describing dynamics of manybody systems that do thermalize stems from the rapid generation of quantum entanglement. Initially simple, non-entangled states, quickly develop non-local correlations; faithfully describing a time-evolved state requires, in general, a number of parameters which grows exponentially with the evolution time. Various efficient numerical methods based on tensor networks have been introduced [18][19][20][21]. Examples of tractable interacting Floquet models, both integrable and thermalizing, have been found and are being actively investigated [22][23][24][25][26].
The goal of this paper is to formulate what we call the influence matrix approach to quantum many-body Floquet dynamics. This approach can be viewed as an extension of the celebrated Feynman-Vernon influence functional approach [27] to interacting Floquet systems. In the original formulation, Feynman and Vernon developed a path-integral description of a quantum-mechanical system coupled to a bath. Typically, the bath consists of physical degrees of freedom that are of different nature than those composing the system itself, such as in the case of two-level systems coupled to a bath of harmonic oscillators [28,29] or particle reservoirs [30]. In our case, in contrast, the bath and the system will be composed of the same physical constituents.
We will consider quantum systems on a lattice, and < l a t e x i t s h a 1 _ b a s e 6 4 = " z j l l p T u Z l i T q 0 A k r H 2 4 j c y + r 9 d 0 = " > A A A C F n i c b V C 7 S g N B F J 3

G e M r a m k z J A g R J O x K Q O 2 C N t p F M A / I h n B 3 c h M H Z x / M 3 B X C k t 5 P 8 C t s t b I T W s L / 8 X d J I V G T 3 U 4 5 7 u u c e L l D R k 2 5 / W w u L S 8 s p q b i 2 / v r G 5 t V 3 Y 2 W 2 a M N Y C G y J U o W 5 7 Y F D J A B s k S W E 7 0 g i + p 7 D l 3 V k f u s e t Z F h c E O j C L s + D A M 5 k A I o l X q F o u s D 3 Q p Q y d W 4 7 C a u k U M f j l w P 9 J S 6 4 8 N e o W R X 7 A n 4 X + L M S I n N U O 8 V v t x + K G I f A x I K j O k 4 d k T d B D R J o X C c d 2 O D E Y g 7 G G I n p Q H 4 a L r J 5 J c x P 4 g N U M g j F w q P h H x 5 0 Y C v j E j 3 0 s n s + R m 3 s v E / 7 x O T I P T b i K D K C Y M R H a I p M L J I S O 0 T E t C 3 p c a i S B L j l w G X I A G I t S
S g x C p G K e t 5 d M + n P n v / 5 L m c c W p V s 6 u q 6 X a + a y Z H N t n R V Z m D j t h N X b J 6 q z B B H t g T + y Z v V i P 1 q v 1 Z r 1 P R x e s 2 c 4 e + w X r 4 x t H + p + f < / l a t e x i t > { ,¯ } < l a t e x i t s h a 1 _ b a s e 6 4 = " + w z i H g M 7 M h + n G v Y 7 g 8 O v K x e y J 6 M = " > A A A C C X i c b V D L S g N B E J z 1 G e M r K p 6 8 D A b B g 4 R d C a i 3 o B e P E c w D s k v o n X T i k N k H M 7 1 C W P I F f o V X P X k T r 3 6 F B / / F T b I H T a x T U d V N V 5 c f K 2 n I t r + s p e W V 1 b X 1 w k Z x c 2 t 7 Z 7 e 0 t 9 8 0 U a I F N k S k I t 3 2 w a C S I T Z I k s J 2 r B E C X 2 H L H 9 5 N p e J T E X 9 v p B A Y M w r 8 b D I A e j D z 3 k T 8 z + s k 1 L / 0 U h n G C W E o J o d I K p w e M k L L r B f k P a m R C C b J k c u Q C 9 B A h F p y E C I T k 6 y o Y t a H M / / 9 I m m e V 5 x q 5 e q u W q 5 d 5 8 0 U 2 B E 7 Z q f M Y R e s x m 5 Z n T W Y Y C l 7 Z i / s 1 X q y 3 q x 3 6 2 M 2 u m T l O w f s D 6 z P H 2 b j m l U = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " z j l l p T u Z l i T q 0 A k r H 2 4 j c y + r 9 d 0 = " > A A A C F n i c b V C 7 S g N B F J 3 1 G e M r a m k z J A g R J O x K Q O 2 C N t p F M A / I h n B 3 c h M H Z x / M 3 B X C k t 5 P 8 C t s t b I T W 1 s L / 8 X d J I V G T 3 U 4 5 1 7 u u c e L l D R k 2 5 / W w u L S 8 s p q b i 2 / v r G 5 t V 3 Y 2 W 2 a M N Y C G y J U o W 5 7 Y F D J A B s k S W E 7 0 g i + p 7 D l 3 V 1 k f u s e t Z F h c E O j C L s + D A M 5 k A I o l X q F o u s D 3 Q p Q y d W 4 7 C a u k U M f j l w P 9 J S 6 4 8 N e o W R X 7 A n 4 X + L M S I n N U O 8 V v t x + K G I f A x I K j O k 4 d k T d B D R J o X C c d 2 O D E Y g 7 G G I n p Q H 4 a L r J 5 J c x P 4 g N U M g j 1 F w q P h H x 5 0 Y C v j E j 3 0 s n s + R m 3 s v E / 7 x O T I P T b i K D K C Y M R H a I p M L J I S O 0 T E t C 3 p c a i S B L j l w G X I A G I t S S g x C p G K e t 5 d M + n P n v / 5 L m c c W p V s 6 u q 6 X a + a y Z H N t n R V Z m D j t h N X b J 6 q z B B H t g T + y Z v V i P 1 q v 1 Z r 1 P R x e s 2 c 4 e + w X r 4 x t H + p + f < / l a t e x i t > Figure 1. a) Cartoon illustration of the influence matrix approach to many-body Floquet dynamics. The many-body quantum state of a spin system becomes increasingly entangled due to periodic local interactions and kicks. The influence matrix I, defined in Eq. (15) below, describes the dynamical effects of all spins k < p on a spin p, in terms of a path integral weight affecting the trajectories {σ,σ} of spin p forward and backward in time. b) In translationally-invariant one-dimensional systems, the influence matrix is the same for all p, which, as illustrated, leads to a self-consistency equation for it.
will be interested in the dynamics of a finite subsystem, treating its complement as a bath. The effect of this bath on the subsystem can be described by the influence functional. Although here we focus on Floquet systems, generalizations to Hamiltonian systems appear possible. More specifically, we will study a class of kicked Floquet systems, which can be viewed as many-body extensions of the celebrated quantum rotor model [31]. In this case, the influence functional becomes discrete and we will therefore refer to it as the influence matrix (IM). The setup and the key idea of the approach are summarized in Fig. 1. We consider a one-dimensional system of quantum spins, or qudits, σ k (with local Hilbert space dimension q). We choose spin p, and treat all degrees of freedom to the left as a bath, as illustrated in Fig. 1a. The starting point of our analysis is the Keldysh pathintegral formulation of time evolution. In this formalism, forward and backward trajectories of spins arise; for spin k they are denoted by σ τ k ,σ τ k , τ = 0, 1, 2, . . . . To capture the effect of the bath on spin p, we trace out the bath degrees of freedom, which gives rise to the influence matrix I({σ τ p ,σ τ p }). It enters as an additional weight into the path integral describing the evolution of the reduced density matrix of spin p. In general, this influence matrix is non-local in time.
For translationally invariant systems, spin p, subject to the IM I({σ τ p ,σ τ p }), should produce exactly the same IM for its right neighbor. Using this observation, we will formulate a self-consistency equation for the IM, pictorially illustrated in Fig. 1b. This equation is a key ingredient of our approach.
The knowledge of the IM allows to determine the local dynamical properties of the system, including all temporal correlation functions. The IM naturally incorporates the initial state, and allows averaging over ensembles of initial states. The self-consistent IM describes how a system acts as a bath on itself. This approach thus gives access to detailed information regarding thermalization time scales and memory time of the bath. As discussed below, it also allows one to analyze the effect of the bath on an impurity spin. In this sense, the self-consistent IM provides a more complete characterization of a manybody system as a bath, compared to spectral properties (such as presence or absence of level repulsion [7]), and the statistics of matrix elements studied both in the context of the eigenstate thermalization [3] and many-body localization [32].
While in general the self-consistency equation for the IM is complicated, it admits exact solutions in special cases. For kicked Ising models (see, e.g., Ref. [24] and references therein) with certain parameter values, we find that the infinite-temperature IM can be obtained exactly. This solution describes a bath that is a perfect dephaser (PD), that is, its effect on a spin at each step of the Floquet evolution is to exactly cancel off-diagonal matrix elements of its reduced density matrix. Phrased differently, the bath measures every spin of the system at each time step. Interestingly, for the case of kicked Ising models, the perfect dephaser class coincides with models that can be recast as dual-unitary circuits introduced recently [24,25]. While dual-unitarity implies PD form of the IM, we note that at present we do not know whether the converse holds.
Perfect dephaser systems serve as remarkably simple quantum baths: when coupled to an impurity spin, the system would act on it as an exactly Markovian bath, with a relaxation rate that depends on the coupling strength. This is remarkable, as Markovianity is usually an approximation which requires the internal dynamics of the bath to be much faster than the quantum system that it measures. Thus, the IM approach allows one to identify quantum systems which act as Markovian baths.
At the mathematical level, the IM approach bears a similarity to a tensor-network numerical method introduced by Bañuls et al. [19] for modelling dynamics of Hamiltonian systems. Building on the previous insights [19], we apply a matrix-product state (MPS) ansatz to construct the IM away from the PD points.
This tool allows us to shed light on the more general structure of the IM in ergodic Floquet systems.
At the PD points, the IM, viewed as a "wave function" in the space of single-spin trajectories, is effectively nonentangled. We further find that away from PD points, the IM "wave function" exhibits slow growth of temporal entanglement, which allows us to analyze the system's dynamics at longer times than those accessible via exact diagonalization. This observation provides a tool for identifying regimes of thermalizing Floquet dynamics that are amenable to efficient MPS-based methods.
To characterize the structure of IMs in ergodic systems detuned away from PD points, we adopt a statisticalmechanics-like description, viewing "quantum" intervals of a spin trajectory (i.e., the intervals where the forward and the backward path differ, σ =σ) as "particles". We study the weights of these particles and their interactions, demonstrating that in thermalizing systems they decay exponentially. We further use this insight to predict how the system thermalizes a slower impurity spin, finding a good agreement with numerical simulations.
The influence matrix approach has several additional attractive features. Perhaps most importantly, it provides a direct, physically intuitive way to describe a many-body system as a quantum bath for its constituent parts, giving access to relevant time scales and various correlation functions, including the Loschmidt echo (see below). Furthermore, it allows one to describe dynamics for different ensembles of initial states and in the thermodynamic limit. As mentioned above, the IM approach admits certain exact solutions that are likely not limited to PD which will be our focus here.
The rest of the paper is organized as follows. In Section II we introduce the IM approach, and formulate the self-consistency equation for the IM. In Section III we will discuss the cases where the IM can be found exactly, and has a perfect dephaser form. Further, Section IV is dedicated to dynamics away from PD points; we introduce and justify the use of MPS-based methods, develop an analytical characterization of the IM structure, and discuss implications for dynamics. Finally, in Section V we will summarize our results, and provide an outlook.

II. INFLUENCE MATRIX FORMULATION
We will start by introducing the models to be considered. We then describe a path integral representation of the dynamics. The correlation functions are expressed in terms of a transfer matrix which acts on the space of single-spin trajectories. We discuss general properties of such transfer matrices, arguing in particular that they have a pseudoprojection property. The eigenvector of the transfer matrix encodes the dynamical properties in the thermodynamic limit, and allows one to find all temporal correlation functions. We interpret the eigenvector as an influence matrix, and formulate a self-consistency equation for it.

A. Model
We consider a class of kicked one-dimensional Floquet systems. At each site of a periodic chain of length L, we place a spin, or qudit, with q basis states, denoted by |σ . During each driving period, a two-spin operatorP j+1/2 acts on all neighboring pairs j, j + 1. This operator represents an Ising-type coupling: it is diagonal in the |σ j , σ j+1 = |σ j ⊗ |σ j+1 basis, and symmetric under exchange j ↔ j + 1 : with φ(σ, σ ) = φ(σ , σ). This is followed by unitary single-spin operatorsŴ j ("kicks"). The combination of the two steps gives the following Floquet operator (evolution operator over one driving period): A graphical representation of the corresponding Floquet evolution is provided in Fig. 2. Here we choose to focus on this class of models for clarity of presentation. The following discussion, however, can be extended to more general Floquet systems, where a shallow quantum circuit is applied periodically in time. For q = 2, this model reduces to the kicked Ising model (KIM), which we will frequently use for illustration purposes.

B. Dual transfer matrix
Next, we use a discrete path-integral representation of the dynamics of the system. An amplitude for the evolution over t periods, from an initial product state |{σ 0 j } = ⊗ j |σ 0 j at τ = 0 to a state |{σ t j } at time τ = t, is given by This amplitude can be expressed by a path integral: The sum runs over all possible trajectories of L spins, with fixed initial and final configurations.
To describe the evolution of the density matrix, we will employ a Keldysh-type formalism. To that end, let us introduce a superoperator R which describes the evolution of the system's density matrix (DM). It is the tensor product of A and its conjugate A * , so that its matrix elements are given by  , the DM at time t is obtained by contracting R with ρ 0 . In particular, R gives the probability of the system's transition from an initial classical state to a final classical state, if we put σ 0 j =σ 0 j and σ t j =σ t j for all j. Further, following Ref. [19], we introduce a dual transfer matrix T j , which will provide a convenient alternative representation of the reduced density matrix evolution, with time and space interchanged (see Fig. ??). Its matrix elements are given by We treat initial and final spin configurations as parameters. The matrix acts on the q 2(t−1) -dimensional space of single-spin forward and backward trajectories at times 1 ≤ τ ≤ t − 1. Then, for a periodic system, Eq. (4) can be expressed via the dual transfer matrices at different sites as follows [33]: In the tensor-network language, this representation corresponds to contracting the network in the space direction, as discussed in Ref. [19] which introduced a new numerical method for Hamiltonian systems based on this idea. As we will see below, this representation has several advantages. It is particularly suited for describing ensembles of initial product states (or, more generally, weakly entangled states). For product states, the initial DM is a tensor product of individual spins DMs ρ 0 Figure 3. a) Graphical illustration of the pseudoprojection property in Eq. (11). The Keldysh path integral can be represented by a folded circuit and hence evaluated sequentially in the space direction, yielding the iterationT of the dual transfer matrix (blue shaded tensor) as in Eq. (9) (cf. Fig. 2). Unitarity of time-evolution allows the contraction of the network using the rules in the framed box on the left. Performing all possible contractions, one finds that for > t the input (left) and output (right) legs belong to disconnected networks, bounded by the upper and lower light cones. These two networks define vectors v| and |u , respectively. Hence, the result can be interpreted as the decomposition in Eq. (11). b) Graphical illustration of the eigenvector equation of the dual transfer matrix, which gives the self-consistency equation for the influence matrix of the system (see below).
In this case, to obtain ρ t we can conveniently contract T j with the corresponding initial DM of spin j. Moreover, if we are interested in the evolution of a given spin p, we should contract the final indices for all other spins, by putting σ t j =σ t j , and summing over them. This new dual transfer matrix, illustrated graphically in Fig. 2, can be expressed using Einstein notation for tensor contraction as Then, the evolution of spin p's DM is generated by the superoperator It is straightforward to adapt this formalism to different boundary conditions. Furthermore, for translationally invariant systems and initial states, all dual transfer matrices are the same,T j =T . To further simplify Eq. (9), we first discuss properties ofT .
Dual transfer matrices constructed as above share the following pseudoprojection property. Unitarity of time evolution implies that Tr(T L ) = 1 for all integer powers L. Thus,T has a single non-vanishing eigenvalue λ = 1. We denote by |u the corresponding right eigenvector, with components u {σ,σ} . Transposing the eigenvector equatioñ T |u = |u , we find that the associated left eigenvector v| has components All other eigenvalues are zero. As a result, iterations ofT eventually produce a projection onto |u . The system size dependence is encoded in the Jordan blocks ofT . Indeed, the strictly linear light-cone effect in such Floquet models graphically illustrated in Fig. 3, implies that the dimensions of Jordan blocks are upper bounded by * = 2t in general (and * = t for infinite-temperature initial ensembles, as in the Figure). Hence, for > * , we havẽ where the normalization is such that v|u = 1. The evolution superoperator R of the density matrix of spin p in an infinite homogeneous system [34], Eq. (9), can be conveniently expressed via the eigenvector |u as Knowledge of R allows us to compute all temporal correlation functions of a single spin. Thus, the problem of describing the dynamics of subsystems reduces to characterizing the properties of the eigenvector of the transfer matrixT . We remind the readers that this transfer matrix depends on the ensemble of initial states, see Eq. (8).

C. Self-consistency equation for the influence matrix
Let us take a closer look at the eigenvalue equation forT , and give it a transparent physical interpretation. We will assume translational invariance. We now denote the components of the eigenvector |u (in the space of trajectories of a spin on the Keldysh contour) via where {σ,σ} denotes a trajectory of a spin, i.e., σ 1 ,σ 1 , . . . , σ t−1 ,σ t−1 . Then, using Eq. (5), we can rewrite the eigenvalue equation as where W is the "folded" kick operator, which is the tensor product ofŴ acting along the forward path s τ , andŴ * acting along the backward paths τ .
Next, note that Eq. (14) has a simple physical origin. We can think of the r.-h.s. as describing the propagation in time of a spin {s,s}, which experiences on-site kicks W, and which is coupled to a neighboring spin {σ,σ}. Finally, the term I({s,s}) describes the effect of all the degrees of freedom to the left of our spin on its evolution. Viewing it now as a functional of the trajectory, rather than a vector, we can interpret it as the Feynman-Vernon influence functional (see Ref. [27] and Appendix A), or, rather, influence matrix, given the discreteness of time. In fact, the influence functional is obtained by tracing out the degrees of freedom to the left of a given spin p (see Fig. 3): This expression describes the effect of the environment, composed of all spins k < p, on the time evolution of spin p. Thus, Eq. (14) states that a spin, subject to an influence matrix I created by a bath to its left, creates an equal influence matrix for its neighboring spin to its right. This self-consistency equation for I is pictorially illustrated in Fig. 1b.
In more abstract terms, the influence matrix I({σ p ,σ p }) can be regarded as the overlap between the forward and the backward propagations of the environment formed by the spins k < p, subject to the trajectories {σ p }, {σ p } of the spin p, respectively: By unitarity of quantum evolution, one has |I({σ,σ})| ≤ 1. We will call trajectories where the forward and backward paths are identicalσ τ p = σ τ p classical trajectories since they are the equivalent of classical field configurations in the Keldysh formalism. For classical trajectories, one has I({σ, σ}) = 1.
Other general properties of the influence matrix can be obtained by extending the analysis of Ref. [27] to a discrete-time evolution, and are reported in Appendix A.
Although the self-consistency equation (14) encodes much of the complexity of a many-body system's dynamics, we will show below that there are cases when it can be solved analytically. Furthermore, we will show that the local relaxation dynamics can be naturally linked to a statistical-mechanics interpretation of the influence matrix elements I({σ,σ}).
Below we will be interested in averaging over the infinite-temperature ensemble.
Thus we put ρ 0 Other initial product states will lead to different transfer matrices. It is also possible to consider entangled initial states, at the expense of increasing the dimensionality of the transfer matrix.

III. PERFECT DEPHASERS
The process of thermalization starting from an initial product state is accompanied by the growth of entanglement between a spin and the rest of the system. In the language of the influence matrix, this corresponds to the suppression of non-classical paths, |I({σ =σ})| < 1. The exact form of this suppression encodes the dynamics of thermalization and the decay of correlation functions in a many-body system. In general, we may expect the IM to be a complicated functional that depends on the precise nature of the path {σ,σ}; parametrizing such a functional requires a number of parameters which is exponential in evolution time.
Surprisingly, there is a class of models for which the exact form of the IM is extremely simple: it vanishes exactly for all non-classical trajectories. This solution of the self-consistency equation (14) has a direct physical interpretation: The environment cancels out all the interference terms. Phrased differently, the environment completely dephases a spin at each evolution step, erasing the off-diagonal elements of its density matrix. Thus, we call such models perfect dephasers (PD).
Below we discuss examples of such solvable points for q = 2 (kicked Ising model of spins 1/2), and generalizations of KIMs to higher spins, q > 2. We find that for this family of models, the perfect dephaser points coincide with the self-dual points introduced by Akila et al. [24], and subsequently studied in Refs. [25,35,36]. Indeed, it is possible to show that dual-unitarity implies PD property, see Refs. [25,37] and Subsection C below. We, however, believe that the converse may not hold. Throughout this Section, we will consider infinite-temperature initial ensembles.
A. Kicked spin-1/2 Ising model First, we will consider the case of spin-1/2, q = 2. We will choose the basis |σ to be the eigenbasis of the z spin projection operator, such that σ = ±1. We will also employ the conventional Pauli matrix notations,σ x ,σ y ,σ z .
As the single-spin kick operatorŴ , we will choose a combination of a rotation around the z axis followed by a rotation around the x axis, Further, the two-spin term that depends on phases φ(σ, σ ), reduces to the Ising interaction with a coupling strength J: Thus, for the q = 2 case, our model is equivalent to the much studied kicked Ising model (KIM), which is known to display a variety of dynamical regimes, depending on the values of parameters h, J, . In particular, for h = 0 this model becomes solvable by mapping onto free fermions via Jordan-Wigner transformation. In general, when all three parameters are non-zero, the model is nonintegrable, and obeys the eigenstate thermalization hypothesis (ETH), and exhibits thermalizing dynamics [38]. Next, we assume that the influence matrix has a perfect dephaser (PD) form, Let us plug this PD influence matrix (IM) into the selfconsistency equation (14). We will see that for some special choices of the system's parameters , J, this IM indeed solves the equation. With this form of I, the summation in the r.-h.s. of Eq. (14) is performed only over classical trajectories, so one has to keep track just of s τ (sinces τ = s τ ). Then, the r.-h.s. can be rewritten using a one-dimensional transfer-matrix, composed from the matrix elements ofŴ and e iJs τ σ τ . The field h drops out and its value can be arbitrary. Then the selfconsistency equation takes the following form: where A = cos 2 sin 2 sin 2 cos 2 , corresponds to the infinite-temperature averaging. For = ±π/4, J = ±π/4, this equation is satisfied. To see this, note that the A matrix projects onto the vector 1 √ 2 (1 1) T , while the expectation value of the matrix B on this vector equals cos(2J(σ τ −σ τ )) = δ σ τστ . For any non-classical configuration, at least one of the factors will be zero. For classical configurations this expression gives 1, and therefore the self-consistency equation (20) is satisfied.

B. Higher spins
We now turn to the case q > 2, and show that extensions of the kicked Ising model also become perfect dephasers for a suitable choice of parameters. We identify PD points using a generalization of the approach outlined in the previous Subsection. As a result, we obtain examples of PDs with q > 2 which coincide with a family of dual-unitary models recently introduced by Gutkin et al. [36]. We emphasize that here our goal is to demonstrate that the IM approach allows one to construct examples of perfect dephasers, rather than to provide a complete classification of such solvable cases; this is left for future work.
Instead of a spin-1/2, we consider a clock variable σ = 1, 2, . . . , q, and an Ising-like spin-spin interaction, φ(σ, s) = Jσs. Let us assume the PD solution of the IM, and derive the suitable model parameters from the self-consistency equation. To that end, we can introduce the generalizations of the transfer matrices A, B from the previous Subsection, Eq. (21), which now become q × q matrices. The self-consistency equation will be satisfied by the PD IM, provided these transfer matrices satisfy the following properties: A αβ = A α β ∀α, β, α , β = 1, . . . , q.
These conditions indeed hold if we fix the parameter J = π/2q, and specify the single-qudit kickŴ to be of the following form, σ |Ŵ q |σ = 1 √ q e i σσ e ihσ , with = π/2q and arbitrary h. Then, exactly the same mechanism as for the KIM with q = 2 is effective at q > 2, yielding the PD influence matrix I({σ,σ}) = t−1 τ =1 δ σ τστ .

C. Relation to dual-unitary circuits
The above examples of perfect dephasers fall into the wider class of dual-unitary circuits, as recently highlighted in Akila et al. [24] and further characterized in Refs. [25,37,[39][40][41]. In fact, it is possible to show that the eigenvector of the dual transfer matrix has the perfect dephaser form in Eq.  2), the depicted contractions in the space direction are allowed. b) In this case, the network contraction proceeds within the light-cone, and the influence matrix reduces to a perfect dephaser form, i.e., a projection on classical paths. The effect on the spins on the right is that of a Markovian bath, such that coherences are cancelled at each evolution step. c) The influence matrix acting on an impurity spin coupled to a perfect dephaser generates strictly Markovian dynamics with a dephasing strength that is tunable by changing the coupling to the spin (see Sec. III D).
U γδ,αβ such that they maintain unitarity when reshaped to dual gatesŨ propagating along space direction, i.e., U βδ,αγ = U γδ,αβ [25]. The proof of this statement can be graphically carried out by tensor contractions, as shown in Fig. 4 with tensor notations adapted to our setting. The parameter values of the perfect dephaser points found in previous Section allow the "horizontal" contractions shown in panel a), which express the dual-unitary property of the gates. The result of these iterative contractions is then shown in panel b), and is nothing but the graphical representation of the PD influence matrix in Eq. (19).

D. A many-body system that is a Markovian bath
The knowledge of the influence matrix provides a complete characterization of a quantum bath and, in particular, gives a tool for describing its effect on dynamics of another quantum system coupled to it. In particular, one can analyze not just the dynamics of a translationally invariant system, but also the dynamics of an impurity immersed in the system. To illustrate this, we now study how a coupling to a perfect dephaser bath affects the evolution of a single spin.
We consider an impurity spin placed at site p, and coupled to its neighbors via a generic Ising coupling. We will see that the system acts on the impurity spin as a mem-oryless, Markovian bath. While the Markovian approximation is commonly employed in a range of problems, it normally relies on the separation of time scales between the system and the bath; in contrast, here it is exact, as a consequence of the PD property.
For simplicity, we focus on the kicked Ising model of spins 1/2, with an impurity spin at the edge, site p. We choose all couplings except those involving spin p to be = J = π/4, while spin p is coupled to spin p − 1 by an Ising coupling with possibly different strengthJ. The kick operatorŴ p may also be different from the π/2 rotation of the other spins. Then, using the PD property, we can derive the influence matrix for the impurity spin: This equation is most easily obtained graphically, noting that the last contraction on the right in Fig. 4b has now non-dual-unitary gates, producing a nontrivial but factorized IM. The latter is represented in Fig. 4c, where red tensors have a parameterJ that differs from that of the orange tensors. The IM in Eq. (25) gives rise to a time-independent superoperator ("quantum map" or "channel") that evolves the local reduced single-spin density matrix, where the dephasing superoperator D damps the offdiagonal entries by a factor cos(2J): Unless the couplingJ equals 0 or π, the damping is present. The corresponding thermalization time is dictated by the leading nontrivial eigenvalue of the superoperator in Eq. (26), which depends on both the spin's autonomous dynamicsŴ p and on the couplingJ. For instance, let us consider a kick operatorŴ p = exp(i˜ σ x ). We define the temporal correlation function whereF is the Floquet operator of the "chain + impurity" and · 0 denotes infinite-temperature averaging. The polarization decay rate γ eff of the impurity spin p, defined by the asymptotics depends on both˜ andJ. In particular, ifJ = π/4, we have e −γ eff = cos(2˜ ), while if˜ = π/4, we have e −γ eff = cos(2J).
The above discussion shows that the effects of perfectly dephasing baths can be computed for arbitrary times, and naturally raises the question about the structural properties of more general quantum ergodic baths possessing nontrivial temporal correlations. We tackle this intriguing question below.

IV. THERMALIZATION AWAY FROM PERFECT DEPHASER POINTS
The discussion above shows that certain fine-tuned Floquet systems can act as perfect dephasing baths for themselves. In other words, these systems can induce complete decoherence of the local degrees of freedom at each time step. The resulting thermalization dynamics is memoryless: i.e., it can be exactly described by a (discrete-time) Lindblad equation [42,43].
This property is surprising, as the Born-Markov approach to open system dynamics typically requires neglecting the bath temporal correlations [44]. This approximation usually relies on the separation of timescales between the system and the bath. While this is suitable in many physical situations involving the interaction of a "slow" particle with distinct, "fast" degrees of freedom (phonons, electromagnetic fields, . . . ), it is generally inadequate for the local dynamics of homogeneous extended quantum many-body systems.
The factorized form of the influence matrix in Eq. (19) of perfect dephasers, remarkably, makes their local dynamics exactly Markovian. Such property is tied with a complete absence of temporal entanglement in the influence matrix. It is thus natural to investigate the structure of the influence matrix in systems detuned from perfect dephaser points. One may expect that such systems would provide more generic examples of thermalizing quantum systems. In this Section, we shall focus on this problem.
In Sec. IV A we will introduce the matrix-product state (MPS) approach, which we adopt as a computational tool throughout this Section. In Sec. IV B we unveil the underlying structure of the influence matrix of generic thermalizing systems away from PD points. We motivate and validate a description of the IM, reminiscent of the statistical mechanics of massive, weakly interacting particles in one dimension, where the role of particles is played by intervals of a Keldysh trajectory with σ =σ. Finally, in Sec. IV C, we show that this approach can be successfully applied to compute the polarization decay rate of a slow impurity spin coupled to an ergodic chain.

A. Low temporal entanglement and matrix-product operator approach
For simplicity, we focus on the case of the spin-1/2 kicked Ising model described by Eq. as a function of detuning δJ = δ . As the system approaches the PD point, this slope goes to zero. We note that the slope remains small in a broad range of detuning, enabling efficient MPS description.
operators (17) and (18). In the following, it will be useful to "fold" the backwards and forwards contour, grouping together each spin σ τ on the forward time branch and its equal-time counterpartσ τ on the backward branch, to form a composite four-dimensional local Hilbert space (cf. Fig. 2). The influence matrix can be interpreted as a "wavefunction" living on a chain of those fourdimensional qudits. In this picture, the influence matrix of perfect dephasers [Eq. (19)] is represented by an exact product state. It is natural to assume that for sufficiently small detuning from the PD point, the influence matrix can be described in the folded picture by a matrixproduct state (MPS) with a moderate bond dimension. We set up a code based on the TeNPy library [45] which applies the dual transfer matrixT repeatedly, starting from a product state. Due to the pseudoprojection property discussed in Sec. II B, it is sufficient to apply the dual transfer matrix t times in order to obtain the influence matrix in the thermodynamic limit. The dual transfer matrix can be expressed as a matrix-product operator (MPO) with bond dimension 4. After each step, the influence matrix MPS is compressed using conventional SVD truncation sweeps. Up to time t = 2 log 2 χ the compression yields no truncation and the results from MPS match exact diagonalization to machine precision.
We take the convergence of the entanglement entropy S of this wavefunction upon increasing the bond dimension as a witness of the quality of our MPS representation. As shown in Fig. 5, the MPS approach allows us to explore the properties of the IMs for larger times t than those accessible via exact diagonalization, as S(t) converges for generic values of the parameters at reasonably low bond dimensions. The results indicate that the growth of S(t) is linear in t with a slope decreasing to zero as the PD point is approached. Accordingly, in the following, we will use the MPS approach described here as a computational tool.
The occurrence of low entanglement entropy in the folded picture is supported by an intuitive argument similar to that in Ref. [46]: In the absence of a longitudinal field h = 0, the system can be described by quasiparticles which move to the left on the forward and to the right on the backward branch. In the unfolded picture, this leads to a strong entanglement between sites on the two branches. At the self dual point, those quasiparticles can only propagate along the light cone edges [25], which means that in the folded picture correlations can only exist between a forward site and its corresponding backward site. Thus, the entanglement entropy of the folded MPS is zero. Detuning from the self dual point introduces a small density of slower quasiparticles, which gives rise to a parametrically slow growth of entanglement between different folded lattice sites.

B. Statistical-mechanics description of the influence matrix
The influence matrix is expected to develop a complex structure in generic systems away from perfect dephaser points, corresponding to the appearance of memory effects, intricate decoherence dynamics and temporal entanglement. To analyze this, it is convenient to introduce some additional formalism. Within the previously introduced folding map (cf. Fig. 2), we perform a discrete Keldysh rotation: We denote the "classical" configurations as and the "quantum" configurations as We further introduce the states Note that these four basis states may be seen as the vectorization (via the folding map) of the four basis op-erators1,σ z ,σ x , iσ y , respectively. With this notation, the initial infinite-temperature density matrices are represented as Figure 6. Illustration of the "sojourn-blip" parametrization in Eq. (38) of a generic spin trajectory on the Keldysh contour. and the perfect dephaser influence matrix assumes the following simple form: General properties of the influence matrix (see Appendix A for details) dictate that |I({σ,σ})| ≤ 1, and I({σ,σ}) = 1 for all classical trajectories, i.e., More generally, as discussed in the Appendix [see in particular Eq. (A5)], the value of I on mixed classical/quantum trajectories depends only on the configuration extending between the first and the last quantum site [47] I c ↑,↓ . . . c ↑,↓ (q ατ . . . q α τ )c ↑,↓ . . . c ↑,↓ = I(q ατ . . . q α τ ).
(36) This property motivates us to view classical Keldysh paths as "vacuum", and to interpret the quantum excursions (σ τστ ) = q ↑,↓ = (↑↓) or (↓↑) of a path {σ,σ} as "particles". In this description, the influence matrix is regarded as a complex statistical weight for a particle configuration, where we refer to S as the influence action. Note that Re S ≥ 0. In the limit of vanishing coupling to the bath, one trivially obtains S({σ,σ}) = 0 for all trajectories, irrespective of their particle content. In the opposite limit of a perfect dephaser, one has I({σ,σ}) = 0 [i.e., Re S({σ,σ}) = ∞] whenever a trajectory has some particles [see Eq. (19)]. In the generic case, a particle configuration will be penalized by a non-zero value of the action, with Re S({σ,σ}) > 0. Since the influence matrix contains all possible information on the dynamical effects of a part of the spin chain on the spins coupled to it, it is clear that the "interactions" encoded in the action S characterize the ergodicity of the quantum dynamics, or lack thereof.
This parameterization is illustrated in Fig. 6. (Here we assumed for simplicity that the initial density matrix is diagonal.) In the following, we will use a combination of analytical insights and numerical computations to show that in quantum ergodic many-body systems, blips behave as a "gas" of massive, short-range interacting particles.

Blip weights
Paths without blips, n = 0, are entirely classical, and correspond to the influence matrix I = 1. The simplest nontrivial configurations have n = 1, i.e., a single blip. Their influence matrix elements, hereby called influence weights, only depend on the internal blip structure, and not on the external sojourns.
To develop some intuition regarding the blip weights, let us first consider a constant blip, (a blip with constant q ↓ can be described in an analogous manner). In this case, it follows from the definition (16) [or, more generally, (A5)] that the influence matrix I reduces to the discrete-time, infinite-temperature version of the so-called Loschmidt echo, i.e., the overlap between two wavefunctions evolved from the same initial state with two slightly different Hamiltonians (see, e.g. Refs. [48,49] for reviews). Denoting by a subscript E the "environment" spins at positions 0 ≤ k < p, we have The two evolution operatorsÛ + E andÛ − E differ by a classical "field" e iφ(σp−1,±1) = e ±iJσp−1 acting on the boundary spin of the environment, due to the spin σ p being ↑ in the forward and ↓ in the backward evolution.
In the limit of weak coupling J → 0 between spins p−1 and p, the two time evolutions in Eq. (42) differ only slightly, and it makes sense to consider L as a measure of the sensitivity of the time-evolution to small perturbations in the Floquet operator. In fact, the continuoustime Loschmidt echo has been introduced as a quantifier of chaotic behavior in quantum systems [50,51]. It has been found that Loschmidt echoes generically exhibit an exponential decay in ergodic systems [52], while in MBL systems it displays a power-law behavior [53].
We are generally interested in the case of a finite, rather than weak, coupling strength J. Nevertheless, it is natural to expect that in ergodic systems weights of constant blips should decay exponentially upon increasing the blip size, similar to the Loschmidt echo. In Fig. 7 we report the results of numerical computations of the influence weight of individual blips with constant q ↑ in the kicked Ising chain, for a range of parameters in a neighborhood of the perfect-dephaser point. In all cases, a clear exponential decay is evident, consistent with our expectations.
Furthermore, the influence weight of non-constant blips may be viewed as a generalized time-dependent Loschmidt echo of the environment, formed by the overlap of two time-evolutions subject to a classical boundary field that differs at all times. It is natural to expect the influence weight of all blips to decay exponentially as a function of the blip duration in generic ergodic models. To show that this exponential decay is a generic feature of long blips, we next prove that the average weightĪ(q) of all blips of duration ∆τ also decays exponentially as (cos(2J)) ∆τ . We start with the self-consistency equation (19) for the IM, summing over all internal configurations of a blip of length ∆τ , which yields: The W transition amplitudes are antisymmetric with respect to the change of sign of the first sojourn. Since such a change does not affect the value of the influence matrix, the contributions of configurations with at least one blip cancel out in the average. This leaves just the contributions from purely classical trajectories, without blips (as in Sec. III D). For those trajectories, the influence matrix is always I({σ, σ}) = 1. Further, note that for all classical trajectories of σ the J dependent term in Eq. (43) is the same and equal (cos(2J)) ∆τ . Thus, the W transition amplitudes can be summed separately, which yields one. As a result, we get We have numerically verified this relation for various values of the model parameters, confirming that average weight of a single blip decays exponentially with its duration. Interestingly, as is evident in Fig. 7, the weight of a blip with a fixed internal structure, exhibits fluctuations around the approximate exponential decay.

Blip interactions
Within the above suggestive analogy between an IM and a statistical-mechanical complex "action" for blips, a long blip q may be pictured as a structured aggregate of particles, the "mass" S(q) of which increases approximately proportionally to its length. The disconnected-blip contribution is isolated by choosing the identity state 1 2 |c+ c+| in the resolution of the identity on every leg, and is highlighted by the red shading. In this case the exterior network equals 1, as it can be completely contracted via the rules in Fig. 3a. In contrast, for all the other contributions, the exterior network is equivalent to a temporal correlation function of traceless local operators, corresponding to states |c− , |q+ , |q− , denoted Oi here. In ergodic dynamics, these correlations generically decay as the temporal separation ∆t1 or ∆t2 increase.
While the weight of long blips is strongly suppressed, this need not be the case for configurations with multiple blips separated by long sojourns. It is natural to define the interaction influence weight of a multiple-blip configuration as the excess influence weight compared to the product of the influence weights of the individual blips: e −Sint(q1,c1,q2,c2,...,qn) = e −S(q1,c1,q2,c2,...,qn)+ n i=1 S(qi) (45) Note that S int depends on the specific configuration of the sojourns c j between the blips. The question of the range of the blip interactions is intimately related to the memory time of the system. It is expected that a weak non-Markovianity of the reduced dynamics can be translated into a fast decay of interactions.
In fact, the decay of the blip interactions may be precisely related to the temporal decay of the correlations of certain local operators, via the following argument. Let us consider a spin trajectory with n blips, {|q i } i=1,...,n , and let us evaluate its IM element by contracting the corresponding tensor network, as shown in the example in Fig. 8  the Figure we average over the configurations of the intermediate sojourns, although this is not necessary for the argument that follows.) To isolate the disconnected blip contributions, we insert a spectral resolution of the identity formed by the states |c + , |c − , |q + , |q − defined in Eqs. (30), (31), on each four-dimensional tensor leg connecting a blip with the rest of the network, as indicated by the arrows in the Figure. In this way, we have formally decomposed the IM element as a sum of O(4 n ) contributions, each given by the product of n + 1 disconnected network contractions: that is, n "small" singleblip networks including |q 1 , . . . , |q n respectively, and the remaining "large" exterior network comprising all sojourns. In this sum, the term where we choose 1 2 |c + c + | on every leg (which is explicitly represented in the Figure) yields exactly the product of the individual blip weights (red shaded), because the exterior network is equivalent to the influence weight of a classical trajectory and thus evaluates to 1 (using the contraction rules in Fig. 3a). For all the other contributions, the large exterior network may be viewed as a (possibly complicated) temporal correlation function of local operatorsÔ 1 , . . . ,Ô n , all with spatial support near the boundary of the chain, and with temporal support near the respective blips. In quantum ergodic systems, such temporal correlations are expected to decay exponentially upon increasing the time separations ∆t 1 , . . . , ∆t n−1 between the operators (i.e., the duration of the intermediate sojourns). Thus, in general, one expects the interaction strength between clusters of blips to vanish for sufficiently long sojourns in between: i.e., if ∆t m → ∞, e −S(q1,c1,q2,c2,...,qn) e −S(q1,c1,...,qm) × e −S(qm+1,cm+1,...,qn) . (46) This scenario is supported by our numerical computations. First, we verify that temporal correlations of local operators decay upon increasing the time separation. For the sake of illustration, we report in Fig. 9 the dynamics of the spin polarization autocorrelation C zz , defined in Eq. (28). The decay occurs for all values of the parameters, as expected in generic ergodic systems. By the above arguments, the generic decay of temporal correlations such as the ones we reported in Fig. 9 justifies the cluster approximation of the blip action. We further verify this directly by inspecting the IM elements. In Fig. 10 we report the influence weight of a pair of blips for the kicked Ising chain, averaged over the configurations of the intermediate sojourns, as a function of their temporal separation, for two examples of blip structures, and for a range of parameters in a neighborhood of the perfect-dephaser point. As the plots clearly show, in all cases the IM approaches a constant plateau equal to the product of the disconnected blip weights.

C. Decay rate of a slow impurity
In the previous Subsection, we analyzed the influence action, finding that blips have a finite range of interactions, which is related to the relaxation time of the system. Next, we analyze the effect of a many-body bath on a "probe" spin embedded a spin chain, whose own dynamics is slower than that of the environment. In this case, thanks to the separation of time scales, the blip "gas" is effectively dilute and can be treated as noninteracting. This will enable us to predict the relaxation rate of the impurity spin.
To make this intuition quantitative, we compute the polarization decay rate of an impurity spin at the edge of a kicked Ising spin chain defined by Eqs. (2), (17), (18), with parameters = π/4 − δ , J = π/4 − δJ detuned from a perfect dephaser point. This impurity spin differs from the other spins of the chain by the smaller strength of its kicks, Our aim is to compute the persistence of the impurity spin polarization, where the autocorrelation function C zz (t) has been previously defined in Eq. (28). We can express this transition probability as a sum over Keldysh trajectories of the impurity spin. Denoting the path variables by σ,σ, we can write the above equation as follows, Switching to the sojourn-blip parameterization in Eq. (38), see Fig. 6, we can rewrite Eq. (49) as follows: Now the sum over paths is arranged as a summation over the number of blips, over their positions in time, and over the internal configurations of individual sojourns and blips. In this equation, we have combined the intrasojourn and intra-blip kick transition amplitudes into the terms where ) is the number of domain walls in the sojourn or blip configuration, and we have defined the total magnetization a blip, (53) The above equations show that for weak kicks |˜ | 1 of the impurity spin, the sum over trajectories will be dominated by terms with a low density of blips. In fact, the kick strength˜ acts as a "chemical potential" for blips. Thus, in most of the relevant blip configurations, we may treat the blips as noninteracting, i.e., to substitute The success of this non-interacting blip approximation (NIBA), discussed at length by Leggett et al. in the context of a two-level system coupled to a bath of independent harmonic oscillators [28], relies on the typical interblip distance being larger than the blip interaction range.
The factorized form of the IM in Eq. (54) allows us to analytically compute the polarization decay rate of the impurity. After this simplification, the influence matrix in Eq. (50) does not depend on the sojourns' internal structure. The summation over all configurations of one sojourn yields c:|c|=∆t Note that the first and last spin are restricted to being σ 0 = σ t = ↑, which gives an extra overall factor of 1/4; moreover, the term n = 0 has both ends fixed and gives Since we are interested in the asymptotic decay rate of P ↑↑ (t), it is convenient to compute its generating function (or a discrete Laplace transform), We split t into intervals according to Eq. (40), and, for each n, we perform the sums over the durations of blips and sojourns. For the blips, we define Putting everything together, we find The long-time behavior of the polarization is governed by the asymptotics ofP ↑↑ (λ) for small argument. In particular, an exponential decay of the correlation function C zz (t) in Eq. (28) with a rate γ eff gives:  For small˜ , this reduces, in the lowest non-vanishing order in˜ , to This equation expresses the decay rate of a slow impurity spin coupled to a spin chain detuned from the PD point, which is a non-Markovian bath. In the limit of a perfect dephasing bath with | | = |J| = π/4, where the blips are completely suppressed, the above equation reduces to the Markovian dynamics, as in Eq. (29). In fact, the first term in Eq. (63) represents the decay rate in this limit. This contribution is purely classical, as the environment completely suppresses interference between all pairs of quantum trajectories. The second term represents instead the noninteracting-blip contribution. This correction can be of either sign, due to the complex blip weights.
Blip interactions generate higher-order corrections to γ eff . Leaving their detailed analysis for future work, we note that such corrections can be systematically taken into account via a renormalization procedure, in which one progressively includes connected contributions of clusters of two, three, etc blips into I(∆τ ). Such a scheme can be truncated as long as the series (50) is dominated by terms with a low density n/t of blips. For small˜ , this holds, since the blip density scales as O(˜ 2 ). The first correction arising from blip interactions is suppressed as r˜ 2 , where r is an effective parameter describing the range of two-blip interactions. Note that the value of r may depend on the structure of the blips involved and model parameters, see e.g. Fig. 10. This will lead to a correction of the order O(r˜ 4 ) to the relaxation rate γ eff in Eq. (63). For fixed˜ , this correction becomes increasingly important as one detunes the system from a perfect dephaser point, effectively enhancing the blip interaction range r.
We tested the predictions of Eq. (63) against our numerical computations, finding a good agreement in a broad range of model parameters. The comparison is shown in Fig. 11. In the left panel, the dynamics of the probe spin polarization autocorrelation C zz (t) is plotted for increasing detuning δ = δJ of the chain from the perfect dephaser point, where the behavior is exactly exponential, C zz (t) = cos t (2˜ ). As shown, the decay of C zz (t) remains approximately exponential even for sizeable detunings. In the right panel, we compare the measured decay rate γ with the prediction of the NIBA in Eq. (62). For small detunings the approximation is excellent. A small discrepancy appears for larger detunings, which we attribute to the neglected contribution of blip interactions. The latter could be accounted for via the renormalization procedure briefly sketched above. We note that we found the range of quantitative validity of the NIBA to be sensitive to the value of the integrabilitybreaking parameter h. This can be attributed to the fact that, as we saw in the previous subsection, changing h tunes the range of blip interactions.
Finally, we briefly comment on the difference between the discrete spin dynamics studied above, and continuous evolution of a two-level system coupled to a bath of harmonic oscillators [28]. In the"strong decoherence limit", when blips are completely suppressed, the dynamics in the two cases is qualitatively different. In the continuous setting, a complete suppression of quantum interference resulting from strong interaction with the environment freezes the spin state and it cannot relax to equilibrium. This phenomenon is known as the quantum Zeno effect [54]. In the discrete case, blips are suppressed for˜ = π/4. Then, in contrast, the spin thermalizes over just one time step, as discussed in Sec. III D.

V. SUMMARY AND OUTLOOK
We have developed an approach to analyzing highly non-equilibrium dynamics of isolated many-body systems, inspired by the Feynman-Vernon influence functional formalism [27]. Focusing on a class of interacting Floquet models, we formulated the self-consistency equation for the influence matrix, and demonstrated that it can be analyzed using complementary analytical and numerical considerations in a whole range of model parameters. Therefore, both "solvable" perfect dephaser baths that are Markovian, and more generic non-Markovian models can be characterized within the IM approach. As an example, we analyzed the effect of a generic manybody system on an impurity spin coupled to it.
The influence matrix provides a novel probe for characterizing quantum chaos and its absence in many-body systems. Compared to spectral probes, such as level statistics, the advantage of the IM is that it provides detailed information regarding the memory time scales and temporal correlation functions. IMs also contain the information expressed by other previously studied dynamical probes, such as the Loschmidt echo, and the statistics of matrix elements of local operators. In future work, it will be interesting to develop a connection between the properties of the IM and the eigenstate thermalization hypothesis.
We have found a family of ergodic models (in particular, in the vicinity of the perfect dephaser points), where MPS-based methods are efficient for computing the IM. This stems from viewing the IM as a "wave function", whose temporal entanglement, as we showed, grows parametrically slowly as a function of the detuning from PD points. One may expect that this property would hold for a broad family of models that thermalize quickly, enabling an efficient MPS treatment of their dynamics. This insight complements previous works on Hamiltonian systems [19,46].
Our work also suggests several promising future directions. First, it appears that PD circuits can be found in higher dimensions. Second, it would be interesting to analyze the precise relation between PD and dual-unitary circuits. As we discussed (see also Refs. [25,37]), dualunitarity implies that a system is a PD, but we do not know whether the converse holds. Other promising generalizations of our approach include Hamiltonian systems and their dynamics at a finite, rather than infinite, tem-perature.
In future work, we plan to further analyze the interacting blip gas, and the effects of blip interactions on the impurity spin dynamics, as briefly discussed in the last Section. We envision that, detuning from PD points, the self-consistency equation for the IM may be solved perturbatively in the detuning. This, along with applying the IM approach to non-ergodic systems [55], will likely lead to a more complete characterization of nonequilibrium quantum matter.

VI. ACKNOWLEDGMENTS
This work was supported by the Swiss National Science Foundation. We thank Soonwon Choi and Ehud Altman for inspiring discussions.

Appendix A: General properties of influence matrices
In this Appendix we review the basic properties of Feynman-Vernon's influence functional, adapted from Ref. [27] to our discrete-time dynamics. For simplicity, we focus on two-level systems (q = 2); the generalization to the case of q > 2 is straightforward.
We consider a quantum spin (S) interacting with an environment (E). In the main text discussion, S is a spin at position p in a spin chain, and E is formed by all spins on previous positions k < p; here, however, for the sake of generality we keep the discussion more abstract. We assume the spin-environment interaction to be of the formV int =V S ⊗V E . Without loss of generality, we takê V S = Jσ z , such that U int = e i(Jσ z ⊗V E ) . (A1) The global Floquet operator has thus the form whereÛ S ,Û E are the parts of the Floquet operator acting on the spin and environment separately. We further assume that in the initial state such that the spin and environment are uncorrelated. According to the standard laws of quantum mechanics, the transition probability P σ 0 →σ t for the spin S can be expressed as a discrete path integral [56].