Introduction.— The last quarter century has witnessed out-of-equilibrium quantum matter passing from being a mere theoretical abstraction to the subject of an active and rapidly expanding multidisciplinary research field. To date, some of the most important challenges in the field include understanding the principles governing the emergence of equilibrium from microscopic quantum dynamics [1]–[5] and the dynamics of quantum information [6]–[12], characterising new forms of emergent universality [13], [14] and novel out-of-equilibrium phases of matter [15], [16], and achieving fault-tolerant quantum computation [17], [18].

From a theoretical perspective, the main obstacle to addressing these questions is the lack of general methods for describing interacting quantum many-body systems out of equilibrium. In generic spatial dimension, numerical approaches based on classical computation face unavoidable fundamental limitations [19]–[26], while existing (semi)analytical techniques typically rely on severe and often uncontrolled approximations [27]–[33]. Even in one dimension, where advanced tensor-network methods can access remarkably large systems [34]–[39], the simulation of quantum dynamics remains restricted to relatively short times: matrix-product-state methods are limited by the growth of entanglement entropy, while matrix-product-operator approaches face the infamous entanglement barrier [40], [41].

A recent breakthrough came with the realisation that some relevant non-equilibrium properties can be described efficiently by exchanging the roles of space and time. The key observation is that, upon interchanging space and time, the evaluation of finite-time properties in large systems can be mapped onto the study of stationary properties in finite systems [42], [43]. The price to pay is that the original unitary time-evolution operator is replaced by a generally non-unitary transfer matrix. Although this approach was originally introduced as a numerical technique [44], [45], it was subsequently recognised as a powerful framework for investigating general aspects of quantum many-body dynamics [46]–[48] and, in some cases, for obtaining exact analytical results [43], [49], [50].

A particularly fruitful arena where these ideas found application is that of interacting integrable systems, where one can describe stationary properties through exact analytical formulas [51]–[56], while dynamical properties are generically out of reach for both analytical and numerical approaches (with the exception of the late-time, ‘hydrodynamic’ regime [57]–[61] and the quasiparticle picture for entanglement entropy [6], [11], [62], [63]). In these systems the exchange of space and time, dubbed ‘space-time duality approach’ (SDA), provided new closed form predictions for quantities that were considered beyond the capacity of analytical methods such as the rate of growth of Rényi entanglement entropies [42] and the behaviour of charge fluctuations [64] after quantum quenches from integrable initial states [65]. Further research, however, revealed that the SDA suffered from an ambiguity that limited its applicability to homogeneous quenches and charge fluctuations in homogeneous initial states [66]–[69]. In essence, this ambiguity was caused by the fact that in the ‘dual system’ the momentum of quasiparticles is no longer a monotonic function of the rapidity.

Here we resolve this ambiguity by developing an ab-initio thermodynamic Bethe ansatz (TBA) treatment of systems with non-monotonic quasiparticle momentum. Building on this framework, we derive closed-form analytical expressions describing both entanglement growth and the dynamics of charge fluctuations following general inhomogeneous quenches in integrable systems. We argue that these results are exact: they reproduce all known predictions in the non-interacting limit [66], [70], as well as in the interacting quantum cellular automaton Rule 54 [71]–[73], while displaying excellent agreement with TEBD simulations in the XXZ spin chain. Furthermore, we show that only after resolving the ambiguity according to our prescription do the SDA expressions for Rényi entropies recover, in the von Neumann limit, the predictions of the quasiparticle picture [62], [74]. At the same time, for any Rényi index different from one, the resulting expressions remain genuinely beyond the scope of the quasiparticle framework, mirroring what was previously observed in the homogeneous setting [42].

Setting.— We consider a one-dimensional quantum system prepared in a spatially inhomogeneous non-equilibrium state \(\ket{\psi_0}\) and then let to evolve according to the integrable Hamiltonian \(H\). The initial state is assumed to have short range correlations, and to be integrable according to Ref. [65]. For example, one can think of a ‘bipartitioning protocol’ producing an initial state of the form \(\ket{\psi_0} = \ket{\psi_{\rm L}} \otimes\ket{\psi_{\rm R}}\), where \(\ket{\psi_{\rm L}}\) and \(\ket{\psi_{\rm R}}\) are homogeneous states on the left and right halves of the system, see Fig. 1. The homogeneous limit is obtained by setting \(\ket{\psi_{\rm L}} = \ket{\psi_{\rm R}}\).

Integrability ensures that the Hamiltonian’s spectrum is constructed in terms of single-particle excitations [75] and the equilibrium properties of the system are described by the TBA [51], [52], [76]. In contrast with noninteracting systems, however, the single particle excitations in a given many-body state experience non-trivial two-body interactions. Specifically, in the TBA framework one describes eigenstates \(\ket{\{\rho_{\rm t},\vartheta\}}\) using two functions: the ‘total root density’ \(\rho_{\rm t}(\mu)\), and the ‘filling function’ \(\vartheta(\mu)\) which are related through a set of integral equations which account for the interactions Note50?. The first represents the density of single particle states and the second the portion of these states that are occupied. The rapidity variable \(\mu\) represents a convenient parameterisation of the single particle dispersion relation \((p(\mu),\varepsilon(\mu))\) [77]. It is also useful to introduce the ‘dressed velocity’ \(v^{\rm dr}(\mu)\) representing the velocity of an elementary particle added to the state \(\ket{\{\rho_{\rm t},\vartheta\}}\) [78]. Because of the interactions, \(v^{\rm dr}(\mu)\) differs from the ‘bare velocity’ \(v(\mu)=\varepsilon'(\mu)/p'(\mu)\) of the quasiparticle and it is related to the latter via an integral equation involving \(\rho_{\rm t}(\mu)\) and \(\vartheta(\mu)\) Note50?.

The Generalised Hydrodynamics (GHD) framework [57]–[59] allows one to extend this approach to describe the stationary states reached at late times by local subsystems after inhomogeneous quenches. The quasi-stationary state at the space-time point \((x,t)\) is specified by means of the space-time dependent total root densities and filling functions, \(\{\rho_{{\rm t}}(\mu,x,t),\vartheta(\mu,x,t)\}\), fulfilling (non-linear) continuity equations involving \(v^{\rm dr}(\mu)\) Note50?. In the particular case of a bipartitioning protocol \(\{\rho_{{\rm t}},\vartheta\}\) depend on \((x,t)\) only through the space-time ray \(\zeta =x/t\), where \(x\) is the distance from the junction.

Importantly, despite the fact that they explicitly depend on time, the functions \(\{\rho_{{\rm t}}(\mu,x,t),\vartheta(\mu,x,t)\}\) only describe the state of the system at asymptotically large times, when local subsystems have already reached equilibrium. In this regime the residual time dependence is only due to the inhomogeneity of the system and is accounted for by classical hydrodynamic equations. This ‘equilibrium’ regime is entirely distinct from the genuine ‘non-equilibrium’ regime that we aim to study in this Letter and a priori retains no information on the growth of entanglement or the evolution of correlation functions. Describing the non-equilibrium regime in the presence of interactions is currently an open problem in integrable systems: here we argue that the SDA solves this problem in the case of integrable initial states.

Space-Time Duality.— The main premise of the SDA is that the non-equilibrium regime can be mapped to the stationary regime of a dual theory in which the role of space and time are exchanged. Whenever the latter is integrable (hence the requirement of an integrable initial state) one can then use TBA and GHD to characterise this equilibrium regime and therefore the non-equilibrium regime of the original theory.

