Entanglement spreading in a minimal model of maximal many-body quantum chaos

TL;DR

This paper provides an exact, rigorous method to compute entanglement dynamics in a class of maximally chaotic, periodically driven quantum spin chains, revealing universal features of entanglement spreading in such systems.

Contribution

It introduces a novel duality-based approach to exactly calculate entanglement growth in maximally chaotic quantum models, bridging different theoretical frameworks.

Findings

Exact entanglement entropy evolution in self-dual kicked Ising chains.

Results match predictions from random matrix theory and conformal field theory.

Entanglement spreading saturates the minimal cut bound, consistent with universal models.

Abstract

The spreading of entanglement in out-of-equilibrium quantum systems is currently at the centre of intense interdisciplinary research efforts involving communities with interests ranging from holography to quantum information. Here we provide a constructive and mathematically rigorous method to compute the entanglement dynamics in a class of "maximally chaotic", periodically driven, quantum spin chains. Specifically, we consider the so called "self-dual" kicked Ising chains initialised in a class of separable states and devise a method to compute exactly the time evolution of the entanglement entropies of finite blocks of spins in the thermodynamic limit. Remarkably, these exact results are obtained despite the models considered are maximally chaotic: their spectral correlations are described by the circular orthogonal ensemble of random matrices on all scales. Our results saturate the so called "minimal cut" bound and are in agreement with those found in the contexts of random unitary circuits with infinite-dimensional local Hilbert space and conformal field theory. In particular, they agree with the expectations from both the quasiparticle picture, which accounts for the entanglement spreading in integrable models, and the minimal membrane picture, recently proposed to describe the entanglement growth in generic systems. Based on a novel "duality-based" numerical method, we argue that our results describe the entanglement spreading from any product state at the leading order in time when the model is non-integrable.

Equations

$$H_{\mathrm{KI}}[\mathbf{h};t]=H_{\mathrm{I}}[\mathbf{h}]+\sum _{m=-\infty }^{\infty }\delta (t-m)H_{\mathrm{K}},$$

$$H_{\mathrm{I}}[\mathbf{h}]\equiv J\sum _{j=1}^L{\sigma }_j^z{\sigma }_{j+1}^z+\sum _{j=1}^L{h}_j{\sigma }_j^z,$$

$$H_{\mathrm{K}}\equiv b\sum _{j=1}^L{\sigma }_j^x.$$

$$|J|=|b|=\frac{\pi }{4}.$$

$$J=\frac{\pi }{4},b=-\frac{\pi }{4},$$

$${\sigma }_j^a\mapsto T{\sigma }_j^{a}T={\sigma }_j^{a\ast }.$$

$${\mathcal{B}}_L=\left\{|\mathbf{s}⟩=|s_1,\dots ,s_L⟩,s_j\in \{\pm 1\}:{\sigma }_j^z|\mathbf{s}⟩=s_j|\mathbf{s}⟩\right\}.$$

$$U_{\mathrm{KI}}[\mathbf{h}]=\mathfrak{T}\exp \left[-i{\int }_0^1\mathrm{d}t{H}_{\mathrm{KI}}[\mathbf{h};t]\right]=e^{-i{H}_{\mathrm{K}}}{e}^{-i{H}_{\mathrm{I}}[\mathbf{h}]},$$

$$e^{-i{H}_{\mathrm{I}}[\mathbf{h}]}=\prod _{i\in \mathrm{odd}}{e}^{-i\frac{\pi }{4}{\sigma }_j^z{\sigma }_{j+1}^z}\prod _{i\in \mathrm{even}}{e}^{-i\frac{\pi }{4}{\sigma }_j^z{\sigma }_{j+1}^z}.$$

$$H_{\mathbf{\theta },\mathbf{\varphi }}=-\sum _{j=1}^L{g}_j{\vec{n}}_{{\theta }_j,{\varphi }_j}\cdot {\vec{\sigma }}_j,$$

$${\vec{n}}_{\theta ,\varphi }=\left(\sin \theta \cos \varphi ,\sin \theta \sin \varphi ,\cos \theta \right),$$

$$|{\psi }_{\mathbf{\theta },\mathbf{\varphi }}⟩=\underset{j=1}{\overset{L}{⨂}}\left[\cos \left(\frac{{\theta }_j}{2}\right)|\uparrow ⟩+\sin \left(\frac{{\theta }_j}{2}\right)e^{i{\varphi }_j}|\downarrow ⟩\right].$$

$$|{\psi }_{\mathbf{\theta },\mathbf{\varphi }}(t)⟩=(U_{\mathrm{KI}}[\mathbf{h}]{)}^t|{\psi }_{\mathbf{\theta },\mathbf{\varphi }}⟩.$$

$$\rho_{A}(t)={\rm tr}_{\mathcal{H}_{L-N}}\!\bigl{[}\ket{\psi_{\bm{\theta},\bm{\phi}}(t)}\bra{\psi_{\bm{\theta},\bm{\phi}}(t)}\bigr{]}\,,$$

$$S_A^{(\alpha )}(t)=\frac{1}{1-\alpha }\log \mathrm{tr}\left[({\rho }_A(t){)}^{\alpha }\right],\alpha >0.$$

$$S_A^{(\mathrm{vN})}(t)=\underset{\alpha \to 1}{\lim }{S}_A^{(\alpha )}(t)=-\log \mathrm{tr}\left[{\rho }_A(t)\log {\rho }_A(t)\right],$$

$$\mathcal{T}=\left\{|{\psi }_{\frac{\pi }{2}1,\mathbf{\varphi }}⟩,{\varphi }_j\in [0,2\pi ]\right\},$$

$$\mathcal{L}=\left\{|{\psi }_{\bar{\mathbf{\theta }},\mathbf{\varphi }}⟩,{\bar{\theta }}_j\in \{0,\pi \}\right\},$$

$$|{\psi }_{\bar{\mathbf{\theta }},\mathbf{\varphi }}(1)⟩=U_{\mathrm{KI}}[\mathbf{h}]|{\psi }_{\bar{\mathbf{\theta }},\mathbf{\varphi }}⟩\simeq |{\psi }_{\frac{\pi }{2}1,\bar{\mathbf{\theta }}-\frac{\pi }{2}1}⟩,$$

$$e^{-i{H}_{\mathrm{K}}}|{\psi }_{\bar{\mathbf{\theta }},\mathbf{\varphi }}⟩\simeq |{\psi }_{\frac{\pi }{2}1,\bar{\mathbf{\theta }}-\frac{\pi }{2}1}⟩.$$

$$|{\psi }_{\bar{\mathbf{\theta }},\mathbf{\varphi }}(t)⟩\simeq |{\psi }_{\frac{\pi }{2}1,\bar{\mathbf{\theta }}-\frac{\pi }{2}1}(t-1)⟩,t\ge 1.$$

$$S_A^{(\alpha )}(t){|}_{\bar{\mathbf{\theta }},\mathbf{\varphi }}=S_A^{(\alpha )}(\max (t-1,0)){|}_{\frac{\pi }{2}1,\bar{\mathbf{\theta }}-\frac{\pi }{2}1},t\ge 0,$$

$$\underset{L\to \infty }{\lim }{S}_A^{(\alpha )}(t)=\min (2t,N)\log 2.$$

$$S_A^{(\alpha )}(t)\le {\ell }_{\min }\log q,\forall \alpha ,$$

$$\underset{L\to \infty }{\lim }{S}_A^{(\alpha )}(t)=\min (2t,N)S_{\mathbf{\theta },\mathbf{\varphi }}^{(\alpha )},$$

$$\mathrm{tr}\left[(U_{\mathrm{KI}}[\mathbf{h}]{)}^t\right],t\in \mathbb{N},$$

$$\mathrm{tr}\left[(U_{\mathrm{KI}}[\mathbf{h}]{)}^t\right]=\mathrm{tr}\left[{\tilde{U}}_{\mathrm{KI}}[h_1\tilde{1}]\cdots {\tilde{U}}_{\mathrm{KI}}[h_L\tilde{1}]\right],$$

$$\tilde{1}=(\underset{t}{\underbrace{1,\dots ,1}}),$$

$$\tilde{J}$$

$$\tilde{b}$$

$$\mathrm{tr}\left[({\rho }_A(t){)}^n\right]=$$

$$\sum _{\{{\mathbf{a}}_i\},\{{\mathbf{b}}_i\}}⟨{\psi }_{\mathbf{\theta },\mathbf{\varphi }}|(U_{\mathrm{KI}}[\mathbf{h}]{)}^{-t}|{\mathbf{a}}_1,{\mathbf{b}}_2⟩⟨{\mathbf{a}}_1,{\mathbf{b}}_1|(U_{\mathrm{KI}}[\mathbf{h}]{)}^t|{\psi }_{\mathbf{\theta },\mathbf{\varphi }}⟩$$

Full text

1. I Introduction

2. II Model, quench protocol, and observables of interest

3. III Outline of the results

4. IV Duality Mapping for the Entanglement Entropies

5. V Separating states

6. VI Entanglement spreading from separating states

7. VI.1 Maximal eigenvalues of the transfer matrix

8. VI.2 Entanglement dynamics

9. VII Entanglement spreading from generic states

10. VIII Conclusions

11. A Duality of traces

12. B Simplified transfer matrix for longitudinal separating states

