Linear growth of the entanglement entropy and the Kolmogorov-Sinai rate

The rate of entropy production in a classical dynamical system is characterized by the Kolmogorov-Sinai entropy rate $h_{\mathrm{KS}}$ given by the sum of all positive Lyapunov exponents of the system. We prove a quantum version of this result valid for bosonic systems with unstable quadratic Hamiltonian. The derivation takes into account the case of time-dependent Hamiltonians with Floquet instabilities. We show that the entanglement entropy $S_A$ of a Gaussian state grows linearly for large times in unstable systems, with a rate $\Lambda_A \leq h_{KS}$ determined by the Lyapunov exponents and the choice of the subsystem $A$. We apply our results to the analysis of entanglement production in unstable quadratic potentials and due to periodic quantum quenches in many-body quantum systems. Our results are relevant for quantum field theory, for which we present three applications: a scalar field in a symmetry-breaking potential, parametric resonance during post-inflationary reheating and cosmological perturbations during inflation. Finally, we conjecture that the same rate $\Lambda_A$ appears in the entanglement growth of chaotic quantum systems prepared in a semiclassical state.


Introduction
Entanglement plays a central role in the thermalization of isolated quantum systems [1][2][3]. The paradigmatic setting consists in a Hamiltonian system prepared in a pure state and evolving unitarily, |ψ t = e −iHt |ψ 0 . The objective is to study the thermalization of observables O A belonging to a subalgebra of observables A A which define a bipartition H = H A ⊗ H B of the system in a subsystem A and its complement B. While the von Neumann entropy of the system vanishes at all times, the entropy of the subsystem A, in general does not vanish and has a non-trivial evolution. The origin of this entropy is the entanglement between the degrees of freedom in the subsystem A and its complement. Equilibration in the subsystem A occurs when the entanglement entropy S A (t) approaches an equilibrium value S eq , with thermalization corresponding to S eq given by the thermal entropy.
A generic behavior has been observed for various systems prepared in a state with initially low entanglement entropy, S A (t 0 ) S eq : After a transient which depends on the details of the initial state of the system, the entropy of the subsystem goes through a phase of linear growth, until it saturates to an equilibrium value as described in figure 1. This behavior is observed in the evolution of various isolated quantum systems, in particular in systems that show the signatures of quantum chaos [4][5][6][7][8], in many-body quantum systems [9] and quantum fields [10][11][12] after a quench, and in the thermalization of strongly-interacting quantum field theories studied using holographic methods [13][14][15][16][17]. Understanding the mechanism of this process is of direct relevance for the puzzle of fast thermalization of the quark gluon plasma produced in heavy-ion collisions [18][19][20], in models of black holes as fast scramblers of quantum information [21], and in the study of the quantum nature of space-time [22][23][24][25][26][27]. In particular, being able to predict from first principles the rate of growth Λ A of the entanglement entropy in the phase of linear growth can provide us with crucial information on the time-scale of thermalization.
On the other hand, at the classical level -in Hamiltonian chaotic systems -the coarsegrained entropy S cl (t) shows a behavior similar to the one described in figure 1, with a linear phase which has a known rate of growth h KS [28,29], where h KS is the Kolmogorov-Sinai rate of the system, an information-theoretic quantity that measures the uncertainty remaining on the future state of a system, once an infinitely long past is known. The Kolmogorov-Sinai rate has dimension of time −1 and for regular Hamiltonian systems is given by the sum of the positive Lyapunov exponents λ i of the system [30][31][32][33].
In quantum systems that have a classical chaotic counterpart, a relation between the rate of growth of the entanglement entropy Λ A and the classical Lyapunov exponents λ i is S eq τ eq Initial transient Linear production Saturation phase expected [4][5][6][7][8][34][35][36], despite the fact that Lyapunov exponents are global quantities which probe the phase space of the full system, not just of the subsystem A.
In this paper we investigate the relation between Λ A and the Lyapunov exponents λ i by studying the evolution of Gaussian states in many-body systems and quantum field theories with quadratic time-dependent Hamiltonians. Non-trivial Lyapunov exponents arise in the presence of instabilities and of parametric resonances. In this context we prove that the linear growth of the entanglement entropy S A (t) has a classical counterpart: The entanglement rate Λ A equals the exponential rate of growth of the volume of a cell in the sub phase space of the subsystem A. We then provide an algorithm for computing Λ A in terms of the Lyapunov exponents λ i of the classical system and the choice of subsystem A. The methods developed apply both to quantum systems with finitely many degrees of freedom and to quantum fields in external time-dependent backgrounds when the subsystem is given by a finitely-generated Weyl subalgebra A A .
The paper is organized as follows. In section 2, we present our main result, theorem 1, which determines the asymptotic rate of growth of the entanglement entropy for Gaussian states in systems with quadratic Hamiltonians and establish its relation to the Kolmogorov-Sinai rate, including bounds for non-Gaussian initial states. Our results are then applied in section 3 to the study of entanglement production in three example systems with finitely many degrees of freedom, including unstable potentials and periodic quantum quenches. Section 4 is dedicated to quantum field theories, where we consider again three example systems: a scalar field in a symmetry-breaking potential, parametric resonance during post-inflationary reheating and cosmological perturbations during inflation. We show that our results agree with numerical evaluations of the entanglement entropy for sufficiently large times. In sections 5 and 6 we present technical results required for the derivation of our main result in section 2. In particular, section 6 reflects the structure of our proof for theorem 1. Finally, in section 7 we discuss limitations and possible extensions of our work, and in particular a conjecture on entanglement production in chaotic systems. Moreover, we discuss the relation to linear growth of the entanglement entropy after a generic quantum quench. The paper is supplemented with appendices where we provide a summary of the relevant results in the study of dynamical systems and Lyapunov exponents and of general symplectic techniques for the study of the time-evolution and entanglement and Rényi entropies of Gaussian states.

Results: Linear growth of the entanglement entropy
We state the main result which relates the asymptotic rate of growth of the entanglement entropy of a quantum system to classical instabilities encoded in the Lyapunov exponents of the classical system. Our proof is based on a set of technical results presented in sections 5 and 6.

Entanglement entropy growth, instabilities and the volume exponent
We consider a quadratic bosonic system with N degrees of freedom. We denote linear observables by ξ a = (q 1 , . . . , q N , p 1 , . . . , p N ) and assume canonical commutation relations [q i , q j ] = [p i , p j ] = 0 and [q i , p j ] = iδ ij . These relations can be more compactly phrased by stating [ξ a , ξ b ] = iΩ ab where Ω ab is a symplectic form. The most general quadratic Hamiltonian is given by which is independent from the metric we choose to measure the length. A system decomposition V = A ⊕ B of the classical phase space into subsystem phase spaces A and B induces an equivalent decomposition V * = A * ⊕B * of the dual phase space. Here, we can generalize the notion of Lyapunov exponents to define the subsystem exponent Λ A defined by where V A ⊂ A * is an arbitrary parallelepiped in the subspace A * . The subsystem exponent captures the exponential volume growth of subsystem regions. The volume vol is measured on the subspace where M (t) V A lives, but the subsystem exponent is independent of the global metric on V * one chooses to define the volume form on arbitrary subspaces. We explain the relation between Λ A and λ in section 5, while more technical details are summarized in appendix A.
• Entanglement of Gaussian states It is well-known that a Gaussian bosonic state |ψ can be completely characterized by its expectation value ζ a = ψ|ξ a |ψ and its covariance matrix G ab = ψ|ξ a ξ b + ξ b ξ a |ψ − 2 ζ a ζ b . Recent progress on unifying methods for bosonic and fermionic Gaussian states [37] suggest an equivalent description where G ab is replaced by a linear complex structure J a b = −G ac ω cb with ω being the inverse of Ω. Choosing a system decomposition A⊕B with complementary subsystems A and B allows us to compute the entanglement entropy S A (|ψ ) between them. For a Gaussian state |ψ , this entanglement entropy can be directly computed from J, which we use in section 6 and review in appendix B.
(i) The entanglement entropy is bounded by the Renyi entropy: We define the asymptotic rate of growth of the entanglement entropy as its long-time linear scaling lim t→∞ 1 t S A (U (t)|J 0 , ζ 0 ). We note that quadratic time-dependent Hamiltonians evolve the initial Gaussian state into a Gaussian state, (B.18). In section 6.1 we prove that the entanglement entropy of a Gaussian state is bounded from below by the Rényi entropy R A (U (t)|J 0 , ζ 0 ) and from above by the Rényi entropy plus a state-independent constant, inequality (6.4). Therefore, we have the equality i.e., the asymptotic rate of growth of the entanglement entropy and of the Rényi entropy coincide. (ii) The Renyi entropy is given by a phase space volume: In section 6.2 we prove that the Rényi entropy of a Gaussian state equals the logarithm of the phase space volume of a symplectic cube V A spanning the subsystem A, (6.8).
The volume is measured with respect to the metric induced by the state, (6.5). In the case of the time-dependent Gaussian state U (t)|J 0 , ζ 0 , we can measure the volume with respect to the time-dependent induced metric G t = M (t)G 0 M (t). Equivalently, we can consider the time-dependent symplectic cube M (t)V A and measure its volume with respect to the initial metric G 0 induced by the initial state, (iii) The Renyi entropy grows as regions in phase space are stretched: The subsystem exponent Λ A introduced in Eq. (2.4) and discussed in section 5.3 provides a generalization of the notion of Lyapunov exponents of a classical Hamiltonian system. It involves the choice of a subsystem A, a symplectic dynamics M (t) and a reference metric G 0 , Despite the metric G 0 is needed for the definition, the value of the subsystem exponent Λ A is independent of G 0 for regular Hamiltonian systems. The subsystem exponent can be expressed in terms of the Lyapunov exponents of the system using the algorithm described in theorem 3, (5.32).
Using (i), (ii) and (iii), we find that the asymptotic rate of growth of the entanglement entropy is given by the subsystem exponent Λ A , for all initial Gaussian states, therefore proving the statement of the theorem.
We note that, as the entanglement entropies of complementary subsystems A and B coincide, S A (|ψ ) = S B (|ψ ), also their asymptotic rates of growth have to coincide. Consistency with the statement of the theorem implies that, at the classical level, the subsystem exponents defined in section 5.3 for a symplectic decomposition V = A ⊕ B coincide (2.10) This statement can be proven using the expression (5.32) of the subsystem exponents or more directly using the property det[J] A = det[J] B for the restriction of a complex structure J to complementary symplectic subspaces.

Entanglement and the Kolmogorov-Sinai entropy rate
Theorem 2 (Entanglement growth -generic subsystem). Given a quadratic time-dependent Hamiltonian H(t) with Lyapunov exponents λ i , the long-time behavior of the entanglement entropy of a generic subsystem A is for all initial Gaussian states |J 0 , ζ 0 and all generic subsystems with N A degrees of freedom.
In particular, the rate of growth of the entanglement entropy is bounded from above by the Kolmogorov-Sinai rate h KS , (2.12) The decomposition in two complementary subsystems both with dimension larger than the number of instabilities results in an entanglement growth proportional to the Kolmogorov-Sinai rate, 13) and therefore saturates the bound (2.12).
Proof. The asymptotic rate of growth of the entanglement entropy of a Gaussian state is given by the subsystem exponent Λ A as stated in theorem 1, (2.9). For a generic subsystem, theorem 4 states that the subsystem exponent equals the sum of the 2N A largest Lyapunov exponents, (5.45). Together with Pesin's theorem (5.43), this result implies that the asymptotic rate of growth is bounded from above by the Kolmogorov-Sinai rate of the system, Moreover, the subsystem exponent Λ A equals the Kolmogorov-Sinai rate h KS when its dimension is in the range N I ≤ 2N A ≤ 2N − N I , (5.47). Recalling that N A + N B = N , this range coincides with the requirement that the dimension of each of the two complementary subsystems is larger than the number of instabilities, 2N A ≥ N I and 2N B ≥ N I . In this case the bound (2.12) is saturated.
We note that quantum many-body systems often have only a small finite number of unstable directions N I compared to the number of degrees of freedom of the system, N I N . A generic decomposition in two complementary subsystems that encompass the fractions f A = N A /N and f B = 1 − f A of the full system satisfies (2.13) if the fractions are in the range As a result, in the limit N → ∞ with N I finite, we have that the entanglement growth is proportional to the Kolmogorov-Sinai rate S A (t) ∼ h KS t for all partitions of the system into two complementary subsystems each spanning a finite fraction of the system, except for a set of partitions of measure zero.

