Entanglement Rényi entropies from ballistic fluctuation theory: The free fermionic case

The large-scale behaviour of entanglement entropy in finite-density states, in and out of equilibrium, can be understood using the physical picture of particle pairs. However, the full theoretical origin of this picture is not fully established yet. In this work, we clarify this picture by investigating entanglement entropy using its connection with the large-deviation theory for thermodynamic and hydrodynamic fluctuations. We apply the universal framework of Ballistic Fluctuation Theory (BFT), based the Euler hydrodynamics of the model, to correlation functions of branch-point twist fields , the starting point for computing Rényi entanglement entropies within the replica approach. Focusing on free fermionic systems in order to illustrate the ideas, we show that both the equilibrium behavior and the dynamics of Rényi entanglement entropies can be fully derived from the BFT. In particular, we emphasise that long-range correlations develop after quantum quenches, and accounting for these explain the structure of the entanglement growth. We further show that this growth is related to fluctuations of charge transport, general-ising to quantum quenches the relation between charge fluctuations and entanglement observed earlier. The general ideas we introduce suggest that the large-scale behaviour of entanglement has its origin within hydrodynamic fluctuations.


Introduction
The understanding of entanglement in quantum many-body systems received a considerable boost in the last decades, with the introduction and characterization of many different quantities which "measure" the amount of entanglement in a given quantum state [1][2][3][4]. An important set of such measures are the so-called entanglement Rényi entropies. Given a quantum system described by a density matrix ρ and a subsystem A of the total system, withĀ denoting its complement, consider the associated reduced density matrix ρ A = trĀρ. Then, for any α ∈ + , the α-Rényi entropy is defined as They are good entanglement measures for all pure quantum states, i.e. states of the form ρ = |Ψ〉〈Ψ|. They fully characterise the entanglement spectrum, and an important property is that in the limit

SciPost Physics
Submission particle pairs, but the method is new, and brings out, we believe, important new physics underlying the entanglement entropies. The main two observations are: (1) We obtain an exact relation between the growth of the Rényi entanglement entropies after a so-called integrable [26][27][28], pair-production quench, and static and dynamic "full counting statistics" in the final GGE. Consider N <,> := |v(θ )| <,> ξ/2 dθ ψ † θ ψ θ the conserved quantity giving the total number of "slow" and "fast" fermionic modes ψ θ , with speeds |v(θ )| < ξ/2 and |v(θ )| > ξ/2, respectively, where ξ = x/t is a spacetime ray. Let us denote Consider for simplicity α to be even. Then, as x, t → ∞ with x/t = ξ fixed, the Rényi entanglement entropy on the interval [0, x], at time t after the quench, has asymptotic form: This extends earlier observations of the connection between entanglement entropy and full counting statistics [29][30][31] to non-equilibrium quenches. Our calculations also provide a fundamental explanation of such relations in terms of twist fields and the large-deviation theory for their asymptotic behaviours, which, as far as we know, has not been noticed before.
(2) We give a new exact derivation of the so-called quasi-particle picture in the case of free fermions with generic dispersion relation. Our derivation is completely independent from the other exact result for the Ising model in [15], which was based instead on Toeplitz matrix representation and multidimensional phase methods. In particular, our method makes transparent how it is the simple structure of long-range correlations induced by particle pairs in integrable quenches that allows one to describe both the growth and saturation of entanglement in a simple and universal way in terms of the long-time GGE, as this structure allows the separation of the contributions of fast and slow modes as per (3). The emphasis on the structure of long-range correlations also gives a clear understanding as to why for quenches starting from more complicated states, for instance producing correlated groups of more than two particles, more information about the initial state is needed to describe the entanglement growth; in these case no simple formula exists (as showed for example in [32,33]).
We concentrate on free fermion models for simplicity and in order to most clearly illustrate the method and physics. However, as the method is based on general large-deviation and hydrodynamic principles, it is expected to be much more widely applicable, which we leave for future works. In particular, it suggests that hydrodynamic modes and hydrodynamic projections are the more accurate notions at the root of the large-scale behaviour of the entanglement dynamics, rather than particles and their productions.
The paper is organized as follows. Sec. 2 is an introduction to BFT and its relation to twist fields. In Sec. 3 we review the replica approach and the associated branch-point twist fields,

SciPost Physics
Submission and we discuss the simplifications occurring in the free fermionic case. Sec. 4 is the core of the paper, where we derive an expression for correlation function of twist fields from BFT, both inand out-of-equilibrium, and use them to obtain (3) and recover the known formulas for Rényi and entanglement entropies. A discussion of our method and results is given in Sec. 5. The appendices complement the main text with observations and details of the calculations. In particular, App. A contains remarks on notions of locality and twist fields. App. B contains all the details about the applicability of BFT in the different situations we consider, by explicitly computing correlations, and their long-range behaviour, after a quench from a state with pair structure. Finally, App. C is about the structure of the S-matrix in the α-copy theory.

Ballistic Fluctuation Theory and twist fields
The BFT [24,25], detailed below, is a theory describing the large-scale, ballistic fluctuations. It is expected to apply to a large class of quantum and classical many-body, extensive systems. It applies to generic systems with space-translation invariant dynamics and interaction range that is short enough. It has been developed originally for states that are spacetime stationary and clustering in space, but many of the ideas have been extended to more general situations [34,35].
In this paper, we use the BFT as originally developed in [24,25],and show how, and under which assumptions, it can be applied to states that emerge after quantum quenches as well. Quantum quenches give rise to states that are locally spacetime stationary, but present time-varying long-range correlations, as we will explain below. We will explain how simple ideas based on the principles explained in [24] allow us to nevertheless use the BFT. We mention that long-range correlations can also, in principle, be accounted for directly by using the more sophisticated ballistic macroscopic fluctuation theory (BMFT) [35], which we leave for further studies.

General setting
The main strength of the BFT is that it stipulates that only some emergent properties of the system are required in order to describe the large-scale, ballistic fluctuations: the data of its Euler hydrodynamics. We assume the system of interest to have a certain number N (which in our application to free fermions will be infinite, as the system is integrable) of conserved quantities such that dQ i /d t = 0. They are assumed to be hermitian (in quantum systems), or real (in classical systems). They include the Hamiltonian H = d x h(x, t), which generates time translations. For simplicity of the discussion we assume these to be in involution, [Q i, Q j ] = 0 for all i, j ∈ {1, · · · , N }, however this is not necessary in general. They have associated conservation laws ∂ t q i + ∂ x j i = 0. The observables q i , j i are the corresponding charge density and current, assumed to be "local". In the present paper, locality of q i and j i simply means that Q i has appropriate extensivity properties; we keep the general discussion formal in order to avoid technicalities, but see the remark about locality concepts in the literature in App. A. Within such systems, we focus on states belonging to the manifold of maximal entropy states (MES). Each state is characterized in terms of a vector β = {β 1 , · · · , β N } of "Lagrange multipliers", with N components (there are as many number of components as the number of conserved SciPost Physics Submission quantities Q i ). Given such a vector, the density matrix defining the system reads 4 These include the GGEs studied in integrable systems. We consider the system to be in infinite volume. Below, when no ambiguity occurs, expectation values on such states will be denoted simply as 〈·〉. Importantly, we assume such states to be clustering strongly enough: connected correlations tend to zero quickly enough at large spatial separations, 〈a(x, 0)b( y, 0)〉 c := 〈a(x, 0)b( y, 0)〉 − 〈a(x, 0)〉〈b( y, 0)〉 → 0 (|x − y| → ∞).
Here and below, a(x, t) is a local observable at the position x and evolved to time t.
As per basic principles of statistical mechanics, the averages of densities are generated by the free energy and the mapping 〈q〉 ↔ β (8) is bijective (in appropriate regions of values of 〈q i 〉 and β i ).
As mentioned, states that are spacetime stationary for local observables, but with spacetime varying long-range correlations, arise naturally in quantum quenches even after long times. This is because for thermalisation to happen on large regions of the system takes a long time. In such cases, the state is not described by the MES (5). A precise description of states with long-range correlations is more involved, see e.g. [35,38]. But the concept of MES is still useful in these situations, as, nevertheless, averages of all local observables, or observables supported on regions smaller than the correlation range, still are described by (5). We will explain how to use this fact in order to "avoid" long-range correlations and apply the results of the BFT.

Large deviation theory of currents
Consider some conserved quantity Q = Q i * (for a given i * ∈ {1, · · · , N }), with associated density q = q i * and current j = j i * .
It is instructive to start with a description of the large-deviation theory for extensive charges at equilibrium, before discussing currents. In a given state ρ β , a natural question is to characterize the restriction ∆J(x) = − x 0 d x q(x , 0) of the charge Q to a spatial interval [0, x], and its fluctuations within this interval. Here x denotes the horizontal path from (0, 0) to (x, 0) (and the notation ∆J(x) is adapted to generalising to currents, as done below). The fluctuations are fully characterized by the cumulants of ∆J(x). It is a simple result that the cumulant generating function at large x is given by a difference of specific free energies f (β) of the system, Here and below, we use the notation A(s) B(s) with the meaning that lim s→∞

SciPost Physics
Submission and this generates the scaled cumulants c m , When studying transport, similarly, we are interested in characterizing the total current passing by a given spatial point (e.g., the origin) in a given interval of time [0, t]. One is therefore interested in the total transfer of Q in time t, i.e., ∆J(t) = t 0 d t j(0, t ), where now t denotes the vertical path from (0, 0) to (0, t). As argued in [24] using hydrodynamic principles, the structure parallels closely the equilibrium case in ballistic systems, such as those admitting many conserved charges. At large times t → ∞, a large-deviation principle holds generically for linear scaling with t.
In fact, one can go further and consider the 2-current j = ( j, q), and the integral along a more general path , starting in (0, 0) and ending in (x, t), over its perpendicular component to the path. This defines the following object where d = (d x , d t ) and j ∧ d = jd t − qd x . By current conservation, the integral (12) is in fact independent of the path chosen, and therefore the result only depends on the end-points (0, 0) and (x, t) (for lightness of notation, we keep implicit the dependence on (0, 0)) For example, we may choose to connect the initial and final points of the path via a segment of ray as will be done below. Let us consider the Euclidean distance between the initial and final points, (we do not assume Euclidean spacetime symmetry, this is simply a convenient way of controlling the scale of x and t). As mentioned, in ballistic systems a large-deviation principle holds generically for ∆J(x, t) ∝ . Then, as in the equilibrium case, we can define the scaled cumulant generating function The function F (λ) also depends on the angle γ, which we keep implicit. The coefficients c m (γ) are the scaled cumulants of ∆J(x, t), which are m-point correlation functions of currents and densities integrated over the path . The finiteness of scaled cumulants depends on the asymptotic behavior of density and current correlation functions at large spacetime separations. See App. B for more details.