13. C Proof of Property 1

14. D Simplified commutation relations

15. E Proof of Eqs. (86) and (87)

16. F Proof of Property (84)

17. F.1 Proof of Lemma 1.

I Introduction

Entanglement is arguably the most distinctive feature of quantum mechanics. It generates a special kind of non-local correlations which can be present in quantum states but have no analogues in the classical realm. While its elusive nature puzzled physicists for many years EPR ; Bell , it is currently regarded as a powerful resource for advances both in technological applications and in the theoretical understanding of the physical world. In particular, it plays a crucial role in the study of quantum many-body systems out of equilibrium Entrev1 ; Entrev2 . This is due to two main reasons. First, the growth of entanglement during the non-equilibrium dynamics measures the increasing complexity of a time-evolving quantum state, with immediate implications on the feasibility of tensor network simulations SWVC:PRL ; SWVC:NJP ; PV ; HCTDL ; Dubail . Second, the evolution of the entanglement gives crucial information on how equilibrium statistical mechanics emerges from many-body quantum dynamics. Specifically, it is now understood that the thermodynamic entropy of the statistical ensemble describing local observables at infinte times is a measure of the entanglement accumulated during the time evolution p-18 ; dls-13 ; bam-15 ; Gur14 ; SPR11 .

Moreover, the very way in which the entanglement spreads for finite times appears to be among the most universal aspects of many-body dynamics. Consider for instance an initial separable state, where none of the local constituents is entangled with any other. Switching on spatially local Hamiltonian interactions throughout the system (a procedure called “global quench”), one finds quite generally that the bipartite entanglement between a large connected region of the system and the rest grows linearly in time. This scenario has been observed in a large number of analytical and numerical investigations CC ; DeChiara ; leda-2014 ; FC:exactXY ; ep-08 ; nr-14 ; bhy-17 ; hbmr-17 ; lauchli-2008 ; KH:NonIntEnt ; fc-15 ; buyskikh-2016 ; cotler-2016 ; Alba ; mbpc-17 ; mkz-17 ; FNR:longrangehigherd ; ckt-18 ; BTC:beyondpairs ; BC:beyondpairs2 ; CLM:minimalcut ; LS:CFT ; alba-inh ; BFPC18 ; daley-2012 ; evdcz-18 ; PL:kickedIsing ; ABGH:CFT ; LM:CFT and recently even in cold atomic experiments expS . Known exceptions to this empirical fact are systems exhibiting real space localization DeChiara ; ZPP , confinement KCTC:confinement , and quenched disorder creating weak links Nahum3:disordered . In particular, the logarithmic spreading of entanglement in the presence of many-body localization (MBL) MBL is one of the main defining features of the MBL phase.

The linear growth of entanglement after a global quench has been first observed in the context of $(1+1)$-dimensional conformal field theory (CFT), where it has been explained in terms of an intuitive quasiparticle picture CC . The initial state, which is not an eigenstate of the Hamiltonian, can be thought of as a collection of pairs of oppositely moving quasiparticles. These, in the course of time, spread the entanglement across the system, in a similar way to the one conceived in the historical gedanken experiment by Einstein, Podolsky and Rosen EPR . In this picture, the entanglement between two different portions of the system is given by the number of pairs sharing a particle with each portion. The same idea can be used to explain the entanglement spreading in systems with stable quasiparticle excitations, for instance free FC:exactXY and interacting Alba integrable models. It does not account, however, for the linear growth of entanglement observed in systems with no identifiable quasiparticle content such as generic interacting systems KH:NonIntEnt ; lauchli-2008 ; PL:kickedIsing or holographic CFTs CLM:minimalcut ; LS:CFT ; ABGH:CFT ; LM:CFT .

More recently, a fruitful avenue of research came from the study of the so-called random quantum circuits RandomCircuitsEnt ; r-2017 ; Nahum2 ; Nahum3:disordered ; Nahum4 ; Nahum5 ; CDC1 , where the dynamics is completely random in space and the only constraint is given by the locality of interactions. In this case the linear growth of entanglement can be explained using a “minimal membrane” picture RandomCircuitsEnt ; Nahum4 , which is conjectured to apply, at least at the qualitative level, to generic (non-integrable) clean and noisy systems in any spatial dimension. In essence one quantifies the amount of entanglement between two portions of the system by measuring the surface of the minimal space-time membrane separating the two portions. This picture has been analytically tested in certain limiting regimes of random quantum circuits, specifically assuming that the Hilbert space dimension $q$ per local constituent is large ($q\gg 1$). Results are available both when the dynamics are random also in time RandomCircuitsEnt ; Nahum5 , and when they are periodically driven CDC1 . No analytical result, however, exists on entanglement dynamics in specific non-integrable models with local interactions and small finite $q$ or, in general, for clean systems.

In this paper we fill this gap providing exact results for the entanglement dynamics of quantum-chaotic spin chains with two-dimensional local Hilbert space ($q=2$). Specifically, we consider a family of Floquet-Ising chains which undergo a transition between integrability and ergodicity (or quantum chaos) by turning on a longitudinal magnetic field. The latter may be either spatially homogeneous or arbitrarily spatially modulated. Note that the non-integrability of the model for non vanishing longitudinal magnetic fields has recently been proved by computing exactly the spectral statistics letter . We identify a class of separable initial states, homogeneous or arbitrarily modulated in space, from which the entanglement dynamics can be computed exactly for any non-disjoint bipartition. These results are of high significance for three main reasons. (i) They provide an exact verification of both the linear growth of entanglement and its relaxation to the thermodynamic entropy in concrete quantum-chaotic models. (ii) They provide a general method allowing one to obtain exact results for the non-equilibrium dynamics of many-body quantum systems even in the absence of integrability. (iii) They are valid in both the integrable and the non-integrable case, allowing for a unified interpretation of the entanglement spreading.

The paper is laid out as follows. In Sec. II we present the model considered, define the protocol used to drive it out of equilibrium, and introduce the entanglement measures of interest. In Sec. III we present a comprehensive summary and discussion of our results. In Sec. IV we explain the duality mapping which is the key for our analytical calculations. In Sec. V we identify the classes of initial states leading to an exactly solvable entanglement dynamics. In Sec. VI we sketch the main steps of the calculation. In Sec. VII we present a thorough numerical analysis of the entanglement spreading from generic separable initial states and advocate that our exact results give the leading-order in time description of the non-integrable case. Finally, Sec. VIII contains our conclusions. A number of technical points and proofs are reported in the appendices.

II Model, quench protocol, and observables of interest

The main objective of this paper is to determine a minimal quantum chaotic model note2 , with local interactions and finite local Hilbert space, allowing for an exact determination of the entanglement spreading. A candidate emerging naturally in our quest is the so called kicked Ising model KI_JPA ; KI_Ruelle ; KI_PRE , which describes a classical Ising model in the presence of a longitudinal magnetic field and periodically kicked with a transverse magnetic field. This model is quantum chaotic in the sense that its spectral statistics are described by the circular orthogonal ensemble of random matrices PP-KI , but, as we recently proved letter , at some specific points of its parameter space it allows for exact calculations. This is because at these points, called “self-dual” points (see below), the system acquires a remarkable algebraic structure making of it a maximal scrambler with local interactions.

To be more specific, let us introduce the Hamiltonian of the kicked Ising model. Setting to one the time interval between the kicks we have

[TABLE]

where $\delta(t)$ is the Dirac delta function and we defined

[TABLE]

In these equations $L$ is the volume of the system, the matrices $\sigma_{j}^{a}$, with ${a\in\{x,y,z\}}$, are the Pauli matrices at position $j$, and we use periodic boundary conditions adopting the notation convention $\sigma_{L+1}^{a}\equiv\sigma_{1}^{a}$.

The parameters $J$ and $b$ are, respectively, the Ising coupling and the transverse “kicking” field, while the $L$-component vector $\bm{h}=(h_{1},\ldots,h_{L})$ describes a position dependent longitudinal field. As anticipated before, in this paper we consider some specific points in the parameter space of the model. In particular we focus on the “self-dual” points, specified by the condition

[TABLE]

In Secs. IV, V, and VI we explain how, at these points, a duality symmetry of the model allows for an analytical treatment of the entanglement dynamics. To be specific, from now on we set

[TABLE]

but our results apply to all the four combinations fulfilling (4). The longitudinal magnetic fields ${\bm{h}}$ are instead left arbitrary and are used to switch between the integrable and the non-integrable case. Indeed, for ${\bm{h}=\bm{0}}$ the Hamiltonian (1) is integrable (it can be mapped to a problem of non-interacting fermions), while it is ergodic (non-integrable) for a generic choice of longitudinal fields. In the latter case the only symmetry of (1) is time reversal. This is represented by the anti-unitary operator $T$, defined by its action on the spin variables as follows

[TABLE]

Here $(\cdot)^{*}$ denotes the complex conjugation in the “computational” basis, composed by simultaneous eigenstates of $\{\sigma_{j}^{z}\}$ for all $j$ in $\{1,2,\ldots,L\}$

[TABLE]