Bounds on non-Gaussian initial states
Computing the entanglement entropy growth of non-Gaussian states is a non-trivial problem as efficient tools similar to the ones discussed in (B.45) are not available. Nevertheless upper bounds that generalize theorems 1 and 2 can be established in the case of evolution driven by a quadratic time-dependent Hamiltonian. Let us consider an initial non-Gaussian state |ψ 0 and the unitary evolution U (t) generated by a quadratic time-dependent Hamiltonian of the most general form described in (B.14). The symmetric part of the connected 2-point function at the time t is given by where |ψ t ≡ U (t)|ψ 0 and M a b (t) is the symplectic matrix defined in (5.12). There always exists a mixed Gaussian state ρ 0 which has the same correlation function G ab 0 at the time t = 0 [38]. By construction, the 2-point function of the unitarily evolved Gaussian state U (t)ρ 0 U −1 (t) is the function G ab t of (2.16). Moreover one can show that the entanglement entropy of the non-Gaussian state |ψ t is bounded from above by the entanglement entropy of the mixed Gaussian state having the same 2-point function G ab t , i.e. S A (ρ N G ) ≤ S A (ρ G ) where ρ N G = Tr B (|ψ t ψ t |) is the reduced density matrix of the non-Gaussian state, and is the reduced density matrix of the Gaussian state [36]. The proof is immediate: Recalling that the relative entropy is a positive function [39,40], we have where S(ρ N G ρ G ) is the relative entropy and the term in parenthesis vanishes as the two states have the same correlation function by construction. On the other hand, theorem 1 generalizes to mixed Gaussian states implying the asymptotic growth S A (ρ G ) ∼ Λ A t for the entanglement entropy of a subsystem A. As a result we find the inequality which states that the asymptotic rate of growth of the entanglement entropy of a non-Gaussian state |ψ t which evolves unitarily with a quadratic Hamiltonian is bounded from above by the subsystem exponent Λ A . This result generalizes theorems 1 and 2 to non-Gaussian states and establishes the Kolmogorov-Sinai rate h KS as the upper bound for the asymptotic rate of growth of the entanglement entropy. Preliminary numerical investigations indicate that, under the unitary evolution given by a quadratic time-dependent Hamiltonian, the upper bound Λ A in (2.19) might in fact be saturated by all initial states and not just by Gaussian states [41].

Applications: unstable potentials and periodic quenches
We briefly discuss three examples of simple systems that show a linearly growing entanglement entropy and allow us to test our results.

Particle in a 2d inverted potential
In our first example, we study a simple system consisting of just two degrees of freedom. Despite its simplicity, the example captures the main features of the theorems presented above. It also resembles the system studied in [35] and thereby illustrates how our theorems simplify the involved steps to understand the asymptotic behavior of the entanglement entropy. We consider a system which can be described as a quantum particle with mass m = 1 moving in a plane with coordinates (q 1 , q 2 ) and corresponding momenta (p 1 , p 2 ). The instabilities arise from an inverted harmonic potential V (q 1 , q 2 ) = − λ 2 1 2 q 2 1 − λ 2 2 2 q 2 2 with λ 1 ≥ λ 2 > 0. The Hamiltonian of this system is explicitly given by If we choose the Darboux basis D V = (p 1 , p 2 , q 1 , q 2 ), the matrices h and K = Ωh become The Lyapunov exponents (λ 1 , λ 2 , −λ 2 , −λ 1 ) are given by the eigenvalues of K and the Lyapunov basis D L = ( 1 , 2 , 3 , 4 ) = (Q 1 , Q 2 , P 2 , P 1 ) are the corresponding eigenvectors (3) t Figure 2. Particle in a 2d inverted potential. The plot shows the exact behavior of the Rényi entropy R A (t) (thick) and entanglement entropy S A (t) (thin) in comparison to the predicted asymptotics Λt (dashed). The system is defined in (3.1) with subsystems specified in (3.4-3.6). For the computation, we choose λ 1 = −λ 4 = 2 and λ 2 = −λ 3 = 1/2. The initial state is chosen to be |J 0 with associated metric G 0 ( i , j ) = δ ij . In the case of examples (1) and (2), With these definitions, let us consider the three different choices of subsystem A discussed also in section 5.4: We can study the entanglement entropy for these subsystems numerically where we start with the (entangled) initial state given by G 0 ( i , j ) = δ ij . Figure 2 shows excellent agreement with our predictions. In particular, we also see that the entanglement entropy S A (t) and the Rényi entropy R A (t) only differ by the constant 1 − log(2) if the system is strongly entangled.

Quadratic potential with instabilities
The second example consists in the evolution in a time-independent potential with instabilities. Let us consider a classical system with N degrees of freedom which we parametrize by N conjugate pairs (q i , p i ) of coordinates in phase space. The Hamiltonian H of the system consists of a standard kinetic term and a quadratic potential, The potential is determined by the symmetric matrix V ij with eigenvalues v i . This classical system has 2N Lyapunov exponents λ i determined by the eigenvalues of V ij and given by Positive eigenvalues correspond to stable directions of the potential, lead to oscillatory motion and vanishing Lyapunov exponent. In the presence of negative eigenvalues v i < 0, the system is unstable, the classical motion is unbounded and nearby trajectories in phase space diverge at an exponential rate given by the λ i = +Im( √ v i ). Now we consider the associated quantum system prepared in a Gaussian state and study the behavior of a generic subsystem. If the potential V ij couples the subsystem A with the rest of the system, we expect that the entanglement entropy of the subsystem changes in time. Our theorem states that the entanglement entropy of a generic subsystem A with N A degrees of freedom asymptotically grows at a rate given by the sum of the 2N A largest Lyapunov exponents, (2.11). We give a concrete example: We consider a system with N = 20 degrees of freedom and quadratic potential specified by a N × N real symmetric random matrix V ij . Negative eigenvalues of V ij correspond to unstable directions of the potential and non-vanishing Lyapunov exponents with the following values: shows the growth of the entanglement entropy of an initially un-entangled Gaussian state for generic subsystems of different dimensions. For a one-dimensional subsystem N A = 1, theorem 2 predicts the asymptotic growth S A (t) ∼ (λ 1 + λ 2 ) t. Note that, as the Lyapunov exponents appear in couples (λ, −λ), the asymptotic growth of the rest of the system has the same rate, S B (t) ∼ (λ 1 + . . . + λ 18 ) t = (λ 1 + λ 2 ) t, as expected for the entanglement entropy of a pure state. Clearly, the statement that the entanglement growth is linear in time applies only to generic subsystems and not to subsystems that are aligned to the shape of the potential. As an example, let us call (Q 1 , . . . , Q N , P N , . . . , P 1 ) the Lyapunov basis of the system and consider the subsystem spanned by the canonical couple (Q 1 , P 1 ). Theorem 1 predicts a sublinear rate , which is consistent with the statement that the entropy will stay constant as the subsystem is isolated. Moreover, note that the difference of Lyapunov exponents can also appear in the asymptotic rate of the entanglement growth, for instance the subsystem spanned by (Q 2 , P 2 + P 3 ) has asymptotic rate S A (t) ∼ (λ 2 − λ 3 ) t.  We plot five subsystems with 1 ≤ N A ≤ 5. The entanglement entropy (colored lines) agrees with the predicted asymptotics (dashed lines). Note that for N A ≥ 3, the asymptotic behavior is the same due to the stable Lyapunov exponents λ i = 0 for 6 ≤ i ≤ 15. These are exactly the cases for which we have Λ A = h KS , which means that entropy production rate coincides with the classical Kolmogorov-Sinai entropy rate. This is the quantum version of the example discussed in (5.39). It is important to remark however that these subsystems are non-generic and form a subset of measure zero as discussed in the proof of theorem 4. The random quadratic potential V ij with Lyapunov exponents specified in (3.8) has N I = 5 unstable directions and classical Kolmogorov-Sinai rate h KS = λ 1 + . . . + λ 5 1.94. Theorem 2 states that at the quantum level the long-time behavior of the entanglement entropy is linear with rate h KS for all generic subsystem decompositions such that 2N A ≥ N I and 2N B ≥ N I , i.e., S A (t) ∼ h KS t for generic subsystems of dimension 3 ≤ N A ≤ 18. Figure 3 shows the numerical evolution of the entanglement entropy and the cases N A = 3, N A = 4, N A = 5 exhibit a linear growth with rate given by the Kolmogorov-Sinai rate as predicted.
A remarkable feature of the predictions stated in theorems 1 and 2 is that the asymptotic growth of the entanglement entropy is determined by the subsystem and completely independent from the choice of initial Gaussian state. This feature is a consequence at the quantum level of the fact that the Lyapunov exponents of a classical system are independent of the choice of metric used to measure the distance between trajectories. In the language of complex structures J 0 that specify the initial Gaussian state |J 0 , ζ 0 , the Lyapunov exponents are given by the eigenvalues of the matrix L = lim t→∞ 1 2t log M −1 (t) J 0 M (t) defined in (5.19) and are independent of J 0 . Figure 3 shows only initial states with vanishing entanglement entropy, but the two theorems 1 and 2 apply to all Gaussian states, even to ones that have large initial entanglement entropy. Clearly the theorem applies only to the asymptotic behavior of the entanglement entropy. In fact we could take as initial state the time-reversal of the Gaussian state used in figure 3 at late times. In this case the entanglement entropy would initially decrease, reach a minimum and eventually start growing linearly as predicted by the theorems on the asymptotic growth.

Periodic quantum quenches in a harmonic lattice
As a third example of a system that displays a linear growth of the entanglement entropy, we discuss the case of a harmonic lattice subject to periodic quantum quenches. The Hamiltonian of the system is which describes the dynamics of a one-dimensional chain of N bosons with nearest-neighbor coupling κ and boundary conditions q N +1 = q 1 , p N +1 = p 1 . The one-particle oscillation frequency Ω(t) is periodically switched between the values Ω 0 ± ε with period 2 T 0 , Ω(t + 2 T 0 ) = Ω(t) and ε Ω 0 . (3.11) The system is prepared in the ground state of the instantaneous Hamiltonian H(0) at the time t = 0 and then let evolve unitarily. The state of the system at stroboscopic times t n = 2 n T 0 which are multiples if the period 2 T 0 can be obtained by computing the Floquet Hamiltonian of the system, H F : where H 1 = H(T 0 ) and H 2 = H(2 T 0 ). We note that the Floquet Hamiltonian H F is quadratic and time-independent, but is not of the standard form consisting of the sum of a kinetic and a potential term as in the case of (3.7). This general form is taken into account in theorems 1 and 2 which apply to all quadratic Hamiltonians. The classical system described by (3.9) shows instabilities when small perturbations from the equilibrium configuration are amplified via the mechanism of parametric resonance. Floquet theory [42,43] provides the tools for the description of the dynamics driven by an Hamiltonian which is periodic in time, as is the case for (3.9). The eigenvalues of Figure 4. Instability bands in a periodic quantum quench. We sketch the instability region of Floquet exponents µ with positive real part (red: Re(µ) > 0.03, blue: Re(µ) > 0.05) for the system described in (3.13) as a function of ε and p k . We chose the values T 0 π, Ω 0 0.3 and κ 0.3. Furthermore, we indicate the discrete momenta p k for a periodic chain with N 20 and ε 0.735, such that there are two modes with unstable Floquet exponents, namely p ±4 = ± 2π 5 .
the symplectic evolution matrix (5.12) evaluated at a period M (2 T 0 ) come in quadruplets (e +2 T 0 µ , e −2 T 0 µ , e +2 T 0 µ * , e −2 T 0 µ * ) where the complex numbers µ are the Floquet exponents of the system. The stability of the system is measured by the real part of the Floquet exponents which coincide with the Lyapunov exponents, λ = Re(µ).
The Lyapunov exponents of the system (3.9) can be easily determined. In Fourier transformed variables 1 Q k , P k with k = 0, ±1, ±2, . . . , ±(N − 1)/2, the Hamiltonian takes the form with ω k (t) ≡ Ω 2 (t) + 4κ sin 2 (p k /2) and p k ≡ 2πk N . (3.14) In particular, the speed of sound of the mode of momentum p k switches periodically between the values v k (T 0 ) and v k (2 T 0 ), with v k (t) ≡ ∂ω k (t)/∂p k . As Fourier modes with different |k| are decoupled, we can analyze the stability of the system mode by mode. The classical evolution of the coupled modes (Q k , Q −k , P −k , P k ) is given by given by a symplectic block of the symplectic matrix M (2 T 0 ) defined in (5.12), and ω 1 = ω(T 0 ) and ω 2 = ω(2 T 0 ). The Lyapunov exponents of the system are the real parts of the Floquet exponents, i.e.
Analytic expressions of λ k can be found assuming that the periodic perturbation is small ε Ω 0 and the mode is at or near a parametric resonance. Defining δω ≡ ω 2 − ω 1 , ω 0 ≡ ω 1 +ω 2 2 , with δω ω 0 , we find that the system is in parametric resonance when the average frequency ω 0 of the mode is an half-integer multiple of the frequency of the perturbation, i.e., At the parametric resonance, the positive Lyapunov exponents of the system are given by if n odd, if n even.
For a finite perturbation, the stability of the system can be determined numerically. Figure 4 shows which modes p k are unstable for a given finite value of the perturbation parameter ε.
In the example we have N = 20, Ω 0 0.3, κ 0.3 and T 0 π. For ε 0.735 we have two unstable modes with k = ±4 and Lyapunov exponents λ +4 ±0.046, λ −4 ±0.046. The Kolmogorov-Sinai rate of the system is h KS = 0.092. Figure 5 shows the growth of the entanglement entropy of a subsystem and the relation to h KS .

Applications: quantum field theory
In section 2 we presented our main results for a bosonic quantum system with N degrees of freedom. These results can be extended with minor modifications to the case of a bosonic quantum field. In particular, the formulation of theorems 1 and 2 in terms of complex structures J is motivated by and tailored to applications to quantum field theory in curved spacetimes [44][45][46].