F (λ) in MES: biased measure, flow equation
From the explanations above, at γ = π/2 (in the spatial direction) we have that F (λ) = −∆ f (λ), explicitly given in terms of thermodynamic quantities (cf. (9)). What is the corresponding quantity for finite values of γ (i.e., involving the time direction)? It turns out that the answer is completely given using the data of the Euler hydrodynamics of the model. The Euler hydrodynamics controls the motion and correlations of the many-body system at large scales of space and time. For our purposes, it is sufficient to recall that it is completely fixed by the flux jacobian matrix, defined as (using the bijection (8)) In integrable systems, the flux Jacobian is in fact an infinite dimensional matrix, or more precisely an integral operator [39].
From the definition it is clear that A i j is basis-dependent. Its fundamental information is contained in its spectrum v eff k k=1,··· ,N , which is composed of eigenvalues if N is finite, and which admits a continuum in integrable systems (where N = ∞). The spectrum is interpreted as the set of effective velocities, or "generalized sound velocities," associated to normal modes of the hydrodynamics [39].
In order to compute F (λ) in (16), the idea is to bias the measure ρ β defining the MES (where expectation values are considered), by multiplying it by the exponential of the integrated 2-current ∆J(x, t), as appears in (16). The state thus becomes λ−dependent, and it is in fact a MES, which we write as ρ λ;β . Associated to the segment of path from (0, 0) to (x, t), with tan γ = x/t, one can derive a flow equation for ρ λ;β , within the space of MES, starting from ρ 0;β = ρ β . Conveniently, this can be written as a flow equation for the Lagrange multipliers β(λ; γ) themselves as follows [24] ∂ where we explicitly introduced the dependence on λ as well as on the path (through γ) in all quantities.
The main assumption underlying the validity of (20) is that of "strong enough" clustering along the space-time ray of velocity x/t = tan γ. Specifically, spatial clustering (6) is not enough: one needs vanishing of correlation functions of perpendicular currents j ⊥ (x, t) := j(x, t) ∧ d , the integrand in (12), when the distance along the ray (0, 0) → (x, t) goes to infinity,

SciPost Physics
Submission and similar requirements for all multipoint functions of perpendicular currents. The vanishing must be fast enough to make integrals defining the cumulants rapidly converging. A crucial remark for our work, concerning the requirement of clustering, is Remark 3.3 of [24]. (12) is independent of the path , and only depends on the endpoints (0, 0) (which we have kept implicit) and (x, t). One may therefore hope to apply the general result of the BFT for each path element and obtain, in place of (20), the expression for any path (with a well-defined large-scale limit → ∞). As explained in [24,Rem 3.3], this is expected to be correct if and only if there are no strong correlations between the perpendicular currents amongst different points on the path. Eq. (22) is the main result from BFT which will be applied below to expectation values of observables (specifically, of twist fields) both at equilibrium in states described by GGEs (in Sec. 4.1), and out-of-equilibrium in states emerging after quenches (in Sec. 4.2 and the following ones). In the latter case, this can be done by choosing a path that "avoids" the long-range correlations that such states present, so that expectation values in the pre-quench state can be approximated with the ones in the corresponding long-time GGE (where correlations can be shown to decay fast enough).

Remark (clustering and the ballistic large-deviation principle).
If there is no path that satisfies strong clustering, then typically the ballistic large deviation principle is broken, and the SCGF resulting from the BFT may be infinite or zero. Exponential clustering, which is sufficiently strong, is expected to hold on rays away from the fluid velocities, tan γ = v eff i ∀i, in generic equilibrium states. In GGEs of integrable systems at nonzero entropy density, the spectrum of the flux Jacobian contains a continuum, and clustering is in fact a power-law, with power 1/ 2 for the perpendicular currents, Eq. (21), for all rays within this continuum, see App. B.1 for the case of free fermions. This is still strong enough for the BFT results to hold, as confirmed numerically for the hard rods [25]. By contrast, in non-integrable systems, where the flux Jacobian has a discrete spectrum, clustering is too weak on rays along any of the eigenvalues of the spectrum of the flux Jacobian (fluid velocties), tan γ = v eff i for some i. In these directions, the ballistic large-deviation principle is broken. See the discussion in [24]. In general, at zero temperatures when there is no gap, weak power-law clustering is seen and the ballistic large-deviation principle is also broken. When the breaking of the large-deviation principle happens as we change a parameter (a (generalised) temperature, a coupling), this can be seen as a "dynamical phase transition".

Application of the BFT to twist fields
One of the most important application of the BFT is to two-point correlation functions of so-called "twist fields". This is useful, because, as explained in the introduction, twist fields are probably the most efficient way of studying entanglement entropy, the main object of this work.
There are a number of ways of defining twist fields, and we will discuss two natural approaches. The first is natural in the context of the large-deviation theory as recalled above and based on the explicit knowledge of extensive conserved quantities; it applies to classical and quantum systems alike. The second is a more abstract formulation that does not require the explicit knowledge of extensive conserved quantities, but that is better adapted to quantum systems.

Submission
Consider as above an extensive conserved quantity Q. Recall that Q has associated density and current q(x, t) and j(x, t). It is convenient to define the height field ϕ(x, t) via the relations This ensures that the continuity equation ∂ t q + ∂ x j = 0 is automatically satisfied. The height field is unbounded (because the charge is extensive), and we note in particular that differences grow linearly, The height field may be written as an integral over half space 5 , One finds that the boundary term at infinity is constant in time, ϕ(∞, t) = ϕ(∞) (which can be chosen to vanish). Further, it is clear that the result is independent from the choice of path in spacetime thanks to the conservation laws, and thus the height field ϕ(x, t) only depends on the point (x, t). This justifies the notation. One may also choose a different direction for the half-space integral, the difference being encoded within Exponentials of ϕ(x, t), that is are known in general as "twist fields". Because of the expression (25), they are not "local" in the naïve sense, but are usually referred to as "semilocal" in the literature (see App. A for an overview); this is made clearer below using exchange relations. As an immediate result of the large-deviation analysis, using with (24), we get the leading exponential behavior of the two-point correlations of twist fields as That is, if the ballistic large-derivation principle holds for the charge Q, then the associated twist field shows an exponential behaviour at large space-time separations. Because e λϕ is not bounded for λ ∈ , T λ should be referred to as an "unbounded twist field". In quantum models, because of the lack of commutativity of observables, and a fortiori in quantum field theory (QFT), because, additionally, of the necessary renormalisation procedure applied to the twist fields, the relation (29) is not strictly valid. However, corrections do not affect the leading exponential behaviour of the two-point function, as argued for the XX quantum chain in [41].

SciPost Physics Submission
A second viewpoint on twist fields is as follows. A twist field T (x, t) is in general an operator associated to a symmetry transformation that "acts locally enough". This is a transformation a(x, t) →ã(x, t) of the model, that maps local observables to local observables, and that preserves the operator algebra; one usually also requires that it preserves the Hamiltonian density, h(x, t) = h(x, t). The property of an operator T (x, t) that makes it a twist field associated to such a symmetry, is the following equal-time exchange relation: for every local observable a( y, t). Corrections at larges distances | y − x| should be small enough, for instance exponentially decaying. This equal-time exchange relation is known in the literature as expressing the "semilocality" of the twist field 6 T (x, t).
In general, symmetry transformations act via unitary operators U asã(x, t) = U a(x, t)U −1 . For instance, in quantum spin chains, one often can write U = x∈ U x where U x acts non-trivially on a neighbourhood of x (possibly up to exponentially decaying corrections); this is indeed local enough. In such instances, one can simple set How does the exchange-relation formulation (31) connect with the height-field formulation (28) of twist fields? This is from the general principle that every extensive conserved quantity gives rise to a continuous, one-parameter unitary group of symmetry transformations that act locally enough. That is, given an extensive Q (recall that it is assumed to be hermitian), we may consider the symmetry transformationã (x, t) = e iηQ a(x, t)e −iηQ (33) for a real parameter η ∈ . It is simple to see that this acts locally enough, as described above. Indeed, by conservation we can replace Q by Q(t) in (33), and by locality of the density q(x, t), we have that , a(0, t)] approaches its limit as L → ∞ quickly enough, and gives in the limit a local observable supported at x = 0 at time t. Therefore, by using the Baker-Campbell-Hausdorff formula, at least for all η small enough, e iηQ a(x, t)e −iηQ gives rise to a local observable at (x, t). Using the same principles, it is then a simple matter to show that (31) holds with the choice (25) for the height field, and the identification in (28) Because e −iηϕ is bounded, we will referred to these as "bounded twist fields". Bounded twist fields are the ones usually considered in the literature. In particular, because they are bounded, their two-point functions should decay in spacetime, One may in fact argue that, viceversa, to every symmetry transformation that acts locally enough, we can associate an extensive conserved quantity. In quantum field theory, Noether's theorem shows

SciPost Physics
Submission that, for continuous symmetry groups, we indeed have U = e iηQ for some conserved quantity Q associated to a conserved 2-current; again this is local enough. In general, for any local enough transformations, one can identify, formally, an extensive charge Q with the operator −i log U (thus taking (33) with η = 1), a formal construction that, we expect, could be used fruitfully within the BFT. In the present paper, we will consider a twist field associated to a discrete symmetry transformation, but this will be embedded within a continuous symmetry group thanks to the free-fermion structure, thus the charge Q will be explicit.
Finally, we observe that applying the BFT for the bounded twist fields, using the identification (34), requires an analytic continuation of λ in the BFT formulae, to purely imaginary values −iη. This is a subtle aspect, as for purely imaginary values of λ, the modified state by the flow equation is not strictly a MES (because the resulting linear functional on the algebra of observables is not necessarily positive). We believe that if fluid velocities are well separated, as is typically the case in non-integrable systems, then the analytic continuation can be obtained meaningfully by keeping the sign of the eigenvalues constant in the flow equation (19) (as the analytic continuation will not "see" the jumps in eigenvalues), and integrating the flow in the complex λ-plane. In free fermion models, the analytic continuation can be performed directly on the explicit result for F (λ), as done in [41] and in the next section. We will show below that the BFT indeed predicts decay of correlation functions in this case. We leave the discussion for interacting integrable systems to future works. See App. A for a brief discussion of notions of locality and twist fields.