In this paper we interchangeably use ${s=+1\equiv\,\uparrow}$ and ${s=-1\equiv\,\downarrow}$ to designate eigenvalues of Pauli matrices. We stress that, as proven in Ref. letter , the self-dual kicked Ising model is ergodic for any $\bm{h}\neq 0$. This means that, due to its special structure, the model never displays Floquet-MBL MBL ; Papic ; Roderich . Note that the special structure of the self-dual model has recently been used to design a multiparty entanglement generation algorithm NSM:multipartyentanglement .

The Floquet operator generated by (1) reads as

[TABLE]

where ${\mathfrak{T}}\!\exp\left[\cdot\right]$ denotes a time-ordered exponential. The time-evlution generated by (13) admits a straightforward local 2-qubit quantum circuit representation. This is observed by writing

[TABLE]

To drive the system out of equilibrium we consider a global quantum quench protocol: we initialise the system in the ground state of a short-range Hamiltonian and suddenly, say at $t=0$, we start evolving with (1). In particular, here we consider the ground states $\ket{\psi_{\bm{\theta},\bm{\phi}}}$ of the following family of local, non-interacting, magnetization Hamiltonians

[TABLE]

where ${\bm{\theta}=(\theta_{1},\ldots,\theta_{L})}$ and ${\bm{\phi}=(\phi_{1},\ldots,\phi_{L})}$ are $L$-component vectors with components ${\theta_{j}\in[0,\pi]}$ and ${\phi_{j}\in[0,2\pi]}$, while $g_{j}>0$ is arbitrary, and

[TABLE]

is the radial unit vector in the three-dimensional space.

The ground states $\ket{\psi_{\bm{\theta},\bm{\phi}}}$ of (10) are separable (i.e. they have zero entanglement): the spin at site $j$ points in the direction $\vec{n}_{\theta_{j},\phi_{j}}$ of its Bloch’s sphere. Namely, the states $\ket{\psi_{\bm{\theta},\bm{\phi}}}$ are explicitly written as

[TABLE]

After $t$ periods of the Floquet evolution the state of the system then reads

[TABLE]

We stress that this protocol is constructive and simple enough to be realizable experimentally, for instance in the context of cold atoms coldatoms ; coldatoms1 ; coldatoms2 .

In this work, the dynamics of the system are characterised by studying the time evolution of the entanglement between a contiguous subset of $N$ spins ${A=\{1,2,\ldots N\}}$ and the rest of the system, see Fig. 1. This is encoded in the density matrix of the system reduced to the subsystem $A$, defined as

[TABLE]

where ${\cal H}_{L-N}$ is the Hilbert space associated with the complement $A^{\rm c}=\{N+1,\ldots L\}$ of $A$. The entanglement content of $\rho_{A}(t)$ is quantified by the Rényi entropies $S^{(\alpha)}_{A}(t)$, also called entanglement entropies. These are a one-parameter family of functionals of $\rho_{A}(t)$ defined as follows

[TABLE]

A particularly important member of this family is the von Neumann entropy

[TABLE]

which is the most used measure of bipartite entanglement for pure states Entrev1 .

In summary: in this paper we study the time evolution generated by (1) of the Rényi entropies (15) for a system initialised in the states $\ket{\psi_{\bm{\theta},\bm{\phi}}}$. As we shall see, our analytical results apply in the thermodynamic limit ${L\to\infty}$. We stress that, in contrast with what we did elsewhere letter , we do not introduce any averaging on the longitudinal magnetic fields. The thermodynamic limit is taken for fixed initial state and time-evolving Hamiltonian.

III Outline of the results

Our main result consists in finding two specific but physically interesting subclasses of the states (12) for which the time evolution of all Rényi entropies in the thermodynamic limit can be found exactly, for any configuration of magnetic fields $\{h_{j}\}$ and subsystem size $N$. These special classes of states are defined as

[TABLE]

where $\bm{1}$ denotes a vector of length $L$ with all entries equal to $1$ note3 . We respectively name them “transverse separating states” and “longitudinal separating states”, while we generically call “separating state” a state belonging to either $\mathcal{T}$ or $\mathcal{L}$. These states are called “transverse” and “longitudinal” because they are respectively eigenstates of the operators $\cos\phi_{j}\sigma^{x}_{j}+\sin\phi_{j}\sigma^{y}_{j}$ and $\sigma^{z}_{j}$ for all $j$s. Therefore, they can be thought of as configurations of a magnet where the spins lie on the $x$-$y$ plane (“transverse plane”) or along the $z$ axis (“longitudinal axis”). The adjective “separating” refers to their key mathematical property and it is thoroughly explained in Sec. V. We stress that the property of being separating is more restrictive that being just separable: all the states $\ket{\psi_{\bm{\theta},\bm{\phi}}}$ are separable but only a subset of them are separating. Specific instances of separating states, which are most relevant from the experimental point of view, are the ground states of the two parts in the Floquet protocol. For example, when $J>0$ and $|h_{j}|<J$ the ground state of $H_{\rm I}$ is $\ket{\psi_{\pi\bm{1},\bm{0}}}\in{\cal L}$, while when $b>0$, the ground state of $H_{\rm K}$ is $\ket{\psi_{\frac{\pi}{2}\bm{1},\pi\bm{1}}}\in{\cal T}$.

To simplify the analysis of the results it is useful to note that the time evolution of the states in ${\cal L}$ can be related to that of those in $\mathcal{T}$. This is easily seen by means of the following identity

[TABLE]

where $\simeq$ denotes equality up to a global phase. This identity is proven by observing that, since the states in $\cal L$ are eigenstates of $\sigma^{z}_{j}$, they are also eigenstates of $H_{\rm I}[\bm{h}]$. Therefore the application of $e^{-iH_{\rm I}[\bm{h}]}$ only changes $\ket{\psi_{\bar{\bm{\theta}},\bm{\phi}}}$ by a global phase. Moreover, an explicit calculation shows that

[TABLE]

Eq. (19) means that the first time step of evolution for states in $\cal L$ keeps them in a separable form, and hence does not change the entanglement, but turns them into states in $\cal T$. An immediate consequence of Eq. (19) is

[TABLE]

Considering the entanglement entropy we then have

[TABLE]

where we explicitly reported the initial state dependence and used $S^{(\alpha)}_{A}(0)=0$. By virtue of this simple argument we can restrict our attention to the states in $\cal T$, the time evolution of the entropy for the states in $\cal L$ is found using Eq. (22).

The time evolution of entanglement entropies from states in ${\cal T}$ (in the thermodynamic limit) can be explicitly determined by means of the “duality method” developed in Secs. IV–VI. The result reads as

[TABLE]

This result is indisputably remarkable: when evolving from separating states all entanglement entropies take the same universal form, independent of the fields $h_{j}$ and of details of the initial states. In particular, at fixed $N$, the entanglement entropies display a linear growth in time up to a non-analytic saturation point where they become constant, see Fig. 2.

The independence of $n$ of the result means that the spectra of the reduced density matrices $\rho_{A}(t)$ are flat. In other words the reduced density matrices have exactly ${2^{\min(2t,N)}}$ eigenvalues equal to ${2^{-\min(2t,N)}}$ while the others vanish.

A form like (23) for the evolution of the entanglement entropies has been found in a number of different physical settings, both in closed and periodically driven systems. Examples range from conformal invariant systems CC ; LS:CFT , to non-integrable closed systems KH:NonIntEnt , from random in time Nahum2 ; r-2017 ; RandomCircuitsEnt , to periodically driven CDC1 random unitary circuits. As opposed to all these cases, however, our result does not hold only at the leading order for large $t$ and $N$. It holds with no corrections for any $t$ and $N$. This gives further evidence for the special status of the self-dual kicked Ising model as minimal solvable model for quantum chaos.

Another interesting feature of our result is that is saturates the “minimal cut” bound CLM:minimalcut . The bound states that, when evolving from a product state,

[TABLE]

where $q$ is the dimension of the local Hilbert space (2 in our case) and $\ell_{\rm min}$ is the minimal number of links intersected by a cut separating the region A from the rest of the system in a local quantum circuit representation of $U_{\rm KI}^{t}$. The fact that the bound is saturated means that entanglement spreads with the maximal possible speed allowed by the range of the Hamiltonian and the dimension of the local Hilbert space. This can be seen in a more physical way by noting that (23) implies that, for $t\leq N/2$, at each time step two more spins of the block $A$ become maximally entangled with the rest of the system. Following Ref. LS:CFT , this can be imagined as an entanglement wave propagating into the block $A$ from the two boundaries. The fact that our exact result (23) saturates the bound (24) also means that it agrees with the “minimal membrane” picture recently put forward in Ref. RandomCircuitsEnt , where a coarse grained version of the minimal cut has been proposed to describe the leading-order-in-time features of the entanglement dynamics in generic systems. Interestingly, however, our system also contains an integrable point, namely $\bm{h}=0$. At this point our result agrees with the the quasiparticle picture of Ref. CC , because in our case all quasiparticles move at the same, maximal, speed ($|v|=1$). We note that, at the integrable point, the result (23) has also been found in Ref. MLS:multipartiteentanglement , for the evolution of the von Neumann entanglement from the separating state $\ket{\psi_{{\frac{\pi}{2}}\bm{1},\bm{0}}}$.

If the initial state is not separating the problem is not amenable to an analytical treatment. In Sec. VII, however, performing a thorough numerical analysis we argue that the entanglement spreading from all the states (12) is still described by (23) at the leading order in time, provided that the system is away from the integrable point $\bm{h}=\bm{0}$. Note that this conjecture is physically very reasonable: since the system is quantum ergodic it is natural to expect the entanglement entropies to eventually become state-independent. On the contrary, in the integrable case our numerical results are consistent with