Definition of a subsystem and the algebraic approach
The presence of infinitely many degrees of freedom in quantum field theory has two immediate consequences which are relevant for our analysis [47]: i) the existence of unitarily inequivalent representations of the algebra of observables, ii) the lack of a factorization of the Hilbert space into a tensor product over local factors.
The algebraic approach to quantum field theory -together with the language of complex structures -provides a natural setting for discussing both aspects and formulating the analysis of the growth of the entanglement entropy of a subsystem in quantum field theory.
At the classical level, the phase space V of a free scalar field has coordinates ξ a = (ϕ( x), π( x)) with x a point on a Cauchy slice Σ. We adopt abstract indices and use the symbol ω ab for the symplectic form on the infinite-dimensional vector space V . Carrying out the rigorous construction of the infinite-dimensional phase space requires the choice of a positive definite metric g ab compatible with the symplectic form ω ab , such that V arises as the completion with respect to this metric. Contracting ω ab with the inverse metric G ab gives rise to the complex structure J a b = −G ac ω cb : V → V . Given a reference complex structure J 0 and a symplectic transformation M , we can define a transformed complex structure J M = M −1 J 0 M . The transformation M is said to belong to the restricted symplectic group if the commutator A = [J 0 , J M ] is a Hilbert-Schmidt operator, i.e. tr(A † A) < +∞ [48][49][50].
At the quantum level, the choice of a complex structure J 0 defines a Gaussian state |J 0 which can be used as vacuum for building a Fock representation of the algebra of observables [44][45][46]. Representations built over Fock vacua |J 0 and |J M are unitarily equivalent if and only if the symplectic transformation M belongs to the restricted symplectic group described above. When interpreted in terms of particle excitations, the state |J M describes a superposition of particle pairs over the vacuum |J 0 . A symplectic transformation M which does not belong to the restricted group corresponds to a Bogoliubov transformation that produces an infinite number of particles [51].
Gaussian states and quadratic time-dependent Hamiltonians appear in the description of particle production in the early universe [52,53], Hawking radiation in black hole evaporation [54], in the Schwinger effect [55,56], in the dynamical Casimir effect [57,58] and more generally in all cases where the free quantum field evolves in a time-dependent background. It is known that, for some time-dependent backgrounds, the time-evolution -which, at the classical level, is encoded in a symplectic transformation M -cannot be implemented as a unitary operator in a Fock space at the quantum level [59]. Nevertheless the correlation functions in the quantum theory are still well-defined in terms of a complex structure J 0 and a symplectic transformation M as described in (B.13) [60]. The algebraic approach focuses on correlation functions and does not involve the construction of a Fock space. It provides sufficient structure for defining the (abstract) state of the system and computing the evolution of the entanglement entropy of a subsystem, despite the potential lack of a standard unitary implementation of the time evolution in a Fock space, (i).
The second aspect which needs some clarification regards the definition of a subsystem in quantum field theory, (ii). It is a well-known fact about the ground state of a quantum field that the entanglement entropy of a region of space is divergent and -when an ultraviolet cutoff is introduced -it scales as the area of the boundary of the region [61][62][63]. The divergence of the geometric entanglement entropy has an algebraic origin: The local subalgebra of observables associated to a region in space is of type III, i.e. it does not identify a factorization of the Fock space in a tensor product of Hilbert spaces [47,64]. Three standard strategies to address this issue are: (a) a modification the ultraviolet behavior of the theory, for instance introducing a lattice cut-off [61,62,65], or (b) computing the mutual information between a region and a carved version of its complement so to introduce a "safety corridor" [64,66,67], or (c) focusing on the excess entropy of a state with respect to the one of the ground state [67,68]. Here we illustrate a different strategy: We focus on the entanglement entropy of a subsystem with a finite number N A of degrees of freedom. The geometric entanglement entropy which captures infinitely many degrees of freedom can be recovered in the limit of increasingly large subsystems [69].
A simple example of a subsystem with a single degree of freedom, N A = 1, is provided by a linear smearing of the fields against given test functions f ( x) and g( x): (4.1) The observablesφ f andπ f generate a Weyl algebra A A of type I which, as in section 6, induces a factorization of the Hilbert space into H = H A ⊗ H B . The symplectic structure Ω A of the subsystem can be computed from the commutator, Given a Gaussian state |J of the quantum field, the symmetrized correlation function restricted to the subsystem A is A ) cb and its eigenvalues ±ν determine the entanglement entropy of the subsystem A through (B.44).
We illustrate our result on three paradigmatic cases in quantum field theory where the time-dependence of the entanglement entropy of a subsystem can be computed and our results on the linear growth can be tested.

Dynamics of symmetry breaking and the inverted quadratic potential
We consider a scalar field ϕ(x) which goes through a symmetry breaking transition in real time [70,71]. A simple model is described by the action 2 with a quartic potential, The quadratic coupling α(t) is chosen so that, for t > 0, a minimum of the potential breaks the symmetry ϕ → −ϕ. We set The system is initially prepared in the ground state at t < 0 and then let evolve. For small quartic coupling ε and short time, the evolution is described perturbatively by a tachyonic instability: At the onset of the symmetry-breaking transition, the scalar field evolves as if it was free and had a negative mass-squared, −µ 2 . We focus on this initial phase.
We set ε = 0 and study the free evolution governed by a quadratic Hamiltonian which transitions from a stable phase to an unstable phase. It is useful to adopt Fourier transformed canonical variables so that the canonical commutation relations read [ϕ( k), π( k )] = i(2π) 3 δ 3 ( k + k ) and the symplectic structure in these coordinates is For t < 0, the Hamiltonian is and the system is stable. The ground state is the Gaussian state |J 0 with correlation functions The complex structure of the ground state is therefore For t > 0, that is, after the transition from the stable to the unstable phase, the Hamiltonian governing the free evolution is given by Modes with k 2 smaller than µ 2 are unstable and have Lyapunov exponents which come in pairs ±λ( k), with (4.14) As a result, infrared modes are unstable and the largest Lyapunov exponent λ(0) = µ is associated to the homogeneous mode. The classical evolution generated by the unstable Hamiltonian (4.13) is given by the As a result, in the quantum theory, the correlation functions at the time t for the system initially prepared in the Gaussian state |J 0 are given by the evolved complex structure J t ( k, k ), which defines a Gaussian state |J t at the time t.
Let us consider a measuring device which probes the field and its momentum only in a neighborhood of the point x = 0 with a linear size R. This device can be modeled by a Gaussian smearing function f ( x). The subsystem A defined by such measurements is encoded in the subalgebra of observables A A generated bŷ This subsystem has N A = 1 bosonic degrees of freedom. Fluctuations of the observablesφ A andπ A at the time t are encoded in the correlation functions of the subsystem . Similarly, the restricted complex structure is the 2 × 2 matrix The eigenvalues of i [J t ] A come in pairs ±ν(t) and the entanglement entropy of the subsystem A is given by where S(ν) is the function (B.44). The predicted asymptotic rate of growth of the entanglement entropy of a subsystem with N A = 1 is given by the subsystem exponent Λ A = 2µ, which is the sum of the two largest Lyapunov exponents, i.e.,  Figure 6. Symmetry breaking and the inverted quadratic potential. We compute the entanglement entropy S A (t) numerically for the time evolution with the unstable Hamiltonian (4.13). We set µ = 1.
The subsystem A is defined in (4.17) with R = 1.
A numerical plot of the entanglement entropy as a function of time, together with the predicted rate of growth, is shown in figure 6.
We note that, as the positive Lyapunov exponents of the system appear in a continuous band λ( k) = µ 2 − k 2 , the prediction for the asymptotic growth of the entanglement entropy of a subsystem with N A degrees of freedom is simply S A (t) ∼ 2N A µ t. In the case of a subsystem with infinitely many degrees of freedom, it is useful to induce an infrared cutoff, for instance a cubic volume V = L 3 . The boundary conditions induce a quantization of the momentum k = 2π L n x , 2π L n y , 2π L n z which splits the degeneracy of the Lyapunov exponents and results in a discrete sequence λ( k). We can now define the number N I of unstable degrees of freedom of the system. In the limit L 2π µ we find A generic subsystem which probes infinitely many degrees of freedom, as in the case of the geometric entanglement entropy of a region of space, would probe all the unstable degrees of freedom of the system. As a result the asymptotic growth of the entanglement entropy is expected to be given by where h KS is the Kolmogorov-Sinai rate per unit volume, As a result, for a generic subsystem with infinitely many degrees of freedom, the quantity h KS describes the asymptotic behavior of the entanglement entropy per unit space-time volume.
In this analysis we assumed that the quartic term 1 4! ε ϕ 4 is not present in the potential and the evolution is simply described by a quadratic Hamiltonian with instabilities. In section 7.1 we discuss when this approximation is expected to be valid.

Preheating and parametric resonance
The simplest model of parametric resonance in quantum field theory is described by the Hamiltonian which has a quadratic potential that oscillates in time with period 2π/M 0 . For small values of the amplitude of oscillation, we have a narrow resonance band around the frequency of the perturbation, Modes with momentum | k| ≈ M 0 are parametrically amplified. The Lyapunov exponents of the system can be determined via Floquet analysis as we already did in section 3.3. The canonical subsystem spanned by ϕ( k), π(− k) with k in the band (4.27) has Lyapunov exponents ±λ( k) with Given an initial Gaussian state -for instance the Minkowski vacuum -the evolution of the correlation functions of the system can be computed analytically in terms of Mathieu functions. In figure 7 we show the time-evolution of the entanglement entropy of a subsystem A defined by a subalgebra of observables A A generated by the linear observables (4.17). At the classical level, the subsystem exponent Λ A can be easily computed: It is given by the sum of the two largest Lyapunov exponents of the systems, which are degenerate in value and correspond to modes exactly at the resonance | k| = M 0 . Therefore we have Preheating and parametric resonance. We compute the entanglement entropy S A (t) numerically. The subsystem is defined in (4.17), the time evolution is governed by the Hamiltonian presented in (4.25) and we start in the Minkowski vacuum state. We set R = 1, M 0 = 1 and v 2 0 = 0.1 leading to the entanglement production rate Λ A = 0.05. The quantity S A (0) is the entanglement entropy of the Minkowski vacuum and the linear production phase is reached after an initial transient.
As predicted by theorem 1, the entanglement entropy of the subsystem initially prepared in a Gaussian state |J 0 is observed to grow as S A (t) ∼ Λ A t.
At the end of cosmological inflation, the inflaton oscillates coherently around the minimum of its potential. Such oscillations excite the vacuum of matter fields via the phenomenon of parametric resonance. This phase of explosive non-thermal particle production is called preheating and is followed by a thermalization phase which provides the initial conditions for Big Bang Nucleosynthesis. A simple model of preheating consists in a coupling V = 1 2 g 2 Φ 2 ϕ 2 between the inflaton Φ and a field ϕ which serves as proxy for Standard Model fields. Coherent oscillations of the inflaton, Φ( x, t) = Φ 0 sin(M 0 t), result in an effective dynamics for the matter field described by an Hamiltonian of the form (4.25) with a coupling constant v 2 0 = g 2 Φ 2 0 . At the beginning of the oscillatory phase, the state of matter can be assumed to be the vacuum |J 0 because of the dilution effect of the inflationary phase. Its evolution in the preheating phase results in a Gaussian state |J t which is far from equilibrium: The Floquet instability of the Hamiltonian results in an explosive production of particles with momenta in the resonance band. For a given subalgebra of observables A A , such as the one discussed in (4.17), theorem 1 predicts a linear growth of the entropy with a rate given by the subsystem exponents Λ A . Phenomenologically, the relevant choice of subalgebra of observables or coarse graining of the system is dictated by the interaction with its environment. The interaction of the produced particles and the expansion of the universe have the effect of reducing the efficiency of the resonance and eventually lead to a thermal-equilibrium radiation-dominated phase, with an expected entropy profile qualitatively similar to the one illustrated in figure 3.
A similar preheating phenomenon is discussed in the context of relativistic heavy-ion collisions where, in the late stages of the evolution of a quark-gluon plasma, a chirally symmetric state rolls down and oscillates around to the minimum of the effective chiral potential [79]. Coherent pion excitations are described by a quark condensate Φ = qq and φ = q σq and a O(4) linear sigma model with φ a = (Φ, φ ) with explicit symmetry breaking. The action of the system is with parameters g ≈ 20, f π ≈ 90 MeV and m π ≈ 140 MeV. The coherent field φ a is initially in a chirally symmetric state. As it rolls down the potential and oscillates around the minimum φ a ≈ (f π , 0 ), a squeezed state of coherent pion pair excitations is produced via parametric resonance. In this phase, the entanglement entropy of a generic subsystem A is predicted to grow with a rate given by the subsystem exponent Λ A . As discussed in [18], the linear entropy growth is expected to be bounded by the Kolmogorov-Sinai rate of the system.
Time-dependent Hamiltonians of the form (4.25) appear also in the description of stimulated quasi-particle production in cold atomic Bose gases [81,82]. A periodic modulation of the external potential that traps the gas induces a response in the condensed portion of the gas which acts as a time-dependent background for quasi-particles. The study of entanglement entropy growth in cold atomic Bose gases is of particular relevance because of current experiments which can probe the non-separability of phonon pair creation [83,84].