Explicit expression of F (λ) in free fermionic theories
Up to now, the theory was general, and all equations correctly give the ballistic part of large deviations in general systems with the properties detailed in Secs. 2.1 and 2.3. When considering integrable systems, Eq. (20) gets further simplified. In fact, using the theory of generalized hydrodynamics (GHD) [42,43], an explicit expression of the hydrodynamic quantities 〈q〉 λ , 〈 j〉 λ along the flow can be worked out. We now overview the simplified expression for F (λ) in the special case of free fermions in the continuum, arguably the easieast among 1D integrable systems, which is our focus in this paper (the corresponding results in the case of generic interacting integrable models can be found in [25]).
In order to keep the structure general, we simply assume that a fermionic, complex field ψ(x, t) exists with interactions that are quadratic and short-range. Its Fourier modes are denoted ψ θ , with anti-commutation relation {ψ † θ , ψ θ } = δ(θ −θ ). Here θ represents the momentum, which we assume takes values in for simplicity (for quantum chains, this would be a bounded subset instead, but the general ideas are not affected). We also denote the dispersion relation as E(θ ), which we assume is strictly convex and symmetric E(θ ) = E(−θ ). Thus, under canonical normalisation, As it is integrable, the model possesses an infinite number of conserved quantities. A "scattering basis" for these is given by Q θ = ψ † θ ψ θ , θ ∈ . Strictly speaking, the Q θ 's are not linearly extensive, but (for generic dispersion relation) any extensive conserved quantity can be obtained by a suitable "linear combination", or basis decomposition, Q i = dθ h i (θ )Q θ . Here, h i (θ ) is the one-particle eigenvalue of the extensive charge Q i . Examples are the number of fermions N = dθ Q θ , the total momentum P = dθ θQ θ , and the total energy H = dθ E(θ )Q θ .

A typical GGE (5) takes the form ρ
is the generalised Boltzmann weight in the particle basis . For example, for a thermal state, we have In general, the physical meaning of the Lagrange multipliers β i depends on the choice of the set of charges Q i , i.e. on the choice of the set of one-particle eigenvalues h i (θ ). But there is no need to choose any particular infinite set of charges Q i , or to write explcitly w(θ ) in a basis decomposition w(θ ) = i β i h i (θ ). The function w(θ ) fixes the GGE, and only few basic requirements constrain w(θ ) for ρ w to be a valid GGE (we will ask that it be positive and grow sufficiently fast as |θ | → ∞). We note that the conserved charge densities take the standard form and that, in a system of length L with periodic boundary conditions, we have 〈Q θ 〉 = L 2π n(θ ). In our calculations, we will assume that n(θ ) has an analytic extension in a neighbourhood of , and that n(θ ) → 0 as |θ | → ∞.
Consider the large-deviation problem for the charge Q = Q i * , with one-particle eigenvalue where v(θ ) = d E(θ )/dθ is the group velocity 7 . The function f (ε) is the fermionic free energy function (the free energy density per distance and per unit rapidity θ ), and the function ε λ (θ ; γ) is the Boltzmann weight along the flow (19) in the particle basis. Its initial condition is ε 0 (θ ; γ) = w(θ ), and the corresponding flow equation (which simply follows from (19) in terms of β i ) simplifies to which is explicitly solved as As mentioned above, in order to apply the BFT to bounded twist fields associated to symmetry transformations, one needs to perform an analytic continuation in λ to the purely imaginary direction λ = −iη, η ∈ . In free Fermion systems this is simple to do, as the above formulae can be directly analytically continued. The resulting F (λ) possesses, in general, both a real and an imaginary part. The real part describes the exponential decay of the two-point correlation functions of twist fields, while the imaginary part describes oscillations. In the following, we will not discuss oscillations, as their full description would require a more in-depth analysis; we will concentrate on the exponential decay, hence the real part of F (λ).
Correlation functions of twist fields are expected to be decaying at large spacetime distances. It is simple to show from (38) that indeed 8 One important remark is that the only information required about the current ∆J(x, t) whose SCGF is taken, is the one-particle eigenvalue h(θ ) of the corresponding total charge Q. Thus, the BFT predicts that only a limited amount of information about the twist field is required in order to evaluate the exponential asymptotic of its two-point function. Note that this is true also in the interacting case.

Entanglement and branch-point twist fields
In this section, we recall how entanglement entropies can be computed using a certain type of twist fields, called branch-point twist fields, associated to permutation symmetries. We then recall that, in free fermionic theories, these can be re-written in terms of U(1) twist fields. This will be used in the next section in order to apply the BFT to the calculation of entanglement entropies.

Replicas and branch-point twist fields
Within the replica method, in order to compute entanglement entropies (cf. Eqs. (1)-(2)) in a given theory, one re-writes the quantity trρ α A in terms of an appropriate expectation value in the replica model. This is the model composed of α independent, commuting copies of the original model (α ∈ ). For a one-dimensional system in a state with density matrix ρ, and with the subsystem A being a single interval, e.g., A = [x 1 , x 2 ], it is a simple matter to show [7,23] that trρ α A is exactly identified with the two-point function of branch-point twist fields, The expectation value on the r.h.s. is computed in the density matrix ρ ⊗α = ⊗ α i=1 ρ i , where ρ i is the original density matrix, on copy i. Branch-point twist fields in the replica theory are twist fields associated to the symmetry under replica cyclic permutations of order α (which generate the group Z α ). They take the product form (32), involving on-site copy-permutation operators 9 [23]: Here, denoting by a i (x) observables lying on (that is, acting nontrivially only on) copy i ∈ {1, 2, . . . , α} and position x, and identifying a α+1 (x) ≡ a 1 (x), the permutation unitary is defined by This implies the equal-time exchange relations (see (31)) From Eq. (42), Rényi entanglement entropies can be simply obtained via Eqs. (1)-(2). In the context of (1+1)-dimensional QFT, exchange relations of the form (45), (46) give the most appropriate formulation for working definitions of the branch-point twist field. It is in this context that they were first introduced [7], as a way of evaluating partition functions on branched surfaces, taking inspiration from [6].
We note that the action of branch-point twist fields can be diagonalized by going to the Fourier basis in the replica index (we choose anti-periodic boundary conditions in replica space, see below), and similarly forT α (x, t). In the next subsection we will use a similar construction, albeit in a different basis of the replica model. As will be explained in Sec. 4, for our purposes, the most general object we need to consider is the two-point correlation functions at different spacetime points.

Enhanced symmetry in free fermions: from Z α to U(α)
Consider the special case of free fermions, see Sec. 2.5. In this case, because of the quadratic nature of free fermion Hamiltonians, the Z α symmetry of the replicated theory turns out to be embedded into the larger symmetry group U(α), which accounts for not only permutation of replicas, but also rotations amongst them. Thus, the branch-point twist field is a twist field associated to a particular symmetry transformation, part of a continuous symmetry group. As explained in Sec. 2.4, using Noether's theorem, this then allows one to write an explicit extensive charge associated to the twist field [7]. The U(α) symmetry is most clearly expressed in a different basis of the replica theory than that used above, obtained by "fermionising" the replica theory. By the basic construction of the replica theory, different replicas commute with each other. However, in order to extract the symmetry U(α), one needs fermions in different copies to anti-commute. One simply defines the replica theory by asserting that fermion fields anti-commute. This is of course natural, but changes the action of the branch-point twist field (the exchange relations (45), (46)) by introducing an extra minus sign, as worked out in [7]. From now on, we denote

SciPost Physics
Submission the Dirac fermion on the i-the copy in this new basis; . We now recall the main arguments of [7] in order to obtain a useful form of the branch-point twist field. In the new basis, the U(α) symmetry is explicitly a linear action on the fermions, in its fundamental representation. Most importantly, in this basis, the cyclic permutation ψ i → ψ i+1 is a particular element of U(α), which is in fact an element of a U(1) subgroup. The Fourier transform (47) can be performed in this basis, This diagonalises that U(1) subgroup; the action of the twist field is then diagonalised. In fact, it turns out that the anti-commuting basis is also the one that guarantees that the Fourier transform operation keeps the S-matrix diagonal, see Appendix C. Thus both the charge associated to the twist field, and the S-matrix, are diagonal in terms of the particles corresponding to ψ p -this is at the root of the simplification. The fermion fields ψ p (x, t) after Fourier transform are still independent free fermions with canonical anti-commutation relations. Each Fourier sector admits an independent U(1) symmetry, and, as shown in [7], the branch-point twist field can be written as a product of U(1) twist fields acting nontrivially on each Fourier sector. Because of the extra minus sign in the twist field action, it is simpler to concentrate on the case of α even (the full dependence on α is obtained by analytic continuation), in which case the product goes over the following values of momenta: Specifically, it is found that [7] T α = p∈I α with τ α p being a U(1) twist field acting non-trivially only on ψ p (as a phase), (cf. (48)). The decomposition (53) allows us to factorise the branch-point twist field two-point functions into products of U(1) twist field two-point functions. This however only holds if the state can be likewise factorised. This is nontrivial: the state ρ ⊗α is naturally factorised in copy space, but not necessarily in the Fourier-copy space. It is a simple matter to verify that if ρ satisfies Wick theorem, then ρ ⊗α also factorises as a tensor product of states ρ in Fourier-copy space; this is because such states are completely determined by fermion two-point functions, which stay diagonal in Fouriercopy space. Therefore, we have, in Wick-theorem states ρ, Note how on the right-hand side, each factor is evaluated in the state ρ for the fermion ψ 2q−1 .

Submission
In the following we are going to apply the BFT machinery to each correlation function of the U(1) twist fields. The crucial fact that makes it simple is that, for any given p, τ α p (x, t) is the (bounded) twist field associated to the U(1)-charge with explicit expressions as exponential of half-space integrals of charge densities, as per Eq. (28): Q p acts on the single-particle basis as (54)). With Q = Q p , the twist field τ α p is identified with τ α p = T −i in the notation of (34) (that is, with η = 1), acting on the sector p. Recall that the action of the charge on the single-particle basis is all we need to know in order to apply the BFT (cf. (38) and (40)).