[TABLE]

where $S^{(\alpha)}_{\bm{\theta},\bm{\phi}}\leq\log 2$ is an initial-state-dependent constant. In the numerical analysis that led to these results, a crucial role has been played by a duality-based numerical approach (see Sec. VII) that allows us to treat the system in the infinite volume limit. Supplemented with some analytical information it can reach $t=17$ Floquet periods of evolution even in cases where the entanglement grows at the maximal speed. Even if based on the “duality-method”, this approach does not rely on the special algebraic structure arising at the self dual points (4) and it is applicable in the whole parameter space of the kicked Ising model.

Another marked difference between the integrable and non-integrable case is observed when the system is confined in a finite volume $L$. Indeed, in this setting the integrable system displays finite-size related recurrences when $t\sim L$, while these recurrences are absent (or at least negligible) in the non-integrable case. These results respectively agree with the predictions of the quasiparticle and the minimal membrane pictures and are also consistent with the numerical analysis of the entanglement spreading in the kicked Ising model carried out in Ref. PL:kickedIsing .

IV Duality Mapping for the Entanglement Entropies

In Ref. Guhr the authors pointed out that the traces of integer powers of the Floquet operator (8) enjoy a useful space-time “duality symmetry”, which can be demonstrated as follows. First note that

[TABLE]

can be represented as a partition function of a classical Ising model on a $t\times L$ lattice, where $U_{\rm KI}[\bm{h}]$ acts as transfer matrix in time, see Fig. 3 for a pictorial representation and Appendix A for the explicit expression.

Second, observe that, due to the short-range nature of the couplings in (2) and (3), the same quantity can also be written in terms of a transfer matrix defined on a lattice of $t$ sites and propagating in space (see Appendix A). Namely we have

[TABLE]

where “tilded” bold symbols denote vectors of $t$ components, in particular

[TABLE]

has all entries equal to $1$, and $\tilde{U}_{\rm KI}[\tilde{\bm{h}}]$ is the transfer matrix in space, also called “dual” transfer matrix. It turns out that $\tilde{U}_{\rm KI}[\tilde{\bm{h}}]$ has the same form as the Floquet operator (8) (with $L$ replaced by $t$ in (2) and (3)) where the longitudinal magnetic field vector is given by $\tilde{\bm{h}}$, while Ising coupling $\tilde{J}$ and the transverse field $\tilde{b}$ are given by the following functions of $J$ and $b$

[TABLE]

Since $\tilde{J}$ and $\tilde{b}$ are generically complex, the transfer matrix $\tilde{U}_{\rm KI}[\tilde{\bm{h}}]$ is generically not unitary. The dual couplings become real only when the model is at one of the self-dual points (4).

In Ref. letter we showed that such duality symmetry can be used to compute non-trivial observables, considering the example of the disorder averaged spectral form factor. In that case, even if the quantity cannot be written in terms of a transfer matrix in time, it can still be written in terms of a transfer matrix in space. This allowed us to perform an analytical calculation. The unitarity of the matrix $\tilde{U}_{\rm KI}[\tilde{\bm{h}}]$, however, proved itself to be a necessary requirement for the analytical approach to be feasible. This clarifies the special status of the self dual points (4): they are the only points of the parameter space where this duality mapping leads to an analytic solution.

Here we develop a similar duality mapping for the calculation of the entanglement entropies, or, more precisely, of the traces of integer powers of the reduced density matrix $\rho_{A}(t)$. We will see that also ${\rm tr}\left[(\rho_{A}(t))^{n}\right]$ can be written as the trace of a power of an appropriate transfer matrix in time. In Secs. V and VI we then show that, at the self dual points and for the special initial states (17) and (18), such a trace can be analytically evaluated.

Considering ${\rm tr}\left[(\rho_{A}(t))^{n}\right]$ and using the definitions (13) and (14) we find

[TABLE]

where $\ket{\bm{a}_{i},\bm{b}_{j}}=\ket{\bm{a}_{i}}\otimes\ket{\bm{b}_{j}}$, $\ket{\bm{a}_{i}}\in\mathcal{B}_{N}$, $\ket{\bm{b}_{i}}\in\mathcal{B}_{L-N}$, for $i,j\in\{1,2,\ldots,n\}$. Here we denoted by $\mathcal{B}_{j}$ the computational basis of

[TABLE]

An explicit expression of ${\cal B}_{j}$ is obtained replacing $L$ by $j$ in the expression (7).

Eq. (31) allows one to interpret the trace of the $n$-th power of the reduced density matrix as the partition function of a classical statistical mechanical model on a multi-sheeted two-dimensional lattice, see Fig. 4 for a pictorial representation in the case $n=3$. To see it more explicitly we consider a single “building block”

[TABLE]

and show that it is equivalent to the partition function of a classical Ising model (with complex weights) on a ${t\times L}$ lattice with periodic boundary conditions in space and fixed boundary conditions in time.

This is seen in two steps. First, we insert $t$ resolutions of the identity operator in the basis (7) into (33), obtaining

[TABLE]

Then we evaluate the matrix elements

[TABLE]

and

[TABLE]

Here, to find the last equation we set $r_{L+1}=r_{1}$ and we used the identity

[TABLE]

to treat the “kick” part of the Floquet operator $U_{\rm KI}[\bm{h}]$. Putting all together we have

[TABLE]

which, as promised, is the partition function of the classical Ising model on a two-dimensional cylinder.

Representing in this way each of the $2n$ building blocks in (31) and summing over $\{\bm{a}_{j},\bm{b}_{j}\}$, one connects together the $2n$ different cylinders obtaining the announced multi-sheeted lattice. Explicitly we have

[TABLE]

where ${\rm sgn}(x)$ is the sign function (we adopted the convention ${{\rm sgn}(0)=1}$), ${{\rm mod}(m,n)=m\,{\rm mod}\,n}$ is the mod-function, and we introduced a new index ${\nu\in\{1,2,\ldots,2n\}}$ such that strings ${\bm{s}_{\nu,\tau}}$ with $\nu\leq n$ belong to terms in (31) with forward time evolution, while those with $\nu>n$ belong to terms in (31) with backward time evolution.

The second line of (39) is obtained by explicitly summing over $\{\bm{a}_{i}\}$ and $\{\bm{b}_{i}\}$ with the help of the identity

[TABLE]

We see that this line forces the configurations of spins in the subchains $A$ and $A^{c}$ on the edges of different cylinders to be the same. These “frozen” configurations are represented by coloured strips in Fig. 4.

To proceed, it is useful to introduce the tensor product space $\mathcal{H}_{t}^{\otimes 2n}$, composed of $2n$ copies of $\mathcal{H}_{t}$, which is the space where the dual Floquet operator $\tilde{U}_{\rm KI}[\tilde{\bm{h}}]$ acts. More formally

[TABLE]

Then, we define the operators $\mathbb{T}_{\theta,\phi}[h]$ and $\mathbb{R}_{\theta,\phi}[h]$ on $\mathcal{H}_{t}^{\otimes 2n}$ through their matrix elements in the computational basis

[TABLE]

and

[TABLE]

where the first subscript labels spin variables in the different copies of $\mathcal{H}_{t}$ composing $\mathcal{H}_{t}^{\otimes 2n}$.

Using the above matrix elements it is immediate to see that the expression (39) can directly be rewritten as a trace (on $\mathcal{H}_{t}^{\otimes 2n}$) of products of $\mathbb{T}_{\theta,\phi}[h]$ and $\mathbb{R}_{\theta,\phi}[h]$, namely

[TABLE]

where we defined on ordered product of non-commuting operators $\{\mathbb{O}_{j}\}$ as

[TABLE]

The rewriting achieved by (44) is pictorially represented in Fig. 5, again in the case $n=3$.

In upcoming analysis it will be useful to think of $\mathcal{H}_{t}^{\otimes 2n}$ as a tensor product of two copies of $\mathcal{H}_{nt}$, grouping together the first and the last $n$ copies of $\mathcal{H}_{t}$, see Fig. 6. Namely we write each element of the basis of $\mathcal{H}_{t}^{\otimes 2n}$ in the following way

[TABLE]

We call these two copies of $\mathcal{H}_{nt}$ the “positive-time” and “negative-time” spaces respectively, as the components of $\mathbb{T}_{\theta,\phi}[h]$ and $\mathbb{R}_{\theta,\phi}[h]$ acting on those spaces come from terms in (31) respectively propagating forward and backward in time.

It is useful to note that $\mathbb{T}_{\theta,\phi}[h]$ and $\mathbb{R}_{\theta,\phi}[h]$ are the same up to a cyclic permutation of the copies of $\mathcal{H}_{t}$ composing the negative-time space (i.e. a cyclic permutation of the second row of Fig. 6), namely

[TABLE]

where we defined

[TABLE]

Here $P_{(\nu,\tau),(\nu-1,\tau)}$ is an elementary transposition

[TABLE]