Cosmological perturbations and slow-roll inflation
During slow-roll inflation, quantum perturbations of the metric and the inflaton field are stretched and squeezed. We illustrate this phenomenon -together with the associated growth of the entanglement entropy -using a simple model consisting of a minimally-coupled massless scalar field in a cosmological spacetime. The action of the system is where with a metric g µν that defines the line element of a Friedmann-Lemaître-Robertson-Walker spacetime, The evolution of the field in the cosmic time t is generated by the time-dependent Hamiltonian where k is the comoving momentum and π( k) = a(t) 3 dϕ( k)/dt. During slow-roll inflation, the Hubble rate changes slowly in time. To illustrate the analysis of the stability of the system, here we model this quasi-de Sitter phase with a de Sitter scale factor, The canonical subsystem spanned by ϕ( k), π(− k) with comoving momentum k is not a regular Hamiltonian system because of exponential collinearity (see appendix A.3). In fact the angle between the two Lyapunov vectors 1 ( k) and 2 ( k) approaches 0 as e −H 0 t . As a result, the Lyapunov exponents of the mode k do not have to be opposite in sign. In fact they are found to be given by where H 0 is the Hubble rate. In the quantum theory we consider an initial state at the time t 0 → −∞ given by the Bunch-Davies vacuum. The correlation functions of this state at the time t can be determined in closed form and are given by from which we can read the complex structure J t ( k, k ).
We analyze the entanglement growth of a subsystem spanned by a linear smearing of the field and the momentum. In order to guarantee that the dispersion of the linear observables are finite, we consider a smearing of the form is the comoving Laplacian and the Gaussian smearing is over a region of comoving size R. The eigenvalues of the restricted complex structure [iJ t ] A come in pairs ±ν(t) and are given by The entanglement entropy S A (t) = S(ν(t)) of the subsystem is plotted in figure 8 and for long time, i.e. for large number of efoldings, grows linearly as S A (t) ∼ H 0 t. This is exactly . Quantum field in de Sitter space. We plot the entanglement entropy S A (t) = S(ν(t)) from (4.37) associated to the subsystem described in (4.36). In the phase of linear growth, the entropy is observed to grow with rate given by the Hubble rate H 0 as predicted by theorem 1. We set H 0 = 1 and R = 1 in this plot. the asymptotic growth predicted by theorem 1 written in terms of the subsystem exponent Λ A = H 0 . 3 We note that previous studies of the growth of the entanglement entropy of cosmological perturbations focus on the ( k, − k) subsystem [85][86][87][88]. On the other hand the results presented here apply to all subsystems defined by smeared fields.

Proof, part I: classical ingredients
In this section and the subsequent section we collect and prove results used in the proof of the main theorem presented in section 2.
We consider a classical dynamical system with N degrees of freedom. We assume that the system has a Hamiltonian dynamics defined in a linear phase space [89]. We also restrict attention to quadratic time-dependent Hamiltonians. In this case we discuss the notions of stability, Lyapunov exponents and the growth of the volume of subsystems [32,33].

Linear phase space and quadratic time-dependent Hamiltonians
We consider a system with N degrees of freedom described by a linear phase space V = R 2N . Phase space observables O are smooth functions of 2N real variables denoted ξ a , The space of observables is equipped with a Lie algebra structure defined by the Poisson brackets where Ω ab is a nondegenerate antisymmetric matrix. In this paper we mostly focus on linear observables v = v a ξ a and quadratic observables O = 1 2 h ab ξ a ξ b . We call V * the vector space formed by all linear observables, and denote by v a the elements of V * and by w a the elements of V . The restriction of the Poisson brackets to the space of linear observables is A Darboux basis 4 (also called symplectic basis) of phase space functions consists of a set The notions of symplectic vector space and symplectic transformations play a central role in the description of the system. A symplectic structure on V is an antisymmetric nondegenerate bilinear map ω ab : V × V → R. It provides a canonical map from V to V * given by v a = ω ab v b . The couple (V, ω ab ) defines a symplectic vector space. The inverse of the symplectic structure, denoted Ω ab , is the antisymmetric bilinear map defined by Ω ac ω cb = δ a b and is a symplectic structure on V * . The space V * of linear observables on a linear phase space, equipped with the bilinear map Ω ab describing the restriction of the Poisson brackets to V * as in (5.3), is a symplectic vector space. In a Darboux basis, the symplectic structure Ω ab and its inverse ω ab take the 2N × 2N matrix form Ω = (Ω ab ) and ω = (ω ab ), Technically, DV is a basis of the dual phase space V * , but we refrained from bloating our notation by writing DV * . All Darboux bases D in this paper will live in the dual phase space V * . can be read as a concrete representation of the Kronecker delta δ a b = ξ a c ξ c b = ξ a b when we read both indices as abstract indices. However, when referring to an explicit basis (q1, . . . , qN , p1, . . . , pN ), we will use lower indices to match standard conventions.

The linear symplectic group Sp
In matrix form we have M ΩM = Ω. Note that the inverse of a symplectic matrix is given by M −1 = ΩM ω. The matrices M a b can be interpreted as linear maps either on V or on V * , and preserve the corresponding symplectic structures.
The dynamics of a Hamiltonian system is prescribed by a Hamilton function H(t) that we allow to be time-dependent. The Hamilton equations of motion of an observable O arė In particular, for the linear observables ξ a we havė In this paper we focus on time-dependent quadratic Hamiltonians, i.e. phase space functions of the form In this case the Hamilton equations simplify to the linear equatioṅ where the matrix K a b (t) is defined in terms of the quadratic term in the Hamiltonian by The solution of this equation provides the time evolution of the linear observable ξ a , The matrix M a b (t) solves the differential equationṀ a b (t) = K a c (t)M c b (t) with the identity as initial condition, and can be expressed as a time-ordered exponential, As the time evolution preserves the Poisson brackets, the matrix M a b (t) belongs to the linear and it vanishes if the linear term f a (t)ξ a is not present in the Hamiltonian.
A simple example of time-dependent quadratic Hamiltonian of the form (5.8) is given by a system of coupled harmonic oscillators with time-dependent couplings or driven by external forces. Another important example arises in the description of the lowest-order expansion of the evolution of a time-independent non-linear system around a classical solution chosen as background. In this case the time-dependence of the effective Hamiltonian arises from the background classical solution.

Linear stability and Lyapunov exponents
To characterize the linear stability of a dynamical system we consider a small perturbation δξ a (t) of a classical solution ξ a 0 (t) that satisfies the Hamilton equations. Substituting ξ a (t) = ξ a 0 (t) + δξ a (t) into (5.7) and expanding at linear order in the perturbation we find the linear equation with K a b (t) ≡ Ω ac ∂ c ∂ b H(t)| ξ 0 the stability matrix of the classical solution ξ a 0 (t). For the quadratic Hamiltonian (5.8) the stability matrix is simply given by K a b (t) = Ω ac h cb (t). As a result, the time evolution of the perturbation is given by where the symplectic matrix M (t) is given by (5.12). In order to measure the separation of two configurations in phase space we introduce a metric g ab , i.e. a positive definite symmetric bilinear, and define the norm ||δξ|| ≡ g ab δξ a δξ b . The exponential rate of separation of two sufficiently close classical solutions is given by the Lyapunov exponent λ δξ defined as We note that the Lyapunov exponent λ δξ is independent from the choice of metric g ab used to measure the distance between the classical trajectories ξ a 0 (t) and ξ a 0 (t) + δξ a (t). See appendix A.2 for a proof of this statement.
It is also useful to define Lyapunov exponents of linear observables (δξ) = a δξ a that probe the perturbation δξ a (t) and live in the dual phase space a ∈ V * . From the time evolution equation (5.15), we can read off that (t) evolves as In order to define a norm = G ab a b , we use the inverse metric G ab , such that G ac g cb = δ a b . Again, the Lyapunov exponent will be independent of the metric that we choose. The metric g ab , used above to define Lyapunov exponents, is said to be compatible with the symplectic structure ω ab with inverse Ω ab if the matrix J a b ≡ Ω ac g cb is symplectic, J a c J b d Ω cd = Ω ab , and squares to minus the identity J a c J c b = −δ a b . In this case, J a b defines a complex structure. The inverse metric G ab is then compatible with the symplectic structure Ω ab in the dual space. A compatible metric g ab allows us to define the limiting matrix L a b , that characterizes the long-time stability of the system. 6 Provided that the Hamiltonian system is regular in the sense of appendix A.3, Lyapunov exponents exist for all linear observables a and are given by the eigenvalues of the limiting matrix L a b . As the matrix L a b is symmetric and belongs to the symplectic algebra sp(2N ), its eigenvalues are real and come in pairs with opposite sign (λ, −λ). The Lyapunov spectrum consists of the ordered Lyapunov exponents given by Note that D L is not unique because it depends on our choice of metric G ab , but subsequent results will be independent of this choice [90,91].
We will discuss examples of unstable quadratic systems in section 2. The prototypical example is the inverted harmonic oscillator, with a potential V that is unbounded from below. The Lyapunov exponents of the system are related to the unstable directions of the potential. Another example is provided by periodically driven systems, i.e., systems with a quadratic time-dependent Hamiltonian of the form (5.8) with periodic coefficients h ab (t + T ) = h ab (t). In this case instabilities appear due to the phenomenon of parametric resonance [89]. The real part of the Floquet exponents of the system coincide with the notion of Lyapunov exponents described above.

Subsystems and the subsystem exponent Λ A
The partition of a Hamiltonian system in two complementary subsystems corresponds to a decomposition of the phase space V and its dual V * into direct sums This set of observables provides us with a Darboux basis of linear observables A * that only probe the degrees of freedom in A and it can be completed to a Darboux basis of the full system by introducing a Darboux basis of B * , Given a Darboux basis θ r of A * and its dual basis ϑ r of A satisfying θ r a ϑ a s = δ r s , we can restrict tensors to the subsystem by appropriate contractions. Most importantly, we will consider the restriction [J] A of the complex structure J a b and [G] A of a metric G ab : Note that, as θ r a is a Darboux basis, the restriction of the symplectic structure Ω ab is still a symplectic structure, [Ω] A = Ω A . On the other hand, the restriction [J] A of a complex structure J is not in general again a complex structure, because it does not necessarily satisfy We give some examples of subsystems. Consider for instance a linear chain of N oscillators, with the oscillator at site i having canonical coordinates (q i , p i ). A first example of subsystem A corresponds to the subset of observables (q i , p i ) with i = 1, . . . , N A associated to a geometric decomposition of the chain in two complementary intervals. A second example of subsystem is given by a subset of normal-mode observables (φ k ,π k ) with k = 1, . . . , N A corresponding to the long-wavelength perturbations of the system. A third example is provided by a detector that makes measurement of the localized observables (Q, P ) only, with p i probing only average properties of a localized subset of oscillators. Each example shows that the choice of a subsystem A corresponds to a coarse graining of the system that preserves the symplectic structure of the accessible observables.
Given a subsystem A we introduce a new notion of characteristic exponent Λ A that generalizes the notion of Lyapunov exponents of the system. A Darboux basis D A = (θ 1 , . . . , θ 2N A ) of the subsystem defines a symplectic cube Given the metric G ab , we can compute the volume Vol G (V A ) of the symplectic cube V A as the square root of the determinant of the 2N A × 2N A Gramian matrix (θ r a G ab θ s b ), We will be interested in how this volume changes under time evolution. Let us recall that the action on V * is given by the transpose M (t) a b = M b a (t). If we evolve the symplectic cube V A with M , we have We define the subsystem exponent Λ A as the exponential rate of growth of the volume of the subsystem measured with respect to the metric G ab , For a regular Hamiltonian system this limit exists and is independent of the metric g ab , see appendix A.3. Note that the exponent of the full system Λ V vanishes because [M (t)] V = M (t) and the determinant of a symplectic matrix is equal to one. The vanishing of Λ V is a special case of the Liouville theorem. We say that the subsystem A is unstable under time evolution if it has a positive exponent Λ A .

Relation of the exponent Λ A to Lyapunov exponents
We now show how to compute the subsystem exponent Λ A once the Lyapunov spectrum of the system is known. The result is stated and proven below.

