Entanglement of harmonic systems in squeezed states

The entanglement entropy of a free scalar field in its ground state is dominated by an area law term. It is noteworthy, however, that the study of entanglement in scalar field theory has not advanced far beyond the ground state. In this paper, we extend the study of entanglement of harmonic systems, which include free scalar field theory as a continuum limit, to the case of the most general Gaussian states, namely the squeezed states. We find the eigenstates and the spectrum of the reduced density matrix and we calculate the entanglement entropy. We show that our method is equivalent to the correlation matrix method. Finally, we apply our method to free scalar field theory in 1+1 dimensions and show that, for very squeezed states, the entanglement entropy is dominated by a volume term, unlike the ground-state case. Even though the state of the system is time-dependent in a non-trivial manner, this volume term is time-independent. We expect this behaviour to hold in higher dimensions as well, as it emerges in a large-squeezing expansion of the entanglement entropy for a general harmonic system.


Introduction
It is well-known that entanglement entropy in scalar field theory in its ground state is dominated by an area law. This peculiar feature resembles the famous property of black hole entropy, giving rise to simple, but very fundamental, questions: Can black hole entropy be attributed to quantum entanglement? Can gravity be described as a statistical entropic force, attributed to quantum statistics due to entanglement [1][2][3]? Such a description of gravity is at least not contradictory to holographic duality. It has been shown that Einstein's equations in the bulk, or at least the linearized Einstein's equations around the pure AdS geometry, emerge as a holographic realization of the first law of entanglement thermodynamics [4,5], i.e. of the relation which is a trivial identity for any quantum system, and, thus, in the context of the holographic duality, trivially holds for the boundary theory.
The study of entanglement in scalar field theory was initiated long ago. In 1986, Bombelli et al. argued that entanglement entropy in scalar field theory should obey an area law [6]. The actual numerical calculation was performed a few years later by Srednicki [7], who showed that entanglement entropy in free massless scalar field theory in its ground state does obey an area law. For a review of entanglement entropy calculations in field theory, the reader may consult [8][9][10]. A basic element of Srednicki's calculation is the discretization of the degrees of freedom of the scalar field theory on a lattice of spherical shells. This discretization gives rise to a Hamiltonian of an infinite, but countable number of degrees of freedom, i.e. a textbook quantum mechanics harmonic system. Srednicki's calculation is further based on two facts that apply to the specific ground state scenario: 1. The reduced density matrix can be calculated explicitly in coordinate representation via the performance of Gaussian integrals, due to the fact that the ground state of a harmonic system is a Gaussian state.
2. The reduced density matrix is also a Gaussian kernel in coordinate representation and its eigenstates and eigenvalues can be found explicitly. They resemble the eigenstates of an effective harmonic system, with as many degrees of freedom as those that have not been traced out. However, this effective harmonic system does not lie in its ground state, but in a mixed state. More specifically, each normal mode of the effective system appears to lie in a thermal state. However, the reduced system is not strictly thermal, as each mode is characterized by a distinct temperature.
Since the above hold specifically for the ground state of a harmonic system, the generalization of Srednicki's techniques to the study of entanglement in more general states presents a high level of difficulty. However, both facts listed above still hold for a massive field theory. In this case, the inverse of the mass can be used as a perturbative parameter to bypass the numerical part of Srednicki's calculation [11]. Interestingly, both facts are still true in the case that the overall quantum harmonic system lies in a thermal state [12,13]. In such a case, the quantity that is proportional to the area of the entangling surface is not the entanglement entropy, but rather the mutual information of the system. Finally, both facts are also true when the harmonic system lies in a coherent state. It turns out that the entanglement entropy in such a case is identical to that when the system lies in its ground state [14]. Interestingly, the time evolution of the reduced density matrix is unitary and is described by an effective quadratic Hamiltonian with explicit time-dependence [15].
In a seminal paper [16], Page showed that, in an arbitrary quantum state, the entanglement entropy is close to maximal. The possible maximum is proportional to the number of degrees of freedom of the smaller of the two subsystems. The proportionality constant depends on the dimensionality of the Hilbert space of a single degree of freedom. This bound has been connected to the Page curve obeyed by the entropy of black hole radiation [17]. In scalar field theory, the dimensionality of the Hilbert space of a local degree of freedom is infinite, rendering the above bound also infinite. However, it is still proportional to the number of degrees of freedom, implying that we should expect the entanglement entropy to be proportional to the volume of the smaller subsystem, and not the area (see also [18]). In this sense, the states of scalar field theory that have been investigated, namely the ground state and coherent states, cannot be considered as arbitrary quantum states, but as states with special entanglement characteristics. Therefore, it is worth investigating how Page's argument applies to more general quantum states of scalar field theory. Srednicki's method appears to be more easily generalizable to Gaussian states, indicating as a more promising direction the study of the squeezed states. In a similar manner, random Gaussian states have been studied in fermionic systems [19], and it was shown that indeed the mean entanglement entropy approaches that of the Page curve.
In the case of the most general Gaussian state of the overall harmonic system, i.e. a squeezed state, the first fact that facilitates Srednicki's calculation still holds, namely the reduced density matrix can be found explicitly via Gaussian integrals. However, as we will discuss in what follows, the second fact does not apply, making the calculation of entanglement entropy a more complicated task. This problem has been studied at a more abstract level [20,21]. In this work, we develop a method for the calculation of the spectrum of the reduced density matrix, which is the direct generalization of Srednicki's method [7]. Then, we apply this method to the system of free massless scalar field theory in 1+1 dimensions. Furthermore, we study the dependence of entanglement entropy on squeezing and we quantify the non-unitarity of the time evolution of the reduced density matrix.
The structure of the paper is as follows. In section 2 we review basic facts about the coherent and squeezed states of the quantum harmonic oscillator. In section 3 we study the special case of two coupled harmonic oscillators, which can be solved directly with Srednicki's method, and study the dependence of entanglement on squeezing. In section 4 we generalize Srednicki's method to harmonic systems with an arbitrary number of degrees of freedom lying in a squeezed state. In section 5 we develop an expansion for the entanglement entropy for very squeezed states. In section 6 we apply our method to free massless scalar field theory in 1+1 dimensions. In section 7 we discuss our results. In order not to distract from the main part of the analysis, several secondary aspects, as well as technical details of the calculations, have been relegated to a series of appendices. In Appendix A we discuss the time evolution of the reduced density matrix. Appendix B contains the analysis, with all technical details, of how squeezing affects entanglement in the case of the two oscillators. Appendix C presents an algebraic construction of the eigenstates of the reduced density matrix. Appendix D discusses the special solvable case of a harmonic system in a very specific squeezed state, which demonstrates technicalities of our method and provides a consistency check of the expansion of section 5. In Appendix E the validity of our method is checked in a special example that can be solved in a different way, via the extension of the entanglement entropy to the family of Rényi entanglement entropies. In appendix F we develop an expansion for small squeezing parameters.

Gaussian Solutions of the Simple Quantum Harmonic Oscillator
In this section we define the basic terminology in relation to the Gaussian solutions of the simple quantum harmonic oscillator and remind the reader of some basic properties that they possess. The general Gaussian solutions are characterized as "squeezed" states. A special subclass of those with enhanced properties are the so-called "coherent" states.
Let us consider the simple quantum harmonic oscillator with Hamiltonian H =p 2 2m + mω 2x2 2 . (2.1) The time-dependent Schrödinger equation possesses several Gaussian solutions. The most well-known class of such solutions consists of the so-called coherent states, whose wavefunction reads where
4. The ground state is the special case of coherent state with X 0 = 0.
If the Schrödinger equation is solved with an initial condition identified with a Gaussian state of minimal and balanced uncertainties for position and momentum, then its solution is necessarily a coherent state. On the contrary, if it is solved with an initial condition which is a Gaussian state that does not obey both these conditions, i.e. either it is not a minimal uncertainty state, or the uncertainties of position and momentum are not balanced, then its solution is a squeezed state. A squeezed state is still Gaussian at all times, like a coherent state, but it is more general: The functions x 0 (t) and p 0 (t) are given by equations (2.4) and (2.5), as in the case of the coherent states. The parameter z is called the squeezing parameter. The squeezed states retain some, but not all, properties of the coherent states. Namely: 1. The mean position and momentum follow a classical orbit, exactly like in the case of coherent states.
2. They are not minimal uncertainty states at all times. More specifically: It follows that they are minimal uncertainty states exactly four times during each period of the corresponding classical harmonic oscillator, namely at times t min = t 0 + T 8 + n T 4 , n ∈ Z, (2.13) where T is the period of the oscillator, i.e. T = 2π/ω. At these instants w, defined in (2.10), is real.
3. Even at the instants that the squeezed states are minimal uncertainty states, the uncertainties of position and momentum are not balanced. In general ∆x = 2mRe (w) , ∆p = mRe (w) 2 1 + Im (w) Re (w) 2 , (2.14) implying that, at the instants that the squeezed state is a minimal uncertainty state, we have ∆x min = 2mωe ±z , ∆p min = mωe ±z 2 .
(2.15) 4. The quadratic part of the exponent of the Gaussian wavefunction is not real. 5. The coherent states are special cases of squeezed states with z = 0.
Because of the properties of the coherent and squeezed states, one can conceive the coherent states as the closest to classical states of the quantum harmonic oscillator, whereas the squeezed states as the next to closest. For this reason, the study of entanglement in systems lying in squeezed states presents a certain interest, since it reveals the behaviour of the system in states which are not very close to classicality. For example, it is known that the spectrum of the reduced density matrix for an arbitrary coherent state is identical to that for the ground state. This does not hold for the squeezed states.
As we will discuss in what follows, the study of entanglement in squeezed states of the overall system is much more difficult than the study of entanglement in coherent states, because the quadratic part of the exponent of the Gaussian state is not real. Although this fact does not complicate the explicit calculation of the reduced density matrix via the use of Gaussian integrals, the calculation of its spectrum is much more involved.

The Special Case of Two Oscillators
Let us consider the case of two coupled oscillators. From now on for simplicity we use units where = 1. Without loss of generality we consider that the mass of each oscillator is equal to one. Furthermore, for simplicity of the presentation of this toy case, we consider identical self couplings of the two oscillators. The Hamiltonian of the system is In terms of the canonical coordinates the Hamiltonian assumes the form where ω + = √ k 0 and ω − = √ k 0 + 2k 1 . As expected, the normal modes are decoupled.

The Reduced Density Matrix
We consider the overall system lying in a squeezed state; by that we mean that each normal mode is described by a wavefunction of the form of equation (2.9). It follows that, at any given time, the state of the two-oscillator system can be written as Obviously, we demand that Re (w + ) > 0 and Re (w − ) > 0, so that the wavefunction is normalizable. The system's density matrix is given by We would like to trace out the second oscillator and find the reduced density matrix of the first one. In order to do so, we need to express the density matrix in terms of the original coordinates x 1 and x 2 . It is convenient to define in a similar manner to (3.2). The reduced density matrix assumes the form where we defined The reduced density matrix for the first oscillator, ρ 1 (x 1 ; x 1 ) = dx 2 ρ (x 1 , x 2 ; x 1 , x 2 ), can be easily calculated via Gaussian integrals. We only need to complete the square for the coordinate that is integrated. After some simple algebra we find The reduced density matrix is appropriately normalized, i.e. Trρ 1 = 1. Notice that it does not depend on the parameter p 02 at all. It is easy to show that As a result, for any w + and w − , as long as their real parts are positive, which is required for the normalizability of the wavefunction of the overall system, Re (γ) is always positive and larger than |β|.