With $\nu,\nu^{\prime}\in\{1,\ldots,n\}$ and $\tau,\tau^{\prime}\in\{1\ldots,t\}$. The matrix $\sigma^{a}_{\nu,\tau}$ acts as the Pauli matrix $\sigma^{a}$, $a\in\{x,y,z\}$, at the site $\tau=1,\ldots,t$ of the $\nu$-th copy of $\mathcal{H}_{t}$ in $\mathcal{H}_{nt}$, i.e.

[TABLE]

Note that the property $P_{(\nu,\tau),(\mu,\sigma)}=P_{(\nu,\tau),(\mu,\sigma)}^{-1}$ implies $\mathbb{P}^{\dagger}=\mathbb{P}^{-1}$.

Writing (42) in matrix form we have that the transfer matrix is a simple tensor product of single-copy transfer matrices

[TABLE]

The matrix $\mathbb{T}^{(\nu)}_{\theta,\phi}[h]$ acts non-trivially only on the $\nu$-th copy of $\mathcal{H}_{t}$ in both the positive-time and negative-time spaces ($\nu$-th column of Fig. 6), and it is explicitly written as

[TABLE]

where we introduced the Hermitian matrix $\mathbb{B}^{a}_{\nu,\tau}[\theta]$, the projector $\mathbb{G}^{a}_{\nu,\tau}$, and the unitary matrix $\mathbb{U}_{\phi}^{(\nu)}[h]$ defined as follows

[TABLE]

and finally

[TABLE]

Note that since

[TABLE]

the order in the product (51) is irrelevant.

Putting everything together we have

[TABLE]

This equation accomplishes the duality mapping of the entanglement entropies: we wrote the entanglement entropies in terms of the trace of products of an appropriate transfer matrix in space.

Before continuing with the evaluation of (60) two comments are in order. First we note that the mapping described can be performed also when $J$ and $b$ in (2) and (3) do not fulfil the self-duality condition (4). For generic $J$ and $b$ we obtain that the entropy is still given by (60) but the transfer matrix $\mathbb{T}_{\theta,\phi}[h]$ is modified in two ways. (i) the matrix $\mathbb{U}_{\phi}^{(\nu)}[h]$ is not unitary anymore. The Ising coupling in (56) replaced by $\tilde{J}$ (cf. (29)) and the transverse fields in (55) (the coefficients of $iM^{x}_{\nu}$ in the positive and negative time copy) are respectively replaced by $\tilde{b}$ and $\tilde{b}^{*}$ (cf. (30)). (ii) the projector $\mathbb{G}^{z}_{\nu,t}$ in (60) is replaced by

[TABLE]

As we will see in the next section these changes are enough to hinder the analytical evaluation of (60), however, the duality approach can still be useful for perturbative calculations or numerical approaches (see Sec. VII).

We also observe that when the initial state is in the class $\cal L$ (cf. (18)), namely when

[TABLE]

the expression (60) can be further simplified by effectively reducing the dimension of the space where the trace acts. This is explicitly shown in Appendix B. The final result is again of the form (60) with the replacement

[TABLE]

Here the barred operators have exactly the same form as the non-barred ones (respectively (51) and (48)) but act on $\mathcal{H}_{t-1}^{\otimes 2n}$ instead of $\mathcal{H}_{t}^{\otimes 2n}$. Note that this is nothing but a restatement of Property (22).

V Separating states

Our goal is to use Equation (60) to determine $S^{(n)}_{A}(t)$ in the thermodynamic limit. To do that, however, we need some information on the Jordan normal form of the matrix $\mathbb{T}_{\theta,\phi}[h]$. Indeed, since the matrix is not normal, it is not guaranteed to be (and it is generically not) diagonalisable.

As proven in Appendix C, the form (51)–(52) of the transfer matrix has some simple but useful consequences on its Jordan normal form. Specifically we have

Property 1**.**

The following facts hold

• (i)

$|\lambda_{j}|\leq\lambda_{\rm max}\equiv(1+|\cos\theta|)^{n}$, $\forall\lambda_{j}\in{\rm Spec}\left[\mathbb{T}_{\theta,\phi}[h]\right].$

• (ii)

If an eigenvalue $\lambda$ of $\mathbb{T}_{\theta,\phi}[h]$ fulfils $|\lambda|=\lambda_{\rm max}$ then

• a.

$\lambda$* has trivial Jordan blocks (its geometric and algebraic multiplicities coincide).*

• b.

the associated left eigenvector $\bra{A}$ satisfies

[TABLE]

where ${\rm Spec}\left[A\right]$ denotes the spectrum of the matrix $A$.

Property 1 introduces the crucial simplification of this work. If the maximal eigenvalues of $\mathbb{T}_{\theta,\phi}[h]$ saturate the bounds at point $(i)$ the problem of finding the maximal eigenvalues of the transfer matrix is separated into three much easier ones, consisting of finding eigenvalues and eigenvectors of simple hermitian and unitary matrices.

The bound at point $(i)$, however, cannot be always saturated. To see this, let us consider some constrains on the structure of the matrix $\mathbb{T}_{\theta,\phi}[h]$ coming from the identity (44). These are most easily found by considering the translational invariant case

[TABLE]

Setting $N=0$ in (44) we have

[TABLE]

where in the second step we used that the state (13) is pure. This relation implies that the eigenvalues of $\mathbb{T}_{\theta,\phi}[h]$ are all 0 but one, which is equal to 1. Moreover, the Jordan block corresponding to the eigenvalue 1 is one-dimensional, while the eigenvalue 0 might have (and does have!) a highly nontrivial Jordan structure. More formally:

[TABLE]

From the conditions (69) it follows that the bound at point $(i)$ of Property 1 can be saturated only when $\lambda_{\rm max}=1$. Note that the cases for which $\lambda_{\rm max}=1$ include ${\theta={\pi}/2}$, but also $\theta=0,\pi$. Indeed in the latter case the matrix $\mathbb{T}_{\theta,\phi}[h]$ can be replaced by $\bar{\mathbb{T}}_{{\pi}/2,\pi/2-\theta}[h]$ (see (63) and Appendix B). In other words, the requirement $\lambda_{\rm max}=1$ selects the two classes of states ${\cal T}$ and ${\cal L}$ introduced in Sec. III. This clarifies the meaning of their name. We called them “separating” states because if the initial state is one of them the problem of finding the maximal eigenvalues of the transfer matrix (and the corresponding eigenvectors) can be separated. In the upcoming section we explicitly solve the separated problem (64)–(66) for $\lambda_{\max}=1$, and, incidentally, we also verify that it has no solution for $\lambda_{\max}\neq 1$.

Finally, we note that away from the self dual points (4) the conditions (69) still hold and a property similar to Property 1 is still valid. In that case, however, the bound can never be saturated and no separation can be performed. This makes the problem analytically intractable, at least in an exact fashion.

VI Entanglement spreading from separating states

Here we explicitly solve the entanglement evolution from separating states. In particular in Sec. VI.1 we solve the separated problem (64)–(66) for $\lambda_{\max}=1$ and in Sec. VI.2 we evaluate (60). To be concrete we focus on initial states in the class $\mathcal{T}$, the result for states in the class $\mathcal{L}$ is obtained using (22).

VI.1 Maximal eigenvalues of the transfer matrix

Our strategy is to determine the maximal eigenvalues of $\mathbb{T}_{\pi/2,\phi}[h]$ and the associated eigenvectors, by searching for all the vectors fulfilling (64)–(66) with $\lambda_{\max}=1$. To simplify our analysis we make two observations. First, we note that

[TABLE]

so that (64) becomes trivial. Second, we note that all $\mathbb{G}^{z}_{\nu,t}$ and $\mathbb{U}_{\phi}^{(\nu)}$ commute for different $\nu$s so we can look for simultaneous eigenvectors. The problem is then reduced to finding all vectors $\bra{A}$ fulfilling

[TABLE]

To solve these equations it is convenient to introduce the following one-to-one vector-to-operator mapping (cf. Ref. letter ) $\bra{A}\leftrightarrow A$:

[TABLE]

where $\{\bra{k}\}$ is a basis of $\mathcal{H}_{nt}$ and $(\cdot)^{*}$ denotes complex conjugation in the computational basis $\mathcal{B}_{nt}$, such that

[TABLE]

for any operator $O$. Using the mapping (73), Eqs. (71)–(72) are directly rewritten in operatorial form as follows

[TABLE]

for some $\alpha_{\nu}\in\mathbb{R}$ and all $\nu\in\{1,\ldots,n\}$. In this formulation our goal is to find all independent linear operators $A$ over ${\cal H}_{nt}$ solving the commutation relations (75)–(76). As shown in Appendix D, these commutation relations are equivalent to

[TABLE]

for all $a\in\{x,y,z\}$, $\tau\in\{1,\ldots,t\}$, $\nu\in\{1,\ldots,n\}$. Namely, they are equivalent to requiring that $A$ commutes with the entire algebra of observables in $\mathcal{H}_{nt}$. Since the latter is irreducible, Shur’s Lemma implies that the unique (up to multiplicative factors) solution to (77) is given by

[TABLE]

We then find that the eigenvalue of $\mathbb{T}_{\pi/2,\phi}[h]$ with maximal magnitude is $1$ and corresponds to the unique left eigenvector

[TABLE]