31)
We refer to the 2N columns of T as t i .
3. Find the first 2N A linearly independent 7 columns t i of T which we can label by t i k with k ranging from 1 to The subsystem exponent Λ A is then given by the sum over the 2N A Lyapunov exponents λ i k , where the index i k is defined above. 7 Here, we mean that ti cannot be expressed as a linear combination of the vectors ( t1, . . . , ti−1) standing to the left in the matrix T .
Proof. The rectangular matrix T in (5.31) allows us to express the elements of the Darboux basis D A of the subsystem in terms of the Lyapunov basis, θ r = 2N T r i i . Denoting the columns of T by t i we can select the first 2N A linearly independent columns in the ordered set ( t 1 , . . . , t 2N ). We label them t i k and organize them in the 2N A × 2N A square matrix U , Due to their linear independence, the inverse U −1 exists and turns T into an upper triangular matrixT of the form 0 · · · 0 1 * * · · · · · · · · · · · · · · · · · · * 0 · · · · · · · · · 0 1 * * · · · · · · · · · · · · * . . . . . . . . . . . . 0 · · · · · · · · · · · · · · · · · · 0 1 * * · · · * where the * represents an unspecified value. Acting with U −1 on the left and the right-hand side of (5.31) we findθ Moreover the 2N A vectorsθ k are linearly independent and form a (generally non-symplectic) basis of A * . Therefore the cube M (t) V A is given by a time-independent linear transformation of the one spanned by M (t)θ k . In the limit t → ∞ the volume of the subsystem scales as if there are no directions that become collinear in an exponentially fast way under time evolution. As the exponential collinearity is excluded by the assumption of regularity (see appendix A.3), the subsystem exponent is given by (5.32).
Corollary. The subsystem exponent is non-negative, Proof. Let us denote the elements of the Lyapunov basis by i = Q i for i ≤ N and i = P 2N +1−i for i > N so that ( 1 , . . . , 2N ) = (Q 1 , . . . , Q N , P N , . . . , P 1 ). Each vectorṽ i k with i k > N consists of a linear superposition of momenta P i only, as follows from (5.35). As a result, to span a symplectic subspace, for each suchṽ i k there has to be aṽ i k with i k ≤ 2N + 1 − i k so to contain the conjugate position Q i in the linear superposition. Therefore, negative Lyapunov exponents λ i k with i k > N are paired with positive Lyapunov exponents λ i k , resulting in a sum of non-negative terms λ i k + λ i k ≥ 0 in (5.32). subsystem phase space Figure 9. Subsystem exponents due to phase space stretching. We illustrate the statement of theorem 3. We start with a symplectic cube V A ⊂ A * in the subsystem and time-evolve it to the deformed cube M (t)V A that is dominantly stretched into the 2N A directions of M (t) i k with Lyapunov exponents λ i k . Consequently, the logarithm of its metric volume behaves as log Vol( In generic situations, λ i k are just the 2N A largest Lyapunov exponents, as explained in theorem 4. The quantity log Vol(M (t)V A ) is related to the entanglement entropy S A as explained in section 6.2.
We illustrate this result with some examples of subsystems and the associated exponents. Consider a system with N = 2 degrees of freedom, Lyapunov spectrum (λ 1 , λ 2 , −λ 2 , −λ 1 ) (5. 37) and Lyapunov basis D L = ( 1 , 2 , 3 , 4 ) = (Q 1 , Q 2 , P 2 , P 1 ). A subsystem A with N A = 1 degree of freedom can be identified by specifying a canonical couple (φ, π). Here we give three examples: In particular the subsystem given by a single couples (Q i , P i ) has vanishing subsystem exponent Λ A = 0. Note also that the difference of positive Lyapunov exponents can appear as in example (5.39). We will reconsider these examples in section 3.1 and relate the subsystem exponents Λ A to the production of entanglement entropy.

Relation of the exponent Λ A to the Kolmogorov-Sinai entropy rate
From an information-theory perspective, the Hamiltonian evolution of a dynamical system with sensitive dependence on initial conditions produces entropy. This is because two initial conditions that are indistinguishable at a fixed resolution will evolve into distinguishable states after a finite time. The Kolmogorov-Sinai entropy rate provides a quantitative characterization of this behavior: It measures the uncertainty remaining on the next state of a system, if an infinitely long past is known [30][31][32][33]. It is defined as follows.
We decompose the phase space V into cells (C 1 , . . . , C n ) belonging to a partition P. Given a sampling time ∆t, we can compute the probability µ(C 1 , . . . , C n ) that a trajectory starting in cell C 1 will successively go through C 2 , C 3 and so on. The entropy per unit time with respect to such a given partition is given by Shannon's formula, The quantity h KS is a global invariant of the system and it provides a quantitative characterization of the notion of deterministic chaos in a Hamiltonian system. A positive Lyapunov exponent corresponds to the exponential divergence in time of some initially nearby trajectories. This phenomenon results in the unpredictability of the evolution at finite resolution, and therefore contributes to h KS . Pesin's theorem [90,92] states that, for Hamiltonian dynamical systems, the Kolmogorov-Sinai entropy rate is equal to the sum over all the positive Lyapunov exponents of the system. Let us call N I ≤ N the number of non-vanishing positive Lyapunov exponents. Using the ordering (5.20) of the Lyapunov spectrum, we have This formula, together with (5.32), clearly shows that the Kolmogorov-Sinai entropy rate provides an upper bound to the characteristic exponent Λ A of a subsystem, In the following we discuss when this inequality is saturated and show that, for a large class of system decompositions, the characteristic exponent Λ A equals the rate h KS . The following theorem is instrumental.
Theorem 4 (Subsystem exponent -generic subsystem). The subsystem exponent of a generic subsystem A of dimension N A is given by the sum of the first 2N A Lyapunov exponents, a special case of (5.32),

45)
This behavior holds for all subsystems A ∈ V , except for a set of measure zero.
Proof. The space of 2N A -dimensional symplectic subspaces of V has the structure of a differentiable manifold and is called the symplectic Grassmannian SpGr(2N A , V ). Let us consider the set of points on this manifold where the generic asymptotics (5.45) does not apply. The statement of the theorem is that this set forms a lower dimensional submanifold. All standard measures on differentiable manifolds will therefore assign a measure zero to this subset. By applying theorem 3, (5.32), we find Λ A generic = 2N A i=1 λ i whenever the first 2N A columns of the transfer matrix T are linearly independent. Let us therefore analyze for how many system decompositions this does not hold. The space of 2N A -dimensional symplectic subspaces SpGr(2N A , V ) can be identified with the space of transformation matrices such that the restricted symplectic form [Ω] A is non-degenerate, modulo GL(2N A ): Let us now compare this to the space of subspaces for which the subsystem exponent Λ A is not given by the sum over the first 2N A largest Lyapunov exponents. For this to happen, it is a necessary condition that the first 2N A columns of the transfer matrix are linearly dependent. This space is (4N N A −1)-dimensional which we still need to quotient by GL(2N A ). Therefore, the subset of spaces of subsystems A, for which we find Λ A = Λ A generic , has a dimension of at most 4N A (N − N A ) − 1. This is a set of measure zero with respect to any standard measure on SpGr(2N A , V ) because it lies in a lower dimensional submanifold.
This behavior was conjectured by Zurek and Paz in [4] and later discussed by Asplund and Berenstein [35] and ourselves [36].
Of the three examples discussed at the end of section 5.3, only the one-dimensional subsystem (φ, π) with Λ A = λ 1 + λ 2 is generic, (5.38). Note that most numerical algorithms for the computation of the Lyapunov exponents of a dynamical system start with the computation of the exponential rate of expansion of the volume of a subsystem [93]. Lyapunov exponents are computed by taking the difference between the exponential rate of expansion of subsystems of different dimension. The efficiency of these algorithms relies on the generic behavior discussed above. Now we investigate when the subsystem exponent equals the rate h KS assuming that the subsystem is generic, (5.45).
In a stable Hamiltonian system, all Lyapunov exponents vanish. The system becomes unstable as soon as a single Lyapunov exponent turns positive. We call N I the number of non-vanishing positive Lyapunov exponents. Pesin's formula for the Kolmogorov-Sinai entropy rate then reads h KS = N I i=1 λ i . On the other hand the characteristic exponent of a generic subsystem A of dimension N A in the range N I ≤ 2N A ≤ 2N − N I is given by Λ A generic = N I i=1 λ i . Therefore, we have the equality that identifies subsystems that saturates the inequality (5.44). Note that unstable many-body systems often have a number of unstable directions N I that is much smaller that the number of degrees of freedom of the system, N I N . A generic subsystem that encompasses a fraction f = N A /N of the full system satisfies (5.47) if the fraction is in the range In particular, in the limit N → ∞ with N I finite, we have Λ A = h KS for all partitions of the system into two complementary subsystems each spanning a finite fraction f of the system, except for a set of partitions of measure zero.

Proof, part II: quantum ingredients
The main result, theorem 1, is proven in three steps which heavily rely on the three ingredient presented in the following subsections. We show how the Renyi entropy provides bounds for the entanglement entropy, we explain how the Renyi entropy can be understood as the logarithm of the volume of a region in the dual phase space and finally, we derive the time evolution of the Renyi entropy as the volume deformation of this region under the classical symplectic flow.

Upper and lower bounds on the entanglement entropy
We recall that there are different entanglement measures that quantify the amount of correlations in a state |ψ with respect to some system decomposition into subsystems A and B.
Beside the entanglement entropy S A (|ψ ), we have the class of Renyi entropies defined by . For a Gaussian state |J, η labeled by a complex structure J, all these entropies can be computed directly from the eigenvalues ±iν i of [J] A , the complex structure restricted to the subsystem A. If we take the positive value ν i of each eigenvalue pair, the Renyi entropy 8 R A = R (2) A and the entanglement entropy S A are given by

2)
8 From now on, we will refer to the Renyi entropy of order 2 as the Renyi entropy. holding for ν ≥ 1 as shown in figure 10. An immediate consequence of this inequality is that the entanglement entropy of a Gaussian state is bounded from below by the Rényi entropy and from above by the Rényi entropy plus a state-independent constant, This means that the Renyi entropy R A determines a corridor for the entanglement entropy S A . This implies immediately that both of them will grow asymptotically with the same rate which we use for the main result of this paper.

Renyi entropy as phase space volume
A Gaussian state |J, η equips the dual phase space V * with a metric G ab defined by (B.9), which is really just the covariance matrix of the state. The complex structure J can be expressed in terms of the metric J a b = −G ac ω cb . (6.5) Furthermore, the restriction of the complex structure to the subsystem A can be written in matrix form as a product of the symplectic Ω A and the restriction of the metric, In a Darboux basis, where we have det ω A = 1 and det[G] A > 0, we find that the determinant of the restriction of the complex structure can be expressed in terms of the phase space volume Vol G (V A ) of a symplectic cube V A (spanned by a Darboux basis and with symplectic volume 1) measured with respect to the induced metric.
As a result, the Rényi entropy of a Gaussian state is given by the logarithm of the phase space volume of a symplectic cube D A defining the subsystem, measured with respect to the metric G ab defined by the state, Note that the symplectic cube V V associated to a Darboux basis of the full system satisfies Vol G (V V ) = 1 and therefore the Rényi entropy vanishes R V (|J, ζ ). On the other hand, the restriction to a subsystem A can result in a larger volume Vol G (V A ) ≥ 1 and a non-vanishing Rényi entropy.

Entanglement entropy growth as phase space stretching
Let us consider a one-parameter family of Gaussian states |J t , η t = U (t)|J 0 , η 0 generated under time evolution of some quadratic Hamiltonian. In particular, we have J t = M (t)J 0 M −1 (t) where M (t) : V → V is the classical Hamiltonian flow on phase space. We call G t the timedependent metric associated with J t , and G 0 the initial metric associated with J 0 . The evolution of the Rényi entropy is given by (6.8), where the volume is now measured with respect to the time-varying metric: In this formula, the symplectic cube V A is kept fixed while the metric evolves. However, the same volume is obtained if we let the symplectic cube evolve according to M (t)V A for a fixed metric G 0 . Hence, we can compute The symplectic basis of the subsystem A is stretched by the Hamiltonian flow M (t) : V * → V * on the dual phase space, and the variation in its volume determines the evolution of the Rényi entropy.
Since the absolute difference between the entanglement entropy and the Rényi entropy of a Gaussian states is bounded by a state independent constant, we have that: i.e., the asymptotic rate of growth of the entanglement entropy and of the Rényi entropy coincide. This allows us to compute the asymptotic rate of growth of the entanglement entropy from (6.10) as: 12) in terms of the stretching of the symplectic cube under time-evolution.

Discussion
We discuss the role of interactions in the saturation phase of the entanglement growth, explain the relation to results on quantum quenches, present a conjecture on entanglement and chaos, and summarize our results.

Interactions and the production of non-Gaussianities
Quadratic Hamiltonians appear naturally in the analysis of small perturbations around equilibrium configurations, both stable and unstable. Let us consider for instance a dynamical system with a "Mexican hat" potential and an initial Gaussian state which is sufficiently peaked at the top of the potential. For short times the evolution of this initial state is well described by a perturbative quadratic Hamiltonian with unstable directions as in (3.7). As a result, if the scales of the problem are sufficiently separated, the entanglement entropy of a subsystem will show an intermediate linear growth with rate Λ A , followed by a non-Gaussian phase. In particular the linear growth driven by the perturbative instability stops when the spread of the state starts to probe the bottom of the potential and interactions become nonnegligible. As the full Hamiltonian of the system is stable at the non-perturbative level, the entanglement entropy of the subsystem is bounded from above by the entanglement entropy S eq of the thermal state with the same energy as the initial Gaussian state. In the presence of an equilibration and a thermalization mechanism, the entropy S eq provides also the saturation value as shown in figure 1. Quadratic Hamiltonians appear also in the analysis of small perturbations of classical solutions. In this case the perturbative Hamiltonian inherits the time-dependence of the classical solution. For instance if the classical solution is periodic in time, then it provides a time-dependent background for the perturbations which leads to a perturbative Hamiltonian that is periodic in time. The stroboscopic dynamics of the system can be analyzed with the same methods discussed for the Hamiltonian (3.9). In particular, in the presence of parametric resonances, the Floquet exponents of the system determine the subsystem exponent Λ A and the growth of the entanglement entropy as described in theorems 1 and 2. When the conditions (2.13) for the subsystem are satisfied, the rate of growth is given by the classical Kolmogorov-Sinai entropy rate as discussed also in [35]. After the initial phase of linear growth, two distinct phenomena render the parametric resonance inefficient and lead to a saturation phase. The first phenomenon is dephasing: large perturbations are not harmonic; their period depends on the amplitude of the oscillation and, when the period is driven far from resonance, the periodic background cannot pump energy efficiently into the perturbation. The second phenomenon is backreaction: clearly the linear entropy growth is accompanied by the production of a large number of excitations that at some point start to interact and backreact, thus leading to a saturation phase in which non-Gaussianities cannot be neglected. The phase of linear growth manifests itself only if the typical scales of the problem are sufficiently separated. A preliminary numerical investigation of the effect of interactions and non-Gaussianities on the entanglement growth can be found in [41].