Entanglement entropies from BFT
We arrived to a rewriting of the two-point function of the branch-point twist fields as product of two-point functions of U(1) twist fields, Eq. (55). From there, using BFT, all such components can be accessed via Eq. (30) specified to the twist fields τ α p , which in the notation of Eq. (30) is identified with T −i , with the r.h.s. evaluated via Eq. (38), thus exploiting the free fermionic nature of the problem. Considering different choices of the points in spacetime where the global fields T α ,T α are located, we are able to access Rényi entropies both at equilibrium and after a quench, as we are now going to discuss.

Rényi entropies of a finite interval in a GGE (charge fluctuations in space)
We start by considering the α−Rényi entropy of a finite interval A = [0, x] within a generic GGE ρ w uniquely defined by the function w(θ ) (see Sec. 2.5). This means that we are interested in the following two-point function From the BFT perspective, this is obtained by focusing on the purely spatial direction, namely, we consider an "horizontal path" by setting γ = π/2 in (38) (and we take h(θ ) = h p ). Each two-point function of U(1) twist fields in (55) reads Then we consider the product in Eq. (55), which turns into a sum in the exponent, i.e.,

Submission
We may further evaluate those sums, by considering separately the part which depends and the part which does not depend on p (equivalently q, q , Eq. (53)). The latter is trivial and simply gives a contribution to For the remaining part, let us start by focusing on half of the sum, the terms from q = 1 to α/2, in (61).
where we used the Taylor expansion log(1 + x) = ∞ r=1 (−1) r+1 x r /r (which converges for w > 0), and we performed the sum over q. Next, we want to perform the sum in r in the r.h.s. of (63). To do that, we substitute the values of z and w first: where now we should consider separately three cases: 1. r = αm for integer m: in this case r is even (as α is even), and we have 2. r even but r = αm for any integer m: in this case each term of the sum (64) The sum of the terms for q = −α/2 + 1 to 0 in (61) give exactly the complex conjugate of this result. Thus we get Putting everything together, F α (−i) in (61) can be written as Finally, it is a matter of simple algebra to show that, in terms of the occupation function n(θ ) (37), we get

SciPost Physics
Submission where we defined The α−Rényi entropy is finally given by which coincides with the results obtained in [11,13] (there in the more general context of interacting integrable models).

Long-range correlations due to correlated particle pairs in homogeneous global quenches
We now review the main concepts underlying quantum quenches, restricting to "integrable" pairproducing initial states, and we explain how long-range correlations develop after such quenches. A quantum quench is an initial value problem for the many-body system where the initial state is the ground state of a different Hamiltonian than that used for the time evolution. Typically, one imagines a sudden change of parameter, for instance of the mass parameter. In integrable models, certain quenches are known to be of "integrable" type [26][27][28]. In these cases, the initial state can be written explicitly in terms of the scattering states (or Bethe ansatz states) of the post-quench, evolution Hamiltonian, as a so-called "squeezed state": for some (θ -dependent) factor K θ ,−θ , with N θ denoting a normalization constant, and |0〉 being the ground state of the post-quench Hamiltonian. The squeezed state is generically a finite-density state, where the energy (of the post-quench Hamiltonian) is extensive with the system size. See App. B.2 for a discussion of such integrable initial states in free fermion models. We will use later the fact there is a Bogolioubov transformation of the fermionic mode operators (a transformation between the post-quench and pre-quench fermions), such that the squeezed state satisfies (is defined by) After a long time in a quench problem, the state locally approaches a GGE. In integrable quenches, there is a well-known relation between the squeezed-state representation of the initial state, and the long-time GGE (see e.g. [44]). The statement of convergence to a GGE pertains only to local operators, or operators supported on finite intervals (that do not grow with time): The limit in (76) is expected to be valid everywhere in space. The relation between initial state and long-time GGE in free fermions can be worked out explicitly (see (150))

SciPost Physics
Submission Namely, we see that the map from squeezed states to GGEs is in fact one-to-one. Because of this one-to-one correspondence, it is clear that it is sufficient to know the longtime GGE in order to know the full behaviour of correlation functions in spacetime, as the GGE fixes the initial state uniquely (naturally under the condition that it be a pair-producing squeezed state). However, this relation can be relatively complex. Indeed, the statement of generalised thermalisation -that a GGE is reached -is true, generically, only on finite regions of space (see App. B.3). As is typically the case out of equilibrium, on large regions, say regions that grow linearly with the time after the quench, the state might not correctly be described by a GGE; and this may even be true for all times! Instead, the state may admit long-range spatial correlations, for instance correlations that have a large weight on distances that grow linearly with the time. These are not present in GGEs; recall that, as discussed above, GGEs typically have correlations that decay quickly enough in space (see App. B.1). Thus, even at long times, there may remain effects of the initial state that are not described by a GGE.
In the case of a squeezed state, such long-range correlations indeed exist, as we show in App. B.3 (and also in App. B.5 for "single-mode densities and currents", introduced below in Sec. 4.4). Their interpretation is that they are due to production of correlated pairs of oppositemomentum particles by the quench protocol. These particle pairs carry correlations to large distances as they separate. We evaluate explicitly these long-range correlations for conserved densities and currents in App. B.3 (and App. B.5). For instance, we find, for the charge density ballistic-scale correlations, in accordance with the picture according to which particle pairs are emitted at all velocities admitted by the dispersion relation.
In order to evaluate the Rényi entanglement entropy, as is clear from the calculation for GGEs in Sec. 4.1, we must evaluate the large-deviation theory for fluctuations of charges on large regions of space, and / or, as we will see below, fluctuations of current on large intervals of time. The BFT allows us to do that. However, as mentioned, the BFT requires no long-range correlations along the path in (12), as the flow equation (cf. (19)) assumes that the state along the path is a GGE. Long-range correlations may break the assumption that scaled cumulants are evaluated within a GGE -the effects of long-range correlations on the 2nd cumulants is illustrated in App. B.3. Below, we take them into account by choosing appropriately the path in order to avoid such correlations! Thus, the knowledge of where such correlations exist, and the knowledge of the long-time GGE, is sufficient.
The fact that there exists a path that avoids long-range correlations explains why, in pairproduction states, the full behaviour of the Rényi entropies can be written in a simple and universal way in terms of the long-time GGE. The same is not true when considering non-integrable quenches, namely quenches from more complicated states where groups of more than two correlated particles are emitted. Without the constraint of pair-production, the one-to-one correspondence between the state and the GGE is lost, and long-range correlations carry additional information not present in the GGE. Then, from our approach, we see that a universal description in terms of the long-time GGE is lost because such an "avoiding path" does not exist in general anymore.

Rényi entropies of half system after a quench (current fluctuations in time)
We now turn to the calculation of the α−Rényi entropy of a semi-infinite interval A = [0, ∞) after a global homogeneous quantum quench, at long times t → ∞. This is obtained from the

K I N G ' S C O L L E G E LO N D O N P A O L A R U G G I E R O
• : (for time-dependence) • Correlations and breaking of LDT : further info • BFT ("vertical" path) : • Rényi entropy of semi-infinite system after a quench: branch-point twist field one-point function (73), As expressed in (76), one-point functions of local observables converge to averages within GGEs. However, as we discussed, twist-fields are "semi-local" observables; from the point (0, t) emanates a branch cut, which is sensitive to the state where it passes. The branch cut can be taken on the horizontal half-line {(x, t) : x ∈ [0, ∞)} going from (0, t) to (∞, t), as done in the explicit construction of the field in Sec. 3. As explained in Sec. 4.2, and analysed in App. B.3, along this half-line, there exist long-range correlations due to coherent particle pairs emitted by the initial state. This prevents us from applying the BFT along this path (see Fig.1 (left)).
Using path independence of twist fields correlation functions, we can deform the path, between its initial and final points, in a way to avoid such correlations. Specifically, we choose the piecewise linear path joining the points (0, t) → (0, 0) → (∞, 0). This is shown in Fig. 1 (right). We note that as the final point is at spatial infinity, it can be displaced to time 0 -this in fact implements the correct physics of the entanglement entropy due to the single boundary at x = 0. Then, we may represent the one-point function as where the factors T α (0, t)T α (0, 0) represent the segment of path (0, t) → (0, 0), and the factor T α (0 + , 0), the segment (0, 0) → (∞, 0). This is valid as an asymptotic relation for large t, where the UV singularity due to the proximity of the fieldsT α (0, 0) and T α (0 + , 0) (which occurs because of renormalisation effects) is neglected. We simplify the expression (80) in two steps. First, we note that the segment of path (0, 0) → (∞, 0) does not provide any contribution to the result. This is because we may re-write the branch-point twist field T α (0 + , 0) as is done in SciPost Physics

Submission
Sec. 3, but in the basis of the before-quench canonical free fermions of the replica theory,ψ i (x, 0), Eq. (74). The exchange relations (45) (here at t = 0) hold for any field a i (x, 0), and in particular hold for a i (x, 0) =ψ i (x, 0). Therefore, by the same arguments, we obtain a decomposition as in (53), instead of (57). By (75), we haveψ i (x, 0)|Ψ j 〉 = 0 for all i, j, so it is clear that Hence Second, we note that along the path (0, t) → (0, 0), generic observables do not have longrange correlations coming from pair productions: correlations of generic observables approach those within the final GGE fast enough, in such a way that corrections due to the quench give only sublinear corrections to cumulants of time-integrated fermion bilinears (such as conserved densities and currents). This is because particle pairs always create correlations between points at separate spatial coordinates: it is not possible to create two co-propagating fermions, with the same momentum (here they would be with vanishing momentum, as the total momentum has to be zero). This fact is discussed in App. B.3; the discussion there is for a single copy, but it extends immediately to the α-copy state |Ψ α 〉. Therefore, rewriting the branch-point twist fields in terms of U(1) τ α p with branches in the time direction, using (53), and expanding in cumulants of U(1) currents, we see that on the segment (0, t) → (0, 0), for the purpose of the BFT, the state is correctly described by a GGE 10 .
An important consequence of these arguments is that we may evaluate the twist field one-point function after a quench, as an equal-space, different-time two-point function within the GGE representing the final state, Eq. (77): This is valid at long times, and omits small-time effects that occur before generalised thermalisation (which do not affect the asymptotic regime we look at). In order to apply BFT to the r.h.s. of (86) a last observation is needed. As the GGE is a Wicktheorem state, we can use (55), thus we are interested in the separate two-point functions of U(1) twist fields τ α p with branches in the time direction. It turns out that, as emphasised in Sec. 2.3, the currents j p (0, t ), t ∈ [0, t] in GGEs have time-correlations that decay fast enough so as to give only linearly growing cumulants: this is what allows the application of the BFT (see App. B.1 for

SciPost Physics
Submission a full discussion). We note that this is not true in general of other observables: in GGEs, generic fermion bilinears have cumulants that grow faster than linearly with time. But we are intersted in the currents only.
We are now in position to apply standard BFT. This amounts to repeating the same calculation as above in Sec. 4.1, but now in the purely temporal direction. We use the general formula (38) for the "vertical" path connecting initial and final point by choosing γ = 0, and similarly get (note that the path is in opposite direction as that of formula (38), and thus we must take h(θ ) = −h p ) where we used sgn(v(θ )) = sgn(θ ). Again, after performing the product all two-points functions of U(1) twist fields, we get Using H α (θ ) as defined in (71), the α−Rényi entropy reads This is the result obtained both from exact calculation in [15] and within the quasi-particle picture in [19,20].