where we used the computational basis, omitted complex conjugation as the basis is real, and we included the normalisation factor ${\sqrt{{\rm tr}[\mathbbm{1}]}=2^{nt/2}}$. Note that the unique right eigenvector of $\mathbb{T}_{\pi/2,\phi}[h]$ associated to ${\lambda=1}$ is given by $\ket{\mathbbm{1}}=(\bra{\mathbbm{1}})^{\dagger}$, as it can be directly verified.

We also observe that since we just proved that (65) and (66) have $\bra{A}=\bra{\mathbbm{1}}$ as only solution, and, moreover

[TABLE]

the separated problem (64)–(66) has no solution for ${\theta\neq\pi/2}$.

VI.2 Entanglement dynamics

Our next step is to use the eigenvectors determined above to compute the entanglement dynamics. First we note that the eigenvector $\ket{\mathbbm{1}}$ is independent of $\phi$ and $h$. Moreover, $\ket{\mathbbm{1}}$ is orthogonal to all left generalised eigenvectors corresponding to the eigenvalues [math] of $\mathbb{T}_{\pi/2,\phi}[h]$ for all $\phi$ and $h$. These two facts imply

[TABLE]

where we introduced

[TABLE]

The relation (81) can be used to find the slope of the linear growth of the entanglement entropy. Indeed, taking $N$ to infinity we have

[TABLE]

The simple structure of $\mathbb{T}_{{\pi}/{2},\phi}[h]$, however, allows us to progress further and evaluate (81) exactly for each $N$.

This can be done by making use of the following remarkable identity

[TABLE]

where $\lfloor\cdot\rfloor$ denotes the floor function. Here we adopted the convention

[TABLE]

and, to lighten the notation, from now on we assume that a product $\prod\cdots$ only picks a single factor on its right unless several terms are grouped within a square bracket $[\cdots]$.

The identity (84) is proven in Appendix F using the explicit form of $\mathbb{T}_{{\pi}/{2},\phi}[h]$ and the following useful properties of the state (82)

[TABLE]

where $O_{\nu}$ acts non trivially, as the unitary operator $O$, only on the $\nu$-th copy of $\mathcal{H}_{t}$ in $\mathcal{H}_{nt}$, i.e.

[TABLE]

These properties are proven in Appendix E.

A striking consequence of (84) is that the entanglement entropies evolving from separating states are completely independent of the configuration of longitudinal magnetic fields $\{h_{j}\}$ and of the initial-state angles $\{\phi_{i}\}$. For instance, this means that the same result is obtained in the integrable and in the non-integrable case, with or without disorder.

The evaluation of the r.h.s. of (84) is now straightforward. First, we note that in the computational basis (46) of $\mathcal{H}_{t}^{\otimes 2n}$ we have

[TABLE]

where the matrix elements of $\mathbb{G}^{x}_{\nu,\tau}$ and $\mathbb{G}^{z}_{\nu,\tau}$ are computed by repeated use of

[TABLE]

Then, we plug (89) into the r.h.s. of Eq. (84). For $t>\lfloor\!N/2\!\rfloor$ we find

[TABLE]

Proceeding analogously for $t\leq\lfloor N/2\rfloor$ we have

[TABLE]

Therefore, we finally obtain that for initial states in the class $\mathcal{T}$, all entanglement entropies $S^{(n)}_{A}(t)$ with ${n=2,3,\ldots}$ are exactly given by Eq. (23). This, however, implies that

[TABLE]

where $2^{-\min(2t,N)}$ has multiplicity $2^{\min(2t,N)}$, while [math] has multiplicity $2^{N-\min(2t,N)}$. As a consequence, the result (23) holds for all $S^{(\alpha)}_{A}(t)$ with any real positive $\alpha$.

VII Entanglement spreading from generic states

The exact results derived in the previous sections have three remarkable features. (i) The entropies do not depend at all on the longitudinal magnetic fields. In particular, they are not affected by whether or not the system is integrable; (ii) The entropies grow at the maximal speed allowed by the range of the Hamiltonian and the dimension of the local Hilbert space (they saturate the minimal cut bound (24)); (iii) At each fixed time $t$ all entanglement entropies coincide, signalling a flat entanglement spectrum, i.e., that all non-zero eigenvalues of the density matrix reduced to the block $A$ are equal.

It is interesting to wonder whether these are general features of the entanglement spreading in the self-dual kicked Ising chain or, instead, they are special properties of separating initial states. In other words, it is interesting to ask whether the entanglement dynamics from separating states is an exceptional case or, even though special, it can be used to model the generic behaviour. To this aim, in this section we consider the entanglement spreading from generic product states (12) which are not separating. In this case, as pointed out above, we are unable to address the problem in a fully analytical fashion and we resort to a numerical analysis.

From the physical point of view it is easy to see that the most convenient time regimes to examine possible modifications of the features (i)–(iii) are very different. Indeed, for ${\bm{h}\neq\bm{0}}$ the system is ergodic, and any finite subsystem is expected to relax to the infinite temperature state irrespectively of the initial conditions. This means all entropies are expected to saturate to the universal value $N\log 2$. On the other hand, for ${\bm{h}=\bm{0}}$ the system is integrable and finite subsystems relax to generalised Gibbs ensembles EF:GGE ; VR:GGE . We then expect the stationary values of the entropies to retain some memory of the initial configuration. To highlight the difference between integrable and non-integrable systems it is then convenient to focus on the “saturation regime” ${t\sim N}$, where the entropies become stationary. The generic relaxation to the infinite temperature state, however, also means that to see some dependence of the entanglement spectrum on $\bm{h}$, or on the initial state, one has to stay away from the saturation regime and focus on the “growth regime”, ${t\ll N}$. The latter is obviously also the regime of interest to study variations in the speed of entanglement growth.

The saturation regime can be easily accessed by a “direct” numerical approach. Namely, we consider a finite volume $L$ and determine the time evolving state by means of the efficient time-propagation algorithm described in the supplemental material of Ref. letter . The entanglement entropies are found by computing and diagonalising the reduced density matrices $\rho_{A}(t)$ (cf. Eq. (14)) and using Eq. (15). Note that a similar numerical analysis, in the case of the von Neumann entropy, has been performed in Ref. PL:kickedIsing .

Some representative examples of our results are reported in Fig. 7. First of all we see that the qualitative behaviour of the entanglement entropies is the same as for separating states, both in the homogeneous (translationally-invariant) and in the inhomogeneous case. The entropy grows in an approximately linear fashion until it saturates to a value proportional to the subsystem size. There is, however, a clear qualitative difference emerging between the integrable case and the generic one: in the generic case the entropies always saturate to $N\log 2$ (minus the expected correction due to a finite $N/L$ Page ; Lev ; KH:NonIntEnt ), while this does not happen at the integrable point. In particular, in the inset of Fig. 8 we report the evolution of $S_{A}^{(2)}(t)$ for several homogeneous non-separating initial states evolving under the integrable kicked Ising Hamiltonian. We see that, in contrast to the generic case, the saturation values depend on the initial state.

Interestingly, we see that the evolution of the entropies shows very different finite-size effects in the integrable and non-integrable cases. In the former the entropies start to decrease at times larger than $(L-N)/2$, while in the latter they remain constant once they reached the saturation values. These behaviours respectively agree with the predictions of the quasiparticle and the minimal-membrane picture. Indeed, for $L>2N$ the surface of the membrane is not affected by the system being finite. On the contrary, the quasiparticle picture predicts oscillations of the entropies due to quasiparticles traversing the entire system and going back to their initial positions. In particular, if the initial state is homogeneous, using that in our case the quasiparticles have all unit speed (and taking, for convenience, $L$ even), we find that the quasiparticle-picture prediction is $L/2$-periodic and, for $t\in\{0,1,\ldots,L/2\}$, it reads as

[TABLE]

where $S^{(\alpha)}_{\theta,\phi}\leq\log 2$ is a ($N$- and $L$- independent) constant. This prediction holds in the asymptotic limit $t,N\to\infty$ with fixed $t/N$, but, as shown in the main panel of Fig. 8, it is in fair agreement with our numerical results already for $N=11$. Note that, even if the system is free, $S^{(\alpha)}_{\theta,\phi}$ cannot be generically computed analytically. Indeed, for generic values of $\bm{\theta}$ and $\bm{\phi}$ the states are not Gaussian in terms of the time evolving fermions and this makes the problem analytically untreatable. Moreover, since the dispersion is linear, the usual arguments about Gaussification do not apply Gaussification1 ; Gaussification2 . Interestingly, not even separating states are always Gaussian: transverse separating states are Gaussian only for $\phi_{i}=0,\pi$ notegaussian .

A natural question is what happens to the finite-size oscillations when the integrability is weakly broken. This is investigated in Fig. 9, which compares the behaviour of the von Neumann entropy (cf. (16)) for increasing values of the (homogeneous) longitudinal magnetic field. We see that the finite-size oscillations become damped and disappear at large enough times. This can be interpreted as a sign of the decay of the quasiparticles. Consistently, the decay speed increases with the magnitude of the longitudinal magnetic field. Moreover, for any fixed $\bm{h}\neq 0$ the peaks are observed to decay when the volume of the system increases.

Finally, we note that Fig. 7 also contains some information on the speed of entanglement growth. Indeed, we see that the time evolution of the entropies depends (although weakly for ${\bm{h}\neq\bm{0}}$) on the configuration of magnetic fields, indicating that the feature (ii) is lost at short times. It is, however, very hard to make any definitive statement based on Fig. 7. The direct numerical approach allows us to access the linear growth regime only for very short times and it is impossible to exclude that (ii) remains as an asymptotic feature of the entanglement dynamics. A similar argument holds regarding the entanglement spectrum.