Relation to linear growth in quantum quenches
Quantum quenches lead also to a phase of linear growth of the entanglement entropy, followed by a saturation phase. This phenomenon has been studied extensively in free field theories [10][11][12] and in many-body quantum systems [9,[94][95][96][97][98][99]. Despite the similarities in the behavior of the entanglement entropy, the mechanism behind this phenomenon is distinct from the one discussed in this article.
A standard example of global quantum quench is provided by an harmonic lattice similar to the one discussed in section 3.3. The Hamiltonian of the system consists of two terms, the first term is ultralocal, while the second encodes the coupling of first neighboring oscillators. A quantum quench consists in preparing the system in the ground state |ψ 0 of the Hamiltonian H 0 with vanishing coupling κ = 0. At the time t = 0 the coupling is instantaneously switched on, and the state is let to evolve unitarily, |ψ(t) = e iHκ t |ψ 0 . The entanglement entropy of a local subset of the lattice evolves in a way similar to the one depicted in Fig. 1. Instabilities play no role in this phenomenon. In fact the Hamiltonian (7.1) is stable for κ > −Ω 2 0 /4, as it can be seen from Eq. (3.14). The relevant mechanism for the linear growth of the entanglement entropy in this quantum quench is not instabilities, but transport instead. The quench results in the local production of quasi-particles which travel at a finite speed. The phase of linear growth can be understood to be the result of the entanglement produced by the free propagation of entangled couples of quasi-particles, with the entanglement production rate determined by their propagation speed [10][11][12]99].
Interactions or coupling between many degrees of freedom, together with propagation of quasi-particles, play a key role in the phenomenon of entanglement growth in quantum quenches. On the other hand, the phenomenon studied in this paper relies on the existence of instabilities of some modes of a quantum system, as discussed in section 3 and 4. The difference between the two phenomena is easily illustrated by the case of bosonic systems and Gaussian states for which formula (B.44) holds, [100][101][102][103] In the case of quantum quenches, the number of entangled pairs N e with fixed weight ν i grows linearly in time until saturation. On the other hand, in the presence of instabilities, unstable modes have weight ν i which grows exponentially in time until saturation, therefore leading to an entanglement growth of the form depicted in Fig. 1. While the linear regime for quantum quenches can only be seen for a sufficiently large number of degrees of freedom in the system, linear growth due to instabilities can already occur for a system with two degrees of freedom and a single instability [41]. As a result, despite the intriguing similarity, the two phenomena are distinct.

A conjecture on entanglement, chaos and thermalization times
There is an intimate relationship between chaos, thermalization and entanglement [1][2][3][104][105][106]. Here we discuss how semiclassical methods allow us to estimate the rate of growth of the entanglement entropy in the early phase of the thermalization process. Let us consider a classical Hamiltonian system with linear phase space (R 2N , Ω) and a Hamiltonian H which does not depend on time so that, as a result, the energy of the system is conserved. We assume also that the Hamiltonian is bounded from below and, at fixed energy, trajectories in phase space are bounded. This classical system displays a chaotic behavior if its Kolmogorov-Sinai entropy rate is non-vanishing, i.e. h KS > 0 with h KS defined by (5.42), [107][108][109]. We are interested in the process of thermalization in the associated quantum system with Hamiltonian H. We argue that the Kolmogorov-Sinai entropy rate h KS studied in this paper plays a central role in determining the relevant time scale in the process of quantum thermalization. An isolated quantum system thermalizes when observables that probe only part of the system cannot distinguish a pure state from a thermal state. More precisely, let us consider a pure state and a thermal state with the same energy, 3) The requirement that they have the same energy fixes the temperature β −1 of the thermal state, i.e., E = ψ t |H|ψ t = Tr(Hσ). Now we consider a subsystem A and the subalgebra of bounded observables 9 O A in A. We say that the subsystem A thermalizes if all bounded observables O A attain a thermal expectation value, i.e., This condition can be formulated in terms of entanglement between the subsystem A and its complement B. Let us define the restricted states 10 Thermalization in A is a measure of how distinguishable is the restricted states ρ A (t) from the restricted thermal state σ A . The relative entropy [39,40], provides a measure of such distinguishability. In fact, using the inequalities S(ρ σ) ≥ 1 2 ρ − σ 2 together with the Schwarz inequality σ ≥ Tr Oσ / O , one can prove the inequality which holds for all bounded observables in A. Therefore, proving S(ρ A (t) σ A ) → 0 as t → ∞ provides a proof of thermalization in A. Now, the relative entropy can be expressed in turn as the sum of two terms, The first term is the difference between the equilibrium entropy S eq A (E) = −Tr A (σ A log σ A ) and the entanglement entropy S A (t) of the subsystem. The second term measures energy flow between the subsystem A and its complement, as measured by the effective HamiltonianH A of the subsystem defined in Eq. (7.5). At equilibrium, both terms vanish independently. This paper and the following conjecture focus on the evolution of the first term, i.e., the growth and saturation of the entanglement entropy S A (t).
When a subsystem thermalizes, the entanglement entropy S A (t) approaches the equilibrium value S eq A (E). The eigenstate thermalization hypothesis (ETH) [104][105][106] provides a sufficient condition for such subsystem thermalization to occur. In a chaotic quantum system with local interactions one observes that energy eigenstate, H|E n = E n |E n , in the bulk of the energy spectrum have a non-trivial entanglement structure: their restriction to a local subsystem results in a thermal state of the form of Eq. (7.5), i.e., Tr B (|E n E n |) ≈ σ A (E n ). As a result, the restriction ρ A (t) of a pure state |ψ t = n e −iEnt c n |E n with support in a narrow band of energy E is also well approximated by a thermal state when averaged over time, i.e. 1 T T 0 ρ A (t)dt ≈ σ A (E) for large T . This condition is sufficient to prove thermalization in average, but it does not provide a time-scale for the thermalization process.
We propose a conjecture which complements previous arguments to the quantum thermalization of subsystems [1][2][3][104][105][106]. The conjecture applies to semiclassical states and provides a time-scale for subsystem thermalization: Given an initial state |ψ 0 = n c n |E n peaked on a classical configuration of energy E = ψ 0 |H|ψ 0 with E large compared to the energy gap of the system, and a local subsystem A such that its initial entanglement entropy is small compared to the thermal entropy at the same energy, S A (|ψ 0 ) S eq A (E), the time-evolution of the entropy S A (t) ≡ S A (e −iHt |ψ 0 ) displays a linear phase S A (t) ∼ Λ A (E) t before saturating to the plateau at S eq A (E) as described in figure 1. The rate Λ A (E) can be computed from the classical chaotic dynamics of the Hamiltonian H on the energy-shell E. The rate is given by the subsystem exponent discussed in section 5.3 and, apart from its energy, it is largely independent of the initial state |ψ 0 . In particular, the subsystem exponent sets the time-scale of subsystem thermalization, The conjecture is based on theorem 1 presented in section 2, together with semiclassical arguments. Let us consider a classical solution (q cl i (t), p cl i (t)) with energy H(q cl i (t), p cl i (t)) = E. At the leading order in a semiclassical expansion, the evolution of a perturbation ξ a = (q cl i (t) + δq i , p cl i (t) + δp i ) of the classical solution is governed by the perturbative Hamiltonian where δξ a = (δq i , δp i ) and . (7.10) The Lyapunov exponents of a non-perturbative chaotic system with Hamiltonian H can be computed directly from the perturbative Hamiltonian H pert (t), which is quadratic timedependent and for which our theorem 2 applies. In fact, because of ergodicity of a chaotic system, all trajectories on the same energy-shell E (except a set of measure zero) have the same Lyapunov exponents. 11 Moreover, under standard assumptions of regularity [90], the Kolmogorov-Sinai rate h KS (E) on the shell of energy E is given by Pesin's formula (5.43) in terms of the positive Lyapunov exponents λ i (E). We consider now a symplectic subsystem (A, Ω A ) and define its subsystem exponent Λ A (E) as in section 5.3. This is also the rate of growth of the entanglement entropy derived assuming a Gaussian state and a quadratic Hamiltonian in theorem 3. The conjecture extends this result to a full non-quadratic system with bounded and chaotic motion, within the regime of validity of the semiclassical expansion. The inequality Λ A (E) ≤ h KS (E) provides an upper bound on the rate of entanglement growth during the linear phase. Clearly, the linear phase ends when the semiclassical approximation breaks down, i.e. when the spread of the wavefunction is so large that higher-order terms in the expansion H = E + H pert (t) + . . . cannot be neglected. An estimate of this time is provided by τ eq ∼ S eq A (E)/h KS (E) which measures the ratio between the accessible volume in phase space and the rate of growth of the phase space volume occupied by the perturbation.
The conjectured behavior of the entanglement entropy S A (t) depicted in figure 1 is expected to manifest itself only in the regime where the semiclassical approximation holds. This conjecture can be tested on a model system such as the one described by the Hamiltonian This is a well-studied model which appears in the analysis of the homogeneous sector of Yang-Mills gauge theory [110,111]. Its Lyapunov exponents are known to scale with the energy as λ i (E) ∼ E 1 4 and its equilibrium entropy, estimated as the log of the phase space volume at fixed energy, scales as S eq (E) ∼ log E. As a result, for a semiclassical initial state of energy E we expect our conjecture to apply: the entanglement entropy of a subsystems such as (x, p x ) is expected to initially grow linearly with a rate ∼ E 1 4 and then saturate in a time τ eq ∼ E −1/4 log E. This behavior can in principle be tested via numerical investigations. The numerical analysis involves the unitary evolution of a pure state under a chaotic quantum Hamiltonian, which is beyond the scope of the present paper. A preliminary numerical analysis of the growth of the entanglement entropy in interacting systems prepared in a semiclassical state can be found in [41].
We note that the conjecture is expected to apply only to initial states which are semiclassical, i.e. states with average energy much larger than the energy gap, small spread in energy and, in general, small spread around a point in phase space. On the other hand, when the energy E of the initial state is comparable to the energy gap of the Hamiltonian, classical orbits of that energy have an action comparable to and there is no reason to expect that they provide a useful tool for predicting the behavior of the entanglement entropy in the linear regime of figure 1. In fact, recent results from quantum field theories with a gravity duals [112] show that -at low energy -the rate of growth of the entanglement entropy is bounded from above by the energy of the subsystem divided by and therefore deviates from the semiclassical prediction [113].

Summary
We studied the relationship between entropy production in classical dynamical systems and the growth of the entanglement entropy in their quantum analogue in the semiclassical regime. Most importantly, we found that in both cases the production rates are given by a sum over Lyapunov exponents λ i characterizing stable and unstable phase space directions. For classical systems, there is a standard notion of rate, the Kolmogorov-Sinai entropy rate given by the sum over all positive Lyapunov exponents. We have shown that for the associated quantum system and a subsystem A, the production rate of the entanglement entropy S A (t) ∼ Λ A t is given by a subsystem exponent Λ A . We have shown that Λ A ≤ h KS and found that this inequality is saturated for sufficiently large subsystems. Moreover we found that the rate Λ A is independent of the initial state of the system and -except for a set of measure zero of subsystems -it depends on the choice of subsystem A only via its classical dimension N A , i.e.
where λ i are the 2N A largest Lyapunov exponents of the system. Our rigorous derivation of this result is based on the assumption of unstable quadratic Hamiltonian and Gaussian initial state. The derivation takes into account the case of time-dependent Hamiltonians with Floquet instabilities. The derivation of the main theorem proving S A (t) ∼ Λ A t consists of three parts. First, the subsystem exponent Λ A is introduced at the classical level as a generalization of Lyapunov exponents and defined to encode the exponential rate of growth of the volume of a symplectic cube in a subsystem under Hamiltonian evolution. Second, the time evolution of a Gaussian initial state through an unstable quadratic Hamiltonian is conveniently encoded in terms of complex structures or equivalently phase space metrics and their classical Hamiltonian flow. Third, the evolution of entanglement entropy is shown to be asymptotically the same as the one of the Renyi entropy which can then be shown to grow with the rate of the subsystem exponent Λ A . We interpret the exponent Λ A as a quantum analogue of the Kolmogorov-Sinai entropy rate of a given subsystem A.
The predicted linear growth of the entanglement entropy shows up in a wide range of physical systems such as unstable quadratic potentials, periodic quantum quenches in manybody quantum systems and instabilities in quantum field theory models. We presented three examples of the latter where entanglement is produced through different mechanisms, namely unstable modes due to a symmetry-breaking instability, parametric resonance in models of post-inflationary reheating, and cosmological perturbations in an inflationary spacetime.
We believe that our results are also relevant in the context of thermalization of isolated quantum systems. A subsystem of a chaotic quantum system is expected to thermalize with equilibrium entropy S eq (E) determined by the average energy E of the initial state. In the semiclassical regime we conjecture that the time-scale of this equilibration process is τ eq ∼ S eq (E)/Λ A (E) where Λ A (E) is the subsystem exponent of the energy-shell E.
We thank also Marcos Rigol and Ranjan Modak for discussions on extensions of the presented results to the case of non-Gaussian states. EB thanks Berndt Müller for inspiring conversations on the Kolmogorov-Sinai entropy rate. EB acknowledges also extensive discussions with Renaud Parentani and Bei-Lok Hu which took place during the 22 nd Peyresq Physics workshop. LH thanks Pavlo Bulanchuk for a suggestion leading to the geometric representation of the Rényi entropy. We thank Hal Haggard, Carlo Rovelli and Matteo Smerlak for multiple discussions and feedback at various stages of this project. The work of EB is supported by the NSF grant PHY-1404204. LH is supported by a Frymoyer fellowship. NY acknowledges support from CNPq, Brazil and from the NSF grant PHY-1505411. This research was supported in part by the Perimeter Institute for Theoretical Physics.