Single-mode and pair-mode twist fields
We have discussed in Sec. 2.5 the conserved quantities Q θ = ψ † θ ψ θ , forming a "scattering" or continuous basis for the extensive conserved quantities of the free fermion model. In Sec. 3.2, we discussed the replica model with α copies, and the U(1) charges Q p , which are just the integration Q p = h p dθ Q θ ,p (with h p = πp α ) over all momenta θ of the continuous basis Q θ ,p = ψ † θ ,p ψ θ ,p in the Fourier-copy p. There, we have also discussed the twist fields τ α p associated to these charges, which turned out to be useful in the computation of the Rényi entanglement entropies in Subsections 4.1 and 4.3. A natural extension of these constructions is to the twist fields associated to each conserved quantity Q θ ,p . As we will see, these are indeed useful in evaluating the behaviour of Rényi entanglement entropies for intervals that grow linearly with time.
In order to simplify the notation, we consider a single copy of the fermion, and the scattering basis Q θ ; the discussion immediately adapts to the Fourier-copy p.
In the study of the ballistic behaviours of many-body systems, and in particular in the BFT, it is essential that the conserved charge Q considered be extensive -scale linearly with the volume (typically one requires 〈Q 2 〉 c ∝ L [37,45]). The charges Q θ are not extensive. However, as they form a continuous basis, integrals on small θ -intervals are extensive; thus it is better to define, for ε > 0 as small as desired, These act as

SciPost Physics
Submission hence have one-particle eigenvalues We show in App. B.4 that such Q θ are indeed extensive in GGEs, and we evaluate explicitly their associated densities and currents q θ (x, t) and j θ (x, t), From this, one can immediately construct the associated twist field and, for its correlation functions, apply the corresponding BFT based on the one-particle eigenvalue (92).
In fact, we are interested in studying the squeezed state (73). It is clear that this state factorises into momentum intervals as follows: where and we write the ground state in a naturally factorised way as |0〉 = θ ∈( + 1 2 )ε |0 |θ | 〉. Likewise, we will consider the pair-mode charges Q |θ | = Q θ + Q −θ and the associated densities Both act trivially (as zero) on |Ψ |θ | 〉 if θ = θ (θ , θ ∈ ( + 1 2 )ε). From these, we get the pair-mode twist fields which acts trivially (as the identity) on |Ψ |θ | 〉 if θ = θ . These are still U(1) twist fields, for the sub-U(1) symmetry acting on the tensor factor of modes within [θ − ε/2, θ + ε/2]. Note in particular that the global U(1) twist field τ(x, t) associated to the total charge Q = dθ ψ † θ ψ θ = d x ψ † (x)ψ(x) can be factorised as and that, by factorisation of its action on the state, we have Clearly, as the pair-mode twist fields act trivially on other tensor factors in the state, we may also write, more simply, SciPost Physics

Rényi entropies of an interval after a quench (fluctuations of single-mode densities and currents)
We finally extend the result for the entanglement growth after a quench to a finite, but ballistically growing interval A = [0, x], with x = ξt. To this aim, we should consider the following two-point correlation function in a squeezed state Eq. (73) (or Eq. (79) in the replicated theory): The idea is the same as that used in Sec. 4.3, that we need to deform the integration path in such a way that, everywhere along the path, all points remain uncorrelated (on large scales), thus enabling us to apply BFT. The choice of the path will now depend on the values of ξ, and in fact, we will need re-write the two-point function as a product of two-point functions of pair-mode twists fields, and to choose different paths for each such two-point function.
It will simplify the discussion to already re-write the two-point function in terms of U(1) twist field. As the squeezed state is a Wick-theorem state, we can directly use (55): where |Ψ p 〉 is the squeezed state |Ψ〉 for the fermions ψ p (x), ψ † p (x) on Fourier-copy space p. We start by considering the two asymptotic regimes: • At short times (more precisely in the limit ξ → ∞ of the scaled, long-time asymptotic behaviour), entangled particle pairs coming out from the initial state will correlate points within the original integration path. To apply BFT, then, we deform the straight path (0, t) → (x, t) to the piece-wise straight path (0, t) → (0, 0) → (x, 0) → (x, t), made of three segments (see Fig. 2 (left)). By the same arguments as in Sec. 4.3, the space-like segment will not contribute to the entanglement growth, and the time-like segments will give separated, identical contributions given by the long-time GGE. We are thus left with the contribution of the two, separate time-like segments. The fact that the segments do not correlate with each other is thanks to the assumption that the GGE satisfies n(θ ) → 0 as |θ | → ∞ (that is, the density of pairs produced at large momenta tends to zero), as is discussed in App. B.3.
• At long enough times (either the limit ξ → 0 of the scaled, long-time asymptotic behaviour, or the long-time limit followed by the long-distance scaling), the particles generated from the initial state do not correlate points within the path (0, t) → (x, t): cumulants scaled by the distance x do not receive contributions from such particle pairs. Corrections terms to the GGE values of cumulants can only come from pairs of particles at infinitesimally small momenta, and, it turns out, such corrections become zero when the total number of correlated pairs on the interval [0, x] tend to zero. As there is at most a finite density of pairs produced per unit momenta, there remain no pairs on infinitesimally small momentum intervals 11 . Thus the asymptotic behaviour is that obtained from the long-time GGE. This is discussed in App. B.3. Particle pairs of finite momenta would, of course, correlate points between the path segments (0, t) → (0, 0) and (x, 0) → (x, t) (also discussed in App. B.3),

K I N G ' S C O L L E G E LO N D O N P A O L A R U G G I E R O
• : • Asymptotic regimes : "small" and "large" time ( ) • Rényi entropies of finite interval after a quench : thus we must avoid the piece-wise straight path. Therefore the correct way to use BFT is by using the original path (see Fig. 2

(right).
So we arrive to the following asymptotic results for x, t → ∞: t)). This leads to Namely, at short times (but much larger than microscopic times), the growth is described by the path in the purely temporal direction (89), and at long times, the system goes, uniformly as a function of the velocity x/t → 0, to the equilibrium GGE and there the result is the one of the purely spatial path (72). These are, however, only asymptotic results in ξ, within the scaled regime x, t ∝ . It turns out that we can access all values of ξ = x/t within this regime, by using similar arguments, but now for the single-mode twist fields introduced in Sec. 4.4 (in fact, we need the pair-mode twist fields). Effectively, using these, we will be able to take into account that the meaning of "short" and "long" time depend directly on the speed of the travelling particles v(θ ) = E (θ ).
We start with the decomposition of global U(1) twist fields into pair-mode twist fields (99), which we write in the replica model for each Fourier-copy p, and with single-particle charge eigen-

SciPost Physics
Submission value h p = πp α instead of 1 in (99) (as done in (57)) and q |θ |,p (x, t) has the form (173) times h p . Thus, by factorisation of two-point functions (101), we re-write (103) as Having made this re-writing, the analysis now follows that of the ξ → ∞ and ξ → 0 limits made above: there is an exact parallel for each individual two-point function with the only difference that it is not necessary to take the asymptotic limit in ξ. For each θ (and each p), the factor |Ψ |θ |,p 〉 of the full state |Ψ p 〉, on which τ α |θ |,q act non-trivially, correlates points (recall that θ ∈ ( + 1 2 )ε). The analysis of single-mode correlations is made in App. B.4. Therefore, for ξ > 2v(θ + ε/2) correlations occur on the horizontal path (0, t) → (x, t), but no correlations occur on (0, t) → (0, 0) → (x, 0) → (x, t) (note that, again, the segment of path (0, 0) → (x, 0) does not contribute). Thus we must choose the latter path ( Fig. 2 (left)). On the contrary, ξ < 2v(θ − ε/2), correlation occur between the segment of paths (0, t) → (0, 0) and (x, 0) → (x, t), but not on the horizontal path (0, t) → (x, t). Thus we must choose the latter (Fig. 2 (right)). In making these right choices, the correlation functions of pair-mode twist fields tend to their values in the long-time GGE, Note how in the first line, it is a four-point function that appears. Now again, in order to apply the BFT, we need to consider the correlations between twist fields within GGE. We already argued in Sec. 4.3 (with supporting calculations in App. B.3) that no strong correlation occurs between local operators at equal times and different points in space, thus on the second line of (109) we may apply the BFT. We also argued that no strong correlation occurs between current operators at equal space and different times, and in fact this also holds for single mode currents. However, in order to simplify the first line of (109), we need to address correlations between currents on the path segments (0, t) → (0, 0) and (x, 0) → (x, t). In general, local observables have strong correlations at different space-time points, due to hydrodynamic

SciPost Physics
Submission modes propagating in space-time. As we do not make any strong assumption about the dispersion relation, all hydrodynamic velocities occur, hence correlations occur between generic local observables on these separate path segments. However, single-mode currents only produce hydrodynamic modes at the corresponding velocities; supporting calculations are found in App. B.5. Thus, as long as ξ > v(θ + ε/2), no correlation occurs between these paths for the single mode currents j ±θ ,p . As v(θ ) > 0, then ξ > 2v(θ + ε/2) ⇒ ξ > v(θ + ε/2), hence on the first line simplifies and we have where the BFT can be applied for all two-point functions.
For the two-point functions with spatial separation 〈τ α |θ |,p (0, 0)τ α |θ |,p (x, 0)〉 ρ w , we can use the analysis of (60) made in Sec. 4.1, where we only have to replace, on the right-hand side of (60) inside the θ -integral, the constant one-particle eigenvalue h p by the piece-wise constant function Thus the same analysis goes through, but with the integral restricted to Likewise for the two-point functions with temporal separation 〈τ α |θ |,p (0, t)τ α |θ |,p (0, 0)〉 ρ w , with the analysis of Sec. 4.3. Putting the results together, we obtain For the case of x/t within this excluded region, we do not have explicit results, but the scaled cumulants still are finite (as one can see by doing a calculation similar to App. B.3, for instance). Thus, the result may be deemed valid as well within this excluded region, up to an error of order ε.
Taking the product over θ 's as per (108), and as this holds for all ε > 0, we can take the limit ε → 0 and we obtain This is in full agreement with the quasiparticle picture [19,20].