A useful way to circumvent this problem is offered by the expression for the Rényi entropies resulting from the duality mapping, namely Eq. (60) (a similar duality-based approach has also been proposed in the context of tensor networks BHVC ). This becomes particularly convenient in the translational invariant case

[TABLE]

Indeed, in this case one can use the general constraints (69) to take analytically the limit of infinite $L$ and $N$ and focus on the growth regime. Specifically, for $n=2,3,\ldots$, we find

[TABLE]

where $\bra{M_{\rm L}}$ and $\ket{M_{\rm R}}$ are respectively the left and right eigenstates of $\mathbb{T}_{\theta,\phi}[h]$ corresponding to the eigenvalue $1$. The constraints (69) imply that these vectors exist and are unique. Note that the simplification (96) cannot be generically performed in the inhomogeneous case, since transfer matrices with different $h,\theta$ and $\phi$ have different left and right eigenvectors.

The numerical evaluation of (96) is achieved in two steps. First one has to determine the left and right eigenvectors and then to evaluate the matrix element. Finding eigenvectors is particularly convenient due to the tensor product structure of the transfer matrix (cf. (51)). Indeed, we can search eigenvectors of the form

[TABLE]

where ${\ket{A_{\rm R}},\ket{A_{\rm L}}\in{\mathcal{H}_{t}\otimes\mathcal{H}_{t}}}$. In the notation of Fig. 6, this means that we can effectively work in the Hilbert space of the $\nu$-th copy in both the positive-time and negative-time spaces (corresponding to the $\nu$-th column of Fig. 6). The eigenvectors are efficiently determined by means of a simple “power method”: one starts from a random vector and finds $\ket{A_{\rm R}}$ by repeated application of $\mathbb{T}_{\theta,\phi}[h]$. The left eigenvector $\bra{A_{\rm L}}$ is then determined by using $\ket{A_{\rm L}}=e^{-i\frac{\pi}{8}M^{x}_{1}}\otimes e^{i\frac{\pi}{8}M^{x}_{1}}\ket{A_{\rm R}}$. Proceeding in this way the eigenvector can be determined in ${\cal O}(tm2^{2t})$ operations, where $m$ is the number of iterations of the power method notepowermethod .

The form (97) is also convenient for evaluating the matrix element in (96). Indeed, after a straightforward calculation we find

[TABLE]

where $A_{\rm L,\rm R}$ are the $2^{t}\times 2^{t}$ matrices corresponding to the vectors $\bra{A_{\rm R,L}}$ through the vector-to-operator mapping (73) (performed for ${n=1}$). Note that for transverse separating initial states we have

[TABLE]

but for generic initial states these matrices become non-trivial. The r.h.s. of Eq. (98) can be numerically evaluated for integer $n\geq 2$, whereas Rényi entropies with more general index $\alpha>0$ can be found by analytically continuing (98) in $n$ and diagonalising $A_{\rm R}^{\dagger}A^{\phantom{{\dagger}}}_{\rm L}$ numerically to compute the powers. Evaluating  (98) has complexity $\propto 2^{3t}$ and is the bottleneck of the numerical procedure, meaning that we were able to reach up to ${t_{\rm max}=17}$.

We stress that, since the constraints (69) hold also away from the self dual points, this procedure can be used to study the entanglement spreading in the entire parameter space of the kicked Ising model. Note that close enough to the self-dual points (and to separating initial states) an analytical, perturbative, analysis is also possible. These aspects, however, go beyond the scope of the present manuscript and will be investigated the course of future research. Here we focus on the self-dual points (5) and use this duality-based numerical approach to effectively investigate the fate of the features (ii) and (iii) when the system is initialised in a generic state (12).

Representative examples of our numerical results are reported in Figs. 10 and 11. We see that, consistently with the results in Fig. 7, at short times both (ii) and (iii) are violated. The entropies grow in an approximately linear fashion but the the slope appears to depend on the initial state and on the longitudinal magnetic field $h$ (see Fig. 10). Moreover, different Rényi entropies have different slopes (see Fig. 11). Crucially, however, a more refined analysis suggests that these deviations vanish for large times. To show this we proceed as follows. First we introduce the the “instantaneous” slopes

[TABLE]

Plotting the instantaneous slopes as functions of $1/t$ we see that they become approximately linear at large enough times. We then perform a linear fit and extrapolate the result to $t=\infty$. In the integrable case this procedure gives results consistent with the quasiparticle picture prediction (cf. Eq. (94)), namely

[TABLE]

Instead, in the non-integrable case the results are consistent with

[TABLE]

This is observed for any initial state (12), for any non-vanishing longitudinal magnetic field, and for any Rényi index $\alpha$, as demonstrated in the insets of Figs. 10 and 11. The only seemingly exceptional cases are observed when the system is very close to the integrable point (see, e.g., the red curve in the inset of Fig. 10). This is, however, straightforwardly explained as a prethermalization effect MK:pret ; BEGR:pret ; Pretrev . For small enough longitudinal fields there is an initial transient, before the quasiparticles decay, in which the observables follow the integrable predictions. After this transient, however, the entropies are expected to follow the non-integrable curves, reaching the asymptotic value (102) for the slope. For instance, this is consistent with the behaviour of the red curve in the inset of Fig. 10. Interestingly, since $S^{(\alpha)}_{\theta,\phi}=\log 2$ for separating states, this effect is not observed in our exact result (23). From the physical point of view (102) is very natural. Since the system is ergodic the initial state dependence is washed away at large enough times, and the entropies behave as if they would be evolving from separating states. As expected, this does not happen in the integrable case.

In conclusion, the numerical results presented in this section support in the following general picture. The evolution of the entanglement entropies from a generic initial state (12) in the thermodynamic limit differs from that from separating states (cf. (23)), and depends explicitly on the initial state, the longitudinal magnetic fields, and the Rényi numbers. In the scaling limit $t,N\to\infty$, however, all these dependences are washed away: the entropies collapse to the prediction (23) if the system is non-integrable and to the thermodynamic limit of (94) if the system is integrable. In other words, our exact result (23) serves as an asymptotic description of the entanglement spreading in the non-integrable kicked Ising model. This picture is also supported by the numerical results of Ref. PL:kickedIsing , which found that the evolution of the von Neumann entropy averaged over all separable initial states, is consistent with (23). Finally, we remark that this section also demonstrated that extracting universal information on the behaviour of the entanglement entropies from the numerics is extremely hard, especially in the ergodic case where the entanglement growth is exceptionally fast. This highlights even more the practical importance of our exact result (23).

VIII Conclusions

We have developed a constructive and mathematically rigorous approach for computing the dynamics of bipartite entanglement in a class of “maximally scrambling”, locally interacting, chaotic spin chains. Specifically, we considered the so called “self-dual” kicked Ising spin chains, where the integrability is broken by switching on an external longitudinal magnetic field. We prepared the system in class of ground states of simple local Hamiltonians and determined exactly the dynamics of all Rényi entropies of finite blocks of spins of arbitrary size. The results presented are non-perturbative, no kind of averaging is involved, and, most importantly, they hold in the presence of longitudinal magnetic fields with arbitrary spatial dependence. It is remarkable that such an explicit exact calculation can be performed for a specific non-integrable many-body system.

Our result shows that in the thermodynamic limit the Rényi entropies of finite blocks of spins are independent of the longitudinal magnetic field at all times. Moreover, they obey universal scaling laws that can be predicted both by means of the quasiparticle picture of Ref. CC , put forward for integrable models, and of the minimal membrane picture of Ref. RandomCircuitsEnt , propounded for generic systems.

Using our novel rigorous approach, we also developed a numerical procedure for studying the entanglement spreading from generic product initial states. A thorough numerical analysis suggests that, away from the integrable point, our exact result continues to describe the entanglement spreading at the leading order in time. On the contrary, in the integrable case the entanglement production is generically renormalised by an initial-state-dependent multiplicative coefficient. Further qualitative differences between the integrable and the non-integrable case emerge for finite systems. In particular, we showed numerically that there are recurrences in the integrable case, which are absent in the non-integrable one. We stress that these differences are correctly accounted for by the quasiparticle and minimal membrane pictures, which disagree for finite sizes.

Our analytical method can be used to highlight qualitative differences in the entanglement spreading of integrable and non-integrable systems directly in the thermodynamic limit. To do that one could follow Refs. ABGH:CFT ; LM:CFT and consider the bipartite entanglement of disjoint blocks. Our preliminary results suggest that the scaling forms produced in the two cases are indeed different and respectively agree with the predictions of quasiparticle and membrane pictures. Another possible direction is to perturb the kicked Ising spin chains away from the “self-dual” points, where the predictions of the two pictures disagree also for the entanglement of a single block. This could be tested within our approach by using perturbation theory.

More generally, we expect that our method would allow for explicit calculations similar to the ones presented also for other measures of correlations and dynamical complexity, such as operator space entanglement entropy and out-of-time order correlators.