A Dynamical systems and Lyapunov exponents
We summarize relevant properties of Lyapunov exponents in the context of Hamiltonian systems. In particular, we make precise the notion of regular Lyapunov system.

A.1 Linear Hamiltonian systems
We consider a 2N -dimensional linear phase space V with symplectic form ω. A timedependent Hamiltonian H is a smooth map The equations of motion are given bẏ where Ω ab satisfies ω ac Ω bc = δ b a and (dH) b (t) is the gradient of H at time t. The solution of these equations can be conveniently described by a flow This map is a diffeomorphism that preserves the symplectic form ω, namely the push-forward satisfies (Φ t ) * ω = ω. For a given point ξ 0 ∈ V , the push-forward (Φ t ) * maps a tangent vector δξ ∈ T ξ 0 V to the tangent vector (Φ t ) * δξ ∈ T Φt(ξ 0 ) V . Due to the linearity of V , we can identify the tangent spaces at all points with V itself. Formally, we have an isomorphism φ ξ : V → T ξ V that maps v ∈ V to the tangent vector φ ξ v ∈ T ξ V that acts on a function f : V → R as φ ξ v(f ) = d dt f (ξ +tv). Using φ ξ , we can to define the linear map M ξ 0 (t) : that corresponds to the above push-forward after we identify T ξ 0 V and T (Φt)ξ 0 V with V .
In the special case, where the Hamiltonian H is given by a linear quadratic function for every t, the Hamiltonian flow is given by an inhomogeneous symplectic transformation (M (t), η(t)) via Φ t ξ 0 = M (t)ξ 0 + η(t) , (A.6) whose differential is given by M ξ 0 (t) = M (t), independent of ξ 0 . This follows from the fact that a quadratic Hamiltonian gives rise to linear and homogeneous equations of motion. The symplectic group element M (t) is formally given by the time-ordered exponential where the generator K(t) is an element of the symplectic Lie algebra sp(2N ).
In the general case, where H is not quadratic, we can still find the time-dependent generator K ξ 0 (t) a b = Ω ac h ξ 0 (t) cb . In this case, however, the generator also depends on the initial ξ 0 . To find h ξ 0 (t) cb , we just need to Taylor expand the Hamiltonian H(t) along the trajectory ξ(t) = Φ t ξ 0 which amounts to finding its Hessian At this point, we understand that the difference between the special (quadratic) and general case (arbitrary Hamiltonian) in regards of the linear map M ξ 0 (t) are the following: • If the Hamiltonian is quadratic or affine quadratic, the linear map M (t) describing the push-forward of the Hamiltonian flow Φ t is independent of the starting point ξ 0 and completely characterized by the quadratic part h(t) ab of H(t).
• For a more general Hamiltonian, we can still compute its quadratic part h ξ 0 (t) ab as Hessian of H(t) along the trajectory ξ(t). This means h ξ 0 (t) ab depends on the initial condition ξ 0 and the corresponding solution ξ(t) with ξ(t) = ξ 0 . In particular, the quadratic map M ξ 0 (t) will differ for different initial conditions ξ 0 .
The linear symplectic map M ξ 0 (t) contains all the information about how two sufficiently close trajectories converge or diverge. This behavior will be captured in the so called Lyapunov exponents.
In order to define Lyapunov exponents, we need to equip phase space V with a positive definite metric g ab that gives rise to a norm δξ = g ab δξ a δξ b . Equivalently, we can use the inverse metric G ab to define the norm = G ab a b on the dual phase space V * . We will show that Lyapunov exponents are actually independent of the specific choice of a positive metric. In order to show this, it is useful to have the following theorem at hand. Proposition 1. Given a finite dimensional, real vector space V and two distinct positive metrics g andg, we can compute the following two values which allow us to relate norms and angles measured by the different metrics: • Norm inequality Given a vector v ∈ V , its norm v g with respect tog is related to v g via: • Angle inequality Given an angleψ between two vectors measured with respect tog, it is related to the angle ψ measured with respect to g via the following inequality: This inequality can be simplified to the slightly weaker version given by: (A.12) • Volume inequality Given the d-volume Volg(V A ) of some region V A in an arbitrary d-dimensional subspace A ⊂ V measured by the metricg, it is related to the d-volume Vol g (V A ) measured by g via the following inequality: If consider the same equations for the dual phase space V * with the replacements g → G and g →G, all inequalities hold if we replace a → 1/b and b → 1/b.
Proof. Let us prove the different inequalities: • Norm inequality Let us take two different norms induced by the two positive metrics g andg. Over a finite dimensional vector space V the set S = {v ∈ V | v g = 1} is compact. This means that the continuous function v g will take a minimal and maximum value on S: (A.14) Linearity of the induced norm implies than the inequality that we wanted to prove: • Angle inequality Let us choose a two-dimensional plane P ⊂ V . On this plane, we have the restricted metrics g| P andg| P . The two are related by a linear map D : P → P with where D is not unique. We can always choose it to be diagonalizable with ordered eigenvalues d i and eigenvectors e i . At this point, we can identify the inner product with respect tog as the one with respect to g after having acted with D on the vectors. This implies a ≤ d i ≤ b to not violate the norm inequality. Let us choose two unit vectors v, w ∈ P that form an angle ψ with respect to g and whose angle bisector lies at an angle of φ to e 1 : v = cos(φ + ψ/2)e 1 + sin(φ + ψ/2)e 2 (A.17) w = cos(φ − ψ/2)e 1 + sin(φ − ψ/2)e 2 (A. 18) We can compute the deformed angleψ(ψ, φ) from the deformed vectors Dv = d 1 cos(φ + ψ/2)e 1 + d 2 sin(φ + ψ/2)e 2 (A.19) by using the arctangent rules with respect to g based on v, w g = Dv, Dw g : By taking the derivative with respect to φ, we can find the minimum and maximum of this function for fixed ψ. The minimum is at φ = 0 and the maximum at φ = π/2 (recall that we chose d 2 > d 1 ). Evaluatingψ(ψ, φ) at these values leads to the inequality This interval becomes maximal when d 1 /d 2 is as small as possible, but for a given metric g, we have d 1 /d 2 ∈ [a/b, 1] for any plane P ∈ V . Thus, we find the following bound which holds in general. For small angles, we can Taylor expand this and find (A.24) • Volume inequality: If we use a metric to measure the volume of some region V A ⊂ A, we use the Lebesgue measure in R d by identifying with A with R d by choosing an orthonormal basis in A. For two metrics g andg, linearity implies that there exists a unique number c, such that Volg(V A ) = cVol g (V A ) holds for any region V A ⊂ A. In order to bound this constant, we can use the norm inequality to show that the d-dimensional unit ball This implies a d ≤ c ≤ b d which leads to the volume inequality we wanted to prove.
If we replace V → V * and accordingly g → G andg →G, we can run exactly the same arguments, but we need to compute This follows from the fact that a and b are the smallest and largest eigenvalue of the linear map (Gg) a b = G acg cb . Under above replacement, we need to consider its inverse map (gG) a b = g acG cb whose smallest and largest eigenvalues are therefore 1/b and 1/a, respectively.

A.2 Lyapunov exponents
In what follows, we restrict ourselves to quadratic systems where M (t) is independent of the initial value ξ 0 . This generalizes to non-quadratic systems by replacing M (t) by M ξ 0 (t). In this case, Lyapunov exponents and vectors depend on the specific trajectory ξ(t) = Φ t (ξ 0 ).
Definition 1 (Lyapunov exponent). Given a linear Hamiltonian flow M (t) : V → V and a vector δξ ∈ V , we define the Lyapunov exponent λ δξ as the limit provided it exists. This definition is independent of the positive definite metric g that induces the norm · . Analogously, we define the Lyapunov exponent of a dual vector ∈ V * as the limit provided it exists. Here, the definition is independent of the inverse metric G.
Proof. We need to prove the independence of this definition from the chosen norm · g induced by some metric g. We can use the norm inequality (A.10) to show M (t) δξ g = c t M (t) δξ g with factor c t ∈ [a, b]. Let λ p be the Lyapunov exponent of δξ ∈ V with respect to the norm · g . We can now computẽ where the second term vanishes because c t is a bounded function. For dual Lyapunov vectors ∈ V * , we can use the same arguments where only our bounds for c t change to c t ∈ [1/b, 1/a].
To characterize the Lyapunov exponents of all vectors in a 2N -dimensional vector space, it is sufficient to select a representative sample of 2N vectors. Such a basis is called Lyapunov basis and is defined as follows. Proof. The construction of a Lyapunov basis as eigenvectors of the limiting matrix L is an important part of Oseledets multiplicative ergodic theorem. A comprehensible proof with further details can be found [90]. The fact that the eigenvectors can always be chosen to form a Darboux basis follows from the fact that L a b is an element of the symplectic algebra sp (2N, R).
When restricting to a subsystem A ⊂ V , it is natural to ask what is the Lyapunov spectrum of the subsystem.
Definition 3 (Subsystem Lyapunov basis and spectrum). Given the linear flow M (t) and a symplectic subspace A ⊂ V , we define the subsystem Lyapunov basis of A as the 2N A vectors ( 1 A , . . . , 2N A A ) with associated subsystem Lyapunov spectrum The subsystem Lyapunov spectrum does in general not consist of conjugated pairs (λ, −λ). Moreover, it is important to emphasize that the Lyapunov spectrum of A is defined as those Lyapunov exponents of linear observables θ ∈ A * , rather than of perturbations δξ ∈ A, because the two are not the same.
The following proposition explains in detail how one can compute the subsystem Lyapunov basis and spectrum when the Lyapunov basis and spectrum of the full system is known.
We refer to the 2N columns of T as t i .
3. Find the first 2N A linearly independent 12 columns t i of T which we can label by t i k with k ranging from 1 to 2N A . The result is a map k → i k ∈ (1, . . . , 2N ) with i k+1 > i k .
The subsystem Lyapunov spectrum is given by (λ A 1 , . . . , λ A 2N A ) with λ A k = λ i k and the subsystem Lyapunov basis is given by Proof. The rectangular matrix T in (A.32) allows us to express the elements of the Darboux basis D A of the subsystem in terms of the Lyapunov basis, θ r = 2N T r i i . Denoting the columns of T by t i we can select the first 2N A linearly independent columns in the ordered set ( t 1 , . . . , t 2N ). We label them t i k and organize them in the 2N A × 2N A square matrix U , Due to their linear independence, the inverse U −1 exists and turns T into an upper triangular matrixT of the form 0 · · · 0 1 * * · · · · · · · · · · · · · · · · · · * 0 · · · · · · · · · 0 1 * * · · · · · · · · · · · · * . . . . . . . . . . . . 0 · · · · · · · · · · · · · · · · · · 0 1 * * · · · * where the * represents an unspecified value. Acting with U −1 on the left and the right-hand side of (A.32) and acting on θ k , we find Given an arbitrary vector θ = its Lyapunov exponent is clearly given by the λ A k = λ i k where k is the smalles i ≥ 1, for which c i is non-zero.
In our geometric representations of the Rényi entropy, we are interested in how the volume of some initial region changes under the Hamiltonian flow M (t). Due to the linearity of M (t), we can restrict ourselves to studying the time-dependent volume of parallelepipeds laying in some subspace A ⊂ V . The evolution will in general evolve this parallelepiped out of A.
Definition 4 (Subsystem exponent). Given the linear flow M (t) and a symplectic subspace A ⊂ V of dimension 2N A , we can define the subsystem exponent as the limit provided it exists. The set V A ⊂ A is an arbitrary parallelepiped spanning all dimensions of A. This definition is independent of the metric that one uses to measure the volume and independent of the choice of parallelepiped V A .
Proof. We need to prove the independence of this definition from the choice of positive definite metric G. We can use the volume inequality (A.13) which ensures that for a different metric (A.39) where the second term vanishes because c t is a bounded function.