Submission
Finally, we note that the relation (3) between this formula for Rényi entanglement entropy growth, and the static and dynamic fluctuations, is directly obtained from the above discussion, by identifying and using the explicit expressions of pair-mode densities and currents (176) in App. B.4, and taking the limit ε → 0.

Discussion and conclusion
In this paper we have studied the Rényi entanglement entropy in GGEs and after quenches from integrable (pair-production) states in free fermion theories. Although this has been relatively well studied in the literature, most results were based on specific ways of writing the Rényi entanglement entropy using the free fermion structure (e.g. in terms of determinants), and on the idea of entanglement due to engangled pairs produced by the quench [19,20]. A first-principle derivation that generalises beyond free fermions was still largely missing, while it is known that the simple quasi-particle picture fails for α-Rényi entanglement entropies (with α = 1) in interacting models [22]. We have proposed a new approach based on twist-field correlation functions and hydrodynamic fluctuations. This uses the hydrodynamic theory for free fermions, which is a special case of generalised hydrodynamics (GHD), and the ballistic fluctuation theory (BFT), which relates the exponential decay of twist-field correlation functions to hydrodynamic large-deviation theory. Crucially, in order to have a full understanding of the quench dynamics, we have introduced a new concept: that of single-mode twist fields. These are twist fields associated to the quasi-local charge counting the number of fermions within a small interval of momentum; or more generally twist fields "acting" on the quasi-locality sector of observables supported on a small momentum interval. The approach is potentially more general and more fundamental, as hydrodynamics, the BFT and single-mode twist fields -twist fields associated to individual hydrodynamic modes -are applicable much beyond free fermions. Perhaps most interestingly, it reveals the new physics of thermodynamic and hydrodynamic fluctuations behind the behaviour of the Rényi entanglement entropies.
Three important concepts are brought forward: • Entanglement is deeply connected to fluctuations, and the large-scale behaviour of entanglement, both static and dynamic, is controlled by large-deviation and hydrodynamic principles.
• Hydrodynamic modes and projections onto such modes are more accurate and general notions which replace the idea of particle-pair productions used to understand entanglement dynamics in integrable systems.

Submission
• The fact that the entanglement growth in quenches from "integrable" states can be written as a simple and universal function of the generalised Gibbs ensemble (GGE) reached at long times comes from the particularly simple structure of the long-range correlations that such states present.
To elaborate on these concepts: first we have confirmed that the Rényi entanglement entropy in GGEs is controlled by thermodynamic fluctuations, and related to (a simple analytic continuation of) a difference of thermodynamic free energies. This is in agreement with the observations, made earlier [29][30][31], that the large-deviation theory for charge fluctuations is closely related to the entanglement entropy. Here, this relation appears naturally from completely general concepts: branch-point twist fields and the BFT. Using these, in fact, the conclusion is pushed further: we show that the growth of Rényi entanglement entropy after a quench is controlled by hydrodynamic current fluctuations, and related to a dynamical free energy associated to the large-deviation theory for charge transport, as fully encoded in Eq. (3) in the introduction (we mention that some qualitative arguments in this direction were already present in [46]). The relation between Rényi entanglement entropy and charge fluctuations is a general aspect of quadratic theories (not only free fermions, as our approach could also be generalized to free bosons), as in such theories, the branch-point twist field can be written as a product of U(1) twist fields, which are then associated, by the BFT, to the large-deviation theory of U(1) charge fluctuations.
Second, the methods we have developed show that the notion of quasi-particle used in integrable systems to explain the behaviour of entanglement, should in fact be replaced by that of hydrodynamic mode. Indeed, the BFT, which describes the exponential behaviour of twist field correlation functions, is purely based on the Euler hydrodynamic data of the microscopic model. In free fermion models, and in integrable models, it turns out that hydrodynamic modes are in one-to-one correspondence with quasi-particles (see e.g. the review [47]); and in particular, the single-mode twist fields we have introduced, are twist fields associated with such hydrodynamic modes. But beyond these situations, hydrodynamic modes are the more general objects at play in the large-scale dynamics of many-body systems.
Third, our calculations explain why it is possible to specify the growth and saturation of entanglement after quenches from pair-production states in a simple way in terms of the long-time GGE. This is not based on the conventional physical picture of entanglement produced by entangled pairs of particles. But rather, it is based on the study of long-range correlations that such pairs give rise to after quenches. It has not been appreciated until now that quenches give rise to long-range spatial correlations , of the type found recently in non-equilibrium, long-wavelength states [35,38]. These long-range correlations are generically seen by observables supported on regions of space that are large enough; specifically, with a ballistic scaling of the region's length x with respect to the time t since the quench. Thus, for such observables, the state is not a GGE. This is important, as twist fields are semi-local with respect to the fermions, thus the semi-locality branch is affected by such correlations.
The fact that the long-time GGE can be used to describe not only the saturation of the Rényi entanglement entropies but also their growth in a simple way, is because of the particular structure of long-range correlations in integrable pair-production quenches. Indeed, in such quenches, for every choice of x/t, one can always choose a path in space-time which avoids all long-range correlations. This follows from a simple geometric analysis of trajectories in space-time. Intuitively, long-range correlations occur between positions in space-time where pairs of correlated particles lie, a picture that we fully support by simple calculations of correlation functions of conserved densities and currents and their asymptotic behaviours. Once the path avoids long-range correlation,

SciPost Physics
Submission it only perceives the long-time GGE.
Note that it is obvious that the entanglement growth can be described purely in terms of the final GGE (with no further information from the initial state needed), because of the one-to-one correspondence between initial squeezed state (pair-production state) and GGEs, see Eq. (77) (so no more information about the initial state is present at all). The structure of long-range correlations however allow a simple description, in terms of fluctuations within the long-time GGE. In general, the one-to-one correspondence is lost in quenches from more complicated states, and the universality of the entanglement growth is also lost [32,33].
Importantly, we find that the branch-point twist field can be decomposed into tensor factorsthe single-mode twist fields -that act on each small momentum interval. Each small momentum interval is associated with a family of quasi-local operators, with respect to which branch-point twist fields can be defined 12 . Each such twist field is only semi-local with respect to the fermions pertaining to the single momentum interval, and, pairing intervals of opposite momenta, allowed us to separate the two-point function of branch-point twist field (used to evaluate the Rényi entanglement entropy) into a product of contributions on different momentum intervals. For each factor, the semi-locality branch can be chosen in order to avoid long-range correlations corresponding to the pair it perceives. We provided extensive calculations of correlation functions of single-mode densities and currents that support this physical picture.
From these considerations we fully reproduced the long space-time dynamics of the Rényi entanglement entropy that had been obtained by pair-entanglement argument.
Many extensions of our work are possible. Most importantly, we believe the derivation we have provided, and many of the conclusions, can be extended to interacting integrable models, and potentially to interacting non-integrable models.
In interacting integrable models, it would be interesting to reproduce, and provide a better understanding of, the recent result [22]. This was obtained using crossing symmetry of relativistic quantum field theory, in order to relate Rényi entanglement entropy growth in time to the linear scaling of Rényi entanglement entropy in space, much like one can evaluate currents by crossing from conserved densities [43]. However, as we have argued, the more general understanding of time-extensive behaviours is via hydrodynamic modes. Thus, it is likely that hydrodynamic ideas will provide a more first-principle derivation.
Technically, in interacting models, it is not possible to factorise the branch-point twist fields, in such a simple way as in [7], into U(1) twist fields, a trick that we have used here. We believe this difficulty can be surmounted as follows. First, it is still possible to diagonalise the twist field action, at the price of making the resulting S-matrix non-diagonal, see App. C. The branch-point twist field is then associated with a symmetry that has diagonal action on the new asymptotic particles, and hydrodynamic modes can be constructed from these particles (by constructing the corresponding nested thermodynamic Bethe ansatz). Having twist fields associated to charges that are diagonal in the particle basis, the results of the BFT for generic intergable models can in principle be applied.
It is also immediate that single-mode twist fields exist as well in interacting integrable models, which would be interesting to study, independently form their applications to entanglement.
Another avenue is to use the ballistic macroscopic fluctuation (BMFT) theory developed recently [35]. This is a more general construction which does not involve a flow on GGEs (by constrast to the BFT). This would allow us to apply the principles introduced here -the relation

SciPost Physics
Submission between hydrodynamic fluctuations and entanglement entropy using twist fields -beyond homogeneous quantum quenches, and beyond the simple particle-pair quenches, as the BMFT is applicable to inhomogeneous situations and for generic long-rance correlation structures. In particular, it would be interesting to account for the long-range correlations not by choosing paths that avoid them, but by evaluating directly their influence on the fluctuations. This would be important, as initial states that are inhomogeneous generically produce long-range correlations [35,38] that are not necessarily of particle-pair type. This could also access quenches where multiple-particle processes are involved [32,33].
Beyond integrability, perhaps the main results are those based on a notion of "surface tension" underlying entanglement entropy in chaotic systems [48]. It would be interesting to apply our methods, based on hydrodynamics, to such situations.
A natural further extension is to introduce the effects of diffusion, in particular using the exact results in integrable models [49], and potentially the effects of dispersion [50].
Finally, there is no reason to restrict ourselves to entanglement entropy, or to quantum systems. The methods we have introduced here are immediately applicable, for instance, to the largedeviation theory of U(1) densities and currents after quenches in interacting models, both at and away from integrability, and both for quantum and classical systems. This would be a very interesting direction to investigate.

A Remarks on notions of locality and twist fields
There is a lot more that one can say about the twist fields we introduced in Sec. 2.4, as well as the notions of locality we briefly discussed. Here we collect a number of remarks in order to provide a brief and rough guide to this wide subject.

