Growth of R\'enyi Entropies in Interacting Integrable Models and the Breakdown of the Quasiparticle Picture

R\'enyi entropies are conceptually valuable and experimentally relevant generalisations of the celebrated von Neumann entanglement entropy. After a quantum quench in a clean quantum many-body system they generically display a universal linear growth in time followed by saturation. While a finite subsystem is essentially at local equilibrium when the entanglement saturates, it is genuinely out-of-equilibrium in the growth phase. In particular, the slope of the growth carries vital information on the nature of the system's dynamics, and its characterisation is a key objective of current research. Here we show that the slope of R\'enyi entropies can be determined by means of a spacetime duality transformation. In essence, we argue that the slope coincides with the stationary density of entropy of the model obtained by exchanging the roles of space and time. Therefore, very surprisingly, the slope of the entanglement is expressed as an equilibrium quantity. We use this observation to find an explicit exact formula for the slope of R\'enyi entropies in all integrable models treatable by thermodynamic Bethe ansatz and evolving from integrable initial states. Interestingly, this formula can be understood in terms of a quasiparticle picture only in the von Neumann limit.

Rényi entropies are conceptually valuable and experimentally relevant generalisations of the celebrated von Neumann entanglement entropy. After a quantum quench in a clean quantum many-body system they generically display a universal linear growth in time followed by saturation. While a finite subsystem is essentially at local equilibrium when the entanglement saturates, it is genuinely out-of-equilibrium in the growth phase. In particular, the slope of the growth carries vital information on the nature of the system's dynamics, and its characterisation is a key objective of current research. Here we show that the slope of Rényi entropies can be determined by means of a spacetime duality transformation. In essence, we argue that the slope coincides with the stationary density of entropy of the model obtained by exchanging the roles of space and time. Therefore, very surprisingly, the slope of the entanglement is expressed as an equilibrium quantity. We use this observation to find an explicit exact formula for the slope of Rényi entropies in all integrable models treatable by thermodynamic Bethe ansatz and evolving from integrable initial states. Interestingly, this formula can be understood in terms of a quasiparticle picture only in the von Neumann limit.