The origin of this duality is better understood in terms of models with discrete time dynamics, where space and time are treated on similar footings and at the same time avoiding the complications of continuum theories [42], [64]. Consider for instance the one-point function of a one-site operator \(\mathcal{O}\) on the state at time \(t\), which is represented diagrammatically as in Fig. 1 (a). As it is clear from the diagram, this object can be computed by contracting the tensor network horizontally. Namely, one can express it using a suitably defined space transfer matrix \(\mathbb{T}\), i.e.a super-operator acting on the (doubled) temporal lattice. In particular, in the thermodynamic limit we find \[\require{physics} \label{eq:EV} \lim_{L\to\infty}\expval{\mathcal{O}}{\Psi_t} = {\rm tr}[\tilde{\rho}_{{\rm st}, t} \mathcal{O}].\tag{1}\] Here we set \(\tilde{\rho}_{{\rm st}, t}=R L\), where \(R\) and \(L\) are the (left and right) leading eigenoperators of the space transfer matrix (which are unique under mild conditions on the initial state [79]) and \(t\) is the length of the lattice in which \(\tilde{\rho}_{{\rm st}, t}\) is defined. Similarly, considering the half system reduced density matrix, one can show that its Rényi entropies fulfil (for open boundary conditions) [42] \[\label{eq:tracepowers} S_\alpha(t)= \frac{\log\tr[ \rho_{[0,\infty]}(t)^\alpha]}{1-\alpha} =\frac{\log \tr[\tilde{\rho}_{{\rm st}, t}^\alpha]}{1-\alpha},\tag{2}\] where we explicitly took the thermodynamic limit. Finally, if the original system has a \(U(1)\) charge \(Q\), we also have the full counting statistics of the current at \(x=0\) is given by [64], [66] (for open boundary conditions) \[\label{eq:FCS} Z_\beta(t)\mkern-4mu=\mkern-8mu \mel{\Psi_0}{e^{\beta Q_{[0,\infty]}\!(t)}e^{-\beta Q_{[0,\infty]}}}{\Psi_0} \mkern-6mu=\mkern-4mu \tr\mkern-4mu \big[\tilde{\rho}_{{\rm st}, t} e^{\beta \tilde{Q}_{[0,t]}}\mkern-2mu\big],\tag{3}\] where \(\tilde{Q}\) is the charge of the dual system and the subscript denotes the interval over which the charge is computed. Eqs. (1 –3 ) provide a complete description of expectation values, Rényi entropies, and charge fluctuations in the non equilibrium regime. The analogous statements for continuous time follow by taking a Trotter limit (i.e.continuum limit for the time lattice). Note that in models with continuous space-time, similar features can be argued using the topological nature of the line connecting the twist fields (or charge lines) [80].

Interpreting \(\tilde{\rho}_{{\rm st}, t}\) as a stationary state one then conjectures the following correspondence \[\tr[\tilde{\rho}_{{\rm st}, t}^\alpha \mathcal{O}] \overset{ t \leftrightarrow x}{\longleftrightarrow} \tr[\rho_{{\rm st}, x}^\alpha \mathcal{O}], \label{eq:swap}\tag{4}\] where \(\mathcal{O}\) is an arbitrary observable and \(\rho_{{\rm st}, x}\) is the “regular” stationary state of the time evolution reduced to a subsystem of length \(x\) and \(\overset{ t \leftrightarrow x}{\longleftrightarrow}\) indicates the swap of space and time. Note that this swap bears some conceptual similarities with the mirror approaches used in relativistic quantum field theory and string theory [76], [81], [82], although, contrary to the latter approaches, the exchange is performed in real rather than imaginary time.

In practice, Eq. 4 comprises two main steps [42]. Firstly, one assumes that the filling function \(\vartheta(x,t)\) has the same functional form in the original and swapped theories. Secondly, one performs the exchange of space and time through a swap of the dispersion relation, i.e. \[\label{eq:PEswap} p(\mu) \leftrightarrow \varepsilon(\mu).\tag{5}\] This simple mapping, however, leads to a subtle ambiguity: to have a well-defined total root density \(\rho_{\rm t}(\mu)\) the standard TBA treatment requires a strictly monotonic single-particle momentum \(p'(\mu)>0\), therefore the equations are blind to \({\rm sgn}(p'(\mu))\). Since the single particle energy is generically not monotonic the TBA equations for the dual theory are not uniquely defined by the swap in Eq. 5 , as they can involve arbitrary insertions of \({\rm sgn}(\varepsilon'(\mu))\) or of \({\rm sgn}(v^{\rm dr}(\mu))\) [66]. This ambiguity was reduced in Ref. [42] by requiring the SDA predictions to reproduce the exact results of Ref. [72] for the dynamics of a special interacting integrable quantum cellular automaton. This requirement, however, does not distinguish between \({\rm sgn}(\varepsilon'(\mu))\) and \({\rm sgn}(v^{\rm dr}(\mu))\) and, as a result, it only resolves the ambiguity when the latter are the same, i.e.for homogeneous integrable quenches [66]–[68].

This issue can be resolved by deriving TBA expressions for systems with non-monotonic momentum. To this end, we perform a suitable rapidity re-parametrisation \(\mu \to f^{-1}(\mu)=\phi\), such that \(\tilde{p}(\phi) \equiv p(f(\phi))\) is monotonic (see the End Matter for details). This gives a single set of TBA equations involving factors of \(\text{sgn} \,p'(\mu)\) to account for the non-monotonicity of the momentum. Computing the relevant quantities and performing the swap in Eq. 5 this gives the following prediction for the asymptotic slope of Rényi entropies \[\label{eq:slopeSDA} \mkern-8mu \begin{align} s_\alpha &= \lim_{t\to\infty} \frac{S_\alpha(t)}{t} =\frac{1}{1-\alpha}\int \frac{{\rm d}\mu}{2\pi} \varepsilon'(\mu)\mathcal{L}_\alpha(\mu), \\ \mathcal{L}_\alpha(\mu) &= \text{sgn}(\varepsilon'(\mu)) \log\biggl[(1\!-\!\vartheta(\mu)) ^\alpha \mkern-6mu+\mkern-4mu \frac{\vartheta(\mu)^\alpha}{y_\alpha(\mu)^{\text{sgn}(\varepsilon'(\mu))}}\biggr], \end{align} \mkern-8mu\tag{6}\] where the auxiliary function \(y_\alpha\) is defined self consistently through the integral equation \[\label{eq:yEq} \log y_\alpha(\mu) = \int \frac{{\rm d}\mu'}{2\pi}K(\mu',\mu)\mathcal{L}_\alpha(\mu'),\tag{7}\] where \(K(\nu,\mu)\) is the scattering kernel of the system (logarithmic derivative of the scattering matrix). Analogously, we find that for large times the full counting statistics takes a large deviation form \[\label{eq:largedeviationFCS} Z_\beta(t) \simeq e^{t z_\beta},\tag{8}\] with exponent given by \[z_\beta=\lim_{t\to\infty} \frac{1}{t} \log {Z_\beta(t)} = \int \frac{{\rm d}\mu}{2\pi} \varepsilon'(\mu) \mathcal{L}^{(\beta)}_1(\mu), \label{eq:FCSSDA}\tag{9}\] where \(\mathcal{L}^{(\beta)}_1\) takes the same form as \(\mathcal{L}_\alpha(\mu)\) but with a modified auxiliary function to account for the charge insertion \[\log y_{1}^{(\beta)}(\mu) = -\beta \tilde{q}(\mu) + \int {\rm d}\mu' K(\mu',\mu)\mathcal{L}_1^{(\beta)}(\mu') .\] Equation 9 differs from the prediction of ballistic fluctuation theory (BFT) for the full counting statistics on equilibrium states [83]. The difference, however, is very subtle and comes from using in \(\mathcal{L}^{(\beta)}_1\) the sign of \({\rm sgn}(\varepsilon'(\mu))\), i.e.of the bare velocity, instead of that of the dressed velocity (in the state with filling function \(\vartheta(\mu)/((1-\vartheta(\mu))y_{1}^{(\beta)}(\mu)^{v^{\rm dr}(\mu)}+\vartheta(\mu))\) [66]). Note that the BFT expression was argued not to hold after a quench whenever \({\rm sgn}(\varepsilon'(\mu)) \neq {\rm sgn}(v^{\rm dr}(\mu))\), i.e.for non-symmetric or inhomogeneous initial states, because quasiparticle trajectories are generically curved [68]. Instead, Eq. 6 is precisely the expression postulated in Ref. [42]: our results imply that the latter should hold also when \({\rm sgn}(\varepsilon'(\mu)) \neq {\rm sgn}(v^{\rm dr}(\mu))\).