A.1 Unbounded observables and topological charges
As is made clear from the height field definition of twist fields, a twist field may exist as soon as there is an observable, say ϕ(x, t), that is "unbounded": that takes values in a non-compact space which are not bounded by the dynamics. When this happens, the field can grow from a value it has at x, t to another well separated value at x , t . If this growth is linear, then from this we can construct a twist field as done in Sec. 2.4, with a large-deviation theory described by the BFT. This usually happens when there is a non-compact symmetry 13 , such as the symmetry group of the sine-Gordon model if the sine-Gordon field is taken in , or the symmetry group

SciPost Physics
Submission of the real free massless Boson. In such cases, ϕ(∞) − ϕ(−∞) is a "topological charge", it is an extensive conserved quantity with density q(x, t) = ∂ x ϕ(x, t). The vertex operators e −iηϕ(x,t) are the associated twist fields, and the extensive charge ϕ(∞) − ϕ(−∞) should be considered as part of the space of extensive charges Q i used to construct GGEs. It will appear after appropriate quenches via (generalised) thermalisation.

A.2 Descendant twist fields and semilocality sectors
The exchange relations (31) are a good way of characterising twist fields. However, they characterise not a single field, but a family of fields. Indeed, clearly, the identification (34) is not the unique choice satisfying (31). For instance, will also work, for any local observable a(x, t). The choice (34) may be seen as a "highest-weight" twist field, and the above are usually referred to as "descendant twist fields" (these notions make full sense, for instance, in quantum or conformal field theory, using the concept of dimension). All such descendants are in the same "semilocality sector" T defined by the exchange relation (31). One application of the BFT to descendant twist fields is explained in [41] in the context of the XX quantum chain.

A.3 Non-abelian semilocality
Given two "local enough" symmetry transformations σ and σ , that is a(x, t) → σa(x, t) and a(x, t) → σ a(x, t), we can define two semilocality sectors T σ and T σ . In general, if the transformations do not commute, one has that if T ∈ T σ , then σT ∈ T σ•σ •σ −1 . Further, if T ∈ T σ and T ∈ T σ , then the exchange relation takes the form

A.4 Twist fields in the literature
It is difficult to give a full account of the literature on twist fields. As a guidance we mention that twist fields and their semilocality have been discussed extensively in various contexts, including: phase transitions in classical and quantum statistical models [51][52][53] (see the review [54]); vertex operators, Yangians, parafermions and orbifolds in conformal and integrable quantum field theory [55][56][57][58][59]; tau-functions and Painlevé equations [60][61][62][63][64][65][66][67]; and entanglement entropy in quantum field theory and in quantum spin chains [7,23,68,69]. Twist fields have also been considered in higher dimensions [70]. In most works, the focus is on ultra-local "internal" symmetries, that strictly factorise in space, usually part of a symmetry group such as n , U(1), SU(n), permutations, etc. Note that for ultra-local symmetries, large inequalities , can be replaced by ordinary inequalities <, > in (31) for finite distances between the supports of the observables involved. This is the usual way of writing the exchange relations. More recently, twist fields associated to spacetime boost transformations in QFT, that are not ultra-local, have been considered [71]; this is an example where the hamiltonian density is not preserved by the transformation,h(x, t) = h(x, t).

A.5 Concepts of locality in the literature
The concept of "locality" has been discussed widely in the literature, under a variety of definitions. In relativistic QFT, local fields are those that commute at space-like distances with the energymomentum tensor. It is important to remark that under this general definition, local fields include twist fields associated to internal symmetries. This definition can in fact be used in any quantum model, be it a field theory, spin chain or model of interacting particles. In fact, one defines "locality sectors" containing families of local fields that commute with each other at space-like distances; and a distinguished locality sector is that containing the energy density. These considerations are at the basis of orbifold conformal field theory [55].
In spin chains, the most naïve concept is that of operators supported on finitely many sites. In the C * algebra formulation, this is completed with respect to the operator norm [36]; importantly, the property of operators commuting in the limit of large separations is preserved by this completion. Part of these C * algebra elements are the "quasi-local operators" that have been introduced in order to describe generalised thermalisation in integrable models [45]. These are elements of the C * algebra for which one can still define a finite support, but only up to corrections of exponentially decaying norm.
But much like in QFT, twist fields, which are semilocal with respect to generic observables but may be local with respect to some family of observables including the energy density, can also be adjoined to the C * algebra, as is natural to do for instance in the context of Jordan-Wigner transformations [72]. Then, either from the C * algebra, or from some potentially smaller algebra of observables deemed local (for instance with appropriate decay of correlation functions, and which, again, may include twist fields), other completions are possible, and sometimes more physically relevant. For instance, the Gelfand-Naimark-Segal Hilbert space with respect to a given state, and its space of bounded operators, both are usually larger than the C * algebra. Another Hilbert space completion is that based on susceptibilities [37]. This gives the concept of "extensive" quantities, a generalisation of quantities written as sums over space of local operators (or "local densities"). These form a Hilbert space, a priori without any algebraic structure. In particular, it has been shown [37] that extensive conserved quantities are in bijection with "pseudolocal charges", roughly defined by their extensivity property 〈Q 2 〉 c ∝ L in a system of length L [45]. Extensive conserved charges form a complete set which has been rigorously shown to fully describe the linearised Euler hydrodynamics [40], and they may be used to more formally and precisely construct GGEs (addressing convergence issues) [37].
Despite all these studies, the full relation between extensive conserved quantities, twist fields, local symmetry transformations and GGEs is still not fully unexplored.

B Correlations after a quench from an initial state with pair structure
In this appendix we provide supporting argument for the choice of the integration path for evaluating the SCGF made in the main text. We recall that the choice of path is dictated by two main ideas: • The production of pairs of quasi-particles by the initial state by the after-quench dynamics gives rise to long-range correlations: points reached by a pair of quasi-particles with opposite momenta are correlated. These quasi-particles have the interpretation of fluid modes, and such long-range hydrodynamic correlations are akin to those seen in the Ballistic Macroscopic Fluctuation Theory (BMFT) [35,38].

Submission
• BFT [24] is applicable only when multi-point correlation functions of the densities and currents, the integrand in (12), cluster fast enough along the chosen contour. Otherwise, depending on the structure of such correlations, the SCGF and the cumulants may be divergent under ballistic scaling, or the BFT result may receive additional contributions from these correlations which have to be taken into account. The more general BMFT [35] in principle provides the corrections. However, it is simpler to directly apply the BFT by choosing contours that avoid such correlations, leading to the paths shown in Fig. 2.
Below, we explicitly show the presence/absence of such long-range correlations along the different paths considered in the main text.
We use the notations a(x, t) and b(x, t) for free fermionic fields with the usual normalisation, e.g.
Note that in the main text a variety of canonical free fermion fields were defined: the original fields ψ(x, t), the replicated ones parametrised by a copy index ψ i (x, t), and the fields obtained from these by diagonalising in copy space, ψ p (x, t). These all are canonical free fermion fields (independent from each other for different copy number i, or for different diagonalised copy number p). The calculations below therefore apply to any such choice of fields.
As we discuss quenches, in this appendix we use two consecutive letters for the canonical free fermion fields: a(x, t) (for the pre-quench fermion) and b(x, t) (for the post-quench fermion).

B.1 Global U(1) densities and currents and decay of correlations in GGEs
First, recall the definition of the generalised current as a line integral (12) where j(x, t) and q(x, t) are the current and the density associated the U(1) conserved charge Q in (118). We recall also that the above integral only depends upon the end points of the path due to the conservation law relating current and density, ∂ t q(x, t) + ∂ x j(x, t) = 0. In free models with global U(1) symmetry the fermionic Hamiltonian is of the form where E(θ ) is the dispersion relation (recall that E(θ ) = E(−θ ) is a strictly convex function). The local charge density is given by and it can be easily verified by using the fermionic algebra that Q = d x q(x, t) is conserved, [Q, H] = 0. Let us find the associated local current.