A.3 Regular Hamiltonian systems
The central theorem of this paper connects quantum mechanical entanglement with the classical notion of Lyapunov exponents. In order to avoid technical complications, we introduce the class of regular Hamiltonian Lyapunov systems. Most standard Hamiltonian systems that one studies in classical or quantum physics with a finite number of bosonic degrees of freedom fall into this class. (i) All Lyapunov exponents are well defined. This means that for an arbitrary initial condition ξ 0 as well as for every initial separation δξ ∈ T ξ 0 V , the corresponding limit exists.
(ii) All Lyapunov exponents appear in conjugate pairs (λ, −λ), such that the geometric multiplicity of the two conjugate exponents agrees.
In short, condition (i) excludes systems with above-exponential or below-exponential growth, while condition (ii) excludes systems where two or more vectors become exponentially fast collinear under evolution by M ξ 0 (t). Let us give an example for each condition that is not a regular Hamiltonian Lyapunov system. For both examples, we consider a single degree of freedom, such that we can express everything with respect to the Darboux basis D V = (q, p). which has Lyapunov exponents given by λ 1 = 1 and λ 2 = 0. The symplectic volume is still preserved under time evolution because arbitrary initial vectors become exponentially fast collinear, for instance where the angle between the two vectors behaves as regardless of which positive definite metric we use. Clearly, this system does not have two conjugate Lyapunov exponents and does not fall into the class of regular Hamiltonian Lyapunov systems.
The following notion of Lyapunov defect is important to show that for regular Hamiltonian systems the subsystem exponent can be simply computed using theorem 3.
Definition 6 (Subsystem defect). Given the inverse linear flow M (t) and a subsystem A ⊂ V with Lyapunov associated subsystem Lyapunov spectrum (λ A 1 , . . . , λ A 2N A ), we define the Lyapunov defect where V A ⊂ A * is an arbitrary 2N A -dimensional parallelepiped in A. If this limit exists, it is independent of the metric G with which we measure the volume and we have Λ * A ≥ 0.
Proof. The volume of a parallelepiped can be computed from the length of its 2N A sides M (t) i and the (2N A −1) angles ψ i (t), which is the angle between M (t) i and the hyperplane spanned by the vectors M (t) j with j = 1, . . . , i − 1. The time dependent volume is then given by Given two distinct metrics G andG, we can use the angle and length inequalities from above, to find the volume inequality This inequality already insures that the above limit is independent of the chosen metric. Moreover, the explicit expression in (A. 46) shows also that the volume is bounded from above by This implies Λ A ≤ 2N A and thus, Λ * A ≥ 0.
For regular Hamiltonian systems, we can prove the following statements that we will need in the proof of our central theorem of this paper.
Proposition 3. In a regular Hamiltonian system, the Lyapunov defect Λ * A of any subspace A ∈ V vanishes. This implies that the asymptotic behavior of any volume V A ⊂ A ⊂ V is given by where λ A i refers to subsystem Lyapunov spectrum of A.
Proof. Let us recall that there is a special class of metrics on V , for which every symplectic transformation M (t) and thus also M (t) preserves the 2N -dimensional volume. These are all the metrics that give rise to the same volume form as the one induced by the symplectic form. This implies that the asymptotic behavior of every 2N -dimensional region V ⊂ V shows the following asymptotic behavior which holds with respect to all metrics G. From our previous discussion, we also recall that we must have If all Lyapunov exponents λ i come in conjugate pairs with equal multiplicities the sum in this expression vanishes. Thus, we have Λ V = −Λ * V which implies Λ * V = 0 due to Λ V = 0.
At this point, we only need to show that Λ * V = 0 for the full system implies that we also have Λ * A = 0 for all subsystems A ⊂ V . This follows from the fact that we can choose an initial parallelepiped V = V A × V B with well known inequality Let us emphasize that proposition 2 and 3 together provide an alternative full proof of theorem 3, the main result of this paper. Put simply, the subsystem exponent Λ A for regular Hamiltonian systems is just given by the sum over the subsystem spectrum λ A i which can be computed using the procedure explained in theorem 3 or equivalently in proposition 2.

B Gaussian states and quadratic time-dependent Hamiltonians
We review how symplectic methods and complex structures provide a tool for describing Gaussian states and their quantum evolution. These methods are instrumental in the derivation of a relation between symplectic volumes and the asymptotic growth of the entanglement entropy.

B.1 Bosonic quantum systems and the symplectic group
We consider a quantum system with N bosonic degrees of freedom [51]. The Hilbert space H of the system carries a regular representation of the commutation relations [ξ a ,ξ b ] = i Ω ab . (B.1) Here Ω ab is the symplectic structure discussed in section 5.1 and the operatorsξ a can be understood as the quantization of the classical linear observables ξ a with Poisson brackets {ξ a , ξ b } = Ω ab . A Fock representation of the commutation relations (B.1) is obtained by introducing creation and annihilation operators with canonical commutation relations The Hilbert space H is obtained by completing the span of these vectors in the norm induced by the scalar product 0, . . . , 0; D|0, . . . , 0; D = 1. The label D refers to a Darboux basis D = (q i , p i ) of the classical phase space V . It enters in the definition of the representation of the commutation relations (B.1) in the following way. We define position and momentum operatorsq i = q iaξ a ,p i = p iaξ a with Ω ab q ia p j b = δ ij and relate them to the creation and annihilation operators via 13â These relations can be inverted to represent the operatorξ a in terms ofâ i andâ † i , where the symmetric matrix h ab is defined in terms of the generator of a symplectic transformation by M a b = e Ω ac h cb . An immediate consequence of (B.7) is that, for systems with a finite number of degrees of freedom, two Fock space representations associated to different choices of Darboux basis D and D = M D are related by the unitary transformation U (M ). This is a special case of the Stonevon Neumann theorem [114,117]. A second consequence is that classical quadratic observables O = 1 2 h ab ξ a ξ b promoted to operatorsÔ with symmetric (Weyl) ordering have commutation relations that reproduce the classical Poisson brackets, [Ô 1 ,Ô 2 ] = i {O 1 , O 2 }. 14 A third consequence of (B.7) is that the unitary evolution generated by a quadratic Hamiltonian can be fully described in terms of linear symplectic transformations in phase space. This fact plays a major role in the analysis of this paper. 13 Following our index convention, it would be more natural to writeâi = a ibξ b to emphasize their relation to vectors in the complexified phase space V C , but we follow the standard convention of writing creation and annihilations operators asâ † i andâi. 14 This property cannot be extended to higher order observables as shown by the Groenewold-Van Hove no-go theorem [118].
Higher n-point functions are determined by Wick theorem applied to the operatorξ a − ζ a . This property corresponds to the absence of non-Gaussianities: Correlations are completely determined by J and η. Conversely, given the expectation value ζ a and the connected symmetric part G ab of the 2-point correlation function, the Gaussian state |J, ζ is determined by J a b = −G ac ω cb .

B.3 Quadratic time-dependent Hamiltonians
We consider a quadratic time-dependent Hamiltonian H(t), (B.14) The unitary evolution operator solves the Schrödinger equation i ∂ ∂t U (t) = H(t) U (t) and is given by the time-ordered exponential The evolution of the observableξ a is then given by where M a b (t) and η a (t) are defined by the classical Hamiltonian evolution and given in (5.12).
An important property of Gaussian states is that they provide exact solutions of the Schrödinger equation for a time-dependent quadratic Hamiltonian (B.14), i ∂ ∂t |J t , ζ t = H(t)|J t , ζ t . ∂ ∂t ζ t = K(t) ζ t + k(t) , (B.20) with the matrix K a b (t) = (Ω ac h cb (t)) and the vector k a (t) = (Ω ab f b (t)) are defined in terms of the parameters of the quadratic time-dependent Hamiltonian H(t), (B.14). The linear equations for the complex structure J t and the shift η t can be solved as time-ordered series, is the symplectic matrix discussed in (5.12).
Given an initial state |J 0 , ζ 0 and a quadratic time-dependent Hamiltonian H(t), the evolution of the one-point and two-point correlation functions are given by (B.12) and (B.13) with ζ = ζ t and G ab = −J t a c Ω cb .

B.4 Subsystems and the restricted complex structure
We consider a bosonic quantum system consisting of two subsystems A and B with N A and N B degrees of freedom. The Hilbert space of the system decomposes in the tensor product of the Hilbert spaces of the two subsystems, (B.24) The density matrix ρ A of a pure state |ψ ∈ H restricted to the subsystem A is defined by With this definition, the expectation value of any observable in the subsystem A can be computed directly from the density matrix as a trace over the Hilbert space H A , From an operational point of view a subsystem is determined by a subalgebra of observables, i.e. by a restriction of the set of measurements performed on the system. We discuss how the choice of subalgebra of observables identifies the subsystem A, its complement B, and allows us to compute the density matrix of a Gaussian state.
The observables of a bosonic quantum system form a Weyl algebra A V = Weyl(2N, C) generated by linear observablesξ a with commutation relations [ξ a ,ξ b ] = i Ω ab . We define a subsystem with N A degrees of freedom by choosing a subalgebra A A ⊂ A V generated by a set of N A linear observablesθ r , where Ω rs A = Ω ab θ r a θ s b is required to be a symplectic structure on the vector space A = R 2N A , so that the couple (A, Ω A ) is a symplectic vector space. The Hilbert space H A is obtained as a Fock representation of the Weyl algebra A A = Weyl(2N A , C) as discussed in section B.1.
The algebra of observables describing the rest of the system is given by A A , the commutant of A A in A V defined by where Ω kh B = Ω ab Θ k a Θ h b is required to be a symplectic structure on the vector space B * = R 2N B , so that (B * , Ω B ) is a symplectic space. The requirement that B * is the symplectic complement of A * results in the commutation relation [θ r ,Θ k ] = 0. The Hilbert space H B is obtained as a Fock representation of the Weyl algebra A A = A B = Weyl(2N B , C). We call Φ i , Π i a set of canonical observables in B * dual to the Darboux basis D B = (Φ 1 , . . . , Φ N B , Π 1 , . . . , Π N B ).
The subalgebra A A has a trivial center, 15 i.e. A A ∩ A A = 1. As a result the algebra of observables of the systems decomposes in a tensor product over the subsystem A and its complement, A V = A A ⊗ A B , and the Hilbert space of the system decomposes in the tensor product H = H A ⊗ H B . This decomposition reproduces at the quantum level the decomposition of phase space V in two symplectic complements A and B with Darboux basis D V = (D A , D B ).
(B.32) 15 We give an example of subsystem defined by a subalgebra with non-trivial center. Consider a bosonic system with N = 3 degrees of freedom. The algebra AV of observables of the system is generated by elements of the Darboux basis DV = (q1, q2, q3, p1, p2, p3). Let us consider the subalgebra AC generated by (q1, p1, q2). Its commutant is A C = (q3, p3, q2). As a result this subalgebra has a non-trivial center ZC ≡ AC ∩ A C = (1, q2) . In this case the algebra of observables of the system decomposes in AV = ⊕ λ (AC (λ) ⊗ A C (λ) ) and the Hilbert space decomposes in a direct sum of tensor products H = ⊕ λ (HC (λ) ⊗H C (λ) ) where λ is a basis of simultaneous eigenstates of the operators in the center (eigenstates of q2 in this example). Choosing a symplectic subspace as done in (B.28) guaranties that the center of the subalgebra is trivial and the Hilbert space decomposes into a tensor product.
The unitary operator U (ζ A ) generates a shift in A with parameter ζ r A = θ r a ζ a . Note that U (ζ) = U (ζ A ) ⊗ U (ζ B ). The parameters ν i are the positive eigenvalues of the matrix [iJ] A = (iθ r a J a b ϑ b s ) obtained as the restriction to A of iJ,

B.5 Entanglement entropy and Rényi entropy of Gaussian states
Complete knowledge of the state of a system does not imply knowledge of the state of its subsystems. This genuinely quantum-mechanical property is captured by the notion of entanglement entropy. The entanglement entropy S A (|ψ ) of a pure state |ψ restricted to the subsystem A is given by the von Neumann entropy of the reduced state, (B.39) V B A ν 1 ν 2 · · · · · · · · · ν N e · · · · · · · · · ν i = 1 Figure 11. Entanglement structure of Gaussian states. We illustrate the entanglement structure of an arbitrary squeezed vacuum |J with subsystems A and B: We can always find a Darboux frame D V = (D A , D B ), such that only pairs of degrees of freedom are entangled across A and B with squeezing parameters ν i . Every black dot represents a degree of freedom, or equivalently a conjugate variable pair (ϕ i , π i ) appearing as basis vectors in D A or D B , every link represents the entanglement between the two connected degrees of freedom. Note that we take N A ≤ N B and find that only up to N A pairs can be entangled. The remaining N B − N A degrees of freedom in subsystem B do not have a partner in subsystem A leading to squeezing parameters ν i = 0 for i > N A . This is the reason why the maximal number of entangled degrees of freedom is dictated by the smaller of the two subsystems.
To compute the entanglement entropy it is useful to introduce the function Z(β) defined as the trace of the density matrix raised to the power β, The entanglement entropy of a Gaussian state can be computed from Z(β) and expressed in terms of the eigenvalues ν i , [100][101][102][103] S A (|J, ζ ) = This expression plays a central role in the analysis presented in this paper in section 2.