I. INTRODUCTION
The linear growth of entanglement is arguably the most distinctive and pervasive phenomenon observed in the context of quantum many-body dynamics: whenever a clean, locally interacting many-body system is prepared in a low-entangled state and then let to evolve, the entanglement between a compact region and its complement grows linearly in time. The ubiquity of this phenomenon suggests that a universal underlying mechanism is hidden behind the scenes. An astonishing outcome of recent research, however, suggests that this is not the case. Two distinct mechanisms for entanglement growth have been identified depending on the nature of the dynamics.
The first account of linear growth of entanglement has been given in the context of (1+1)-dimensional conformal field theory [1], where it was explained assuming that the entanglement is "spread" throughout the system by pairs of correlated quasiparticles produced by the quench [1]. This intuitive quasiparticle picture has then been extended to quantitatively characterise the dynamics of the standard measure of bipartite entanglement -the von Neumann entanglement entropy or simply entanglement entropy [2][3][4] -in many different kinds of systems with stable quasiparticles, such as free [1,5,6] and interacting integrable models [6][7][8] in a large variety of physical contexts. Few years later, however, the same phenomenology has been observed in systems with no quasiparticles, for instance holographic conformal field theories [9,10], generic interacting systems [11,12], and chaotic quantum circuits [13][14][15][16]. This unexpected ballistic growth of entanglement in chaotic systems has finally been explained through a "minimal membrane" picture [17,18].
In essence, the idea is that in chaotic systems the entanglement between two complementary regions is measured by the tension of the minimal spacetime surface that separates the two. Even though both quasiparticle and minimal-membrane pictures explain the linear growth of entanglement, they predict qualitatively different phenomenology when considering more complicated partitions of the system [10,17,19] or for finite sizes [14,17].
Such a twofold description of entanglement growth, however, has recently been challenged by studies on the dynamics of Rényi entropies. These are a family of seemingly minor variations of the entanglement entropy (see the definition below), which have been shown to provide highly nontrivial universal information about the system [4], for instance on its topological properties [20]. Arguably the most significant point of interest of Rényi entropies, of integer order α ≥ 2, is that they are accessible in present-day experiments [21][22][23][24][25][26][27].
Given the unquestionable success of the quasiparticle picture in quantitatively capturing the evolution of the entanglement entropy for integrable systems, it is very natural to assume that the same picture also describes the evolution of Rényi entropies. Several basic facts support this scenario: (i) it holds for free systems [1,5,28], (ii) there is no clear qualitative difference between the numerically computed Rényi entropies and the entanglement entropy [29,30], (iii) in chaotic systems the membrane picture describes both von Neumann and Rényi entropies [17,18]. The extension of the quasiparticle picture to describe Rényi entropies in the presence of interactions, however, proved to be very challenging [29][30][31][32]. In fact, Ref. [33] showed that no consistent quasiparticle picture can describe the evolution of Rényi entropies in an integrable quantum cellular automaton. A possible explanation of these findings is that, although the quasiparticle picture describes the evolution of entanglement entropy also in the presence of interactions, it fails to describe the growth of Rényi entropies. This would highlight a very unexpected fundamental difference between the two quantities, which complements the accounts of sub-linear growth of Rényi entropies in certain systems with diffusive conservation laws [34][35][36].
Motivated by this question here we investigate the dynamics of Rényi entropies in one-dimensional quantum many-body systems using a radically different approach. Our main idea is to argue that, if the roles of space and time are exchanged, the slope of a given Rényi entropy is mapped to the density of the same entropy in an appropriate steady state. This essentially means that the exchange of space and time -which we dub "spacetime swap" -maps the calculation of a non-equilibrium quantity into that of an equilibrium one.
To demonstrate the validity of aforementioned correspondence under spacetime swap we begin considering locally interacting systems in discrete space time: the so called local quantum circuits. Indeed, as recently pointed out in Ref. [37], in these systems the correspondence can be established rigorously using the spacetime duality method introduced in Ref. [14] (see also [15,33] for further developments). In particular, for dual-unitary circuits [38], where the dynamics from a class of compatible initial states [14,15] are essentially invariant under the exchange of space and time, one has that the slope of a given Rényi entropy slope coincides with the entropy density of the infinite-temperature state -the stationary state of the space evolution.
Then we consider another class of systems where the dynamics are essentially invariant under an appropriate spacetime swap: relativistic quantum field theories. In this case we show that the correspondence holds in the free case. Assuming that it continues to hold for interacting integrable quantum field theories [39][40][41] when evolving from appropriate compatible initial states [42], we arrive at a formula for the slope of Rényi entropies in all such systems.
Finally we extend our result to all interacting integrable models treatable by thermodynamic Bethe ansatz (TBA) [40,[43][44][45][46] and evolving from compatible initial states [47]. A significant physical insight of our result is that, as for the special case of Ref. [33], the dynamics of Rényi entropies cannot be understood in terms of a consistent quasiparticle picture.
The rest of this manuscript is laid out as follows. In Sec. II we define more precisely the setting considered and the quantities relevant for our analysis. In Sec. III we demonstrate the correspondence of slopes and densities of Rényi entropies under spacetime swap in local quantum circuits, where the discreteness of spacetime allows for a rigorous treatment. In Sec. IV we discuss the case of relativistic quantum field theories. In particular, in Sec. IV B we present our exact formula for the slope of all Rényi entropies in interacting integrable quantum field theories. In Sec. V we derive a direct generalisation of this result to describe general TBA-integrable systems, and test it against exact analytical and numerical results. In Sec. VI we consider the implications of our formula for the validity of the quasiparticle picture, and in Sec. VII we discuss a possible extension to the finite-subsystemregime. Finally, Sec. VIII contains our conclusions. A number of technical points are relegated to the appendices.

II. SETTING
In this work we consider a quantum many-body system prepared in a non-equilibrium initial state, |Ψ 0 , which is pure and has low entanglement. At time t = 0 we let the system evolve under its own unitary dynamics, so that the state at time t > 0 is given by where U is the time-evolution operator. As a result of the unitary evolution the state becomes increasingly more entangled as time advances [2]. The entanglement between a finite region A and the rest of the system can be quantified computing the Rényi entropies S (α) where ρ A (t) is the density matrix of the system reduced to the subsystem A. In the limit α → 1 the above expression is reduced to the von Neumann entanglement entropy At times that are short compared to the subsystem size |A|, Rényi entropies are expected to grow linearly (at least in the systems of interest here), while at sufficiently large times they saturate to their "thermodynamic values", i.e., they coincide with the Rényi entropies of the stationary state describing the subsystem A. These two regimes can be respectively characterised by the entanglement slope s α , and the stationary entanglement density d α .
The density d α is defined as the density of Rényi-α entanglement entropy of a finite subsystem of a thermodynamically large system in the t → ∞ limit, i.e., where |A| is the size of the subsystem A and L that of the total system. The final |A| → ∞ limit is taken to remove the boundary effects and focus on the bulk physics of the subsystem A.
The asymptotic slope s α , sometimes also referred to as "entanglement production rate", is defined as the ratio of the Rényi-α entanglement between a large subsystem and the rest, and time where the factor of 2 accounts for the fact that the subsystem A has two edges through which it develops correlations with its complement. In analogy with the density, we take the t → ∞ limit to remove finite-time effects. Our main goal is to establish a formal connection between the slope and the density based on a spacetime swap, i.e., an exchange of space and time.

III. SPACETIME SWAP FOR LOCAL QUANTUM CIRCUITS
A convenient setting for our analysis is that of local quantum circuits. These are models in discrete space where the time evolution occurs in discrete steps through local updates. The discreteness of time-evolution and the locality of interactions imply that these systems are closed under spacetime swap. Namely the "dual system" obtained exchanging the roles of space and time in a quantum circuit is still a quantum circuit, although the unitarity of the time evolution is generically not preserved [48]. This simple observation gives a tool to analyse several properties of quantum circuits by "evolution in space" [14,38,[48][49][50][51]. In particular, as we now discuss, it can be used to write explicit expressions for the entanglement slope that closely resemble those for the density [14,33,52].
For the sake of clarity we consider brickwork circuits -i.e., circuits consisting of two-site gates applied to first even and then odd pairs of neighbouring sites -but with minor modifications the argument can be repeated for any discrete-spacetime model with local unitary interactions.
More specifically we consider a chain of 2L sites, hosting qudits with d internal states, and where the time evolution operator U L is written as the tensor product of L two-site unitary gates U multiplied by the same product shifted by one site, i.e., where Π 2L represents a periodic shift by one site on a chain of 2L sites. At time t the reduced density matrix of a subsystem A is given by where we repeatedly apply the time-evolution operator U to the initial state, and then trace over the rest of the systemĀ. This can be represented graphically as a 2L × 4t tensor network, see Fig. 1.
The same tensor network can be equivalently thought of as resulting from an evolution in space. Indeed, rather than viewing gates as acting on the 2L sites arranged horizontally and propagating upwards/downwards, one can imagine them acting on the 4t sites arranged vertically and propagating rightwards/leftwards. To this end we introduce two different space transfer matrices playing the role of evolution operators in space (cf. Fig. 1): W t , that describes the space evolution in the complement of the subsystem A, and W A,t , that acts both on 4t temporal sites and on the subsystem A. Using these definitions, we can express the reduced density matrix ρ A as the trace of a large power of W t multiplied by W A,t , i.e., Due to the unitarity of the time-evolution (see, e.g., [15,[53][54][55][56]) the transfer matrix W t has a single nondegenerate eigenvalue 1, while all the other eigenvalues are 0. The latter generically correspond to nontrivial Jordan blocks of size smaller or equal to t. Therefore, when L is larger than t, the matrix power can be replaced with a projector on the fixed points |M L,t and |M R,t , i.e., the left and right eigenvectors corresponding to eigenvalue 1 [54][55][56]. Namely where we chose the normalisation such that Eq. (10) shows that the fixed points completely capture the effect ofĀ on the finite subsystem A (cf. Fig. 1). For this reason, and to stress their connection with the Feynman-Vernon influence functional, Ref. [57] (see also [58][59][60][61]) proposed to dub them influence matrices. The fixed points |M L,t and |M R,t can be thought of either as vectors in the space of 4t temporal sites, or as matrices M L,t , M R,t mapping from 2t temporal sites in the bottom half to 2t sites in the top half [62] (see r.h.s. of Fig. 1 for an illustration). The latter perspective makes fixed points convenient to access s α .
To demonstrate this we begin by observing that, for n integer, tr ρ n A can be represented as n copies of conjugate pairs of the time-evolved initial state that are coupled in a staggered fashion. In the section corresponding to A the pairs are connected, while in A the conjugate copy of a pair is connected to the non-conjugate copy of the next pair -see the left panel of Fig. 2 for a pictorial representation. This means that tr [ρ n A (t)] can be expressed in terms of products of n copies of the space transfer matrix W t as Diagrammatic representation of the reduced density matrix. The reduced density matrix ρA is obtained by evolving an initial state in time and tracing over the complement A of the subsystem A (graphically denoted by connecting the top and bottom legs). Alternatively, we can understand it as a trace of a product of powers of the space transfer matrix Wt and the transfer matrix WA,t (both shaded in grey), as given by Eq. (8). In the limit L → ∞, the section of the tensor network corresponding to the rest of the system can be replaced by fixed-points ML,t| and |MR,t of the space transfer matrix Wt (cf. Eq. (10)).

FIG. 2. Schematic illustration of tr ρ 3
A . Green and red rectangles are condensed representation of the green and red half of the time-evolution from Fig. 1. The white connections on the left represent trace over the rest of the system, while the connections on the right correspond to matrix products and then an overall trace. In the limit of large L − |A| the left and right parts of A can be substituted by fixed points M * L,t and MR,t connecting pairs of time-sheets. Analogously, when the subsystem size becomes large, the section in the middle can be replaced with fixed points M † L,t , and M T R,t . Note the additional transpose of fixed points, which is the consequence of connecting the opposite parity of pairs of neighbours. Thus we obtain the right-most diagram, which corresponds precisely to the r.h.s. of Eq. (13).
where η 2n represents a shift for one copy in the space of 2n replicas, and the complex conjugate of the transfer matrix comes from the exchange of conjugate and nonconjugate copies (cf. Fig. 2).
For L − |A| and |A| both larger than t the powers of the transfer matrix can again be replaced by the fixed points. In this way we obtain which is schematically depicted on the r.h.s. of Fig. 2. Thus, we find the following succinct expression for the slope (5) where we analytically continued the matrix power.
On the other hand, the stationary entropy density (4) is expected to coincide with the Rényi entropy of the reduced density matrix of the stationary state, i.e., A comparison between the expressions (14) and (15) shows that the slope can be written as a density of Rényi entropy as follows where we introduced the pseudo density matrix To gain physical intuition on the meaning ofρ st,t let us begin considering M L,t and M R,t . These matrices are by definition the fixed points, or stationary states, of the space evolution (from left to right and right to left respectively). Moreover, as shown in Appendix A, they fulfil implying that they are Hermitian and positive. This fact can be understood by recalling that, although not unitary, the evolution in space is a hybrid quantum evolution [48], i.e., it preserves positivity and Hermiticity. Eq. (18) guarantees thatρ st,t has real, non-negative eigenvalues, and, moreover, the normalisation condition (11) gives Combining all the above facts together we see that, even thoughρ st,t is not Hermitian and hence it cannot be interpreted as a proper quantum mechanical density matrix, it has many properties of reduced density matrices [63]. We note that a correspondence similar to (16) between entanglement slope in the original model and "steadystate entanglement" in the dual model has been recently unveiled in Refs. [37] and [48]. In particular, although considering a slightly different setting (they studied entanglement evolution in the presence of edge decoherence), these references expressed the entanglement dynamics in terms of states in the time lattice that in our language correspond to |M L,t and |M R,t . Here, however, we provided two significant advances. First, we showed how to establish the correspondence when the evolution in the original model is purely unitary. Second, we provided a direct interpretation of |M L,t and |M R,t as fixed points of the space transfer matrix (or influence matrices). This, in turn, allows us to viewρ st,t as a stationary state of the space evolution.
In the special case when right and left fixed points coincide,ρ st,t is Hermitian, positive definite, normalised to one, and invariant under the space evolution. Therefore, it is a stationary density matrix of the space evolution. This happens, for instance, for dual-unitary quantum circuits [38] evolving from compatible -or solvable -initial states [14,15]. Indeed, these systems and states are designed in such a way that the evolution in space is again a unitary brickwork quantum circuit, and, therefore, is equivalent to the time evolution. In particular, we have [14,15] where 1 x is the identity matrix acting on x qudits. Therefore, in this case the pseudo density matrix coincides with the infinite temperature state, i.e., Note that this is indeed the stationary state reached by the subsystem [0, t] of the time lattice under space evolution. Thus, for dual-unitary circuits (16) is rewritten as whered α in (16) is the stationary density of entropy in the dual model. For dual-unitary circuits one can also repeat the above reasoning to compute the entanglement growth between a subsystem of the time lattice and the rest. This givess

IV. SPACETIME SWAP IN RELATIVISTIC QUANTUM FIELD THEORIES
Let us now change the setting and consider quantum field theories in 1+1 dimensions, i.e., generic quantum systems defined in a continuous spacetime. Since in these systems space and time are both continuous they are again closed under spacetime swap. Therefore, we expect that one can again establish a direct correspondence between slope and density of Rényi entropies.
In 1+1 dimensional quantum field theories, however, it is more convenient to perform a slight variation of the spacetime swap. Specifically, instead of directly exchanging space and time here we consider the following analytic continuation which corresponds to an exchange of space and time in the Euclidean formulation of the theory [40]. In the string-theory literature the dual model obtained via the mapping (24) is often referred to as the mirror model [46,64,65].
Our key observation is that if one considers a 1+1 dimensional relativistic invariant quantum field theory, crossing symmetry implies that the mirror model coincides with the original one in the bulk. Therefore, in this setting 1+1 dimensional relativistic quantum field theories play a role similar to the dual-unitary circuits considered in the previous section. This means that, assuming appropriate compatible initial states, one should have where the quantities with the tilde denote the slope and density of Rényi entropy in the mirror model. To be more concrete and gain some intuition we begin by proving (25) for non-interacting, fermionic 1+1 dimensional quantum field theories.
A. Proof of (25) for free theories Let us focus on a non-interacting quantum field theory of fermions in 1+1 dimensions. We consider the quench problem with the system initially prepared in the Gaussian state where ψ † (µ) is a creation operator for fermionic modes of rapidity µ, |0 is the vacuum state for fermions, and K(µ) is an odd function such that Since a state of the form (26) produces pairs of correlated quasiparticles, the asymptotic slope (5) can be computed using the quasiparticle picture [1] Here is the occupation, or filling function, of the free mode with rapidity µ and we used the explicit form of the relativistic dispersion relation where m is the mass of the fermions and we set the speed of light to one. Next, we recall that for a stationary state described by a filling function ϑ(µ) the density of Renyi entropy reads as We then proceed by observing that at the level of rapidities the transformation (24) becomes [46] µ therefore in the mirror model the occupation corresponding to ϑ(µ) reads asθ Plugging this definition into the equation for the slope (28) we have where Γ is the positively oriented contour parametrised by i π 2 − x with x ≤ 0. Comparing (34) with (31) we see that the former can be interpreted as the density of entropy in the mirror model times i. We note that in the mirror model the integration is performed over Γ to keep ϑ(x) well defined and rapidly decaying.
Analogously, substituting the definition (33) into the expression for the density (31) we see that the density of entropy in the original model can be interpreted as −i times the slope of the mirror model

B. Slope of Rényi entropies in interacting integrable theories
Our next step is to use the correspondence (25) to find a prediction for the slope in the presence of interactions. To this end, a convenient setting to consider is that of massive integrable quantum field theories with a single species of excitations. We consider integrable quantum systems because for these systems the stationary density can be computed exactly using the quench action approach [29], while we focus on theories with a single species of particle excitations for simplicity. Note that for this family of integrable quantum field theories, a similar argument based on the mapping (24) has been employed in Ref. [66] to compute the expectation values of currents in stationary states.
In integrable quantum field theories the scattering is elastic and completely factorised. Therefore, it is fully determined by the two-particle scattering matrix S(µ). This is a meromorphic function of the rapidity in the physical strip S = {0 ≤ Im µ ≤ π}, fulfilling unitarity, crossing symmetry, and real analyticity A powerful method to describe the thermodynamics of interacting integrable theories is provided by the thermodynamic Bethe ansatz (TBA) [40]. In essence, with this method one describes stationary macrostates specifying the density of particle excitations that they contain. This is possible because in these systems particle excitations are stable and hence their densities are conserved.
In particular, we can again describe macrostates using the filling function ϑ(µ), which describes the fraction of available states that are occupied by the particles. Now, however, the density of available states, denoted by ρ t (µ), is not a simple Jacobian as in the non-interacting case. Because of the interactions it depends on the filling function through the following integral equation [40] Here we introduced the scattering kernel T (µ), which is given by the logarithmic derivative of the scattering matrix Since we are considering theories evolving from nonequilibrium initial states, the mirror model will have a nontrivial boundary in space. Here we are interested in the case where this boundary does not break the integrability of the theory, so that we can use the result of Ref. [29] for the density of Rényi entropy, therefore, we have to consider initial states generating integrable boundary conditions. These states are well known in the context of quantum integrability, and they are typically referred to as integrable boundary states [42]. Note that integrable boundary states need regularisation to be considered as initial states of quench problems [67][68][69][70][71]. This is essentially due to the fact that they give non-zero weight to configurations involving particles with infinite energy. A regularised integrable boundary state can be thought of as the generalisation of (26). It is obtained by replacing ψ † (µ) with the operators creating the stable excitations [39] and requiring K(µ) to satisfy instead of being odd. The filling function corresponding to such a state is then obtained as the solution to the following integral equation [67,69] ln Note that (38) and (42) imply that ϑ(µ) is even.
Having introduced the necessary formalism we are finally in a position to write an expression for s α . Our starting point is the exact expression for d α in interacting integrable models derived in Ref. [29]. Specialising it to the case at hand we can express it as where the auxiliary function x α (µ) is the solution to the following integral equation Proceeding as in the free case by substituting the definition (33) ofθ(µ) in (44) we find where we definedỹ The latter fulfils the following integral equation Where we used that, because of the crossing symmetry (37) of the scattering matrix, the kernel satisfies the following relation, Using then the correspondence (25) and rewriting everything for the original model we finally find with This concludes our derivation of s α for interacting integrable quantum field theories with diagonal scattering.

V. SLOPE OF RÉNYI ENTROPIES IN GENERIC TBA-INTEGRABLE MODELS
The argument leading to Eq. (50) can be applied to a much larger class of TBA-solvable models. An immediate generalisation is obtained by considering integrable quantum field theories with non-diagonal scattering as the sine-Gordon field theory. Indeed, since integrable boundary states also exist for these systems [42], one can directly repeat the treatment of the previous section.
In fact, the existence of integrable boundary states is not limited to field theories. Also in algebraic-Betheansatz-integrable lattice systems there exist initial states for which the system obtained by exchanging space and time is integrable [47,[72][73][74][75]. This applies most directly to integrable systems with a discrete time evolution [76]. In these systems the time evolution is generated by an integrable transfer matrix and, by taking appropriate initial states [47,[72][73][74][75], one can ensure integrability of the (boundary) transfer matrix in space [77]. The case of lattice systems with continuous time evolution can then be recovered by taking the Trotter limit [78,79],i.e., sending the discrete time-step ∆t to zero, while keeping fixed the real time t = N · ∆t with N being the number of steps.
In light of these facts here we argue that Eq. (50) can be extended to all TBA-integrable systems by a simple generalisation of the TBA description. In particular we have to account for the following modifications.
(i) Generic integrable models feature multiple species of quasiparticles [40,45]. This means that in general quasiparticles are no longer specified only by their rapidity λ ∈ R but one also needs to introduce a discrete species index n ∈ N. Effectively, this means that we have to make the replacement in the arguments of all functions. Naturally, this also means that when integrating over the rapidity, also the sum over all the possible particle species has to be performed where we followed the standard convention of reporting the species index in the subscript. Note also that the integration and summation boundaries depend on the specific model. Finally, to describe scattering among particles of different species, the scattering kernel needs to be generalised, To keep track of rapidity and species index we employ the following compact notation (ii) In generic TBA integrable systems the dispersion relation does not necessarily coincide with the relativistic one, therefore we make the replacement m cosh λ → ε(λ), m sinh λ → p(λ) .
Taking into account the modifications (i) and (ii), Eq. (50) is rewritten as where we introduced the auxiliary function Here (·) denotes a derivative with respect to the real rapidity µ and the subscript + indicates that the integral range is restricted to positive rapidities. To express (57) we implicitly used that, apart from fine tuned cases, integrable initial states produce reflection symmetric rapidity distributions [42,47,[73][74][75].
Once again (57) closely parallels the expression for the density of Rényi entropy in the post quench stationary state described by the rapidity distribution ϑ(µ). Indeed, for a reflection-symmetric ϑ(µ) we have [29] with In what follows we present strong evidence for the validity of (57) by providing four nontrivial consistency checks. In particular, in Sec. V A we show that Eq. (57) reduces to the exact free-fermion result [29] when the interaction kernel vanishes. Next, in Sec. V B we prove that in the limit α → 1 the expression recovers the quasiparticle prediction [7] for the slope of the von Neumann entanglement entropy. In Sec. V C we show that (57) agrees with the exact result of Refs. [32,33] for a specific interacting integrable model treatable by TBA, i.e., the cellular automaton Rule 54 [80]. Finally, in Sec. V D we compare Eq. (57) with exact numerical results for the XXZ spin-1/2 chain. Before that, however, we rewrite the expression (57) in an equivalent form, which is more convenient for parts of the upcoming analysis. To this end, we introduce two quantities that are very convenient in the TBA analysis of integrable systems, namely total density ρ t (λ) and dressed velocity v(λ). The former is the direct generalisation of the density of available states introduced in Sec. IV B and is defined as the solution to the following integral equation The latter is the velocity of quasiparticle excitations in the state described by ϑ(µ) and is determined by This quantity plays a crucial role in the quench dynamics, which was first observed in Ref. [81]. Inserting Eqs. (61) and (62) into Eqs. (59), and (57), we obtain the following equivalent expressions for density and slope

A. Free fermions
Our first check concerns free-fermionic systems. In this case, the absence of interactions permits ab initio calculations [1,5,28,82], which prove the validity of quasiparticle picture also for Rényi-α entropies. Therefore we expect to recover the quasiparticle result in the limit of the vanishing interaction kernel.

B. Von Neumann
To recover the prediction for the von-Neumann entanglement entropy we consider the limit α → 1 of the expressions (64), and (58). We begin by noting that which follows from the observation that in the α → 1 limit the function y 1 (µ) = 1 solves Eq. (58), combined with the standard TBA assumption of uniqueness of solutions. Evaluating the limit α → 1 in Eq. (64) we then have Here we used that the terms containing ∂ α y α (µ) cancel, and introduced s(µ) for the density of the Yang-Yang entropy, As promised, (68) recovers the quasiparticle prediction [7] for the slope of the von Neumann entanglement entropy.
Rule 54 can be understood as a quantum circuit consisting of 3-site local deterministic gates U with the following matrix elements, where we introduced the binary function χ : Time-evolution is given in two distinct time-steps which involve the gates U applied at odd or even triplets of sites (see Fig. 3 for an illustration) We recall that Π 2L is a periodic shift for one site on the chain of 2L sites. The gate U deterministically changes only the middle site depending on the state of the sites on both edges. Therefore all the local operators applied at the same step commute, and the dynamics is indeed one of a local quantum circuit, albeit with a slightly nonstandard geometry (cf. Fig. 3). This allows us to use the ideas from Sec. III to formally express the slope s α . Moreover, as demonstrated in Refs. [32,33], the solvability of the model allows for an exact calculation of the fixed points (cf. Sec. III), and hence of the entanglement slope, for a family of solvable initial states. More concretely, for a quench from the state with ϕ ∈ [0, 2π] and ϑ ∈ [0, 1], the asymptotic slope of the Rényi-α entropy reads as where y α is the only real and positive solution to the following equation To compare the above exact expression with Eq. (57) we need to recall some facts about the TBA description of the model [83,94], and a quench from the state (74) [33].
(i) The TBA description of states relevant for this quench problem involve only two species of particles, left and right movers, labelled by n = ±1.

D. XXZ model
Finally, let us consider the anisotropic spin-1/2 Heisenberg chain given by the Hamiltonian where σ x,y,z j are Pauli matrices acting at site j and ∆ is the anisotropy parameter, while the boundary conditions are assumed to be periodic. For the initial states of the quench protocol we consider the Néel state |Ψ N , and the Majumdar-Ghosh state |Ψ MG , defined as Since these states are integrable [47], we expect our result to apply. Moreover, we can efficiently characterise the late-time stationary state using the quench-action approach [98,99]. Therefore, we are able to explicitly evaluate the prediction given by Eqs. (64) and (58) (see Appendix C for additional details), and compare it with numerical data obtained through the infinite timeevolving block decimation (iTEBD) [100] method. See Appendix D for the details on the implementation. In particular, we evaluate the Rényi entanglement entropies between the two halves of an infinite chain, and then express the instantaneous Rényi slope s α (t), defined as the time-derivative of the Rényi entropy S Since the subsystem in question is half-infinite, we expect our prediction to coincide with the instantaneous slope in the t → ∞ limit, lim t→∞ s α (t) = s α .
This limit, however, cannot be accessed numerically because the linear growth of entanglement after the quench implies exponential growth of computational complexity to simulate the dynamics. Therefore, we have to compare the prediction with finite-time data. We consider the regime ∆ ≥ 1. In particular, since it is well known that for the initial states (80) and (81) the entanglement slope increases when approaching ∆ = 1 from above (see, e.g., Refs. [6,7]), we restrict ourself to ∆ > 1. The numerical results for the quench from the Néel state, and a range of different values of ∆, is shown in Fig. 4. We observe that the data exhibit complicated non-universal dynamics at short times, and then start approaching the asymptotic value. For all values of ∆ and α, the agreement between the finite-time dynamics and the asymptotic prediction are extremely good, considering the fact that the accessible times are relatively short. Note that the seemingly larger deviations seen at ∆ = 2 are due to long-wavelength damped oscillations (observed generically in integrable systems, see, e.g., Refs. [5,7]). Similarly, in Fig. 5 we test the prediction with the data for the quench from the Majumdar-Ghosh state. The numerics again matches the asymptotic slope very well.

VI. IMPLICATIONS FOR THE QUASIPARTICLE PICTURE
Here we generalise an argument presented in Ref. [33] for the case of Rule 54, to argue that the expressions (63) and (64) cannot be interpreted in terms of the quasiparticle picture of Ref. [1]. To explain our reasoning, let us begin by briefly recalling the essential ingredients of the latter.
The quasiparticle picture is based on two basic postulates [1]: (i) the initial state |Ψ 0 produces pairs of correlated (or entangled ) quasiparticles -objects propagating as free classical particles -at every point in space; (ii) at any time t ≥ 0, the entanglement between a given subsystem A and its complementĀ is proportional to the number of correlated pairs shared between the two.
Admitting that, in general, quasiparticles can come in multiple species -labelled by a positive integer nand have a nontrivial dispersion relation -parametrised by a rapidity µ -the two postulates above lead to the following evolution equation for a given Rényi entropy Here we adopted the shorthand notation of Eqs. (55) and (58), and used that, for solvable initial states, the correlated pairs are formed by quasiparticles of the same species and opposite rapidity [7]. Moreover, we denoted by v qp (µ) = v n,qp (µ), the velocity of the quasiparticles of species n and rapidity µ, and by the contribution to the Rényi entropy of a pair of quasiparticles of species n and rapidities ±µ.
To make Eq. (84) truly predictive one needs to specify v qp (µ) and s α (µ). In particular, Ref. [7] showed that one can describe the dynamics of von-Neumann entropy by making the two following assumptions: (i) s 1 (µ) is the density of entanglement entropy (cf. (69)); (ii) v qp (µ) is given by velocity of excitations on the thermodynamic macrostate describing the stationary value of local observables after the quench. The latter is fixed by Eqs. (61) and (62), where ϑ(µ) is the filling function of the relevant stationary state.
Here we do not use these assumptions and, for the moment, we compare (84) to (63) and (64) leaving v qp (µ) and s α (µ) unspecified. In particular, we consider an initial state producing a filling function of the form ϑ n (µ) = δ n,n ϑ µ ∈ [μ − δ,μ + δ], 0 otherwise, with ϑ ≤ 1 and δ 1. In this case, we see that the three equations are compatible for alln andμ only if The crucial observation at this point is that the right hand side of (88) depends nontrivially on α. Therefore, one needs to allow for an α-dependent quasiparticle velocity v qp (µ). At first sight this might seem enough to exclude the applicability of any quasi-particle picture. Indeed it is natural to require that the properties of quasiparticles have to be fixed by initial state and dynamics and cannot depend on the specific observable (e.g. on α).
Here, however, we allow for more flexibility: since Rényi entropies have a non-linear dependence on the state of the system, their stationary values are described by an α-dependent macrostate with filling function [29] ϑ where ϑ(µ) is the filling function (87). One can then wonder whether the velocity of excitations on the αdependent macrostate -obtained by solving Eqs. (61) and (62) with ϑ(µ) -coincides with (88). However, this is the case only in the limit α → 1.
Since the velocity on the r.h.s. of (88) cannot be interpreted as the velocity of the excitations on a physically meaningful macrostate, we conclude that the quasiparticle picture does not describe the dynamics of Rényi entropies, at least at the quantitative level.

VII. BEYOND THE LINEAR GROWTH REGIME
One of the benefits of the quasiparticle picture is that, just using the few assumptions recalled in the previous section, one can quantitatively determine the evolution of entanglement in a wealth of different settings. Essentially, the dynamics of entanglement becomes a problem of one-dimensional kinematics: knowledge of velocities and entanglement contributions of each species of quasiparticles is enough to immediately determine the whole dynamics of the entanglement of a finite subsystem (cf. (84)). The breakdown of the quasi-particle picture for α = 1 completely changes the game. Determining the full curve S A (t) becomes a highly non-trivial task for interacting integrable systems and our results for slope and density do not seem sufficient to achieve it. Here we present some evidence suggesting that, in fact, they might be enough.
To this end we make two minimal assumptions: (i) Each "mode" with quantum number µ evolves independently; (ii) For each mode the entanglement grows with fixed slope 2s α (µ) until it abruptly saturates to d α (µ)|A| (the factor of 2 comes from the fact that the subsystem has two edges). A way to justify the assumption of abrupt saturation is to argue that modes behave as chaotic systems following the membrane picture [17,18].
The two assumptions above lead to the following evolution equation for a given Rényi entropy which we conjecture apply at the leading order for large t and |A|. Testing (90) numerically in a standard interacting integrable model, such as the XXZ spin chain, is very hard (see the simulations for α = 1 in Ref. [7]): one cannot typically access its regime of validity in a sufficiently controlled manner to wash off all sub-leading corrections. However, it is instructive to test it for Rule 54. Indeed, since in that model there is a single mode for all α's, we can directly verify Assumption (ii), without worrying about the subtle effects of the integration over µ. Our numerical results based on tensor network simulations are reported in Fig. 6. We can clearly see that as increases the data approach our conjecture quite neatly, although we cannot exclude a different crossover close to the saturation point.

VIII. CONCLUSIONS
In this paper we investigated the growth of entanglement after quantum quenches in quantum many-body systems by characterising the evolution of the Rényi entropies of a compact subsystem. In the cases of interest here these quantities exhibit a linear growth in time followed by saturation to their "thermodynamic value", i.e., their value in the steady state. We showed that one can generically express the initial slope of a given Rényi entropy as the density of entropy in a particular state of the dual system, i.e., the system obtained swapping the roles of space and time, cf. Eq. (16). The latter state is directly expressed in terms of the fixed points of the space transfer matrix [54][55][56] -also known as influence matrices [57] -, which characterise the evolution of local observables in the thermodynamic limit.
In cases where the dual system can be interpreted as an isolated quantum many-body system -for example for dual-unitary quantum circuits [38] -the slope of a Rényi entropy is given by the density of entropy in the stationary state of the dual system. Crucially, because of crossing symmetry, this is the case also for relativistic quantum field theories, provided that the direct swap of space and time is replaced by an appropriate analytic continuation.
We used this observation to find a closed-form expression for the slope of Rényi entropies in integrable relativistic quantum field theories, going beyond what is currently achievable by known approaches such as form factor expansions [82,[101][102][103]. Moreover, we argued that this expression can be directly extended to all TBA integrable models. The most general form of our formula is reported in Eq. (57). To support the validity of Eq. (57) we showed that it reproduces the only known results for the slopes of Rényi entropies in an interacting integrable system [33], the quasiparticle picture prediction in the von Neumann limit [7], and it reduces to the correct non-interacting limit [5,28]. We also provided stringent numerical checks for several quenches in the XXZ spin-1/2 chain. We then used Eq. (57) to argue that the quasiparticle-picture does not describe the growth of Rényi entropies away from the von Neumann limit, and, finally, we proposed an extension of our formula away from the initial growth regime.
Our results have two significant merits. First, with Eq. (16) we provided a direct relation between growth of entanglement in a given isolated quantum many-body system and the spatial scaling of stationary entanglement in its dual. This complements an analogous relation discovered in Ref. [37], for systems with edge decoherence. Second, with Eq. (57) we solved the long standing open problem of computing the growth of Rényi entropies in interacting integrable models.
The results presented in this paper open many significant directions for future research. A direct question concerns the possibility of devising generalisations of our approach to treat other relevant quantities or describe more general settings. For instance, a recent point of interest in the research on quantum many-body dynamics is to understand how the entanglement is split among different symmetry sectors, with symmetry resolved entropies directly measured in experiments [104]. While the quasiparticle picture has been shown to hold for free systems [105,106], no result is available in the interacting case. An immediate question is then whether one can use the recently obtained explicit results for equilibrium states [107] to generalise our approach and access the full dynamics of symmetry resolved Rényi entropies. At the same time, it is also interesting to wonder whether our approach can be extended to inhomogeneous settings. Indeed, in this case the late-time quasi-stationary regime is still characterised using integrability via the framework of generalised hydrodynamics (GHD) [66,108], and an appropriate modification of the quasiparticle picture correctly characterises the time evolution of the von-Neumann entanglement entropy [109,110].
A second set of questions, instead, stems from our findings on the inapplicability of the quasiparticle picture to describe Rényi entropies in interacting integrable models. Indeed, as touched upon in the introduction, the fact that the entanglement is propagated by quasiparticlesrather than behaving as a membrane in the spacetime -has direct consequences on its phenomenology. These are revealed, for instance, in the qualitative behaviour of the bipartite entanglement between a disjoint region and the rest of the system [10,17,19,109], or in a system of finite size [17,109,111]. Recent studies, however, suggest that this might be the case also concerning multipartite entanglement, which is conveniently characterised by the entanglement negativity and the higher moments of the partial transpose of the reduced density matrix [112][113][114][115][116]. Using the quasiparticle picture, Ref. [117] argued that, after a quantum quench, the logarithmic negativity coincides with half of the Rényi mutual information with α = 1/2 -a similar statement holds for the higher moments [118]. On the other hand, Ref. [52] has recently shown that in general this relation holds only in the early time regime. Therefore, it is interesting to wonder what happens for interacting integrable systems beyond this regime. Another important question relates to the possibility of developing an alternative emergent picture to describe the growth of entanglement in interacting integrable models. Since these systems are ultimately de-fined by the presence of stable quasiparticles at all energy scales, it is reasonable to expect that a more complicated quasiparticle picture -perhaps involving quasiparticles propagating in the multi-replica space -will be able to account for the growth of Rényi entropies.
Finally, we mention that our formula (57) for the slope of Rényi entropies in all TBA-integrable models has not been rigorously proven here, even though the arguments we provided leave little doubts on its validity. Nevertheless, the correspondence that we established between slope in the original model and steady state entropy in the dual model provides an ideal starting point for such a rigorous proof. A direct question for future research is then to devise such a rigorous proof -for instance using the framework of algebraic Bethe ansatz. Besides the major interest that such a proof would have per se, it would also lead to a rigorous validation of the quasiparticle conjecture in the replica limit α → 1 -a problem that has been open since 2005 [1]. m R,t , m L,t as M R,t = m R,t m † R,t , M L,t = m † L,t m L,t . (A2) Note that we only considered the case with the initial state in the product form. With minor modifications, however, the argument can be repeated also for the initial state in the form of a matrix-product state (MPS), as long as the MPS transfer matrix has a unique dominant eigenvector.
Appendix B: Partially decoupled form of (58) and (60) In systems with multiple types of particle species, Eqs. (58) and (60) involve both an integral over rapidities, and an infinite sum over the particle species. However, using standard TBA manipulations [45] the equations can be put in an equivalent form, referred to as the decoupled form, so that each particle species n is only coupled to n + 1, and n − 1, which makes the set of equations simpler to solve.

(C9)
Appendix D: Details on the iTEBD simulations To perform the simulations we first build the MPS representation for the initial state. Both the Néel state and the Majumdar-Ghosh state (cf. (80) and (81)) admit a MPS representation with small bond dimension. Then we perform the dynamics by applying a second order Trotter decomposition of the time-evolution operator. We verified that a Trotter step δt = 0.1 is sufficient to ensure time-converged results. Due to the linear growth of entanglement the bond dimension of the MPS representing the time-evolved state increases exponentially with time. For this reason, at each step of the evolution we perform a truncation of the MPS using singular value decomposition keeping the largest χ max singular values. To monitor the loss of precision, we perform iTEBD simulations with increasing bond dimension χ max up to χ max = 8192 for ∆ = 2, and χ max = 4096 for other values of ∆. We then compare the data with two consecutive values of χ max and only keep the data for which the two simulations agree. This allows us to reach times of the order t < ∼ 15.