Submission
By the conservation law we have (we suppress the time dependence as all fields are the same time t) Integrating with respect to x we find where the x-independent integration constant is chosen in such a way that the result is a local observable 14 . This is the known expression for the current in the case of a quadratic dispersion relation in the continuum. Actually, restricting the integration over momenta in [−π, π] and taking E(k) = cos(k) one reproduces also the current on the lattice; but here we keep θ , k ∈ for simplicity. We now show that in a GGE, the connected correlation functions of densities decay fast enough in space, and the correlation functions of currents decay fast enough in time, in such a way that scaled cumulants are finite, thus making the BFT applicable. The former in fact is valid for all local observables, while the latter only hold for the currents.
For simplicity, here and in the following subsections, we will concentrate solely on two-point correlation functions -although all higher-point functions (and their respective cumulants) should in principle be investigated similarly.
Let 〈·〉 be a GGE. Let us assume that the occupation function n(θ ) characterising the GGE is analytic in a neighbourhod of . Using we have, on the one hand, For x > 0 (resp. x < 0), contour deformation can be performed as θ → θ −iγ (resp. θ → θ +iγ) for γ > 0 small enough, and we see that the resulting integral decays exponentially as |x| → ∞. This implies exponential decay of all two-point connected correlation functions of local observables formed out of sums of products of b(x), b † (x) and their derivatives, including U(1) densities. It also implies linear scaling of cumulants; for instance this would mean where in the last line we have shifted θ → θ − i sign(x − x )γ and used sign(x)x = |x|. This is the correct ballistic growth of the cumulant.
On the other hand, we find This has a stationary phase at θ * such that E (θ * ) = 0; this point is unique by our assumption of strict convexity (and θ * = 0 by symmetry, although we don't make use of this fact in this calculation), so a saddle point analysis gives Therefore, correlation functions of generic local observables o(x, t), o (x, t) formed out of bilinears of creation and annihilation operators have algebraic decay For such decay, cumulants of total time integrals do not grow linearly, thus breaking the large-deviation principle at the basis of the BFT. However, an important remark is that this generic behaviour of fermion bilinears does not hold in the case of currents, t). Indeed, using (123) with x = 0, we see that we must set θ = k = θ * for the long-time limit of the current two-point function. From we realise that . Therefore, the current two-point function decays faster than 1/t; in fact it decays as This guarantees the correct scaling of cumulants

SciPost Physics
Submission and thus the validity of the BFT. A similar argument shows that the current perpendicular to the path -the integrand in (119) -has a similar property along the path, thus guaranteeing that the clustering requirement (21) holds.

B.2 Quench protocol and initial state
In order to describe the quench protocol considered in the main text (for which we obtain predictions on the dynamics of the entanglement entropy) we define the pre-quench and the post-quench fermionic Hamiltonians respectively as where again the fermions satisfy According to the protocol, the system is initialized in the ground state of H 0 and then let evolve with H. This corresponds to changing the whole dispersion relation (not only a parameter as in typical quenches in the literature), see e.g. [44,73].
The two set of fermions are related by a Bogolioubov-type transformation in the following way Imposing the validity of anticommutation relations one gets the following constraints on the func- Note that the first of these is identically satisfied choosing f θ = f −θ and g θ = −g −θ or viceversa. In our analysis we will keep these functions general. The initial state |Ψ〉 is defined as which, for instance, could be a filled Fermi sea so that operators a θ are to be interpreted as creating excitations on top of these. It can be easily shown that in terms of post-quenches quantities this is described by the following squeezed state where |0〉 is the ground state of the post-quench Hamiltonian satisfying b θ |0〉 = 0 (this state has the nice property of being gaussian so that Wick's theorem applies). The function K θ ,θ = −K θ ,θ can be related directly to the functions f θ and g θ appearing in (136) using the fact that |Ψ〉 is annihilated by a θ . In terms of post-quench operators, this condition reads which used in (140) in combination with b θ |0〉 = 0 gives so that the condition is Note in particular that we must have, by the anti-symmetry K θ ,−θ = −K −θ ,θ , Finally, we may evaluate the predicted long-time GGE for the quench simply by evaluating the post-quench conserved quantities in the pre-quench vacuum state |Ψ〉. Inverting (136), we write For this computation, it will be convenient to look at the values of the extensive conserved quantities; hence we take a finite system of length L. This is warranted, as the quench is homogeneous. The above description of the quench stays valid with the discretisation θ ∈ 2π/L, and with the usual canonical anti-commutation relations with regularisation δ(θ − θ ) → δ θ ,θ L 2π . We consider b † θ b θ , and obtain, after some algebra using Eqs. (137), Using 〈Ψ|a † θ a θ |Ψ〉 = 0 and 〈Ψ|a −θ a † −θ |Ψ〉 = L 2π , we get But also, in a GGE with density matrix ρ w , see Sec. 2.5, we have 〈b † θ b θ 〉 = L 2π n(θ ), thus we identify Again using Eqs. (137), we obtain and therefore, from (37) and (148),

Submission
Note that in our analysis of GGEs, we assume that n(θ ), thus |g(θ )| 2 , has an analytic extension to a neighbourhood of . This analyticity property is not true of the function w(θ ), which must have a singularity (e.g. logarithmic) at θ = 0 because of Eq. (144). We also assume that n(θ ) → 0 as |θ | → 0, and thus this must also be true for g(θ ).
Below we report for completeness all the relevant elementary correlation functions of fermionic operators after the quench, and their symmetries. Those will be used in the following subsections, where evaluating current and density correlations, which are bilinears in the fermions (so application of Wick theorem requires only the knowledge of those). In real space we define and similarly for their hermitian conjugates. Going to momentum space, these take the form where in particular |g(θ )| 2 = n(θ ) = 1 − | f (θ )| 2 . Note the following symmetries so that at equal times and also

B.3 Approach to the GGE
In the previous Sec. we have evaluated the GGE ρ w corresponding to the initial state |Ψ〉 simply by evaluating the averages of the mode occupation. Here we analyse a bit more in detail how the GGE is approached in time. We first note that Thus, by Wick's theorem, the only difference between averages in |Ψ〉 and in 〈·〉 ρ w come from the contraction and its complex conjugate. Thus we evaluate 〈Ψ| b(x, t)b(x , t ) |Ψ〉 in three main situations that are important for our analysis: t = t , x = x (for the cumulants of space-integrated conserved

B.4 Single-mode density and currents and decay of correlations in GGEs
In the main text, when studying the dynamics of the entanglement of an interval after a quench (Sec. 4.5), we are interested in the single-mode twist fields. This requires constructing densities and currents not only for the global U(1) charge as done above, but also for the individual conserved quantities b † θ b θ . These conserved quantities are not extensive -they are not integrals of local or quasi-local observables -however, as reviewed in [47] in the more general context of integrable models, they form a scattering basis for such extensive quantities. Thus, integrations over small θ -intervals give extensive conserved quantities. These are the single-mode conserved quantities that we now investigate.
As we are working directly in the thermodynamic limit and, in the continuum of space, the momenta fill the real axis [−∞, ∞]. Let us write this as a union of disjoint intervals centered at equispaced "target" momenta: ∪ ∞ i=−∞ A θ i where A θ = [θ − ε/2, θ + ε/2) and θ i = (i + 1/2)ε. We can write the total charge as

Submission
Clearly each "regularised" (by ε) single-mode charge Q θ is conserved, [Q θ , H] = 0. It is also extensive: in a GGE in a finite volume L, we have 〈Q 2 θ 〉 c ∝ L: θ +ε/2 θ −ε/2 dθ n(θ )(1 − n(θ )) . (170) As mentioned, if we want to write a density in real space for each b † θ b θ , we will get something non local. However, Q θ 's have quasi-local densities. We seek a function f θ (x, y) such that Going to Fourier space, one can show that (see (117)) The corresponding regularised single-mode density, parametrised by the momentum, and one choice of the density (the only hermitian and PT symmetric one), is given by In terms of Fourier modes, this takes the form As [Q θ , H] = 0, the density q θ (x, t) has an associated current satisfying a continuity equation and by a calculation analogous to (123) one finds where ϑ(x) is the Heaviside theta function. This is basically the same as (123) with a restriction on the q integration around the target mode θ . For convenience, in fact we will consider the momentum-pair densities and currents associated to a pair of opposite momenta; these are the ones used in the twist field decomposition (101). We now analyse the behavior of two-point functions on GGEs following the same lines of B.1. This is because, again, the original BFT was defined for expectation values on GGEs and, eventually, we would like to replace with those the expectations on the initial state before the quench. Later, we will study more carefully the approach to the GGE. In a GGE (characterized by Boltzmann factor w(k)), the only surviving elementary correlator is 〈b † k (t)b k (s)〉 = e iE(k)(t−s) δ(k − k )n(k) with n(k) = (1 + e w(k) ) −1 .
Let us take the case η = 1 : the two step functions square to one and the condition for it to vanish is k * (ζ) ≥ ε/2 + θ or k * (ζ) ≤ −ε/2 + θ . Using ζ = ∂ k E(k * (ζ)) = v(k * (ζ)) and the monotonicity of the velocity, we can invert the above relation, leading to a vanishing step functions whenever For η = −1 we obtain the same condition because we are considering ζ ≥ 0. This means that, under the condition (182), the leading term of correlation (181) (coming from the saddle point) vanishes.
Correlations between the two vertical segments, on the contrary, arise only in a tiny cone (of order ε) around the velocity v(θ ). When those are present, the second cumulant resulting from (181) grows as O(T log T ) at large time T (to be contrasted with (133)).
Note that if we consider the special case t = s ( ζ = +∞), then (181) has no saddle point and therefore its decay is exponentially fast.
In the above calculation we assumed x = 0 (ζ = 0), which corresponds at looking at correlations between the two different aforementioned vertical paths (see again Fig. 2, main text). In order to study the decay of current-current correlations along any of such vertical paths, instead, we need to set equal space (and different times). Again, using translational invariance we just set x = 0 in (180). In this case we can repeat the argument in (131) and below to show with the corresponding cumulant scaling linearly in time. Finally, we recall again that the horizontal segment of the path (0, t) → (0, 0) → (x, 0) → (x, t) does not contribute to the cumulants. This observation, together with (183), and when the condition (182) is satisfied, allows to conclude that within a GGE: (i) the SCGF associated to the two-point function of the pair-mode twist fields can be correctly evaluated along this path (i.e., BFT is applicable); (ii) the total SCGF factorises into those associated to the two vertical cuts, namely it is the sum of the two corresponding SCGFs (again we refer to Remark 3.3 in [24]), and, in fact, in the main text, we make use of such factorization.

B.5 Approach to the GGE
Finally in this subsection, following B.3, we study in more detail the approach in time to the corresponding GGE values of the same single-mode densities and currents correlations considered in B.4. In particular we want to understand how cumulants in the GGE's are modified at large time when taking into account correlations coming from the initial state |ψ〉 (cf. Eq. (139)). Below, we denote by 〈·〉 expectation values on such initial state, while 〈·〉 ρ w the ones on the GGE atteined at infinite time. This time we cannot use time-translation invariance as done before when computing the expectation on a GGE. The first part in the above expression is again the GGE contribution alayzed in (181), while the second comes from quasi-particle pairs and gives rise to saddle points in the ballistic regime. It can be checked that the condition for the saddle point contribution not to vanish is exactly the same as for the density, Eq.(186). Therefore, only when the condition (186) does not hold, the contribution of the second term in (187) is subleading wrt to the GGE one, and the associated cumulant is not shifted with respect to the leading GGE value.
Therefore, depending on the value of θ (apart for θ in a region of order ε, which is to trace back to our regularization of the observables) either the correlations of single-mode densities or those of single-mode currents show a fast approach to their GGE value, thus imposing the right path to choose when using BFT.

C S matrix in the α−copy theory
Consider S(θ , θ ) to be the S matrix in the single-copy theory (we consider the diagonal case for simplicity, but everything below can be generalized to non-diagonal S matrices), i.e.
|θ , θ 〉 = S(θ , θ )|θ , θ 〉 (namely, it is the factor we get by exchanging θ , θ in the two-particle state), then one can define in the α−copies theory acting on the state |θ , i; θ , i 〉, where we introduced explicitly the dependence on the copy-index. Importantly, the ± sign in (189) depends on the commutation relations of the fields among the copies. Note that Eq. (189) is the S-matrix associated to α independent copies, which only describes the symmetry of the α-copies theory (i.e., it does not take into account the constraints of the fields in different copies implemented by the twist fields exchange relations (cf. (45)-(46))). However this is enough for our purposes.
(191) This can be written explicitly as a 2α×2α matrix S (α) (θ , θ ) m,n for α ∈ , with row and column indices m = p, p and n = k, k respectively. For example, for α = 2, it takes the form