Finally, we believe that the remarkable algebraic structure unveiled in this work paves the way for the determination of a new class of exactly solvable, maximally chaotic models. Elements of this class can serve as minimal models for characterising the non-equilibrium dynamics in generic systems.

Acknowledgements

We thank Adam Nahum for very valuable comments on the manuscript and Pasquale Calabrese for stimulating discussions. B.B. and T.P. acknowledge the hospitality of the Erwin Schrödinger Institute (ESI), Vienna, where this project has been conceived. The work has been supported by Advanced Grant of European Research Council (ERC), No. 694544 – OMNES and program P1-0402 of Slovenian Research Agency.

Appendix A Duality of traces

Here we explicitly demonstrate the duality relation (27). Writing ${\rm tr}\left[\left(U_{\rm KI}[\bm{h}]\right)^{t}\right]$ in the computational basis $\mathcal{B}_{L}$ (cf. (7)) we have

[TABLE]

Here $s_{\tau,L+j}\equiv s_{\tau,j}$ and in the second step we used the identity

[TABLE]

where

[TABLE]

This expression can be thought of as the partition function of a two-dimensional Ising model with complex couplings on a $L\times t$ periodic lattice. In other words, the r.h.s. of (103) is proportional to the partition function of a classical statistical mechanical model with configuration energy given by

[TABLE]

Reorganising the sum on the r.h.s. of (103) we also have

[TABLE]

where we defined $s_{t+\tau,j}\equiv s_{\tau,j}$. Using again the identity (104) we finally find

[TABLE]

where “tilded” bold symbols denote vectors of $t$ components and we introduced the dual transfer matrix

[TABLE]

with

[TABLE]

Appendix B Simplified transfer matrix for longitudinal separating states

When the initial state is in the class $\cal L$ (cf. (18)), namely when

[TABLE]

the form (60) can be simplified by effectively reducing the dimension of the space where the trace acts. To see this we note that in this case $\mathbb{B}^{z}_{\nu,1}[\theta]$ becomes proportional to a projector:

[TABLE]

so that we have

[TABLE]

where we introduced

[TABLE]

Here the matrix $\bar{\mathbb{U}}_{\phi}^{(\nu)}[h]$ is defined as

[TABLE]

and the barred operators read as

[TABLE]

So the have the same form as (56) and (57) but at fixed $\nu$ they act non-trivially only in the space $\mathcal{H}_{t-1}$ composed of the last $t-1$ sites of $\mathcal{H}_{t}$. In other words $\bar{\mathbb{T}}_{\theta,s}[h]$ has the same form as ${\mathbb{T}}_{\theta,s}[h]$ but acts on $\mathcal{H}_{t-1}^{\otimes 2n}$ instead of $\mathcal{H}_{t}^{\otimes 2n}$. We stress that the trace operations in expression (113) and below are taken in the subspaces where the operators act nontrivially, for example, for the barred operators in ${\cal H}^{\otimes 2n}_{t-1}\cong\mathcal{\cal H}_{2n(t-1)}$. Noting

[TABLE]

we finally find

[TABLE]

Therefore we see that in this case the entropies are given by an expression of the form (60), with $\theta_{j}=\tfrac{\pi}{2}$ and $\phi_{j}=\tfrac{\pi}{2}s_{j}$, but with matrices acting on $\mathcal{H}_{t-1}^{\otimes 2n}$ instead of $\mathcal{H}_{t}^{\otimes 2n}$. Note that for $\theta_{j}=\bar{\theta}_{j}$ the states (12) do not depend on $\phi_{i}$ and this independence is correctly reflected in Eq. (119).

Appendix C Proof of Property 1

In this appendix we provide the proof of Property 1.

Proof.

For each state $\bra{A}$ we have

[TABLE]

where we used that $\mathbb{G}^{z}_{\nu,t}$ is a projector, so its expectation value on a normalised state is smaller or equal to one. Expanding $\mathbb{B}^{z}_{\nu,1}[\theta]^{2}$ we then have

[TABLE]

Since $P_{\nu,1}^{z,\pm}\otimes P_{\nu,1}^{z,\pm}$ are orthogonal projectors we have

[TABLE]

In particular, choosing $\bra{A}$ to be the left eigenstate of $\mathbb{T}_{\theta,\phi}[h]$ corresponding to the eigenvalue $\lambda$ we have

[TABLE]

which proves the first part of the claim.

To prove the point $(ii)$ $a.$ we proceed by reductio ad absurdum. Suppose that the Jordan block of $\lambda$ is non-trivial: let $\bra{A}$ be the eigenvector associated to $\lambda$ and let $\bra{B}$ be the first generalised eigenvector. As it is always possible, we choose $\bra{B}$ to be normalised and orthogonal to $\bra{A}$ (which is also normalised). We then have

[TABLE]

This implies

[TABLE]

which is impossible because it contradicts (122). Point $(ii)$ $b.$ follows by noting that in order to have the equality sign in (122) we must have

[TABLE]

Using now that $\bra{A}$ is a left eigenvector of $\mathbb{T}_{\theta,\phi}[h]$ we have (66). This concludes the proof.

∎

Appendix D Simplified commutation relations

In this appendix we prove the following property.

Property 2**.**

The commutation relations (75)–(76) imply

[TABLE]

Proof.

First of all we note that multiplying (76) on the left and on the right by $e^{-i\frac{\pi}{4}M_{\nu}^{x}}\,e^{ihM_{\nu}^{z}}{U}^{\dagger}_{\nu,\phi}$ we have

[TABLE]

Using the conditions (75), (76), and (129) we see that $A$ commutes with

[TABLE]

Indeed we have

[TABLE]

where in the the first step we used (129), in the second (76), and in the third (75). Using (75) we then have that $A$ also commutes with

[TABLE]

We then have

[TABLE]

Using (133), (76), and (129) we can then conclude the proof by induction.

We will prove that if

[TABLE]

then

[TABLE]

and then proceeding by induction in $\bar{\tau}=t-1,\ldots,1$. The basis of the induction is given by (133), so we just need to prove the inductive step. Assuming (134) and proceeding as in (131) we can show that $A$ commutes also with

[TABLE]

The inductive hypothesis (134) then implies

[TABLE]

Reasoning now as in (131) we then have that $A$ also commutes with

[TABLE]

So we have

[TABLE]

This concludes the proof. ∎

Appendix E Proof of Eqs. (86) and (87)

Let us start by proving (86). First we note

[TABLE]

This is explicitly proven as follows

[TABLE]

Here we used (74) and the fact that $O_{\nu}$ is unitary. Second, we observe that from the definition of $O_{\nu}$ it directly follows

[TABLE]

Combining (143) and (141) we then have

[TABLE]

So we proved (86) for any $O_{\nu}$ acting non trivially, as the unitary operator $O$, only on the $\nu$-th copy of $\mathcal{H}_{t}$ in $\mathcal{H}_{nt}$. The relation (87) follows immediately by taking the adjoint of

[TABLE]

Appendix F Proof of Property (84)

In this appendix we prove Eq (84).

Proof.

Defining

[TABLE]

we can rewrite the l.h.s. of (84) as follows

[TABLE]

To simplify this expression we proceed as follows. First we commute every possible $\mathbb{J}_{\nu,\tau},\mathbb{Z}_{\nu,\tau},\mathbb{X}_{\nu,\tau}$ to the left by using the following commutation relations

[TABLE]

Then we use

[TABLE]

which follow from (86). Finally, using also

[TABLE]

following from (87), on the rightmost term in the product over $j$ we find

[TABLE]

where we defined

[TABLE]

Using

[TABLE]

we can rewrite (163) as follows

[TABLE]

We now make use the following Lemma, proven in Appendix F.1, to simplify the products of $\mathbb{A}_{\nu,j}$s

Lemma 1**.**

[TABLE]

Using now (156), (159), (160), and (161) we have

[TABLE]

which concludes the proof. ∎

F.1 Proof of Lemma 1.

Here we prove Lemma 1.

Proof.

We proceed by induction in the number of terms in the products of $\mathbb{A}_{\nu,j}$s. First we establish the basis. We begin by computing

[TABLE]

where we used

[TABLE]

We see that (170) agrees with (167) for $n=1$. We then compute

[TABLE]

In the last step we used

[TABLE]

Since (172) agrees with (168) for $n=1$ we successfully established the basis for the inductive procedure. To conclude we need to prove that

• $(i)$

if (167) holds for $n$ then (168) holds for $n$

• $(ii)$

if (168) holds for $n$ then (167) holds for $n+1$.

Let us prove $(i)$. Assuming (167) we have

[TABLE]

where we used that $\mathbb{G}^{x}_{\nu,t}$ commutes with all the terms on the first line to bring it close to $\mathbb{J}_{\nu,t-1}$. Then we employed (171). Moving the projector $\mathbb{G}^{z}_{\nu,t-1}$ on the second line to the left and using multiple times (165) and (171) we have

[TABLE]

which is exactly (168). Let us now prove $(ii)$. Assuming (168) we have

[TABLE]

where we again used that $\mathbb{G}^{x}_{\nu,t}$ commutes with all the terms on the first line to bring it close to $\mathbb{J}_{\nu,t-1}$ and employed (171). Moving now the projector $\mathbb{G}^{z}_{\nu,t-1}$ on the second line to the left and using (165) and (171) we have

[TABLE]

which is exactly (167) for $n+1$. This concludes the proof.

∎