To test our formulas, we then focus on the case where \({\rm sgn}(\varepsilon'(\mu))\) and \({\rm sgn}(v^{\rm dr}(\mu))\) are maximally different, i.e.that of inhomogeneous quenches. As shown in the supplemental material (SM) one can describe inhomogeneous quenches by integrating Eqs. 6 and 9 over each separate fluid cell (each one with filling function \(\vartheta(\mu, x,t)\)) Note50?. The general result can also be specialised to the simpler case of the bipartitioning protocol, which allows for a direct numerical evaluation. In particular, denoting by \(-x\) the position of the bipartition point between \(\ket{\psi_L}\) and \(\ket{\psi_R}\) (cf. Fig. 1) we find that the entropy scales as \[S_{[x,\infty]}^{(\alpha)} (t) = S_{[\zeta t,\infty]}^{(\alpha)} (t) \sim t s_{\alpha}(\zeta),\] where we introduced the ‘ray’ \(\zeta=x/t\).

In the case \(x=0\) the result turns out to only depend on the ray \(\zeta=0\) (see Fig. 1). Namely, we have that \(s_{\alpha}(\zeta=0)\) is described by Eq. 6 with \(\vartheta(\mu)\mapsto\vartheta(\mu,\zeta=0)\). Instead, the solution for different values of \(x\) involves the GHD solution at different rays (see Fig. 1). By means of simple algebraic manipulations the result can be expressed as Note50? \[s_{\alpha}(\zeta) = \frac{1}{1-\alpha}\int \frac{{\rm d}\mu}{2\pi} (\varepsilon'(\mu)-\zeta p'(\mu)) \mathcal{L}_\alpha(\mu,\zeta), \label{eq:Sgeneralzeta}\tag{10}\] where \(\mathcal{L}_\alpha(\mu,\zeta)\) is obtained from the one in Eq. 6 with the replacements \(\vartheta(\mu)\mapsto \vartheta(\mu,\zeta)\) and \(y_\alpha(\mu)\mapsto y_\alpha(\mu, \zeta)\) (that is found by solving Eq. 7 with \(\mathcal{L}_\alpha(\mu)\mapsto \mathcal{L}_\alpha(\mu,\zeta)\)). Note that Eq. 10 is almost a transformation to the co-moving frame of the ray \(\zeta\), i.e.\(\varepsilon'(\mu)\mapsto \varepsilon'(\mu)-\zeta p'(\mu)\), but not exactly as \(\mathcal{L}_\alpha(\mu,\zeta)\) still involves \({\rm sgn}(\varepsilon'(\mu))\). Analogously, we have that the full counting statistics retains the form in Eq. 8 with \(z_\beta\) replaced by \[z_\beta(\zeta) = \int \frac{{\rm d}\mu}{2\pi} (\varepsilon'(\mu)-\zeta p'(\mu)) \mathcal{L}^{(\beta)}_1(\mu,\zeta)\,. \label{eq:FCSgeneralzeta}\tag{11}\]

Consistency Checks.— We begin by considering some relevant limits of Eqs. 10 and 11 . First we note that if \(|x|\geq v_{\rm max} t\), where \(v_{\rm max}\) is the maximal velocity of the quasiparticles, the equations should reduce to the homogeneous case as a consequence of the Lieb-Robinson bound. This is indeed the case: for \(|\zeta| > v_{\rm max}\) one has that \(\vartheta(\mu,\zeta)\) is either equal to the filling function of \(\ket{\psi_{\rm L}}\) or of \(\ket{\psi_{\rm R}}\), implying that it is a symmetric function of the rapidity (both \(\ket{\psi_{\rm L}}\) and \(\ket{\psi_{\rm R}}\) are integrable states [65]). Therefore the term in \(\zeta\) vanishes by antisymmetry of the integrand and one recovers the homogeneous result. For the same reason, one obtains the homogeneous result for all \(\zeta\)s in the special case \(\ket{\psi_L} = \ket{\psi_R}\). Next, we note that Eqs. 10 and 11 reproduce the exact results obtained for inhomogeneous quenches in both free fermions [66], [70] and the (interacting) quantum cellular automaton Rule 54 [72], [73]. Finally, considering \(\alpha\to 1\) we find that Eq. 10 reduces to the quasi-particle picture result for the von-Neumann entropy in interacting integrable models [74]. This is a particularly important test as the result of Ref. [74] accounts for curved quasiparticle trajectories, i.ethe fact that \({\rm sgn}(\varepsilon'(\mu)) \neq {\rm sgn}(v^{\rm dr}(\mu))\).

To test our results numerically we consider the paradigmatic example of the XXZ chain \[H_{\rm XXZ} = \sum_j \sigma_j^x\sigma_{j+1}^x + \sigma_j^y\sigma_{j+1}^y +\Delta \sigma_j^z\sigma_{j+1}^z,\] where \(\sigma_j^{x,y,z}\) act as Pauli matrices at site \(j\) and the anisotropy \(\Delta\) is a real parameter, and study bipartitioning protocols joining states of the (tilted) Néel and (tilted) ferromagnetic families \[\mkern-8mu \ket{N,\theta} = e^{i\frac{\theta}{2}\sum_j \sigma_j^x} \ket{\uparrow\downarrow}^{\otimes \frac{L}{2}}\mkern-4mu,\mkern12mu \ket{F,\theta} = e^{i\frac{\theta}{2}\sum_j \sigma_j^x} \ket{\uparrow}^{\otimes L}\mkern-8mu. \mkern-8mu\] These states are integrable and, for \(\theta\neq 0, \pi\), not symmetric under the \(U(1)\) charge \(\sum_j \sigma_j^z\) which characterises the XXZ chain. A third class of initial states which we consider is the Majumdar-Ghosh (dimer) state, \[\mkern-8mu \ket{D} = \frac{1}{2^{\frac{L}{2}}} \Big(\ket{\uparrow\downarrow} - \ket{\downarrow\uparrow}\Big)^{\otimes \frac{L}{2}}\mkern-4mu,\mkern12mu\] which is integrable and symmetric. All TBA results for homogenous quenches can be found in [84], [85]. Furthermore, in the SM we also consider bipartitions of thermal states at different temperatures, a setup which is commonly considered to study transport [57], [58].

The time evolution is implemented via time-evolving block decimation (TEBD) algorithm [34]–[36], which allows us to reach times of about \(t\sim 15\) in most quenches of interest by describing the states as matrix product states (MPS) with maximum bond dimension \(\chi_{\rm max}=1024\). We are interested in the comparison between the results of this numerical evolution and those of Eqs. 10 and 11 . For reference we also report the results obtained by replacing \({\rm sgn}(\varepsilon'(\mu))\) with \({\rm sgn}(v^{\rm dr}(\mu))\), i.e.deliberately missing the effect of curved quasiparticle trajectories. Since the initial time regimes exhibit strong oscillations which strongly affects the global slopes in Eqs. 6 and 9 within the time regimes accessible to TEBD, we numerically evaluate the instantaneous slopes \({{\rm d} S^{(\alpha)}}/{{\rm d} t}\) and \({{\rm d}\log Z_\beta}/{{\rm d}t}\). The two approaches give the same asymptotic behaviour because of the linear growth.

Considering first the case of Rényi entropies (see a representative example in Fig. 2 and further ones in the SM) we see almost perfect agreement already at short times, which is surprising as Eq. 10 is an asymptotic prediction. Moreover, although for \(\alpha=1,2\) the numerics are not able to distinguish between the predictions using \({\rm sgn}(\varepsilon'(\mu))\) and \({\rm sgn}(v^{\rm dr}(\mu))\), they clearly agree only with the former for \(\alpha\geq 3\), confirming the validity of Eq. 10 also in the non-trivial case where quasiparticle trajectories are curved. Note that the Rényi entropies with \(\alpha\geq 3\) are less dependent on the smaller Schmidt values, which are truncated in the TEBD approach, and give therefore more reliable results.

On the other hand, the numerics is less conclusive for the FCS, as seen in Fig. 3 (additional examples of quenches from both pure and thermal bipartitioned states are reported in the SM). While the numerical data clearly converges to the asymptotic results, strong oscillations make it difficult to distinguish between the two predictions, whose difference is smaller than the oscillations’ amplitude. Similar numerical issues were noted in [66] for the global slope in a homogeneous quench, suggesting that such strong early-time oscillations are a general feature of FCS evaluations.

Conclusions.— In this letter we have resolved from first principles an ambiguity affecting the space-time duality approach to integrable quenches [42], [64], [66] and used it to obtain closed-form predictions for the large-scale entanglement growth and the evolution of charge fluctuations in the presence of both non-trivial interactions and spatial inhomogeneity. We have shown that our predictions recover known exact results, whenever the latter are available, and agree very well with extensive TEBD simulations of the dynamics of the XXZ spin-\(1/2\) chain. Therefore, they seem to provide a hitherto missing exact description of interacting integrable dynamics in the presence of interactions.

Our results prompt two immediate questions for future research. The first is to provide a mathematical proof of our predictions by developing an exact characterisation of the thermodynamic properties the fixed points of the space evolution. The second is whether our approach can be generalised to other relevant probes of the non-equilibrium dynamics which are currently beyond reach, e.g.equal time two-point functions.

TBA with Non-monotonic Momentum.— In this section we obtain a generalisation of TBA equations for a model with non-monotonic momentum. For simplicity, we consider a model with a single species of quasiparticles, as the modifications accounting for multiple species are straightforward.

Within the TBA description thermodynamics is constructed by placing the theory on a ring with periodic boundary conditions, which impose the quantisation condition for the momenta \[z(\mu_i) \equiv \frac{1}{2\pi} p(\mu_i) +(\delta*\rho_p)(\mu_i) = \frac{n_i}{L},\] where \(\delta\) is the scattering phase, \(n_i\) are integer numbers which specify the state, while \(\mu_i\) represents the rapidity variable in the parametrisation which is more convenient for each specific model, and \(\rho_p\) is the TBA particle density.

Typically one considers rapidity variables defined such that the condition \(p'(\mu)>0\) is satisfied and this can be used to define the density of states \(\rho_{\rm t}\) in the thermodynamic limit \[\frac{{\rm d}z(\mu)}{{\rm d}\mu} = \frac{1}{2\pi} p'(\mu) + (K*\rho_p)(\mu) \equiv \rho_{\rm t}(\mu), \label{eq:Bethe-Takahashi}\tag{12}\] where we have denoted by \(K(\mu,\nu)\) the scattering kernel of the system.

The positivity of the driving term guarantees that \({\rho_{\rm t}>0}\), as it is expected for a density. Instead, if the rapidity variable is chosen such that \(p(\mu)\) is not monotonic, Eq. 12 is not a sensible definition for a density of states. This problem can be treated in several ways, for example by considering separately each branch in which \(p(\mu)\) is strictly monotonic, and treating each of such branches as a different quasiparticle species. Here we use a different approach and note that the positivity of momentum can be restored by performing a suitable change of variables \(\mu \to f^{-1}(\mu)\equiv \phi\) with \(f\) injective (hence invertible) but not monotonic, such as a piecewise function given in Fig. 4. Indeed, letting \(\tilde{p}(\phi) = p(f(\phi))\), this implies \[\frac{d\tilde{p}}{d\phi} = {p'(f(\phi))}{f'(\phi)},\] from which we see that \(\tilde{p}'(\phi)\) can be made strictly positive by choosing a coordinate change, such that \(\text{sgn}( p'(f(\phi))) = \text{sgn}(f'(\phi))\). Note that under any rapidity reparametrisation the scattering kernel has to transform in such a way that Eq. 12 make \(\rho_{\rm t}\) transform correctly as a density. This implies the transformation \[K(\phi,\psi)= f'(\phi) K(f(\phi),f(\psi)) \text{sgn}(f(\psi))\,. \label{eq:kerneltransformationlaw}\tag{13}\]

In terms of the new rapidity variable \(\phi\), the TBA evaluation can then proceed following standard treatments. In particular, the partition function can be expressed as \[\begin{align} & \log \mathcal{Z} = \int \frac{d\phi}{2\pi} \tilde{p}'(\phi) \log(1+\tilde{\eta}^{-1}(\phi)), \\ & \log\tilde{\eta}(\phi)=\tilde{\nu} (\phi) + \int \frac{d\phi'}{2\pi} \tilde{K}(\phi',\phi) \log(1+\tilde{\eta}^{-1}(\phi')), \end{align} \label{eq:tbaPhi}\tag{14}\] where \(\tilde{\nu}(\phi)\) is the driving term specifying the TBA in the dual channel, expressed in the rapidity variable \(\phi\). At this point, it is important to recall that the reparametrisation used is not monotonic. Therefore, if one is to rewrite Eq. 14 in terms of the original rapidity variables, the sign of the coordinate change would have to be introduced. Since this sign is constructed precisely to be the same as the sign of \(p'\), this gives \[\begin{align} & \log \mathcal{Z} = \int \frac{d\mu}{2\pi} p'(\mu) \overline{\mathcal{L}}(\mu), \\ & \log\overline{\eta}(\mu)=\tilde{\nu}(\mu) + s(\mu) \int \frac{d\mu'}{2\pi}K(\mu',\mu) \overline{\mathcal{L}}(\mu'), \end{align} \label{eq:tba95dual95theory}\tag{15}\] where we have used the transformation law for the kernel in Eq. 13 together with the transformation of the integration measure, and introduced the shorthand notation \[\overline{\mathcal{L}}(\mu) = s(\mu) \log(1+\overline{\eta}^{-1}(\mu)), \,\,\, s(\mu) = \text{sgn}(p'(\mu)).\] This can be expressed in more convenient form by defining \(\eta(\mu) = \overline{\eta}(\mu)^{s(\mu)}\), which satisfies the TBA equations \[\log \eta(\mu) = s(\mu)\tilde{\nu}(\mu) + \int \frac{d\mu'}{2\pi} K(\mu',\mu) \mathcal{L(\mu')},\] where the transformation implies that \[\mathcal{L}(\mu) = s(\mu) \log(1+\eta^{-s(\mu)}).\] Note that the occupation functions are naturally defined in terms of \(\overline{\eta}(\mu)\) as \(\vartheta(\mu) = (1+\overline{\eta}(\mu) )^{-1}\), and therefore when expressed in terms of \(\eta(\mu)\) it acquires an extra sign factor \[\vartheta(\mu) = (1+\eta(\mu)^{s(\mu)} )^{-1}.\] From the partition function, it is possible to evaluate \(d_\alpha = \log \mathcal{Z}_\alpha - \alpha \log \mathcal{Z}\) where \(\mathcal{Z}_\alpha\) is the partition function obtained by rescaling the driving term \(\nu(\mu)\to \alpha\nu(\mu)\). This can be simplified as \[d_\alpha = \int \frac{d\mu}{2\pi} p'(\mu)s(\mu) \log\biggr[(1-\vartheta(\mu)) ^\alpha + \frac{\vartheta(\mu)^\alpha}{x_\alpha(\mu)^{s(\mu)}}\biggl]\] where \(x_\alpha = \log\eta_\alpha -\alpha \log\eta\) satisfies \[\log x_\alpha(\mu) = \int \frac{d\mu}{2\pi}K(\mu',\mu)s(\mu)\log\biggr[(1-\vartheta(\mu)) ^\alpha + \frac{\vartheta(\mu)^\alpha}{x_\alpha^{s(\mu)}}\biggl]. \nonumber\]

1 Review of the TBA solution of the XXZ model↩︎

Here we present a minimal self-contained introduction to the thermodynamic properties of the XXZ chain, obtained through the Thermodynamic Bethe Ansatz (TBA) approach, of which a complete review can be found in [52]. Since the theory is gapless for \(|\Delta| \leq 1\) and gapped for \(\Delta > 1\), leading to considerably different mathematical structures, we treat the two parameter regimes separately in the following subsections. For a more detailed analysis of the gapped case we refer to references [84], [85], while for the gapless case some more details can be found in the supplementary material of [64].

1.1 \(\Delta > 1\) regime↩︎

The gapped regime is conveniently studied by making use of the parametrisation \(\Delta = \cosh\eta\). The general eigenvalue problem can be solved through the Bethe Ansatz method for any system size. Here, however, we focus on the simplifications arising in the thermodynamic limit.

We consider the system placed on a ring of length \(L\) (different boundary conditions are possible, but the standard TBA analyses are performed with PBC), such that \[L,N \to \infty,N/L \text{ fixed},\] where \(N\) the total number of quasiparticles, which in the XXZ model are given by magnons (spin flips above the fully polarised reference state) and their bound states, which are labeled by a parameter \(j=1,\dots \infty\). In this limit, the macrostate describing the system is fully described in terms of two quantities, namely the total density of states for each bound state \(\rho_{{\rm t},j}\) and the filling function \(\theta_j(\mu)\), which specifies the particle density as \(\rho_j(\mu) = \theta_j(\mu)\rho_{{\rm t},j}\). While in free theories \(\rho_{\rm t}(k) = {1}/{2\pi}\) is completely independent of the filling, in the presence of interactions the two functions are connected through the Bethe-Takahashi equations \[\rho_{{\rm t},j}(\mu) = a_j(\mu) - \sum_k \left(K_{jk}*\rho_k\right)(\mu), \label{eq:BTDelta621}\tag{16}\] where \(a*b\) represents the standard convolution and we have introduced the driving term \(a_j(\mu) = {p_j'(\mu)}/{2\pi}\) and the scattering kernel \(K_{jk}\): these encode the model dependence and in the present context are given by \[a_j(\mu) = \frac{1}{\pi}\frac{\sinh(j\eta)}{\cosh(j\eta)-\cosh(2\mu)}, K_{jk}(\mu) = (1 - \delta_{jk})a_{|j-k|}(\mu) + 2a_{|j-k|+2}(\mu) + \dots + 2a_{j+k-2}(\mu) + a_{j+k}(\mu).\] Equation 16 alone is not enough to characterize both the total and particle densities within the system. In particular, the filling functions are strongly dependent on the specifics of the protocol: for example, if the system is in a GGE \(\sum_n \beta_n Q_n\) these are determined by the generalized TBA equations \[\log \frac{1-\theta_j(\mu)}{\theta_j(\mu)} = \sum_n\beta_n q_{n,j}(\mu)\, - \sum_k( K_{jk}* \log(1-\theta_j))(\mu). \label{eq:GTBA}\tag{17}\] where \(q_{n,j}(\mu)\) are one-particle eigenvalues of the charge operators \(Q_n\) on each bound state. On the other hand, after homogeneous quenches from integrable initial states the the main structure of Eq. 17 is preserved, with the only modification being the replacement of \(\sum_n \beta_n q_{n,j}(\mu)\) with a different driving term which can be extracted from the overlaps between the initial states and the Hamiltonian’s eigenstates [55]. For the standard integrable initial states considered in the literature, namely the tilted Néel state \[\ket{N,\theta} = e^{i\frac{\theta}{2}\sum_j \sigma_j^x} \ket{\uparrow\downarrow}^{\otimes \frac{L}{2}},\] and tilted ferromagnetic state \[\ket{F,\theta} = e^{i\frac{\theta}{2}\sum_j \sigma_j^x} \ket{\uparrow}^{\otimes L},\] and explicit solution of the TBA equations was presented in Ref. [84]. Labelling \(\theta_j(\mu) = {1}/{(1+\eta_j(\mu))}\), the equations determining \(\eta_{n\geq 2}(\mu)\) are the same for the different initial states, \[\eta_{n\geq 2} (\mu) = \frac{\eta_{n-1}(\mu + i \eta/2)\eta_{n-1}(\mu - i \eta/2)}{\eta_{n-2}(\mu)+1}-1\] where \(\eta_0=0\) and the value of \(\eta_1\) is what distinguishes the two states: in particular \[\eta_1(\mu) = -1 + \frac{T_1\left(\mu + i \frac{\eta}{2}\right) T_1\left(\mu - i \frac{\eta}{2}\right)}{\phi\left(\mu + i \frac{\eta}{2}\right)\overline{\phi}\left(\mu - i \frac{\eta}{2}\right)}\] where in the case of the tilted ferromagnet the functions take the form \[\begin{cases} T_1(\mu) = \cos(\mu) \left[ 4 \cosh(\eta) - 2 \cos(2\theta) \sin^2\mu + 3 \cos(2\mu) + 1 \right]\\ \phi(\mu) = 2 \sin^2\theta \sin\mu \cos\left(\mu + i \frac{\eta}{2}\right) \sin\left(\mu - i \frac{\eta}{2}\right) \\ \overline{\phi}(\mu) = 2 \sin^2\theta \sin\mu \cos\left(\mu - i \frac{\eta}{2}\right) \sin\left(\mu + i \frac{\eta}{2}\right) , \end{cases}\] while for the tilted Nèel state \[\begin{cases} T_1(\mu) = -\frac{1}{8} \cot(\mu) \left[ 8 \cosh(\eta) \sin^2(\theta) \sin^2(\mu) - 4 \cosh(2\eta) \right.\\ \left.+ (\cos(2\theta) + 3)(2 \cos(2\mu) - 1) + 2 \sin^2(\theta) \cos(4\mu) \right] \\ \phi(\mu) = \frac{1}{8} \sin(2\mu + i\eta) \left[ 2 \sin^2(\theta) \cos(2\mu - i\eta) + \cos(2\theta) + 3 \right] \\ \bar{\phi}(\mu) = \frac{1}{8} \sin(2\mu - i\eta) \left[ 2 \sin^2(\theta) \cos(2\mu + i\eta) + \cos(2\theta) + 3 \right] . \end{cases}\]

Beyond the description of the macrostate, the TBA framework also allows to characterize the fundamental excitations on top of it, which are responsible for the large space-time scale dynamics. In particular, the presence of interactions is mainly encoded in a renormalization of the energy and momentum of each quasiparticle, through the “dressing” equation \[f_j^{\rm dr}(\mu)= f_j(\mu) - \sum_k K_{jk}*(\theta_k f_k^{dr})(\mu).\] In this light of this, the TBA solution can be interpreted as the thermodynamics arising from a set of free quasiparticles characterized by momentum \(p_j^{\rm dr}\) and energy \(\varepsilon_j^{\rm dr}\), propagating at a velocity \[v^{\rm dr}_j(\mu) = \frac{\partial\varepsilon^{\rm dr}_j(\mu)}{\partial p_j^{\rm dr}(\mu)} = \frac{(\varepsilon')^{\rm dr}(\mu)}{2\pi \rho_{{\rm t},j}(\mu)} \implies 2\pi\, v^{\rm dr}_j(\mu)\rho_{{\rm t},j}(\mu)=(\varepsilon')^{\rm dr}(\mu)\,,\] where we have used that \((p_j')^{dr}=2\pi \rho^t_j\), as implicit in the Bethe-Takahashi equations.

1.2 \(\Delta \leq 1\) at roots of unity↩︎

In the gapless regime, we focus on a specific set of parameter values, for which the TBA treatment admits a considerable simplification, \[\Delta = \cos\eta,\eta = \frac{\pi}{p+1},p\in\mathbb{N}.\] It is possible to show that this parameter choice leads to a truncation in the number of bound states; in particular, the total number of quasiparticle species is \(p+1\). Defining the string parity and string length as \[\nu_j = \begin{cases} 1\text{for } j \leq p \\ -1\text{for } j =p+1 \end{cases} q_j = \begin{cases} j\text{for } j \leq p \\ 1\text{for } j =p+1 \, , \end{cases}\] and introducing the auxiliary function \[\tilde{a}^\nu_n(\mu) = \frac{\nu \sin\frac{n\pi}{p+1}}{\pi\left(\cosh(2\mu) - \nu\cos \frac{\pi n}{p+1}\right)},\] it is possible to define the driving term for the Bethe-Takahashi equations as \[a_j(\mu) = \tilde{a}^{\nu_j}_{q_j}(\mu),\] which gives \[\rho_{{\rm t},j}(\mu) = \nu_j a_j(\mu) - \nu_j\sum_{k=1}^{p+1} (K_{jk} *\rho_j)(\mu).\] In this case the scattering kernel \(K_{jk}\) has the slightly more involved expression \[K_{jk} = (1 - \delta_{q_j, q_k}) \tilde{a}_{|q_j - q_k|}^{\nu_j \nu_k}(\mu) + 2 \left( \tilde{a}_{|q_j - q_k| + 2}^{\nu_j \nu_k}(\mu) + \tilde{a}_{|q_j - q_k| + 4}^{\nu_j \nu_k}(\mu) + \dots + \tilde{a}_{q_j + q_k - 2}^{\nu_j \nu_k}(\mu) \right) + \tilde{a}_{q_j + q_k}^{\nu_j \nu_k}(\mu).\] Again, the filling functions contain the dependence on the initial state determining the quench, and can be determined case by case. Introducing the auxiliary functions \[s(\mu) = \frac{p + 1}{2\pi \cosh \left( (p + 1)\mu \right)}, d(\mu) = \log \left[ \coth \left( \frac{(p + 1)\mu}{2} \right) \right]^2,\] in the case of a quench from the Néel state the eta functions \(\eta_j(\mu)\) are determined as \[\begin{cases} \log \eta_1 &= (1 + \delta_{p,2}) s \star \log(1 + \eta_2) - d(\lambda), \\ \log \eta_n &= s \star \left[ \log(1 + \eta_{n-1}) + (1 + \delta_{p,n+1}) \log(1 + \eta_{n+1}) \right] + (-1)^n d(\lambda), \\ \log \eta_p &= s \star \log(1 + \eta_{p-1}) + d(\lambda) \delta_{p (\text{mod}2),0}\, , \\ \log \eta_{p+1} &= -\log \eta_p \end{cases}\] where \(1<n<p\). Another interesting initial state is the Dimer (or Majumdar-Ghosh) state \[\ket{D} = \frac{1}{\sqrt{2}}\left[ (|\uparrow\downarrow\rangle - |\downarrow\uparrow\rangle) \right]^{\otimes L/2}\] in which the eta functions can be characterised by a similar set of equations \[\begin{cases} \log \eta_1 &= (1 + \delta_{p,2}) s \star \log(1 + \eta_2) - d(\lambda), \\ \log \eta_n &= s \star \left[ \log(1 + \eta_{n-1}) + (1 + \delta_{p,n+1}) \log(1 + \eta_{n+1}) \right] - d(\lambda), \\ \log \eta_p &= s \star \log(1 + \eta_{p-1}) - d(\lambda), \\ \log \eta_{p+1} &= -\log \eta_p \end{cases}.\]

2 Bipartite quenches↩︎

In this section we specify the results of the main text to the typical scenario of the bipartitioning protocol, in which the GHD equation admit an exact solution. For simplicity, we restrict to the case of a theory with a single quasiparticle species, the generalisation to multiple species is straightforward. In inhomogeneous setups, we build a hydrodynamic framework by subdividing the system in mesoscopic fluid cells, and consider Eq. 6 to be valid at the level of each cell, with space-time dependence entering through the occupation functions \(\theta(\mu,x,t)\) defining the state. In the bipartite scenario, the scale invariance of the system implies that the GHD solution for the occupation functions is constant along fixed rays \(\zeta=x/t\), hence giving \(\theta(\mu,x,t) = \theta(\mu,\zeta)\).

We first recall some basic notions of Generalized Hydrodynamics which are needed to proceed (see, e.g., Ref. [60]). Firstly, the occupation functions \(\theta(\mu,\zeta')\) satisfy the convective equation \[\partial_x \theta(\mu,x,t) + v^{\rm dr}(\mu,x,t) \partial_t \theta(\mu,x,t) = 0 \implies v^{\rm dr}_{\zeta}(\mu) \partial_\zeta \theta(\mu,\zeta) = \zeta \partial_\zeta \theta(\mu,\zeta). \label{eq:GHD95occupation}\tag{18}\] Note that an equation of the form like that in Eq. 18 holds for any function of \(\theta(\mu,\zeta)\), such that \(v^{\rm dr}_{\zeta} \partial_\zeta F( \theta(\mu,\zeta)) = \zeta \partial_\zeta F( \theta(\mu,\zeta))\). Moreover, the particle and root densities satisfy continuity equations, reflecting the fact that both are conserved in integrable models, which in terms of the rays are expressed as \[\partial_\zeta (v^{\rm dr}_{\zeta} \rho_{\rm t,\zeta} ) = \zeta \partial_\zeta \rho_{{\rm t},\zeta},\partial_\zeta (v^{\rm dr}_{\zeta} \rho_\zeta ) = \zeta \partial_\zeta \rho_\zeta, \label{eq:GHD95continuity}\tag{19}\] where the rapidity dependence is left implicit. For a bipartite quench, the GHD equations allow for an exact implicit solution for the filling functions in terms of the occupations \(\theta_L(\mu)\) and \(\theta_R(\mu)\) of the two reservoirs \[\theta(\mu,\zeta) = \theta_{\text{L}}(\mu) \Theta\big( v_\zeta^{\text{dr}}( \mu) - \zeta \big) + \theta_{\text{R}}(\mu) \Theta\big( \zeta - v^{\text{dr}}_\zeta(\mu) \big). \label{bipartite95occupation}\tag{20}\] Note that this is an implicit solution, as the dressed velocity depends itself on the value of the filling function. This expression allows, however, for a very efficient numerical evaluation.

We are now in a position to continue with our discussion of SDA results for bipartite quenches. The key idea is that each ray \(\zeta\) contributes a factor which is precisely the density of entanglement at the given ray in the swapped theory, given by 6 specialized to the inhomogeneous values for \(\rho_{\rm t,\zeta}\) and \(\theta_\zeta\) obtained through the GHD solution. This allows us to write \[\begin{align} s_\alpha(\zeta)= \log\mathcal{Z}_\alpha^{(\zeta)} - \alpha \log \mathcal{Z}^{(\zeta)}= \int \frac{d\mu}{2\pi} v^{\rm dr}_{\zeta'}(\mu)\rho_{t,{\zeta'}} \Biggl\{ \theta(\mu,\zeta') \log y_{{\zeta'},\alpha} + s(\mu)\log\left((1- \theta(\mu,\zeta'))^\alpha + \frac{\theta^\alpha(\mu,\zeta')}{y_{{\zeta'},\alpha}^{s(\mu)}}\right) \Biggl\}, \end{align}\] where we have expressed the local entropy slopes in terms of partition functions in the language of [42], [56], see also the End Matter. The leading order solution for the full entropy is obtained by integrating over such contributions in the time interval shown in Fig. 1 of the main text. This gives \[\begin{align} &S^{(\alpha)}_{[x,\infty]}(t) = \frac{1}{1-\alpha}\int_0^t dt' \left(\log \mathcal{Z}^{(\zeta)}_\alpha - \alpha \log \mathcal{Z}^{(\zeta)}\right) = \frac{x}{1-\alpha}\int_{\zeta}^\infty \frac{d\zeta' }{(\zeta') ^2} \left(\log \mathcal{Z}^{(\zeta' )}_\alpha - \alpha \log \mathcal{Z}^{(\zeta' )}\right). \end{align}\] Considering the first term we have \[\begin{align} - x \int_{\zeta}^\infty \frac{d\zeta' }{(\zeta') ^2} \log \mathcal{Z}^{(\zeta' )} &=& x \int_{\frac{x}{t}}^\infty \frac{d\zeta' }{(\zeta') ^2} \int \frac{d\mu}{2\pi} \varepsilon'(\mu) s(\mu) \log(1-\theta(\mu,\zeta' ))\\ &=& x \int d\mu s(\mu) \int_{x/t}^{\infty} \frac{d\zeta' }{(\zeta') ^2} v^{\rm dr}_{\zeta'} \rho_{t,\zeta' } \left\{\log(1-\theta(\mu,\zeta' )) + \theta(\mu,\zeta' ) (\tilde{\nu}_{\zeta}-\log \overline{\eta}_{\zeta'}) \right\}, \label{integraltosolve} \end{align}\tag{21}\] where in the second line we have used the relation between dressed velocity, total density and dressing of the derivative of the energy discussed in the previous section, i.e., \[\frac{\varepsilon'(\mu)}{2\pi} = v^{\rm dr}_{\zeta'} \rho_{\rm t,\zeta' } + \int d\mu K(\mu,\mu') \theta(\mu,\zeta')v^{\rm dr}_{\zeta'} \rho_{\rm t,\zeta' }, \label{eq:dressedmomentum61vrho}\tag{22}\] together with Eq. 15 for the TBA in the dual theory (defined by the driving term \(\tilde{\nu}_\zeta\) which is fixed by the requirement that the occupation functions are the same in the two theories). For notational convenience we will refer to the term in the curly bracket as \(F( \theta(\mu,\zeta'))\). Using the GHD identities discussed above, Eq. 21 can be solved by partial integration, giving \[x \int_{\zeta}^{\infty} \frac{d\zeta' }{(\zeta') ^2} v^{dr}_{\zeta'} \rho_{t,\zeta' } F( \theta(\mu,\zeta')) = -x\left[\frac{1}{\zeta' }v^{\rm dr}_{\zeta'} \rho_{t,\zeta' } F( \theta(\mu,\zeta'))\right]_{\zeta}^\infty + x \int_{\zeta}^{\infty} \frac{d\zeta' }{\zeta' }\partial_{\zeta'} \left( v^{\rm dr}_{\zeta'} \rho_{t,\zeta' } F( \theta(\mu,\zeta'))\right).\] The last term is simplified by making use of Eqs. 18 and 19 , i.e. \[\partial_{\zeta'} \left( v^{\rm dr}_{\zeta'} \rho_{t,\zeta' } F( \theta(\mu,\zeta'))\right) =\zeta' \partial_{\zeta'} (\rho_{t,\zeta' } F( \theta(\mu,\zeta'))),\] allowing us to simplify the integral to \[x \int_{x/t}^{\infty} \frac{d\zeta' }{(\zeta') ^2} v^{\rm dr}_{\zeta'} \rho_{t,\zeta' } F( \theta(\mu,\zeta')) = t(v^{\rm dr}_{\zeta} - \zeta) \rho_{\rm t,\zeta} F( \theta(\mu,\zeta)).\] Finally, putting all together we obtain \[-x \int_{\zeta}^\infty \frac{d\zeta' }{(\zeta') ^2} \log \mathcal{Z}^{(\zeta' )} = t \int d\mu s(\mu) (v^{\rm dr}_{\zeta} - \zeta) \rho_{\rm t,\zeta} F( \theta(\mu,\zeta)) = t \int \frac{d\mu}{2\pi} s(\mu) (\varepsilon'(\mu) - \zeta p'(\mu)) \log(1- \theta(\mu,\zeta)),\] where we have used Eq. 22 and \((p')^{\rm dr}_{\zeta}=2\pi \rho_{\rm t,\zeta}\). Note that an identical discussion can be performed for the term \(\log \mathcal{Z}_\alpha\): in this case, it is possible to view it as the partition function of a model with driving term \(\alpha \tilde{\nu}_{\zeta}\), and, therefore, to define \(\alpha\)-dependent effective velocity and occupation function fulfilling Eqs. 18 and 19 . This allows us to conclude \[(1-\alpha) S^{(\alpha)}_{[x,\infty]}(t) = t \int \frac{d\mu}{2\pi} s(\mu) (\varepsilon'(\mu) - \zeta p'(\mu)) \left(\alpha \log(1- \theta(\mu,\zeta)) - \log(1-\theta^{(\alpha)}_\zeta(\mu))\right),\] where the integrand can be rewritten as \[\alpha \log(1- \theta(\mu,\zeta)) - \log(1-\theta^{(\alpha)}_\zeta(\mu)) = \log\left((1- \theta(\mu,\zeta))^\alpha + \frac{\theta^\alpha_\zeta(\mu)}{y_{\zeta,\alpha}^{s(\mu)}}\right),\] with the auxiliary function \(y_{\zeta,\alpha}\) satisfying a non-trivial ray dependent dressing \[\log y_{\zeta,\alpha} = \int d\mu' K(\mu,\mu') s(\mu') \log\left((1- \theta(\mu,\zeta))^\alpha + \frac{\theta^\alpha_\zeta(\mu)}{y_{\zeta,\alpha}^{s(\mu)}}\right),\] reproducing the form discussed in the main text. Exploiting the latter equation, \(S^{(\alpha)}_{[x,\infty]}(t)\) can be further rewritten in a way that is closer in spirit to the formulation of Ref. [74] \[S^{(\alpha)}_{[x,\infty]}(t) = t \int d\mu\, s(\mu) (v^{\rm dr}_\zeta(\mu) - \zeta)\rho_{\rm t,\zeta }\left\{\log\left((1- \theta(\mu,\zeta))^\alpha + \frac{\theta^\alpha_\zeta(\mu)}{y_{\zeta,\alpha}^{s(\mu)}}\right) + s(\mu) \theta(\mu,\zeta) \log y_{\zeta,\alpha} \right\}.\] Note that the solution for a generic ray \(\zeta\) is essentially the same as the one for the ray \(\zeta=0\) after moving to the comoving frame. Moreover, since in the \(\alpha \to 1\) limit \(\log y_{\zeta,\alpha}^{s(\mu)}\) vanishes just as in the homogeneous case, we see that in the Von Neumann limit the result reduces to \[S_{[x,\infty]}(t) = -t \int d\mu s(\mu) (v^{dr}_\zeta(\mu) - \zeta)\rho_{\rm t,\zeta}\left[(1- \theta(\mu,\zeta))\log(1- \theta(\mu,\zeta)) + \theta(\mu,\zeta) \log \theta(\mu,\zeta)\right].\] which is the quasiparticle-picture formula proposed in Ref. [74].

Note that this also gives the correct prediction in the homogeneous limit: in that case the term in multiplying the global factor \(\zeta\) vanishes (the filling functions are even for integrable states) giving \[S^{(\alpha),\rm hom}_{[x,\infty]}(t) = t \int d\mu s(\mu) v^{dr}_\zeta(\mu)\rho_{\rm t,\zeta}\left\{\log\left((1- \theta(\mu,\zeta))^\alpha + \frac{\theta^\alpha_\zeta(\mu)}{y_{\zeta,\alpha}^{s(\mu)}}\right) + s(\mu) \theta(\mu,\zeta) \log y_{\zeta,\alpha} \right\}\] which corresponds to the result of [42]. As a final comment, we note that for values of \(|\zeta|\) which are larger than the maximum possible value of the dressed velocity, the solution reduces to the homogeneous one with occupation function \(\theta_L(\mu)\) or \(\theta_R(\mu)\) depending on wether \(x<0\) or \(x>0\) respectively. This is a natural implication of Lieb Robinson bounds, which imply that the nontrivial physics takes place only within the light cone.

The discussion regarding the FCS of the current can be carried on similar grounds. First one assigns the contribution in Eq. 9 of the main text to each fluid cell and then integrates over the corresponding time interval. Manipulating the integral by directly repeating the procedure presented here leads to Eq. 11 of the main text.

3 Details on the numerical analysis↩︎

The numerical evaluation is performed implementing the TEBD scheme by means of the iTensor Julia library [39]. Representing the initial state as a Matrix Product State (MPS) the time evolution is approximated via a Trotter-Suzuki decomposition of the time evolution operator, which can be taken as the forward-backward swift \[{\rm{U}}(dt) = \prod_{i=0}^{L} e^{-ih_{i,i+1} dt} \prod_{i=L}^{0} e^{-ih_{i,i-1} dt}\] where each term \(e^{-ih_{i,i+1} dt}\) is considered as a quantum gate acting on two spins. This produces errors of order \(dt^3\) per time step, and allows us to obtain good numerical results for \(dt = 0.01-0.1\). We consider a global system of size \(L=80\) with open boundary conditions, and a subsystem of size \(L_A=40\) placed on the left half of the system.

Truncations are performed to avoid an excessive growth of the MPS. In particular, we use a cutoff of \(10^{-12}\) on the Schmidt values at each step and a maximum bond dimension \(\chi=1024\): these allow us to reach values of \(t\sim 10\) by keeping a good precision in the description of the state, with variations depending on the entanglement slope for each quench.

The cutoff of smaller Schmidt values has the consequence that the lower the Rényi index, the larger will be the cumulative effect of the truncation in time. This implies that larger values of \(\alpha\) allow to minimize such numerical effects, leading to more reliable numerical results. Moreover, it is convenient to consider the instantaneous slope \({\rm d}S_0^{(\alpha)}/{\rm d} t\), rather than the global one \(S_0^{(\alpha)}(t)/t\). This is useful because the global slope is affected by the large spike present in the initial time regime, which eventually gets canceled on time scales that are larger than those accessible in the TN evaluation. The TEBD result is then compared with the analytical predictions obtained in the main text. Representative examples of our results for the slope of Rényi entropies are presented in Figs. 5, 6, and 7, as well as Fig. 2 of the main text. The figures compare the TEBD results with our prediction in Eq. 10 of the main text evaluated at \(\zeta=0\). For reference we also report the prediction obtained by replacing \({\rm sgn}(\varepsilon'(\mu))\) with \({\rm sgn}(v^{\rm dr}(\mu))\) in in Eq. 10 , i.e.deliberately missing the effect of curved quasiparticle trajectories. The results at higher \(\alpha\) clearly show convergence to the prediction of Eq. 10 rather than to the modified one.

The evaluation of the FCS can be performed in a similar manner. In this case the initial state is written as \[e^{-i\beta Q_A(0)} \ket{\psi_0}.\] Since the bond dimension of the operator is just \(\chi_{MPO}=1\), this does not appreciably increase the complexity of the evolution through TEBD.

We perform this analysis both for pure and mixed (thermal) bipartite states obtaining the results reported in Fig. 8, as well as Fig. 3 of the main text. Although the results show convincing convergence to the analytical prediction in Eq. 11 , the precision is not enough to distinguish between the latter and its variant obtained by replacing \({\rm sgn}(\varepsilon'(\mu))\) with \({\rm sgn}(v^{\rm dr}(\mu))\). This happens both because the difference between the two predictions is smaller compared to the one observed for the entropy, and because the FCS exhibits stronger oscillations which do not disappear completely in the time interval accessible through TEBD.

The situation is not solved by focusing on the global slope, which was considered for example in [64], [66], as the strong fluctuations in the early-time cannot be canceled completely for the relatively small values of \(t\) accessible in the simulations. This effect is clearly visible in Fig. 9.

References↩︎