The Spectrum of the Reduced Density Matrix
The coefficient β is real, similarly to the ground or coherent state case. This is enforced by the fact that the reduced density matrix is Hermitian, i.e. ρ 1 (x 1 ; x 1 ) = ρ * 1 (x 1 ; x 1 ). Even though this property of the reduced density matrix is sufficient in order to set the parameter β real in the case of a reduced system with a single degree of freedom, it is not sufficient in the general case. In other words, this is a special property of systems where all oscillators but one are traced out. However, the simple example of the two oscillators exposes an interesting difference to the ground state case or the more general coherent state case: The coefficient γ is complex, namely This is a direct consequence of the fact that the coefficient of the quadratic term w in the exponent of a squeezed state is complex. Nevertheless, the eigenvalues of the reduced density matrix do not depend on the imaginary part of the coefficient γ. Moreover, they do not depend on the parameter p 01 . We may write the reduced density matrix ρ 1 as i.e.ρ 1 is the same as ρ 1 , where we have set the imaginary part of the coefficient γ and the coefficient p 01 equal to zero by hand. Letf (x 1 ) be an eigenstate ofρ 1 with eigenvalue λ, namely Then, the function f ( is an eigenfunction of the reduced density matrix ρ 1 with the same eigenvalue. Therefore, the spectrum of the reduced density matrix does not depend on the imaginary part of γ and the parameter p 01 . In order to specify its spectrum, it suffices to specify the spectrum of the matrixρ 1 . The matrixρ 1 (x 1 ; x 1 ) is of the same form as the reduced density matrix in the case of the ground [7] or coherent states [15]. It is well known that its normalized eigenstates arẽ and H n is the Hermite polynomial of order n. Notice that the parameter α is always real, since Re (γ) > |β|. The corresponding eigenvalues are where The eigenvalues are properly normalized, since obviously ∞ n=0 p n = 1. It directly follows that the normalized eigenstates of the reduced density matrix ρ 1 ( with corresponding eigenvalues given by equation (3.20). The above imply that the reduced density matrix assumes the form As in the case of the ground state, the entanglement entropy is given by The parameter ξ and the spectrum of the reduced density matrix depend on time, unlike the case of the ground state or a coherent state. This implies that the time evolution of the reduced density matrix is non-unitary. However, the simple form of its spectrum allows the description of the non-unitary part of its evolution in terms of a single quantum jump operator. More information on this point is provided in appendix A.2.

Squeezing and Entanglement
We have obtained the spectrum of the reduced density matrix, as well as the entanglement entropy, for a system of two degrees of freedom in a squeezed state. This state is in some sense "less classical" than a coherent state (see discussion in section 2), so it is natural to ask how the entanglement entropy is altered by squeezing. For the case of two simple coupled harmonic oscillators, we study some representative examples. The details are presented in appendix B, while here we only summarize the conclusions. When the overall system is described by a squeezed state, the entanglement entropy depends on time. When only one mode is squeezed, while the other lies in the ground or a coherent state, this dependence results in a periodic oscillation between a maximal and a minimal value, as shown in figure 1. The contribution of the squeezed mode to entanglement grows as the squeezing parameter increases. However, this contribution may act constructively or destructively on the contribution of the mode that lies in its ground state. As a result, the maximal value of entanglement entropy always increases as the squeezing parameter increases, whereas the minimal value in not monotonous with the squeezing parameter. As the latter increases, the minimal value decreases, up to a critical value of the squeezing parameter z 0 = ln ω − ω + , where it vanishes. Further increase of the squeezing parameter results in an increase of the minimal value of the entanglement entropy. These features are depicted in figure 2. The contribution of the squeezed mode is constructive at the instants when the squeezed state is a minimal uncertainty state with maximal position uncertainty, and destructive at the instants when it is a minimal uncertainty state with minimal position uncertainty.
The dependence of the minimal and maximal entanglement entropy on the squeezing parameter is such that the mean entanglement entropy always increases with the squeezing parameter. After some tedious algebra that is included in appendix B.1, it turns out that the mean entanglement entropy is given by equation (B.36). For small squeezing the mean entanglement entropy is quadratic in the squeezing parameter, whereas for large squeezing t S EE it becomes a linear function: where ξ 0 and S 0 are the parameter ξ and the entanglement entropy when the system lies in its ground state, respectively. The mean entanglement entropy as a function of the squeezing parameter is depicted in figure 3. This wave-like addition of the contributions of the two modes to the entanglement entropy persists when both modes are squeezed. As a result, even when both squeezing parameters are large, the minimal entanglement entropy may be small or even vanishing, when the difference of the two squeezing parameters is small. This is depicted in figure 4.
On general grounds, the contribution to entanglement by the two modes is the maximal possible at instants when both modes are minimal uncertainty states, but one of those has maximal position uncertainty and the other has minimal position uncertainty. When Figure 4. The minimal entanglement entropy relatively to S 0 as a function of the squeezing parameters.
they are both minimal uncertainty states and have both either maximal or minimal position uncertainty, squeezing acts competitively and their contribution to entanglement is the minimal possible. When both modes are squeezed, the dependence of the entanglement entropy on time is in general not periodic, as the ratio of the eigenfrequencies of the two modes may be irrational. However, the constructive and destructive addition of the contri-  Unlike the case where a single mode is squeezed, it is not possible to acquire an analytic formula for the mean value of entanglement entropy. However, it appears that the latter is an increasing function of both z ± , as shown in the numerical calculation depicted in figure 6. Indicatively, the small and large squeezing parameter expansions of the mean entanglement Figure 6. The mean entanglement entropy as a function of the squeezing parameters for ω + = 1, 26) in line with the above statement.

General Harmonic System at a Squeezed State
Having studied the system of two coupled harmonic oscillators in the previous section, we can now proceed to study a more general harmonic system with an arbitrary number of degrees of freedom. We consider a system of N coupled quantum harmonic oscillators described by the Hamiltonian where we use the vector notation Similarly to section 3 and without loss of generality, we have assumed that all oscillators have unit mass. The matrix K is symmetric and positive definite, so that it describes an oscillatory system with N degrees of freedom around a stable equilibrium position. There is an orthogonal transformation O, relating the coordinates x i to the normal coordinatesx i , which diagonalizes the matrix K, reducing the system to a set of decoupled harmonic oscillators, one for each normal mode 1 . In other words, The diagonal matrixK contains the squares of the eigenfrequencies of the normal modes It is obviously related to the initial matrix K as We consider states of the system where all normal modes lie in a squeezed state whereW is the diagonal matrix whose diagonal elements are the complex values w i , i.e.
with w i given by the application of formula (2.10) for each normal mode. We define as subsystem 1 the set of n oscillators described by the coordinates x j , where j ≤ n. The N − n oscillators described by coordinates x j , where j > n, comprise the complementary subsystem, which we call subsystem 2. We would like to trace out subsystem 1 in order to find the reduced density matrix for subsystem 2 and the corresponding entanglement entropy.
The bulk of this section is quite technical. For this reason we would like to provide first the reader with the summary of the results and point out the basic differences with respect to the case of the ground state, which has been extensively studied in the literature.

Summary
The reduced density matrix, which describes a subsystem of the overall harmonic system, turns out to be of the form where γ is a complex symmetric matrix, β is a Hermitian matrix and y 2 = x 2 − x 02 . The vector x 2 contains the coordinates of the degrees of freedom that have not been traced out, while the vectors x 02 and p 02 contain the corresponding parameters x 0i and p 0i .
The imaginary part of the matrix γ does not affect the eigenvalues of the reduced density matrix and can be set to zero by hand. The same holds for the parameters in the vectors x 02 and p 02 . Therefore, we are left with the problem of the specification of the eigenvalues of a matrix of the form where γ is a real symmetric matrix and β is a Hermitian matrix. This is very similar to the reduced density matrix in the case that the overall harmonic system lies in its ground state. However,there is one important difference: for the ground state, the matrix β is not Hermitian, but real and symmetric.
In the case of the ground state, the fact that both matrices γ and β are real and symmetric allows the specification of the spectrum of the reduced density matrix in a trivial manner. One has to perform three coordinate transformations. The first one is an orthogonal transformation that diagonalizes the matrix γ. The second is a coordinate rescaling that sets the matrix γ equal to the identity matrix. The last one is an orthogonal transformation that diagonalizes the matrix β. After these transformations, the density matrix is written as the tensor product of N − n matrices of the form describing one degree of freedom each. The coordinatesx i are the coordinates describing the reduced system after the three linear transformations that we performed above, i.e. real linear combinations of the original coordinates. Namely, they equalx i = v T i x 2 , where v i are the normalized eigenvectors of the real symmetric matrix The parametersβ i are the corresponding eigenvalues of the matrixβ. The above density matrix is of the form of equation (3.16), i.e. of the form that we met in the simple case of two coupled harmonic oscillators. We know that its eigenfunctions are given by equation (3.22) and its eigenvalues by equation (3.20). The eigenfunctions are identical to the eigenstates of an effective simple harmonic oscillator with eigenfrequency equal to α i = 1 −β 2 i . This fact, combined with the form of the eigenvalues, implies that the above reduced density matrix is identical to a thermal density matrix describing the effective harmonic oscillator at a temperature Returning to the reduced system, the fact that the reduced density matrix can be factored to matrices of the above form, describing one degree of freedom each, implies that it is identical to the density matrix describing an effective harmonic system in a quasithermal state. The coordinatesx i are the "canonical" coordinates of this harmonic system, whereas the values α i are the corresponding eigenfrequencies. The eigenstates of the system are trivially given bỹ (4.14) The state is quasi-thermal is the sense that each normal mode is in a thermal state, but has its own temperature. Such a state is not unexpected, considering that the normal modes do not interact. In an obvious manner, the spectrum of the reduced density matrix is of the form In our case, the overall system does not lie in its ground state, but rather in a squeezed state. The matrix β is not real; it is Hermitian. We may diagonalize the matrix γ via an orthogonal transformation and even rescale the coordinates in order to set γ equal to the identity matrix. Nevertheless it is not possible to diagonalize the matrix β through another orthogonal transformation.
It follows that the reduced density matrix cannot be factored to the tensor product of matrices, each describing a single degree of freedom. However, we can still perform the first two coordinate transformations and express the reduced density matrix in the form 2 Rather surprisingly, it turns out that the general properties of the eigenstates and eigenvalues of the reduced density matrix remain the same as when β is real and symmetric.
Even though we cannot factor the system, we may search for eigenfunctions of the reduced density matrix that are of similar form to those in the case of the ground state. First, there is a Gaussian "ground" eigenstate where the matrix A satisfies the quadratic equation This equation has many solutions, but only one gives rise to a normalizable Gaussian eigenstate. Second, there are N − n "first excited" eigenstates where A satisfies (4.18) and the vectors v i are the right eigenvectors of the matrix Let us call ξ i the eigenvalue of the matrix Ξ that corresponds to the eigenvector v i . Then, it turns out that the eigenvalue of the eigenstate Ψ 1i of the reduced density matrix is λ 0 ξ i , where λ 0 is the eigenvalue of the "ground" eigenstate (4.17). The vectors v i are in general complex. The matrix Ξ has no specific symmetry property, nevertheless, its eigenvalues are real. If the overall system lay in its ground state, we would upgrade the linear combinations of the coordinates v T i x to Hermite polynomials of those and construct the whole tower of states. In our case, this is not possible. However, testing the function it can be shown that, although it is not an eigenstate, there is always a unique way to add terms of lower order to the polynomial we obtain an eigenstate of the reduced density matrix. The corresponding eigenvalue is It follows that the spectrum of the reduced density matrix in the case of squeezed states has the same form as in the case of the ground state, namely (4.15). The difference is the following: in the case of the ground state, the values of the parameters ξ i are determined by the eigevalues of the matrixβ, defined in equation (4.12), via the formula (4.13). In the case of the squeezed state, the ξ i are the eigenvalues of the matrix Ξ, defined in equation (4.20). In the coherent limit of the squeezed states the two definitions become equivalent.
The matrix Ξ is defined via the matrix A. This introduces an extra difficulty: in order to calculate the matrix A one needs to solve the non-linear matrix equation (4.18) and specify which solution gives rise to normalizable eigenstates of the reduced density matrix. It turns out that there is only a single admissible solution for A and thus for the matrix Ξ. The eigenvalues of the admissible matrix Ξ can be specified as the solutions of the equation that are smaller than 1. The above equation has in general 2 (N − n) solutions that come in pairs of the form (λ, 1/λ). Therefore N − n of its solutions are smaller than 1 and the other N − n are larger than 1. The next four subsections contain all the technical details of the calculation of the eigenfunctions and eigenvalues of the reduced density matrix. The reader who is not interested in these details may skip them and move to the next section. In appendix A.3 we discuss the time evolution of the reduced density matrix. In appendix C we present a systematic iterative method to construct the eigenfunctions of the reduced density matrix based on creation and annihilation operators. In appendix E we present an explicit realization of this calculational process in a toy example. We confirm that the result is in agreement with that obtained through an alternative method using the Rényi entropies, which is not applicable to the general case, but is feasible in this simple example.

The Reduced Density Matrix
In order to trace out the degrees of freedom of subsystem 1, we need to express the state (4.7) in terms of the original coordinates x, The matrix W is a complex symmetric matrix. The density matrix describing the overall system assumes the form We use the block form notation where the matrix A is an n × n matrix, the matrix C is an (N − n) × (N − n) matrix and so on. Notice that the matrices A and C are complex symmetric matrices, whereas the matrix B is not even a square matrix. Using this block form notation, the reduced density matrix describing subsystem 2, , can be easily found via the application of multidimensional Gaussian integrals. It is a matter of simple algebra to show that and Notice that the matrix γ is by definition a complex symmetric matrix; it obeys γ T = γ. On the contrary, the matrix β is by definition a Hermitian matrix; it obeys β † = β. For this reason, in the case of the two oscillators, where the matrices γ and β were numbers, the coefficient γ was complex, whereas the coefficient β was forced to be real, as we saw in section 3.

The Eigenproblem for the Reduced Density Matrix
Similarly to the case of the two oscillators of section 3, the imaginary part of the matrix γ does not affect the eigenvalues of the reduced density matrix. We may write this matrix as Consider the matrixρ 2 which is identical to the reduced density matrix ρ 2 , where we have set by hand the imaginary part of γ and the shifts x 02 and p 02 to zero, i.e.
Then, the function is trivially an eigenfunction of the reduced density matrix ρ 2 with the same eigenvalue. Therefore, it is sufficient to find the spectrum of the simpler matrixρ 2 . We may further simplifyρ 2 via the following linear transformations of the coordinates: 1. A real orthogonal transformation of the coordinates x 2 , which diagonalizes the matrix γ.

2.
A rescaling of the coordinates x 2 , so that the matrix Re (γ) becomes the identity matrix.
Let us denote the coordinates after these two transformations asx 2 . In an obvious manner Then, the matrixρ 2 assumes the form The calculation of its spectrum cannot continue along the same path as in the case of the ground or coherent state. In such a case, the matrices γ and β would be both real and symmetric and so would be the matrixβ. So, we would apply a final real orthogonal transformation, which would diagonalizeβ and effectively factorize the problem to problems of a single degree of freedom, rendering the calculation of the density matrix eigenstates and eigenvalues trivial. This is not possible in our case. The matrixβ is not real and symmetric, but rather it is a Hermitian matrix, and thus it cannot be diagonalized via a real orthogonal transformation.
In the following, we drop the index 2 from the coordinates that describe the degrees of freedom of subsystem 2.

The Eigenstates of the Reduced Density Matrix
In the case that the matrixβ is real and symmetric, we know the form of the eigenstates of the reduced density matrix. They are the Fock space states of an effective system of coupled harmonic oscillators. This is not the case when the matrixβ is Hermitian and not real. The structure of the eigenstates is deformed. However, there are several characteristics that remain invariant and allow the specification of the spectrum of the reduced density matrix.

The "Ground" Eigenstate of the Reduced Density Matrix
First, let us investigate whether there is a "ground" state similar to the case of realβ, i.e. a Gaussian state. Its existence is supported by the fact the the reduced density matrix is also Gaussian. Consider the normalized wavefunction The matrix A is in general a complex symmetric matrix, which needs to be specified so that Ψ 0 is an eigenstate of the reduced density matrix. It is a matter of algebra to show thatρ Performing this integral yields where .
It follows that the Gaussian state (4.36) is indeed an eigenstate of the reduced density matrix, if Then, the corresponding eigenvalue is Since the reduced density matrix is Hermitian, and thus, it has real eigenvalues, the above implies that the matrix I + A has a real determinant, although in general it is a complex symmetric matrix.

The "First Excited" Eigenstates of the Reduced Density Matrix
In the case that the matrixβ is real, there exist eigenstates that are the first excited states of the effective harmonic system. These states are the product of a Gaussian with a linear combination of the coordinates that is the corresponding normal coordinate of the mode that is excited. So let us investigate whether there are eigenstates of the form where v is a constant vector. The algebra is similar to the case of the Gaussian eigenstate: (4.44) In order to perform this integral, we complete the square in the exponent, as in the case of the "ground" eigenstate. This yields It follows that the state (4.43) is indeed an eigenstate of the reduced density matrix, as long as it has the same matrix A as the Gaussian eigenstate, i.e. the solution of equation (4.41), and furthermore the vector v is a right eigenvector of the matrix Notice that the matrix Ξ in general is not symmetric or real or Hermitian. Let v i be the eigenvectors of the matrix Ξ and ξ i the corresponding eigenvalues, i.e.
Then, we have found N − n "first excited" eigenstates of the reduced density matrix. They read and the corresponding eigenvalues are The fact that the density matrix is Hermitian implies that its eigenvalues are real; therefore the eigenvalues ξ i of the matrix Ξ are real, although this matrix has no particular symmetry property. The fact that the density matrix is Hermitian further implies that these eigenstates are orthogonal, although the eigenvectors of the matrix Ξ are not necessarily orthogonal either in the real or the complex sense. It is a matter of simple algebra to show that demanding that the states Ψ 1i are not only orthogonal, but also normalized, yields This implies that the eigenvectors of the matrix Ξ are orthogonal in the complex sense upon the introduction of a real metric, which is equal to the inverse of the real part of the matrix A. In what follows we always choose the eigenvectors of the matrix Ξ to be orthonormal in Obviously, this definition of the vectors v i , combined with the equation (4.50), implies that the appropriate choice for the normalization constant of the "first excited" eigenstates (4.48) is The above also imply that the matrix has the same eigenvalues as Ξ, and its eigenvectors, which are simply It follows that the matrix Ξ is Hermitian.

The Tower of Eigenstates of the Reduced Density Matrix
In the case of realβ the construction of the rest of the tower of eigenstates is trivial, since they constitute the Fock space of an effective harmonic system. In our case of interest, namely that of complexβ, the construction of the whole tower of states is not that simple. Let as consider the "second excited" state The indices i and j may coincide or not. In the case of realβ, we would expect that this is an eigenstate if i = j, whereas, if i = j, the square v T i x 2 should be corrected to the second order Hermite polynomial of v T i x. In the case of complexβ, the above do not hold. Let us study the action of the density matrix on this state: (4.56) We complete the square in the exponent as in the case of the Gaussian eigenstate. This yields The above integral does not vanish either when i = j or when i = j, unlike the case of realβ. Therefore, the state ψ 2ij is not an eigenstate of the reduced density matrix, as the action of the latter on ψ 2ij gives a linear combination of ψ 2ij and the Gaussian "ground" eigenstate Ψ 0 . Nevertheless, it follows that it is trivial to construct an eigenstate of the reduced density matrix by taking an appropriate linear combination of ψ 2ij and Ψ 0 , namely More interestingly though, we do not need to explicitly calculate the coefficient c 0ij in order to specify the eigenvalue of Ψ 2ij . This is the coefficient of ψ 2ij inρ 2 ψ 2ij . The addition of c 0ij Ψ 0 cannot alter this term; it only corrects the subleading terms. Therefore, there is a set of "second excited" eigenstates, the states Ψ 2ij , with corresponding eigenvalues It is not difficult to show that this argument holds inductively. It is always possible to build an eigenstate Ψ {m 1 ,m 2 ,...,mn} (x), whose higher-order term is The eigenfunctions of the matrix ρ 2 can also be constructed algebraically with the use of appropriate creation and annihilation operators. The details of this construction are presented in appendix C.
The definition of the matrix Ξ (4.46), combined with the defining equation of the matrix A (4.41), implies that Using the definition of the matrix Ξ, the trick I = 1 2 (I + A) + 1 2 (I − A) and the above relation, yields This, combined with equation (4.42), implies that Thus, the eigenvalues of the eigenfunctions that we have already discovered are equal to It follows that Therefore, we have discovered all the eigenfunctions of the reduced density matrix, or at least all the eigenfunctions with non-vanishing eigenvalues.
As in the simple case of the two oscillators, the parameters ξ i are time-dependent, implying that the time evolution of the reduced density matrix includes a non-unitary part. This can be described by as many quantum jump operators as the number of degrees of freedom of the smaller subsystem. A discussion of this point is presented in appendix A.3.
Finally, equation (4.68) implies that the entanglement entropy is given by the same formula as in the case of the ground or coherent states, namely

A Comment on the Symmetry Property of Entanglement Entropy
We found above that the spectrum of the reduced density matrix is given by equation (4.68), where ξ i are the eigenvalues of the square matrix Ξ, defined in (4.46), whose dimension is equal to the number of degrees of freedom of the reduced system. It follows that the entanglement entropy is given by the formula (4.70).
It is well known that, when the overall system lies in a pure state, which is the case for the system we study in this work, the entanglement entropy has a symmetry property: the entanglement entropy calculated by the reduced density matrix of subsystem A is identical to the entanglement entropy calculated by the reduced density matrix of the complementary subsystem A C , Within the framework that we performed our calculation of the spectrum of the reduced density matrix, this property may appear peculiar, since the larger of the two subsystems would be characterized by a larger number of parameters ξ i . Actually, the symmetry property does not only require equality of the two entanglement entropies. The two reduced density matrices have identical spectra; the one with the larger dimension has the same eigenvalues as the one with the smaller dimension, plus vanishing eigenvalues.
Given that the spectrum of the reduced density matrix in our case is given by (4.68), these facts imply the following: The matrix Ξ of the subsystem with the larger number of degrees of freedom, namely max (n, N − n), has min (n, N − n) eigenvalues, which are identical to the min (n, N − n) eigenvalues of the matrix Ξ of the subsystem with the smaller number of degrees of freedom. The remaining max (n, N − n) − min (n, N − n) eigenvalues are vanishing.
As long as we are interested in the entanglement entropy or even the spectrum of the reduced density matrix, it is simpler to consider the reduced density matrix of the smaller of the two subsystems. This simplifies numerical calculations, since these are performed with matrices of smaller dimension. More importantly, considering the reduced density matrix for the smaller subsystem eliminates the presence of vanishing eigenvalues in the spectrum of matrices such as β or Ξ, which would render them non-invertible.

The Eigenvalues of the Reduced Density Matrix
The problem of the specification of the spectrum of the matrixρ 2 , and, thus, of the reduced density matrix, has been reduced to the problem of the specification of the eigenvalues of the matrix Ξ, defined in equation (4.46). The matrix Ξ has no specific symmetry property, but nevertheless it has real eigenvalues; we have shown that it is similar to a Hermitian matrix (see equation (4.53)). A problem that appears in this task is that the matrix Ξ is defined in terms of the matrix A, which is a solution of the quadratic equation (4.41). As such, there are many matrices A, and thus matrices Ξ. However, not all of them correspond to normalizable eigenstates of the reduced density matrix. We need to find a systematic way to distinguish which A corresponds to normalizable eigenstates and then calculate the eigenvalues of the matrix Ξ that corresponds to this specific choice.
We write (4.41) in the form of a matrix Riccati equation as The solutions of equations of the form are constructed as follows: We define the matrix Then the original equation is equivalent to the "eigenvalue" problem where Z is a matrix. The first line of this equation implies that Z = M 11 + M 12 W and thus the second line is equivalent to the original Riccati equation (4.73). Given the solutions of the ordinary eigenvalue problem the solutions of (4.73) are where the matrices χ and ψ are constructed using some of the χ j and ψ j as columns, i.e.
In particular when all matrices of (4.73) are square d × d matrices, there are (2d)!/(d!) 2 combinations. Of course not all combinations correspond to valid solutions, since the matrix χ should be invertible. In our case, the matrix M is 2 min (n, N − n) × 2 min (n, N − n) and reads The eigenvalues of the matrix M are specified by the equation Since the determinant is invariant under transposition, the eigenvalues of M come in pairs of the form (λ, 1/λ).
In order to identify admissible solutions for Ξ we have to understand the relation between the eigenvalues of the matrix Ξ and the eigenvalues of the matrix M . It turns out that this relation is quite simple. Recall that the spectrum of the reduced density matrix is given by (4.68). It follows that the eigenvalues of Ξ (and thus the eigenvalues of Ξ T ) should all be not only real and positive, but also smaller than 1.
The eigenvalue problem (4.76) implies that the matrices χ and ψ obey where λ D is a diagonal matrix containing the eigenvalues of M which correspond to the eigenvectors that we used in order to construct the matrices χ and ψ. Equation (4.77) along with (4.83) implies that the matrix Ξ T reads As a direct consequence, the eigenvalues of Ξ coincide with the inverse of the eigenvalues of M which correspond to the eigenvectors that we used in order to construct the matrices χ and ψ. This implies that M has at least min (n, N − n) eigenvalues that are real and larger than 1. Since the eigenvalues of the matrix M come in pairs of the form (λ, 1/λ), it follows that exactly min (n, N − n) eigenvalues of the matrix M are real and larger than 1 and the other min (n, N − n) are real, positive and smaller than 1. This also implies that: • There is a single admissible matrix Ξ. It is constructed using the min (n, N − n) eigenvectors of M which correspond to its eigenvalues that are larger than 1.
• The matrix A, which corresponds to this admissible Ξ, is the only one which gives rise to normalizable eigenstates of the reduced density matrix.
• The eigenvalues of the admissible matrix Ξ are simply the min (n, N − n) solutions of the equation (4.81) that are smaller than 1.
A solvable example that demonstrates the structure of the eigenvalues of the matrix M is presented in appendix D. In this example, although the system lies in a squeezed state, the phases of the modes are selected in a very particular way (they are all equal to zero) and as a result the matrix β is real. Another example, where the matrix β is complex is presented in appendix E, where it is verified that the calculated entanglement entropy is in agreement with a calculation based on the Rényi extension of entanglement entropy.

Large-Squeezing Expansion
In section 3 we studied the system of two coupled oscillators, where we traced out one of them. We showed that in general entanglement entropy increases with squeezing. For small squeezings, the entanglement entropy is a quadratic function of the squeezing parameters, whereas for large squeezings, the entanglement entropy becomes a linear function of the squeezing parameters, see e.g. equation (3.26).
Although it is very difficult to find exact formulae for the case of a general harmonic system with an arbitrary number of degrees of freedom, we would like to find asymptotic expressions in order to check whether the form of dependence of entanglement entropy on squeezing in the system of two oscillators persists in the general harmonic system. In the following we discuss this expansion for large squeezing parameters, which will be relevant for the interpretation of the numerical results of section 6. The expansion for small squeezing parameters is presented in appendix F.
The most important parameter that defines the state of the overall system and consequently the entanglement between subsystems is the coefficient w of the quadratic part of the exponent of the wavefunction of a squeezed mode. We remind the reader that, in the case of the ground state or a coherent state, this coefficient is trivially real and equal to the eigenfrequency of the mode, whereas in the case of a squeezed state, it is complex, depends on the squeezing parameter and is not constant in time, as shown in equation (2.10). Defining = exp(−z), we may rewrite this equation as This can be written as a series in powers of as, i.e. a series forming a large-squeezing expansion for the parameter w.
Notice that this series is convergent only when When taking the limit z → ∞, the above inequality is satisfied for all times except a very small time period centred around the specific instant that the mode in question is a minimal uncertainty state with minimal position uncertainty. For a system containing a large number of modes, the dominant effect to the mean entanglement entropy arises from the bulk of modes, for which this inequality is satisfied. When we study an arbitrary harmonic system with N degrees of freedom, the coefficient w is upgraded to the N × N matrix W , see equation (4.24). The reduced density matrix is expressed directly in terms of the blocks of the matrix W defined in (4.26). We assume that all modes are characterized by large squeezing, i.e. z i 1, and write the mode squeezing parameters as z i = z + ζ i , where z is the mean squeezing parameter. We focus on the case that the parameters ζ i are subleading to the mean parameter z, i.e. ζ i z. We may now define a single parameter = exp(−z), so that the matrix W , as well as its blocks, have an expansion of the form Notice that the imaginary part of W contains only even powers of , whereas the real part of W contains only odd powers of . The leading contribution in is imaginary. In the following, we use the same notation for the expansions of the blocks of W and the matrices A, B and C.
The reduced density matrix is expressed in terms of the matrices γ and β defined in equations (4.28) and (4.29). The above expansion implies that γ and β are given by In the following, we use the notation The important property of the expansions of the matrices γ and β is the fact that the leading contributions are identical, i.e. β (−1) = γ (−1) . Recall that the entanglement entropy is determined by the eigenvalues of the matrix Ξ that corresponds to the normalizable eigenstates of the reduced density matrix. These are identical to the eigenvalues of the matrix M , defined in equation (4.79), which are smaller than 1. This matrix can be written as The matrix M is similar to M and thus, it has the same eigenvalues. It has the expansion (5.12) It is a matter of algebra to show that 11 + M where λ (1) solves the equation Bear in mind that the non-vanishing eigenvalues of the matrix Ξ are as many as the number of degrees of freedom of the smaller subsystem, namely min (n, N − n). The above equation is of order 2 min (n, N − n). We know that the eigenvalues of M come in pairs of the form (λ, 1/λ). It can be shown that the solutions of the above equation also come in pairs of the form λ (1) , −λ (1) . It follows that exactly min (n, N − n) of these solutions are positive. Therefore, the eigenvalues of the matrix Ξ that corresponds to the normalizable eigenstates of the reduced density matrix are In appendix D, the example of a squeezed state with all modes having the same squeezing parameter and vanishing phases is presented. This case is exactly solvable as this specific selection leads to a real matrix β. In this solvable example, the eigenvalues of the matrix Ξ for large squeezing parameters have indeed the form of equation (5.14). This provides a consistency check for our large-squeezing expansion.
To leading order in , the entanglement entropy reads This formula shows that, for large squeezing, the entanglement entropy has a linear dependence on the squeezing parameter. It also shows that the leading term, which depends linearly on the squeezing parameter, is time-independent. The formula is in agreement with (3.26). In that case, z = z + +z − 2 and N − n = 1. Formula (5.16) suggests something very interesting. The leading contribution to the entanglement entropy is proportional to min (n, N − n), i.e. the number of degrees of freedom of the smaller subsystem. It follows that in a continuous harmonic system, like free scalar field theory, in the large squeezing limit, the leading contribution to the entanglement entropy is proportional to the volume of the smaller subsystem. In other words, squeezing generates a violation of the famous area-law property of entanglement entropy. This property apparently holds only when the system lies in a coherent state, which is a closest-to-classical state.

A Field Theory Example
We are particularly interested in the application of the method that we developed in section 4 to the harmonic system of scalar quantum field theory. Our interest is enhanced by the fact that the large squeezing expansion, which we developed in section 5, suggests that the area-law property of entanglement entropy may not persist when the theory lies in a squeezed state.
The calculation of entanglement entropy in scalar field theory in 3+1 dimensions presents several technical difficulties. The usual discretization of the degrees of freedom, which is also employed in the original calculation at the ground state [7], relies on the expansion of the field in spherical harmonic moments. This is obviously a suitable choice when we desire to introduce a spherical entangling surface. However, such a choice makes it difficult to preserve a uniform density of the degrees of freedom. Therefore, an elegant regularization scheme is required, so that both area and volume terms are detectable. On the other hand, a uniform square lattice would solve this problem, but then the entangling surface would not be smooth, giving rise to additional universal terms in the entanglement entropy. For these reasons, we restrict here our attention to scalar field theory in 1+1 dimensions, where these problems do not appear, and leave the study of the 3+1 dimensional system for future work.
The Hamiltonian of a free scalar field in 1+1 dimensions reads We discretize the degrees of freedom, introducing a uniform lattice in space as The discretized Hamiltonian that we obtain reads where we set the boundary conditions ϕ 0 = ϕ N +1 = 0. We introduce this kind of boundary conditions in order to avoid the existence of a zero-frequency mode. This Hamiltonian describes N coupled harmonic oscillators, exactly as studied in section 4. Their Hamiltonian is of the form (4.1), where the couplings matrix K is given by In the following we use a lattice with N = 60. We furthermore consider the case of free massless scalar field theory in 1+1 dimensions, i.e. we assume that µ = 0. In all cases we divide the system in two complementary subsystems; the first one contains the degrees of freedom ϕ j with 1 ≤ j ≤ n, and the second one those with n + 1 ≤ j ≤ N . We indicate the division of the degrees of freedom in these two subsystems by the number n. In our calculations we set a = 1, which is equivalent to measuring time in units of the UV cutoff set by the lattice spacing.
We take advantage of the symmetry property of the entanglement entropy and we always trace out the larger subsystem. This is required in order to apply the method of section 4.6. Furthermore, this speeds up the numerical calculation, since the matrices involved have the smallest possible dimension. Additionally, the required precision of the numerical calculations is achieved more easily when making this choice. The required precision can be high because the local nature of the couplings generates a hierarchy in the eigenvalues of the matrixβ [11]. As a result, an increase in the dimension of the related matrices, not only increases the volume of the required calculations, but also the required precision of them. Indicatively, our calculation, which includes at most 30×30 matrices, requires about 300 significant digits in order to estimate all eigenvalues accurately.
The overall system has N normal modes. The classical motion of the system when the i-th mode is excited is given by where ω i is the frequency of the i-th mode. In the following, we squeeze one or more of these modes and study the entanglement entropy.

Squeezing a Single Mode
Following the example of the toy model of the two coupled oscillators, which we presented in section 3, we first study the system lying in a state where only one normal mode is squeezed; the rest are put in their ground states. In this way, the time evolution of the system is periodic, with period equal to half the period of the corresponding mode, and thus its study is more transparent. Figure 7 shows the entanglement entropy as a function of n for various times. Several cases are presented, which differ with respect to the mode that has been squeezed. The squeezing parameter is always the same. We observe the following: • Squeezing generally increases the entanglement entropy in comparison to that in the ground state. However, there are instants when the entanglement entropy is smaller than that of the ground state for specific values of n. This is more easily visible in the case that the first mode has been squeezed (top left panel of figure 7).
• The increase of the entanglement entropy by squeezing a single mode does not depend strongly on which mode is squeezed, as long as the squeezing parameter is the same.
• The entanglement entropy is oscillating in time with a period half that of the corresponding mode, as expected.
• The oscillation of entanglement entropy with time is generally more intense when a mode with a smaller index has been squeezed.   Figure 7. The entanglement entropy as a function of n for various times when only a single mode has been squeezed with squeezing parameter z = 3.
• The pattern of the amplitude of the oscillation of the entanglement entropy as n varies is interesting. It appears that this pattern is strongly related to the form of the squeezed normal mode.
-The pattern has the form of a stationary wave with several nodes. There are specific values of n, where the amplitude of the oscillation of the entanglement entropy vanishes.
-The nodes of the entanglement entropy oscillation are twice as many as the nodes of the squeezed mode.
-The nodes appear at n where where i is the index of the squeezed mode. The nodes appear at locations where the classical amplitude squared of the oscillation due to the squeezed mode is half of the maximum. Nodes appear at all such locations except for the first and last one. In other words, the k-th node is located at position -The existence of the nodes justifies why the time dependence of the entanglement entropy is suppressed when the squeezed mode is higher.
The dependence of the entanglement entropy on the shape of the squeezed mode is visible on the mean entanglement entropy as well. Figure 8 depicts the mean entanglement entropy as a function of n for various squeezed modes with the same squeezing parameter.
The mean entanglement entropy is always larger than that in the ground state of the system, unlike the entanglement entropy at a given time. We also observe that the mean entanglement entropy is about the same for the vast majority of the modes. Significant differences appear only for the first and last modes. Notice that the curves corresponding to the 11th mode, as well as the 51st mode are almost identical. All intermediate ones are also almost identical, like the one corresponding to the 31st mode, which is also depicted. This can be attributed to the relation between the pattern of entanglement entropy and the shape of the squeezed normal mode that we pointed out above. The classical amplitudes of oscillation of the first and last modes have a strong pattern: there are regions with large and regions with small amplitudes. On the contrary, the classical amplitudes of oscillations for most intermediate modes are more dispersed; thus the similar pattern of the mean entanglement entropy.
The relation between entanglement entropy and the shape of the squeezed mode is also supported by figure 9, which depicts the mean entanglement entropy as a function of the index of the squeezed mode for several divisions of the system in two subsystems, indicated by the integer n. node Entangling Figure 9. The mean entanglement entropy as a function of the index of the squeezed mode for various divisions of the system to two subsystems indicated by the integer n. The squeezed mode has always squeezing parameter z = 3. mean entanglement entropy does not depend strongly on which mode has been squeezed.
The mean entanglement entropy as a function of the order of the squeezed mode i for given n presents as many maxima as n. These maxima are almost equidistant. For example, the mean entanglement entropy for n = 1 has a single maximum around i = 30. This may be attributed to the fact that the amplitude of oscillation A (i) 1 of the first degree of freedom has this kind of dependence on i. Equation (6.5) implies that i.e. indeed A (i) 1 has a single maximum around i = N/2. Similarly we can show that the amplitude of oscillation of the n-th degree of freedom has n almost equidistant maxima, since which are as many as the maxima of the mean entanglement entropy. We studied the dependence of the entanglement entropy on time and on the shape of the squeezed mode. It remains to study its dependence on the squeezing parameter z. Figure 10 depicts the mean entanglement entropy as a function of n for various values of z and for a specific choice of the squeezed mode. We generally observe an increase of the mean entanglement entropy as z increases. For small values of z, this increase changes the shape of the curve as a function of n. Above some critical value of z, a further increase appears to move the curve as a whole. Furthermore, the increase of entanglement entropy appears to be proportional to the increase of z.
In order to clarify this behaviour, we depict in figure 11 the mean entanglement entropy as a function of z for fixed values of n. For small z values, the increase of entanglement entropy is quadratic in z and depends on n. However, after some critical z, this dependence becomes linear and independent of n; all curves have the same slope asymptotically for large z. Furthermore, this slope does not depend on which normal mode is squeezed. In other words, for large z, when only one mode is squeezed where c depends neither on n nor on the order of the squeezed mode. Our discretized version of 1+1 scalar field theory differs from the continuum field theory in three ways: we have introduced a UV cutoff, an IR cutoff and Dirichlet boundary conditions for the normal modes. As we discussed above, it appears that the dependence of the mean entanglement entropy on which mode is squeezed is due to fact that modes are stationary waves, and thus each mode excites the various degrees of freedom with different amplitudes. We expect that this also happens in the continuum limit if we preserve Dirichlet boundary conditions at some specific point, i.e. if we continue defining the theory in a finite segment or in the infinite half-line. If Dirichlet conditions are abandoned, the normal modes would correspond to travelling waves, which are characterized by identical amplitude of oscillation for all degrees of freedom. Therefore, in free scalar field theory defined on the continuous infinite line we expect that the mean entanglement entropy would not depend at all on which mode is squeezed. This should also be a property of the discretized system if periodic boundary conditions are adopted. In this case a mass should be introduced so that a zero-frequency mode is avoided. This investigation is beyond the scope of this work.

Squeezing All Modes
The large squeezing expansion that we presented in section 5 suggests that for large squeezing we should expect that entanglement entropy is dominated by a volume term proportional to the mean squeezing parameter. In the previous subsection, we studied the system of discretized scalar field theory in 1+1 dimensions in a state where only a single mode lies in a squeezed state, whereas all others lie in their ground state. For large squeezing the entanglement entropy is dominated by a term that is proportional to the squeezing parameter. However, this is not a volume term, but rather a constant term. This is not contradictory to the large-squeezing expansion, which requires that the deviation of the squeezing parameter of each mode from the mean is small, while we assumed that only one mode is squeezed.
In order to understand the effect of strong squeezing, in this subsection we study a system in which all modes lie in a squeezed state with the same squeezing parameter. The time evolution of the overall system is much more complicated than in the case when only a single mode has been squeezed. In general, the evolution is not periodic. In figure 12 we show the entanglement entropy at random instants, as a function of the number n that determines the division of the system in two subsystems. We present several cases, which differ in the value of the common squeezing parameter. The initial phases for each normal mode are unimportant. As they change at different rates, even if they are selected to be initially equal, in due time they are more or less random. For this reason, the instants displayed in figure 12 are completely random.
In this figure the black dots depict the entanglement entropy when all modes lie in their ground state. The coloured dots depict the entanglement entropy at the state under study, i.e. when all modes lie in a squeezed state with the same squeezing parameter. Different colors correspond to different instants. Finally, the continuous black line depicts the leading term of the large-squeezing approximation for the entanglement entropy, given by equation Notice that the latter is time-independent and proportional to the volume of the smaller subsystem. We observe the following: • The entanglement entropy generally increases as the squeezing parameter increases.
• The entanglement entropy approaches the large-squeezing approximation formula as the squeezing parameter increases.
• The variations of entanglement entropy with time decrease in comparison to the mean entanglement entropy as the squeezing parameter increases. This is in line with the fact that the leading term of the large-squeezing approximation is time-independent. Figure 13 depicts the mean entanglement entropy as a function of n for various squeezing parameters. The mean has been calculated as the average of 200 random times. The blue dots depict the mean entanglement entropy. The black dots and the black continuous line depict the entanglement entropy in the ground state and the leading term of the largesqueezing expansion, as in figure 12. It is evident that the mean entanglement entropy is dominated by a time-independent volume term when the squeezing parameter z is large. This volume term is proportional to z.
In the continuum limit a → 0, N a → L, i.e. the continuum limit with an IR cutoff equal to 1/L, the large-squeezing expansion suggests that the dominant term of entanglement  entropy is where r is the limit of the product na, i.e. the length of the first subsystem. The leading volume term is UV divergent. If we remove the IR cutoff, defining the theory on the infinite half-line, then In this scenario the one subsystem has length r and is attached to the end of the infinite half-line. The dominant term is not only proportional to the volume of the subsystem, but also time-independent, even though the state of the system has non-trivial time dependence. In order to clarify that the variations of entanglement entropy with time are reduced in comparison to the mean entanglement entropy as the squeezing parameter increases, we calculated the standard deviation of the entanglement entropy at n = 30 for 200 random times, as a function of z. We found that the standard deviation approaches a finite limit as z increases, whereas the mean entanglement entropy increases linearly with z.

Discussion
In this work we studied entanglement in coupled harmonic systems lying in squeezed states. We managed to reduce the problem of the specification of the eigenvalues of the reduced density matrix to a linear eigenproblem, in exactly the same fashion as in the case of a harmonic system lying in its ground state [7] or a coherent state [15].
In an interesting way, the spectrum of the reduced density matrix conserves the same structure as for the ground state or a coherent state. Let the number of degrees of freedom of the considered subsystem be n. Then, the spectrum of the reduced density matrix describing the subsystem is indistinguishable from the spectrum of an effective harmonic system with n degrees of freedom, if each of its normal modes is lying in a thermal state at an appropriate temperature. Notice that this is not a thermal state of the whole effective system; each of its normal modes has a different temperature.
However, there are several important differences compared to the ground or coherent state cases. First, the eigenstates of the reduced density matrix have suffered a non-trivial deformation. Although they can be organized in a similar fashion to the Fock space of an effective harmonic system, there is no real linear combination of the physical degrees of freedom of the considered subsystem which can be identified as a normal coordinate. This is due to the fact that the creation operators associated with this Fock space are a linear combination of positions and momenta which are not conjugate to each other (see appendix C). As a matter of fact it is quite difficult to derive explicit expressions in coordinate representation for the whole set of eigenstates of the reduced density matrix, although it is clear how to construct them iteratively. In other words, in the case of the ground or coherent state of the overall system, the reduced density matrix is separable and it can be written as the tensor product of density matrices that describe a linear combination of the original degrees of freedom each. In the case of the squeezed state, the reduced density matrix is separable, but it is written as a tensor product of matrices that cannot be assigned to any real combination of the original degrees of freedom.
Second, the time evolution of the reduced density matrix, unlike the ground or coherent state case, is non-unitary. However, the non-unitarity is rather elegant. At any instant, the rate of change of the reduced density matrix is determined by an effective Hamiltonian that determines the unitary part of the time evolution, and a single quantum jump operator per degree of freedom of the smaller subsystem that determines the non-unitary part (see appendix A). One should bear in mind that the description of an arbitrary non-unitary evolution of the reduced density matrix requires a number of quantum jump operators that is a quadratic function of the dimension of the Hilbert space of the considered subsystem, which in our case is infinite. In this respect, the non-unitary evolution of the reduced system for squeezed states is surprisingly simple.
Finally, the spectrum of the reduced density matrix, and, thus, the entanglement entropy is in general time-dependent. Although at a given instant the entanglement entropy may be smaller than that at the ground state or a coherent state of the system, the examples that we have investigated suggest that the mean entanglement entropy increases with squeezing. For states in which all the modes are strongly squeezed, the mean entanglement entropy appears to be time-independent and proportional to the absolute value of the squeezing parameter and the number of degrees of freedom of the smaller subsystem (see equation (6.11)). A large-squeezing expansion supports this conclusion (see section 5).
Page has shown that the entanglement entropy of an arbitrary quantum state is close to the maximal possible entanglement entropy [16]. For systems where the local degrees of freedom have a finite-dimensional Hilbert space, this bound on entanglement entropy as a function of the number of degrees of freedom of the considered subsystem has a characteristic concave form. It vanishes when the subsystem is null or coincides with the whole system, while it is maximal when the subsystem contains half of the degrees of freedom of the overall system. In the limit that the total system contains an infinite number of degrees of freedom, this function tends to the union of two linear segments. If n and N denote the number of degrees of freedom of the subsystem and overall system respectively, this curve is approximately S max ∼ min (n, N − n) , (7.1) where the proportionality constant depends on the dimensionality of the Hilbert space of the local degrees of freedom. This curve is of great importance. It has been connected to the Page curve followed by the entropy of black hole radiation [17], which is a critical piece of the information paradox. This relation has also been established within the framework of holographic duality [22,23]. Applying Page's argument to the system of a scalar quantum field theory implies that the entanglement entropy in an arbitrary quantum state should be proportional to the number of degrees of freedom of the smaller subsystem, i.e. proportional to the volume of this subsystem. One has to keep in mind that the Hilbert space of a local degree of freedom in this case is infinite-dimensional, rendering the application of Page's argument in scalar field theory a little hazy. However, assuming that the scaling properties are preserved as we take the limit of the dimension of the Hilbert space to infinity, Page's argument contradicts the seminal results of Bombelli and Srednicki [6,7], which apply to the ground state, as well as their generalizations to coherent states [14,15]. These studies demonstrate that entanglement entropy scales with the area of the subsystem and not its volume. In this sense, the area-law property of entanglement entropy should be considered as a special property of the most classical states of scalar field theory, i.e. of the coherent states.
In order to understand this issue, we applied our method to the system of free massless scalar field theory in 1+1 dimensions. In agreement with the large-squeezing expansion that we developed for an arbitrary harmonic system, we found that states where all modes have been squeezed with large squeezing parameters give rise to entanglement entropy which is dominated by a term proportional to the volume of the smaller subsystem, in agreement with Page's arguments. Furthermore, this volume term is time-independent, although the state of the system has non-trivial time-dependence. This is consistent with a maximal entanglement entropy bound in line with Page. We expect this behaviour to hold for scalar field theory in higher dimensions, which will be the subject of future work.
We can speculate on the consequences of our results for the interpretation of gravity as an entropic force attributed to quantum entanglement statistics. Such an interpretation is supported by holographic calculations [4,5]. However, there are more general arguments that suggest how such an entropic force operates. In 1995 Jacobson argued that the dynamic metric of a theory with two specific properties is subject to Einstein dynamics [1]. The first property is the validity of a first law of thermodynamics. The second one is that the entropies of the horizons are proportional to their area. In other words, the scaling properties of entropy determine the gravitational dynamics. Entanglement is fertile ground in which to realize such a mechanism. The entanglement entropy scales with area (at least when the overall system lies in a coherent state) and also obeys a first law of entanglement thermodynamics with the expectation value of the modular Hamiltonian.
Our investigation suggests that the only Gaussian states that give rise to entanglement entropy dominated by an area law are the minimal uncertainty states, i.e. the coherent states. There is no indication that this property extends beyond the Gaussian states. In this spirit, our results imply that Einstein gravity emerges as a quantum entropic force only when the overall system lies in a closest-to-classical state, i.e. in the ground or a coherent state. When a more "arbitrary" quantum state is considered, the emergent dynamics will be more complicated than Einstein gravity.
As a final comment we point out that the squeezing of quantum states plays an important role in early cosmology. During inflation, a momentum mode of a massless field gets stretched by the rapid expansion. When the mode wavelength becomes larger than the horizon, the scalar fluctuation loses its oscillatory form (it freezes) [24]. After horizon crossing, the field can be viewed as a classical stochastic field, and its quantum expectation value can be replaced by the classical stochastic average. The quantum properties of the field are considered invisible in late-time observations, which focus on classical local quantities [25,26]. However, from a quantum mechanical point of view, the modes of the scalar field evolve from a simple oscillator ground state to an increasingly squeezed state [27]. Quantum entanglement is a purely quantum non-local phenomenon that cannot be encoded in the classical probability distributions. The squeezing of canonical modes increases the entanglement between local degrees of freedom and is expected to increase the entanglement entropy [28][29][30][31][32][33][34][35]. The techniques we presented in this work were employed in [36] in order to compute the entanglement entropy resulting from tracing out local degrees of freedom of a quantum scalar field in an expanding universe. It was shown that the entanglement entropy grows continuously during inflation, as successive modes cross the horizon. The resulting entropy is proportional to the total duration of inflation, and is preserved during a subsequent era of radiation or matter domination. The emergence of a volume term in the entanglement entropy as a result of squeezing was observed in [36] in the context of a toy model in 1+1 dimensions, in agreement with our findings here.

Acknowledgments
The research of D. Katsinis

A.1 General Gaussian States
In section 2 we revisited the Gaussian states of the simple harmonic oscillator, namely the squeezed states and their subclass of coherent states. These states are solutions of a quadratic Hamiltonian with constant coefficients, i.e. the one given by equation (2.1).
More general Gaussian states are of the form (2.9) without the restriction to squeezed states, i.e. without the demand that the time evolution of the parameters appearing in (2.9) be given by (2.4), (2.5), (2.10) and (2.11). Such states emerge as solutions of quadratic Hamiltonians with explicit time dependence. These Hamiltonians can be defined via the introduction of appropriate creation and annihilation operators. Then, a complete tower of solutions can be constructed algebraically.
These states and their corresponding Hamiltonians are interesting for two reasons. First, these Hamiltonians appear as those that evolve the reduced density matrix in harmonic systems lying in a coherent state. In our case, the overall harmonic system lies in a squeezed state and the evolution of the reduced system is not unitary. However, the unitary part of the time evolution is described by such a Hamiltonian. Second, such Hamiltonians, with more general Gaussian time-dependent solutions, appear in several interesting physical systems, like the equations of motion of a scalar field in an expanding background (see e.g. [36]).
Here we review some basic elements regarding these general Gaussian states since we utilize them in the study of the time evolution of the reduced system in section A.2. The reader may consult [15] for more details on the subject.

We define the creation and annihilation operatorŝ
where x 0 , p 0 , w and m are arbitrary functions of time. It is straightforward to verify that they obey the commutation relation The operatorsÂ andÂ † are associated to the tower of states where H n is the Hermite polynomial of order n. The operators act on the states aŝ Writing the above expression in terms of position and momentum operators, we obtain where P, ω 2 and B are given bŷ The Nevertheless, this phase is irrelevant regarding the density matrix, since it is a global unobservable phase. In summary, a Gaussian state with arbitrary time-dependent coefficients is the "vacuum" state of the time-dependent oscillator whose Hamiltonian is given by (A.9).
Imposing that the mass is constant, along with the constraints so that the Hamiltonian (A.9) reduces to that of the simple harmonic oscillator, constraints the general Gaussian state to a squeezed state of the simple harmonic oscillator.

A.2 The Special Case of Two Oscillators
The reduced density matrix describes an open quantum system. As such, it does not undergo unitary time evolution in general. In the case that the overall system lies in its ground state [7], the density matrix, and thus the reduced density matrix are timeindependent. The same holds in the case that the overall system lies in a thermal state [12]. In both these cases, the reduced density matrix undergoes unitary evolution in a trivial manner. When the overall system lies in a coherent state, the reduced system undergoes nontrivial time evolution, which is inherited by the time evolution of the coherent state. However, it turns out that the spectrum of the reduced density matrix does not depend on time. It is identical to the spectrum in the case that the overall system lies in the ground state, and therefore, the reduced density matrix undergoes unitary time evolution. The corresponding Hamiltonian is quadratic, with explicit time dependence [15].
In section 3 we studied the toy case of an overall system with two degrees of freedom lying in a squeezed state. Unlike the coherent state case, it is clear that the spectrum of the reduced density matrix is time-dependent. This is a direct consequence of the fact that the parameters γ and β depend on time, which implies that the parameter ξ and thus the spectrum of the density matrix depend on time. These parameters are determined by the quadratic coefficient of the exponent of the state, i.e. w, which is time-dependent, as shown in equations (3.10) and (3.11). It follows that the time evolution of the reduced density matrix in the case of squeezed states is not unitary. For coherent states the evolution is unitary, as the coefficient w is equal to the eigenfrequency of the corresponding mode, and thus obviously time-independent.
We would like to separate the time evolution of the reduced density matrix into a unitary and a non-unitary component in order to quantify the non-unitarity of the time evolution. In general, a density matrix may evolve in such a way that it preserves its basic properties, i.e. it remains Hermitian, its trace remains equal to 1 and it remains positive definite. For a reduced density matrix the last requirement is a little stricter: it should remain completely positive, i.e. not only the reduced density matrix remains positive, but the overall density matrix as well. These requirements imply that the time evolution of the reduced density matrix is given by a superoperator, i.e.
The differential form of the above equation is the Lindblad equation [37] The Lindblad equation naturally separates the time evolution of the reduced density matrix into a unitary part, given by the first term of the right-hand side of the above equation, and a non-unitary part given by the rest of the right-hand side, which is expressed in terms of the so called Lindblad or quantum jump operators L i . So the task of this section is the identification of appropriate operators H and L i and parameters γ i , so that the reduced density matrix (3.9) obeys equation (A.18). We would like to make clear that we do not use the Lindblad equation in the usual context of Markovian systems. In general, when the complementary system has finite size, the evolution of the reduced system is not given fundamentally by a differential equation, due to memory effects, i.e. information that returns to the reduced system by its environment.
In the limit that the information that escapes from the reduced system to its environment is completely dissipated and never returns, the time evolution of the reduced system is given fundamentally by a differential equation in the form of a Lindblad equation. In such cases the Lindblad operators have a physical meaning that depends on the interaction of the reduced system with its environment, i.e. the Lindblad operators do not depend on the state of the system. In this section, we do not claim that the system is Markovian; in an obvious manner the complementary system is finite and information which escapes from the reduced system may return at a later time. We simply claim that the time evolution of the system is described by a Lindblad equation, which we use to quantify the non-unitarity of the time evolution of the reduced density matrix. The Lindblad operators in general are state-dependent; our result is valid only for squeezed states of the overall system.
In the following we denote which allows the writing of (A.18) aṡ The time evolution of the reduced density matrix (3.16) can be studied taking into account expression (3.23). The latter directly implies thaṫ The eigenstates of the reduced density matrix multiplied with an appropriate overall time-dependent phase are known to be solutions of a time-dependent quadratic oscillator [15], as presented in section A.1, where the phase φ obeysφ We identified the Hermitian operator that generates the unitary part of the time evolution of the reduced density matrix. It remains to identify the Lindblad operators L i . It is natural to assume that they are related to the generalized creation end annihilation operators of the time-dependent oscillator whose Hamiltonian is H eff . These are given by equations (A.1) and (A.2), with the same identification of parameters that we used for the Hamiltonian, namely equations (A.24). It is a matter of algebra to show that Putting everything together, the time evolution of the reduced density matrix is determined by the equatioṅ Therefore the non-unitary time evolution of the reduced density matrix is determined by a single Lindblad operator, which is the generalized creation or annihilation operator, depending on whether the parameter ξ is increasing or decreasing at the given instant. The corresponding coefficient is proportional to the rate of change of the parameter ξ.

A.3 Arbitrary Harmonic System
For an arbitrary harmonic system, similarly to the already studied case of the two coupled harmonic oscillators, the spectrum of the reduced density matrix is in general timedependent. It follows that the time evolution of the reduced density matrix is non-unitary. We would like to separate the time evolution to a unitary and a non-unitary part, as we did in the simple case of the two oscillators, and find the quantum jump operators, which describe the non-unitary part.
Let Ψ {m 1 ,m 2 ,...,mn} (x) be the eigenstates of the reduced density matrix, which we have shown that correspond to eigenvalues (4.63). For notational simplicity, we denote {m 1 , m 2 , . . . , m n } ≡ m. This allows the writing of the reduced density matrix as It follows thaṫ The first term corresponds to unitary evolution of the density matrix. The eigenstates of the reduced density matrix are normalized and orthogonal at any given time. Therefore, their time evolution is necessarily given by a unitary operator, Therefore, where H is the Hermitian matrix which generates the unitary operator U . We define the "creation" and "annihilation" operators A i and A † i , which act on the eigenstates of the reduced density matrix as 3 It is easy to show that The above imply Similarly, The above imply Thus, the time evolution of the reduced density matrix is determined by the equatioṅ (A.41) 3 These operators do not coincide with the creation and annihilation operators that we construct in appendix C. The latter act as in equations (A.33) and (A.34), but on the eigenstates of the matrixρ2 and not on the eigenstates of the reduced density matrix. The operators defined in this section can be constructed from the operators of appendix C after some algebra, which is beyond the scope of this work.
The non-unitary part of the time evolution of the reduced density matrix is described by only as many quantum jump operators as the number of the non-vanishing parameters ξ i . Following the discussion of section 4.5, these are as many as the number of degrees of freedom of the smaller of the two complementary subsystems, namely min (n, N − n). Bearing in mind that the general non-unitary evolution of the system would require an infinite number of quantum jump operators, since the latter are d 2 − 1, where d is the dimension of the Hilbert space of the subsystem 2, the non-unitary part of the time evolution in the case of squeezed states of the system is particularly simple.

B The Dependence of Entanglement on Squeezing in the Case of Two Oscillators
In the case of two oscillators, we have found simple formulae describing the entanglement entropy as a function of the squeezing parameters and time. We can study them in order to acquire an intuition on how squeezing affects entanglement. In this appendix we present all the details of this analysis, whose summary was presented in subsection 3.3. Instead of searching for extrema of the entanglement entropy, it is easier to search for extrema of the ratio r, which is defined as where γ and β are given by equations (3.10) and (3.11). The entanglement entropy is a strictly increasing function of r. It follows that extrema of r correspond to extrema of the entanglement entropy. The ratio r is by definition always larger than 1. Using its definition, we can express the entanglement entropy in a symmetric form, namely This expression is related to the calculation of entanglement entropy in terms of correlation functions (see e.g. [38]). Before presenting our analysis, let us first briefly review the case of the ground/coherent state of the overall system, so we can use it as a basis for comparison. In this case the ratio r is given by It is evident that r 0 and consequently the entanglement entropy depend only on the ratio of the two eigenfrequencies ω + and ω − . Furthermore, it is invariant under the interchange Considering that ρ 0 > 1, the ratio r and thus the entanglement entropy are strictly increasing functions of ρ 0 . The symmetry ω + ↔ ω − obviously implies that, when ρ 0 < 1, ξ 0 is a strictly decreasing function of ρ 0 .

B.1 Squeezing a Single Mode
For simplicity let us first consider the case z − = 0, i.e. we "squeeze" only the symmetric mode. This is a more transparent case, since the time evolution of the parameter ξ, and thus of the entanglement entropy, is periodic with period T + /2 = π/ω + . In this case, the ratio r assumes the simple form which is manifestly positive and larger than 1, as required. This equation directly implies that r is bounded between two values, r ± , It directly follows that the entanglement entropy is bounded between two values S ± that are determined by r ± and equation (B.2). The value r + is obtained at the instants when sin 2ω + t = 1, whereas the value r − is obtained at the instants when sin 2ω + t = −1. At these instants the squeezed state is a minimal uncertainty state, i.e. ∆x + ∆p + = /2 (see equation (2.12)). It follows that, at these instants, the symmetric mode is described by a wavefunction that is indistinguishable from an appropriate coherent state of an effective oscillator with eigenfrequency ω + e ±z + . Recall that the spectrum of the reduced density matrix when the system lies in a coherent state is identical to that when it lies in the ground state. Therefore, at these instants, the symmetric mode is effectively described by the ground state wavefunction of this effective oscillator, at least as long as entanglement is concerned. As is evident from the above discussion, as well as the comparison of equations (B.3) and (B.6), the two bounds on the entanglement entropy in the case of a squeezed symmetric mode are identical to the entanglement entropies of two equivalent effective systems of coupled oscillators at their ground state. The ratios of eigenfrequencies of each of these two equivalent systems are ρ ± = ρ 0 e ±z + . (B.7) We assume that ρ 0 > 1 and z + > 0. It directly follows from equation (B.5) that r + > r − . Therefore, the value of r + determines the maximal value of entanglement entropy, whereas the value of r − determines the minimal value of entanglement entropy.
Recalling the monotonicity of the relation between r and the ratio of eigenfrequencies ρ for a system at its ground state that we analysed above, since ρ 0 > 1 and ρ + > ρ 0 , the maximal entanglement entropy is larger than that at the ground state of the system, S 0 . As long as the minimal entanglement entropy is concerned, there is a change of the qualitative behaviour as the squeezing parameter z + increases. As z + increases from zero to positive values, the ratio ρ − gets smaller, but it remains larger than 1 until the critical squeezing parameter As a result, the minimal entanglement entropy is smaller than S 0 and a decreasing function of the squeezing parameter. For z + = z 0 the minimum entanglement entropy vanishes, i.e. there are instants when the wavefunction of the system is separable. As the squeezing parameter further increases, the ratio ρ − further decreases and it is smaller than 1. As a result, the minimal entanglement entropy becomes an increasing function of the squeezing parameter z + . There is a critical value of the squeezing parameter z + , which results in where the minimal entanglement entropy coincides with S 0 .
For small values of the squeezing parameter, it is not difficult to show that the ratio r and the entanglement entropy perform a sinusoidal oscillation around a mean value, and also calculate this mean value. In particular, one can show that for z + 1 the ratio r assumes the form For z + 1 we obtain Using the inequality x ≥ ln (1 + x) one can show that the coefficient of z 2 + in the small z + expansion of the mean entanglement entropy is manifestly positive. Its minimal value is equal to 1/8, and it is obtained in the limit ξ 0 → 1. This implies that the mean entanglement entropy is increasing with squeezing for small squeezing parameters. 4 In order to perform this calculation the formula 1 2π 2π 0 dt ln (1 + a sin t) = ln is required. 5 Notice that the coefficient diverges in the ω+ = ω− limit. This is an artefact of the order of the limits, i.e. the small z+ limit and the degenerate limit do not commute. To see why this is the case, we turn to the degenerate limit. For ω+ = ω− the parameter ξ and the entanglement entropy simplify a lot. It is straightforward to show that they read ξ degenerate = tanh 2 z+ 4 (B.13) and S degenerate = cosh 2 z+ 4 ln cosh 2 z+ 4 − sinh 2 z+ 4 ln sinh 2 z+ 4 . (B.14) So, the entanglement spectrum of the reduced density matrix and the entanglement entropy change at will.

B.1.1 The Mean Entanglement Entropy
In the special case we are studying, i.e. only the symmetric mode has been squeezed, it is possible to derive an analytic formula for the mean entanglement entropy for an arbitrary value of the squeezing parameter. We begin by rewriting the entanglement entropy, which is given by (B.2), as We use the relatively simple expression (B.5) for r and we split the entanglement entropy into three terms as The first term is constant, soS 1 = S 1 . We are left with the calculation ofS 2 andS 3 . To proceed, notice that (B.20) Thus the absolute value of the coefficient of sin 2ω + t in (B.18) in smaller than one. Therefore,S 2 can be calculated directly by equation (B.12): We are left with the calculation ofS 3 . We expand √ r tanh −1 (1/ √ r) and get The integration can be performed using the formula 1 2π where P k is the Legendre polynomial of order k, to obtain where x and y are given by We separate the k = 0 term, which is constant, to writē We are unable to perform this summation directly. So we introduce a Schwinger parameter and interchange the summation and integration, to obtain The summation can be performed using the generating function of Legendre polynomials to arrive atS or with a trivial change of variablē Since we are calculating a physical quantity, the result should be real for all values of the parameters. To verify this fact, recall that x > 1 and 1 > y > 0. The quantity under the square root gets its minimum value for w = √ y, while we have Thus, the result is manifestly real for any value of the parameters. Setting x = cosh u we obtainS We perform another change of the integration variable to obtain where F and E are the incomplete elliptic integrals of first and second kind respectively, and e −2u is their elliptic modulus. The careful reader would have noticed that initially the sign of u was irrelevant, but we have treated e u and e −u differently. Had we made the opposite choice, we would have ended up with a symmetric formula. Of course, this result can also be obtained using the transformations of the elliptic integrals under the inversion of the elliptic modulus, which is required for consistency. In order to substitute the original parameters we use Gathering all the terms, i.e. the above result, along with (B.17) and (B.21), we obtain where ω 2 This expression for the mean entanglement entropy is an increasing function of the squeezing parameter.
Notice that (B.36) is manifestly invariant under the interchange ω + ↔ ω − along with z + → −z + . Had we performed only one of these transformations, it would require the transformation of the elliptic integrals under the inversion of the elliptic modulus to show that (B.36) is indeed invariant. Equation (B.36) also explains why the small z + and the degenerate ω + = ω − limits do not commute. For ω + = ω − in the small z + limit the quantity inside the absolute value is negative, whereas for finite z + this quantity is positive in the ω + − ω − → 0 limit. On the contrary, in the large |z + | limit this quantity is positive, just like in the ω + − ω − → 0 limit, thus these limits commute. Finally, the elliptic modulus takes any positive value, but it is equal to 1 either for z + = 0 or ω + = ω − .

B.2 Squeezing Both Modes
When both the symmetric and antisymmetric modes are squeezed, the parameter r equals where Φ ± are the phases of the two modes, namely, Φ ± = 2ω ± (t − t 0± ). Without loss of generality, in the following we consider that z + and z − are both positive. The introduction of a negative squeezing parameter is equivalent to a shift of the corresponding phase Φ ± by π. Furthermore, we assume that ω − > ω + . It is a matter of algebra to show that which implies that r is manifestly greater or equal to 1, as required.
Unlike the case where we had squeezed only the symmetric mode, the ratio r and thus the entanglement entropy is not necessarily a periodic function of time. Actually it is periodic if and only if the ratio of the eigenfrequencies of the two normal modes is rational. In such a case, the phases Φ ± follow a one-dimensional closed path in the Φ + Φ − plane. Otherwise, they follow an open trajectory, which in infinite time will cover the whole region of possible (Φ + , Φ − ) pairs, i.e. the two phases get arbitrarily close to any given pair of admissible values at some instant. In this spirit, we search for the extrema of the ratio r, and thus of entanglement entropy, treating the two phases as independent, although actually they are not; they are both given functions of time.
There are stationary points of the ratio r when the phases Φ ± satisfy the equations 6 Notice that there is another mathematical solution, namely However, it is unphysical, since for any values of the parameters it does not correspond to real Φ+ and Φ−.
Whenever the above condition holds, equations (B.40) and (B.41) do have two solutions, both corresponding to the ratio r being equal to 1. Actually, it is evident from equation (B.39) that these are the only values of the phases Φ ± where the ratio r can be equal to 1. Since this is the minimal possible value of r, whenever these solutions exist, they provide the global minimum of the ratio r. Equations (B.42) have always 4 solutions in a trivial manner, namely They correspond to the following values of the ratio r where the symbols s ± take the values ±1.
In an obvious manner r +− is larger than r ++ , r −+ and r −− . Therefore, Φ + = π 2 and Φ − = − π 2 is the position of the global maximum of the ratio r, When z + < ln ω − ω + and z − < ln ω − ω + the smallest of the four r s + s − is r −+ . When the above does not hold and z + > z − the smallest is r −− , whereas when z + < z − the smallest is r ++ . The smallest of the four r s + s − is the global minimum of the ratio r, whenever condition (B.45) does not hold, i.e.
Since the value of the ratio r for the vacuum state of the two oscillators can be written as r 0 = cosh 2 1 2 ln ω − ω + , it turns out that the globally minimal value of the ratio r and thus of the entanglement entropy coincides with that of the vacuum state if The values of the phases Φ ± that correspond to the extrema r ±± , namely Φ ± = ±π/2, are not arbitrary. When a phase takes one of these two values, the wavefunction that describes the corresponding mode is a minimal uncertainty state with maximal or minimal position uncertainty, respectively.
The fact that the ratio r is stationary when the two modes are both minimal uncertainty states suggests that these are the times that contributions of squeezed modes to entanglement are maximal. However, these contributions add when the two modes are in opposite phases. When the two phases are equal, the contributions cancel each other, resulting in weak entanglement. The quantity that receives negative or positive contributions directly equal to the squeezing parameter of each mode at these instants is arccosh √ r. In the above we dealt with the two phases as independent variables. Actually, they are not; they are both functions of time. As time flows, the system follows a specific one dimensional trajectory within the two-dimensional space of the phases of the two modes. The trajectory depends on the ratio of the frequencies of the two modes. Because of this fact, the ratio r may present local minima or maxima with time, which do not coincide with the theoretical minima and maxima r min and r max that we specified above. However, the ratio r is always bound by these values.

C Algebraic Construction of the Reduced Density Matrix Eigenstates
In section 4.4, we showed that the eigenfunctions of the matrixρ 2 , and thus those of the reduced density matrix, form a tower of states, in many aspects similar to the tower of eigenstates of a coupled harmonic system. Actually, we know that when the matrix β is real, the above statement is exact [7,15]. It would be nice to construct creation and annihilation operators which would relate the eigenstates ofρ 2 in the same sense that they relate the eigenstates of a coupled harmonic system.
However, we know that this cannot be that simple. For example, we know that the "second excited" eigenstate of the reduced density matrix, which corresponds to different eigenvectors of the matrix Ξ, cannot be produced by the action of two "ordinary" creation operators on the ground state, as we showed in section 4.4.3. In general, the required operators have to be linear combinations of positions and momenta. They differ though from the "ordinary" ones, as the combination of momenta that appears in one of them cannot be the conjugate momentum of the combination of positions that appears.
Therefore, we search for annihilation operators of the form so that they annihilate the "ground" eigenstate Ψ 0 ofρ 2 (4.36). These operators should act on the "first excited" eigenstates Ψ 1 , which are given by equation (4.48), as Introducing the notation (v ) k = v k , equation (C.2) yields Bearing in mind that the eigenvectors of the matrix Ξ are normalized so that the above equation implies that We proved that if the state Ψ (x) is an eigenstate of the matrixρ 2 with eigenvalue λ, then the state A † i Ψ (x) is also an eigenstate with eigenvalue λξ i . Inductively, this means that the states are normalized eigenstates of the matrixρ 2 with eigenvalues given by equation (4.68). The eigenstates of the reduced density matrix ρ 2 can be trivially found, recalling their relation to the eigenstates ofρ 2 , which is given by equation (4.33).

D A Solvable Example
In this appendix we analyse a solvable example, in order to clarify the properties of the spectrum of the matrix M that we introduced in section 4.6, and furthermore to verify that the asymptotic form of the eigenvalues of the matrix Ξ for large squeezing parameters are indeed of the form that is predicted by the large squeezing expansion developed in section 5. Let us consider the case that all modes lie in a squeezed state with the same squeezing parameter z, and further assume that we study the system at a instant when the phase of the oscillation of all modes is the same and equal to 0. In this case, we have Naturally, when z = 0, W = Ω, as in the usual ground state calculation. It follows that the matrices γ and β assume the form where γ 0 and β 0 are the matrices γ and β in the case of the ground state. The first of the two equations implies that Finally, the above imply that the matrix M , which defined in equation (5.11) and is similar to the matrix M , assumes the form The eigenvalues of this matrix are , (D. 6) whereβ i are the eigenvalues of the matrix γ −1 0 β 0 . These eigenvalues come in pairs of the form (λ, 1/λ), since λ i+ λ i− = 1. The eigenvalues which are smaller than 1 are the λ i+ , and, thus, For |z| 1 we have (D. 8) This formula is in agreement with the large-squeezing expansion developed in section 5.
Similarly, for |z| 1 we obtain The correction to ξ 0 i is always non-negative.

E Entanglement Entropy through Rényi Entropies -a Toy Case
In section 4 we developed a method to calculate the entanglement entropy in harmonic systems with an arbitrary number of degrees of freedom that lie in a squeezed state. It would be nice if we could verify this method via an independent calculation in a non-trivial case.
There is an alternative method to calculate entanglement entropy via the so called entanglement Rényi entropies. In this appendix we use this alternative method as a verifying example in a special case that the reduced system contains two degrees of freedom and the reduced density matrix is complex, but has a specific form.
Rényi entropies constitute a family of entropies which extend the notion of Shannon's entropy. For a probability distribution p i , the Rényi entropy of order a is defined as and recover the entanglement entropy as the limit In the case we study, namely when the overall oscillatory system lies in a squeezed state, we know that the reduced density matrix is of the form where γ is a complex symmetric matrix, β is a Hermitian matrix. The normalization constant is given by c = det Re (γ − β)/π d 1 2 , where d is the number of degrees of freedom of the reduced system.
(E. 24) In the case that the matrixβ is real, e.g. when the overall system lies in the ground state, all matricesγ n andβ n are functions of a single matrix, namely ofβ, and these recursion relations can be solved as recursion relations for numbers. Their solution leads to simple explicit formulas for the entanglement Rényi entropies and the exact same formula for entanglement entropy found in [7]. However, in our case of study the recursion relation forγ n contains the two matricesβ andβ * , which in general do not commute. However, the above formulae can be solved in a special case, where the matrixβ is complex. Assume the case where the reduced system contains two degrees of freedom and β = β 0 I + β 2 σ 2 .

(E.27)
Since the initial condition for the recursion relation forγ n (E.22) isγ 1 = I, it follows that this recursion relation has the trivial solution γ n = a n I, (E. 28) where the coefficients a n obey a n+1 = a n + detβ 1 + a n .
(E. 29) This recursion relation has the solution a n = 1 − detβ 1 + ξ n 1 − ξ n , . (E.41) Using the form of the matricesγ n ,β n andΓ n from equations (E.31), (E.34) and (E.39) and putting everything together yields The Renyi entanglement entropy is defined by equation (E.3). It reads (E.43) The form of the trace Trρ a 2 clearly implies that lim n→1 Trρ n 2 = 1, as expected. It follows that It is a matter of algebra to show that In order to verify that this result is consistent with the general method that we developed in section 4, we need to solve the non-linear eigenvalue equation of the eigenvalues of the matrix M , namely det ξ 2βT − 2ξI +β = 0. (E.47) In our example, the matrixβ is given by equation (E.25). The above equation gives det ξ 2 β 0 − 2ξ + β 0 I + β 2 −ξ 2 + 1 σ 2 = 0. The last equation has four solutions, These indeed form two pairs of solutions that are inverse to each other. Namely ξ 1+ = 1/ξ 2− and ξ 2+ = 1/ξ 1− . The solutions that are smaller than 1 are the solutions ξ 1− and ξ 2− . It follows that the entanglement entropy reads It is a matter of tedious algebra to show that the above expression is identical to equation (E.46).

F Small-Squeezing Expansion
Similarly to the large squeezing expansion that we presented in section 5, we expand the parameter w as a series in the squeezing parameter. This reads The zeroth order term is real, time-independent and equal to the eigenfrequency of the mode. Unlike the case of the large squeezing expansion, all terms in the expansion of w (apart the zeroth order one) contain both a real and an imaginary part. We will use a similar notation to that we used in section 5. Equation (F.1) implies that the matrix W has an expansion of the form where W (i) are in general complex. It follows that its blocks have a similar expansion and the same holds for the matrices γ and β, In all these expansions, the small squeezing parameter z may be the squeezing parameter of a single mode or even a small parameter in terms of which the small squeezing parameters of all modes can be expressed. We would like to perform textbook first order perturbation theory to the spectrum of the matrix M . This would be simpler if the matrix M were Hermitian, at least at zeroth order, so that its eigenvectors are orthogonal. Actually, we can find a matrixM , which is similar to M and Hermitian. This readŝ The matrixM (0) is not only Hermitian but also block-diagonal and its eigenvectors are trivially constructed from the eigenvectors ofβ. Let x i be the eigenvectors of the matrixβ with corresponding eigenvalues equal toβ i , i.e. As expected the eigenvalues come in pairs of the form (λ, 1/λ). Indeed, λ i λ i = 1. The eigenvalues that are smaller than 1 are the ones corresponding to eigenvectors of the first kind, namely the λ i . Indeed, they coincide with the values of the parameters ξ in the original calculation by Srednicki [7]. Notice also that the matrixβ is real and symmetric, and thus its eigenvectors x i are real. It is a matter of algebra to show that Now we can apply perturbation theory to find the eigenvalues of the matrix M at first order. Considering that ξ i = ξ This implies that the entanglement entropy at a given time contains corrections which are first order in z, namely However, this is not the case for the mean entanglement entropy. The quantities ξ (1) i depend linearly on iωe −2iωt , as it results from equation (F.1). However, they are real, therefore they depend linearly on a combination of cos (2ωt) and sin (2ωt). As a result, the mean value of ξ (1) i vanishes and so does the correction of the mean entanglement entropy at first order. It follows that the ground state is a stationary point for the mean entanglement entropy within the space of squeezed states.