# Entanglement in harmonic systems at coherent states

Dimitrios Katsinis1,2,3 and Georgios Pastras2,4

(1Instituto de Física, Universidade de São Paulo, Rua do Matão Travessa 1371,
05508-090 São Paulo, SP, Brazil
2NCSR “Demokritos”, Institute of Nuclear and Particle Physics,
15310 Aghia Paraskevi, Attiki, Greece
3National and Kapodistrian University of Athens,
Department of Physics, 15784 Zografou, Attiki, Greece
4Laboratory for Manufacturing Systems and Automation, Department of Mechanical
Engineering and Aeronautics,
University of Patras, 26110 Patra, Greece
<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

It is well-known that entanglement entropy in field theory at its ground state
is dominated by an area law term, presenting a similarity to the entropy of
black holes. It is interesting to investigate whether this similarity can be
extended by showing that gravitational dynamics emerges from the first law of
entanglement thermodynamics. Answering this question requires the
specification of the modular Hamiltonian. Motivated by the above, we study
entanglement in the toy model of harmonic systems lying at any classicalmost
state, i.e. any coherent state. We specify explicitly the reduced density
matrix and its time-evolution, as well as the modular Hamiltonian.
Interestingly, the time evolution is unitary and we specify the effective
Hamiltonian which generates it. Our results provide the tools to investigate
the similarity between gravity and entanglement in discretized free scalar
field theory in the framework of [1].

###### Contents

  1. 1 Introduction
  2. 2 Generalized Harmonic Oscillator
  3. 3 A Pair of Harmonic Oscillators at a Coherent State
    1. 3.1 The Reduced Density Matrix
    2. 3.2 Time Evolution of the Reduced System
    3. 3.3 The Modular Hamiltonian
  4. 4 System of Harmonic Oscillators at Coherent State
    1. 4.1 The Reduced Density Matrix
    2. 4.2 Time Evolution of the Reduced System
    3. 4.3 The Modular Hamiltonian
  5. 5 Discussion
  6. A Coherent and Squeezed States of the Simple Harmonic Oscillator

## 1 Introduction

In 1993, in a seminal paper, Srednicki [1] showed that the entanglement
entropy in free massless scalar field theory at its ground state is
proportional to the area of the entangling surface. This apparent similarity
of entanglement entropy and black hole entropy gave rise to the question about
whether black hole entropy can be attributed to quantum entanglement between
the interior and the exterior of the black hole. Such a scenario would favour
a description of gravity as a quantum entropic force [2, 3].

The latter hypothesis is greatly supported by the holographic duality. In any
composite quantum system, the first law of entanglement thermodynamics
trivially holds. This relates the variation of entanglement entropy and the
mean value of the modular Hamiltonian111The modular Hamiltonian is defined as
the logarithm of the density matrix $H_{\mathrm{mod}}=-\ln\rho$ under a
variation of the state of the composite system

$\delta S_{\mathrm{EE}}=\delta\left<H_{\mathrm{mod}}\right>.$ (1.1)

This relation is true for the holographic realizations of the entanglement
entropy and the expected value of the modular Hamiltonian for spherical
entangling surfaces in the bulk, if and only if the small variation of the
state in the boundary corresponds to a small variation of the geometry in the
bulk, which obeys the linearized Einstein equation [4]. Therefore, the
holographically emergent gravity is at least not contradictory to a quantum
entropic gravitational force scenario.

Further investigation of the relation between gravity and quantum statistics
attributed to quantum entanglement requires the specification of the modular
Hamiltonian for non-trivial states of the underlying theory. This is an
extremely difficult task, since it requires the explicit calculation of the
reduced density matrix and not only its spectrum. In this manner, the Replica
trick and holographic calculations of relevant quantities, such as the
entanglement entropy, are disadvantageous, since they naturally lead to the
calculation of traces of integer powers of the reduced density matrix. This
data are adequate to reproduce in principle the spectrum of the latter but not
its eigenstates. On the other hand, the simplistic, brute force approach of
[1], although it applies to the non-realistic, non-interacting integrable toy
model of free scalar field theory, it has the advantage that it directly
calculates the reduced density matrix, and, thus, it can lead to the explicit
calculation of the modular Hamiltonian operator.

In this work, we provide the tools to extend the study of the similarity
between gravity and entanglement in the framework of harmonic systems, as
initiated in [1], calculating the modular Hamiltonian, as an operator, which
is a function of the positions and momenta of the local degrees of freedom,
and specifying the time evolution of the reduced system. This task is
performed for a general classicalmost state of the composite system, i.e. a
coherent state. In section 2, we develop the algebraic solution to a
generalized oscillator. This system is described by a Hamiltonian, which is a
quadratic function of both position and momentum, but with explicit time
dependence. This generalized oscillator is going to be the basic building
block of the modular Hamiltonian. In section 3 we find the modular Hamiltonian
and specify the time evolution of the reduced system in the simplest problem
of two coupled oscillators. In section 4, we extend the results of section 3
to an arbitrary harmonic system. Finally, in section 5 we discuss our results.
There is also an appendix, containing an application of the generalized
harmonic oscillator that we developed in section 2, in the specification of
the most general Gaussian states of the simple harmonic oscillator.

## 2 Generalized Harmonic Oscillator

In this section we generalize the usual algebraic solution of the simple
harmonic oscillator to construct a more general quadratic Hamiltonian, which
has explicit time-dependence. This problem has been considered a long time ago
[5, 6] from a completely different perspective. The time-dependent Schrödinger
equation has a tower of solutions, which resembles the Fock space of the
simple harmonic oscillator. More specifically, the lowest state of this tower
is a time-dependent Gaussian state. The entanglement in harmonic systems at
coherent states, which are in fact time-dependent Gaussian states, is
naturally described with the use of this tool.

The problem can be treated algebraically in a similar manner to the simple
harmonic oscillator [7, 8, 9]. Let us recall the main idea behind the
algebraic solution of the simple harmonic oscillator. The ground state is a
Gaussian state, with equal uncertainties in position of momentum. It is a
trivial fact that any Gaussian state is annihilated by an appropriate linear
combination of the position and momentum operators. This linear combination is
what we call the annihilation operator $a$; its Hermitian conjugate
$a^{\dagger}$ is the so-called creation operator. It turns out that the
commutator of these two operators is equal to the identity operator and
furthermore the Hamiltonian $H$ of the harmonic oscillator assumes the form
$\hbar\omega\left(a^{\dagger}a+1/2\right)$. This implies the well-known
commutation relations $\left[H,a\right]=\hbar\omega
a,\left[H,a^{\dagger}\right]=\hbar\omega a^{\dagger}$, which provide the
spectrum and eigenstates of the harmonic oscillator.

We want to generalize this procedure, but instead of the ground state of the
harmonic oscillator, we will consider the most general Gaussian state with
time-dependent parameters. For this reason, in our construction, we consider
as annihilation and creation operators $A$ and $A^{\dagger}$, the most general
linear combinations of the position and momentum operators, which satisfy the
commutation relation

$\left[\hat{A},\hat{A}^{\dagger}\right]=1.$ (2.1)

It is trivial to show that these operators read

$\displaystyle\hat{A}$
$\displaystyle:=\frac{1}{\sqrt{2m\hbar\omega^{R}}}\left[i\left(\hat{p}-p_{0}\right)+m\left(\omega^{R}+i\omega^{I}\right)\left(\hat{x}-x_{0}\right)\right],$
(2.2) $\displaystyle\hat{A}^{\dagger}$
$\displaystyle:=\frac{1}{\sqrt{2m\hbar\omega^{R}}}\left[-i\left(\hat{p}-p_{0}\right)+m\left(\omega^{R}-i\omega^{I}\right)\left(\hat{x}-x_{0}\right)\right],$
(2.3)

where $x_{0}$, $p_{0}$, $\omega^{R}$ and $\omega^{I}$ are arbitrary functions
of time. For the sake of simplicity, we assume that the parameter $m$ is
constant; notice that this parameter appears always in the form of the product
$m\omega^{R}$ or $m\omega^{I}$, and, thus its time dependence can be absorbed
in the time dependence of the parameters $\omega^{R}$ and $\omega^{I}$. One
could also add a global time-dependent phase to these operators, but this is
physically indifferent.

It is straightforward to show that the normalized state, which is annihilated
by $\hat{A}$, reads

$\left\langle
x\middle|\Phi_{0}\right\rangle=\left(\frac{m\omega^{R}}{\pi\hbar}\right)^{1/4}\exp\left[-\frac{m\left(\omega^{R}+i\omega^{I}\right)}{2\hbar}\left(x-x_{0}\right)^{2}+i\frac{p_{0}}{\hbar}\left(x-x_{0}\right)-\frac{i}{2}\varphi\right],$
(2.4)

where $\varphi$ is in general an arbitrary function of time.

In a trivial manner, the states

$\left|\Phi_{n}\right\rangle:=\left(\hat{A}^{\dagger}\right)^{n}\left|\Phi_{0}\right\rangle$
(2.5)

are eigenfunctions of the operator

$\hat{\tilde{H}}:=\hbar\omega^{R}\hat{A}^{\dagger}\hat{A},$ (2.6)

with eigenvalues equal to

$\hat{\tilde{H}}\left|\Phi_{n}\right\rangle=n\hbar\omega^{R}\left|\Phi_{n}\right\rangle.$
(2.7)

When expressed in terms of the position and momentum operators, the operator
$\hat{\tilde{H}}$ assumes the form

$\hat{\tilde{H}}=\frac{1}{2m}\left(\hat{p}-p_{0}\right)^{2}+\frac{1}{2}\omega^{I}\left[\left(\hat{x}-x_{0}\right)\left(\hat{p}-p_{0}\right)+\left(\hat{p}-p_{0}\right)\left(\hat{x}-x_{0}\right)\right]\\\
+\frac{1}{2}m\left[\left(\omega^{R}\right)^{2}+\left(\omega^{I}\right)^{2}\right]\left(\hat{x}-x_{0}\right)^{2}+\frac{\omega^{R}}{2}.$
(2.8)

Since the operators $A$ and $A^{\dagger}$ contain an explicit time dependence,
the operator $\hat{\tilde{H}}$ cannot emerge as the Hamiltonian in a time-
independent Schrödinger equation; The states $\Psi_{n}$ are not solutions of
the time-dependent Schrödinger equation
$i\hbar\partial_{t}\Psi=\hat{\tilde{H}}\Psi$.

In order to upgrade the states $\Phi_{n}$ to solutions of a time-dependent
Schrödinger equation, we need to find a time-dependent Hamiltonian $\hat{H}$,
which obeys specific commutation relations with the creation and annihilation
operators (2.2) and (2.3), namely,

$\displaystyle
i\hbar\dot{\hat{A}}=\left[\hat{H},\hat{A}\right]+\hbar\omega^{R}\hat{A},\qquad
i\hbar\dot{\hat{A}}^{\dagger}=\left[\hat{H},\hat{A}^{\dagger}\right]-\hbar\omega^{R}\hat{A}^{\dagger}.$
(2.9)

It is a matter of algebra to show that this Hamiltonian reads222Notice that
had we considered a time-dependent mass $m$, this Hamiltonian would read
$\hat{H}=\hbar\Bigg{[}\left(\omega^{R}+\frac{\partial_{t}{\left(m\omega^{I}\right)}}{2m\omega^{R}}\right)\left(\hat{A}^{\dagger}\hat{A}+\frac{1}{2}\right)+\frac{\partial_{t}{\left(m\left(\omega^{I}+i\omega^{R}\right)\right)}}{4m\omega^{R}}\hat{A}^{2}+\frac{\partial_{t}{\left(m\left(\omega^{I}-i\omega^{R}\right)\right)}}{4m\omega^{R}}\hat{A}^{\dagger
2}\\\
-\frac{\dot{p}_{0}+m\left(\omega^{I}-i\omega^{R}\right)\dot{x}_{0}}{\sqrt{2m\hbar\omega^{R}}}\hat{A}^{\dagger}-\frac{\dot{p}_{0}+m\left(\omega^{I}+i\omega^{R}\right)\dot{x}_{0}}{\sqrt{2m\hbar\omega^{R}}}\hat{A}+\frac{p_{0}\dot{x}_{0}}{\hbar}\Bigg{]}.$
(2.10)

$\hat{H}=\hbar\Bigg{[}\left(\omega^{R}+\frac{\dot{\omega}^{I}}{2\omega^{R}}\right)\left(\hat{A}^{\dagger}\hat{A}+\frac{1}{2}\right)+\frac{\dot{\omega}^{I}+i\dot{\omega}^{R}}{4\omega^{R}}\hat{A}^{2}+\frac{\dot{\omega}^{I}-i\dot{\omega}^{R}}{4\omega^{R}}\hat{A}^{\dagger
2}\\\
-\frac{\dot{p}_{0}+m\left(\omega^{I}-i\omega^{R}\right)\dot{x}_{0}}{\sqrt{2m\hbar\omega^{R}}}\hat{A}^{\dagger}-\frac{\dot{p}_{0}+m\left(\omega^{I}+i\omega^{R}\right)\dot{x}_{0}}{\sqrt{2m\hbar\omega^{R}}}\hat{A}+\frac{p_{0}\dot{x}_{0}}{\hbar}\Bigg{]}.$
(2.11)

Expressing this Hamiltonian in terms of the position and momentum operators
yields

$\hat{H}=\frac{1}{2m}\left(\hat{p}-p_{0}+m\dot{x}_{0}\right)^{2}+\frac{1}{2}\mathcal{A}\left[\left(\hat{x}-x_{0}\right)\hat{p}+\hat{p}\left(\hat{x}-x_{0}\right)\right]+\frac{1}{2}m\tilde{\Omega}^{2}\left(\hat{x}-x_{0}\right)^{2}\\\
-\left(p_{0}\mathcal{A}+\dot{p}_{0}\right)\left(\hat{x}-x_{0}\right)+\dot{x}_{0}\left(p_{0}-\frac{1}{2}m\dot{x}_{0}\right),$
(2.12)

where $\tilde{\Omega}^{2}$ and $\mathcal{A}$ are given by

$\tilde{\Omega}^{2}=\left(\omega^{R}\right)^{2}+\dot{\omega}^{I}-\left(\omega^{I}\right)^{2}+2\omega^{I}\mathcal{A},\qquad\mathcal{A}=\omega^{I}-\frac{\dot{\omega}^{R}}{2\omega^{R}}.$
(2.13)

The Gaussian state (2.4) does not solve the Schrödinger equation

$i\hbar\partial_{t}\left|\Psi\right\rangle=\hat{H}\left|\Psi\right\rangle$
(2.14)

for any function $\varphi\left(t\right)$. This function is determined by the
Schrödinger equation. Plugging (2.16) into (2.14) yields

$\dot{\varphi}=\omega^{R}.$ (2.15)

The simplicity of this result is due to the specific choice of the $c$-number
terms that we made in the definition of the Hamiltonian (2.11). We define this
solution of Schrödinger equation as the state $\left|\Psi_{0}\right\rangle$,

$\left\langle
x\middle|\Psi_{0}\right\rangle=\left(\frac{m\omega^{R}}{\pi\hbar}\right)^{1/4}\exp\left[-\frac{m\left(\omega^{R}+i\omega^{I}\right)}{2\hbar}\left(x-x_{0}\right)^{2}+i\frac{p_{0}}{\hbar}\left(x-x_{0}\right)-\frac{i}{2}\phi\right],$
(2.16)

where

$\phi:=\int dt\omega^{R}.$ (2.17)

Since we have a solution of the Schrödinger equation (2.14), we can construct
an infinite tower of solutions algebraically, in a similar manner to the
simple harmonic oscillator. Assume that $\left|\Psi\right\rangle$ is a
solution of the Schrödinger equation (2.14). Then, the state

$\left|\Psi^{\prime}\right\rangle=e^{-i\phi}\hat{A}^{\dagger}\left|\Psi\right\rangle,$
(2.18)

where $\phi$ is given by (2.17), is another solution of the Schrödinger
equation. In particular, in view of equation (2.9) it follows that

$i\hbar\partial_{t}\left|\Psi^{\prime}\right\rangle=\hbar\omega^{R}e^{-i\phi}\hat{A}^{\dagger}\left|\Psi^{\prime}\right\rangle+i\hbar
e^{-i\phi}\dot{\hat{A}}^{\dagger}\left|\Psi^{\prime}\right\rangle+e^{-i\phi}\hat{A}^{\dagger}\hat{H}\left|\Psi^{\prime}\right\rangle=\hat{H}\left|\Psi^{\prime}\right\rangle.$
(2.19)

Similarly to the case of the simple harmonic oscillator, we define the states
$\Psi_{n}$ as,

$\left|\Psi_{n}\right\rangle=\frac{e^{-in\phi}}{\sqrt{n!}}\left(\hat{A}^{\dagger}\right)^{n}\left|\Psi_{0}\right\rangle.$
(2.20)

These states obey

$\hat{A}\left|\Psi_{n}\right\rangle=e^{-i\phi}\sqrt{n}\left|\Psi_{n-1}\right\rangle,\qquad\hat{A}^{\dagger}\left|\Psi_{n}\right\rangle=e^{i\phi}\sqrt{n+1}\left|\Psi_{n+1}\right\rangle$
(2.21)

and

$\hat{A}^{\dagger}\hat{A}\left|\Psi_{n}\right\rangle=n\left|\Psi_{n}\right\rangle.$
(2.22)

Recalling the form of the position representation of the wavefunctions of the
simple harmonic oscillator, it is not difficult to show that, the position
representation of these solutions reads

$\left\langle
x\middle|\Psi_{n}\right\rangle=\frac{1}{\sqrt{2^{n}n!}}\left(\frac{m\omega^{R}}{\pi\hbar}\right)^{1/4}H_{n}\left(\sqrt{\frac{m\omega^{R}}{\hbar}}\left(x-x_{0}\right)\right)\exp\left[-\frac{m\omega^{R}}{2\hbar}\left(x-x_{0}\right)^{2}\right]\\\
\times\exp\left[-i\frac{m\omega^{I}}{2\hbar}\left(x-x_{0}\right)^{2}+i\frac{p_{0}}{\hbar}\left(x-x_{0}\right)-i\left(n+\frac{1}{2}\right)\phi\right],$
(2.23)

where $H_{n}$ is the Hermite polynomial of order $n$. This justifies the
choice of $c$-number terms of the Hamiltonian (2.11). They have been
appropriately defined so that the function $\phi$ plays an analogous role to
$\omega t$ in the simple harmonic oscillator.

At any given time instant, the wave-functions (2.23) resemble the usual Fock
space states, i.e.,

$\psi_{n}\left(x\right)=\frac{1}{\sqrt{2^{n}n!}}\left(\frac{m\omega^{R}}{\pi\hbar}\right)^{1/4}H_{n}\left(\sqrt{\frac{m\omega^{R}}{\hbar}}x\right)\exp\left[-\frac{m\omega^{R}}{2\hbar}x^{2}\right]$
(2.24)

of the simple harmonic oscillator with eigenfrequency $\omega^{R}$, but they
differ in the following three aspects:

  * •

they have their argument shifted by $x_{0}$

  * •

they are multiplied by a position dependent phase,

$\varphi\left(x\right)=-\frac{m\omega^{I}}{2\hbar}\left(x-x_{0}\right)^{2}+\frac{p_{0}}{\hbar}\left(x-x_{0}\right),$
(2.25)

which is _identical for all states_

  * •

they are multiplied by _a position independent phase_

$\phi_{n}=-\left(n+\frac{1}{2}\right)\phi,$ (2.26)

which depends on the state.

That is to say,

$\left\langle
x\middle|\Psi_{n}\right\rangle=\psi_{n}\left(x-x_{0}\right)\exp\left[i\varphi\left(x\right)\right]\exp\left(i\phi_{n}\right).$
(2.27)

It follows that although these states are time-dependent, at any given time
instant, they form an orthonormal basis for square integrable wavefunctions.

Orthonormality is a direct consequence of the orthonormality of the Fock space
states,

$\begin{split}\left\langle\Psi_{m}\middle|\Psi_{n}\right\rangle&=\exp\left(i\phi_{n}-i\phi_{m}\right)\int
dx\bar{\psi}_{m}\left(x-x_{0}\right)\exp\left[-i\varphi\left(x\right)\right]\psi_{n}\left(x-x_{0}\right)\exp\left[i\varphi\left(x\right)\right]\\\
&=\exp\left(i\phi_{n}-i\phi_{m}\right)\int
dx\bar{\psi}_{m}\left(x-x_{0}\right)\psi_{n}\left(x-x_{0}\right)\\\
&=\exp\left(i\phi_{n}-i\phi_{m}\right)\int
dx\bar{\psi}_{m}\left(x\right)\psi_{n}\left(x\right)\\\
&=\exp\left(i\phi_{n}-i\phi_{m}\right)\delta_{mn}=\delta_{mn}.\end{split}$
(2.28)

Regarding the completeness of the states, assume that we want to expand some
square integrable function $f\left(x\right)$ as a linear combination of
$\Psi_{n}$ at $t=t_{0}$. It is a well-known fact that the Fock space states
$\psi_{n}\left(x\right)$ form a basis for square integrable wavefunctions. A
shift in the argument of a wavefunction does not alter the fact that it is
square integrable. The same holds for multiplication of a wavefunction with a
pure phase. Therefore, we can always expand the wavefunction
$f\left(x+x_{0}\right)\exp\left[-i\varphi\left(x+x_{0}\right)\right]$ in the
Fock space states,

$f\left(x+x_{0}\right)\exp\left[-i\varphi\left(x+x_{0}\right)\right]=\sum_{n=0}^{\infty}c_{n}\psi_{n}\left(x\right).$
(2.29)

This directly implies that

$\begin{split}f\left(x\right)&=\sum_{n=0}^{\infty}c_{n}\psi_{n}\left(x-x_{0}\right)\exp\left[i\varphi\left(x\right)\right]\\\
&=\sum_{n=0}^{\infty}\left[c_{n}\exp\left(-i\phi_{n}\right)\right]\psi_{n}\left(x-x_{0}\right)\exp\left[i\varphi\left(x\right)\right]\exp\left(i\phi_{n}\right)\\\
&=\sum_{n=0}^{\infty}\left[c_{n}\exp\left(-i\phi_{n}\right)\right]\left\langle
x\middle|\Psi_{n}\right\rangle,\end{split}$ (2.30)

i.e. any square integrable function can be expanded as a linear combination of
the states $\left|\Psi_{n}\right\rangle$. Summing up, _at any time instant,
the states ( 2.23) form an orthonormal basis for square integrable
wavefunctions_.

In the form of a bilinear kernel, the completeness of these wave-functions
reads

$\left\langle
x\middle|\hat{I}\middle|x^{\prime}\right\rangle=\sum_{n=0}^{\infty}\left\langle
x\middle|\Psi_{n}\right\rangle\left\langle\Psi_{n}\middle|x^{\prime}\right\rangle$
(2.31)

and implies that

$\left\langle
x\middle|\hat{A}^{\dagger}\hat{A}\middle|x^{\prime}\right\rangle=\sum_{n=0}^{\infty}n\left\langle
x\middle|\Psi_{n}\right\rangle\left\langle\Psi_{n}\middle|x^{\prime}\right\rangle.$
(2.32)

Before, we conclude this section, we would like to make an important comment.
One can easily observe in the form of the wavefunctions (2.23) that they
actually do not depend separately on $m$ and $\omega^{R}+i\omega^{I}$, but
only on their product. Therefore, there is a whole family of Hamiltonians, one
for each value of $m$, which have exactly the same solutions, and, thus, they
are completely equivalent, as long as the time evolution of the system is
concerned. These Hamiltonians may look quite different. Consider for example
the case that $\omega^{R}$, $\omega^{I}$, $x_{0}$ and $p_{0}$ have been chosen
so that $\mathcal{A}=0$ and $\dot{x}_{0}=p_{0}$. Then, if we write down the
Hamiltonian with $m=1$, there will be no linear term in $\hat{p}$. However, if
we write it down with any other value of $m$, such a term will be present.

## 3 A Pair of Harmonic Oscillators at a Coherent State

Entanglement entropy of coherent states was studied in the framework of
quantum field theory in [10] using the formalism of [11]. It turns out that
the entanglement entropy is identical to that in the case of the ground state.
This is somehow expected, since the coherent states, including the ground
state, are the classicalmost states of the system. These findings are in line
to the analysis of the next two sections.

### 3.1 The Reduced Density Matrix

Before we proceed to study more complex systems, we first consider a system of
two coupled harmonic oscillators. As subsystem $A$ we consider the harmonic
oscillator, which is described by the coordinate $x_{2}$. The complementary
system $A^{C}$ is consequently identical to the other oscillator, which is
described by the coordinate $x_{1}$. The corresponding canonical momenta are
named $p_{2}$ and $p_{1}$ respectively. Without loss of generality, the mass
of both oscillators is considered to be equal to one. The Hamiltonian of this
system is given by

$H=\frac{1}{2}\left[{p_{1}^{2}+p_{2}^{2}+{k_{0}}\left({x_{1}^{2}+x_{2}^{2}}\right)+{k_{1}}{\left({{x_{1}}-{x_{2}}}\right)}^{2}}\right].$
(3.1)

We define the canonical coordinates as usual

$x_{\pm}:=\frac{x_{1}\pm x_{2}}{\sqrt{2}},\quad p_{\pm}:=\frac{p_{1}\pm
p_{2}}{\sqrt{2}}.$ (3.2)

Then, the Hamiltonian assumes the form

$H=\frac{1}{2}\left({p_{+}^{2}+p_{-}^{2}+\omega_{+}^{2}x_{+}^{2}+\omega_{-}^{2}x_{-}^{2}}\right),$
(3.3)

where $\omega_{+}=\sqrt{k_{0}}$ and $\omega_{-}=\sqrt{k_{0}+2k_{1}}$. This
Hamiltonian describes two decoupled oscillators, one for each normal mode of
the system.

We consider states of the composite system, where _all_ normal modes lie at
coherent states. The wave function of a coherent state of an oscillator with
eigenfrequency $\omega$ reads (see e.g. Appendix A)

$\left\langle
x\middle|\Psi\left(\lambda\right)\right\rangle=\left(\frac{\omega}{\pi}\right)^{1/4}\exp\left[-\frac{\omega}{2}\left(x-\lambda\right)^{2}+\frac{\omega}{4}\left(\lambda^{2}-\left|\lambda\right|^{2}\right)-i\frac{\omega}{2}t\right].$
(3.4)

The parameter $\lambda$ evolves with time as

$\lambda=\lambda_{0}e^{-i\omega t}$ (3.5)

and determines the mean values of position and momentum

$\left\langle x\right\rangle=\mathrm{Re}\lambda,\quad\left\langle
p\right\rangle=\omega\mathrm{Im}\lambda.$ (3.6)

It directly follows that the wavefunction of the composite system when both
normal modes lie at a coherent state reads

$\left\langle
x_{+},x_{-}\middle|\Psi\left(\lambda_{+},\lambda_{-}\right)\right\rangle\\\
=\left(\frac{\omega_{+}\omega_{-}}{\pi^{2}}\right)^{1/4}\exp\left[-\frac{\omega_{+}}{2}\left(x_{+}-\lambda_{+}\right)^{2}+\frac{\omega_{+}}{4}\left(\lambda_{+}^{2}-\left|\lambda_{+}\right|^{2}\right)-i\frac{\omega_{+}}{2}t\right]\\\
\times\exp\left[-\frac{\omega_{-}}{2}\left(x_{-}-\lambda_{-}\right)^{2}+\frac{\omega_{-}}{4}\left(\lambda_{-}^{2}-\left|\lambda_{-}\right|^{2}\right)-i\frac{\omega_{-}}{2}t\right].$
(3.7)

The density matrix of the composite system reads

$\left\langle
x_{+},x_{-}\middle|\hat{\rho}\left(\lambda_{+},\lambda_{-}\right)\middle|x_{+}^{\prime},x_{-}^{\prime}\right\rangle\\\
=\left(\frac{\omega_{+}\omega_{-}}{\pi^{2}}\right)^{1/2}\exp\left[-\frac{\omega_{+}}{2}\left(x_{+}-\lambda_{+}\right)^{2}-\frac{\omega_{+}}{2}\left(x^{\prime}_{+}-\bar{\lambda}_{+}\right)^{2}+\frac{\omega_{+}}{4}\left(\lambda_{+}-\bar{\lambda}_{+}\right)^{2}\right]\\\
\times\exp\left[-\frac{\omega_{-}}{2}\left(x_{-}-\lambda_{-}\right)^{2}-\frac{\omega_{-}}{2}\left(x^{\prime}_{-}-\bar{\lambda}_{-}\right)^{2}+\frac{\omega_{-}}{4}\left(\lambda_{-}-\bar{\lambda}_{-}\right)^{2}\right].$
(3.8)

In order to obtain the reduced density matrix of the second oscillator, which
comprises the subsystem $A$, we have to trace out the degree of freedom
corresponding to the first harmonic oscillator,

$\left\langle
x_{2}\middle|\hat{\rho}_{2}\middle|x_{2}^{\prime}\right\rangle=\int
dx_{1}\left\langle
x_{1},x_{2}\middle|\hat{\rho}\middle|x_{1},x_{2}^{\prime}\right\rangle.$ (3.9)

This is a simple task of performing a Gaussian integral, but first we need to
revert to the original coordinates $x_{1}$ and $x_{2}$ from the normal mode
coordinates $x_{\pm}$. It is also convenient to define the moduli
$\lambda_{1}$ and $\lambda_{2}$ so that

$\lambda_{\pm}=\frac{1}{\sqrt{2}}\left(\lambda_{1}\pm\lambda_{2}\right),$
(3.10)

in direct analogy to the definition of the normal mode coordinates (3.2). Bear
in mind that the “local” coherent state parameters $\lambda_{1}$ and
$\lambda_{2}$ do not have the simple time evolution of those of the normal
modes. The reduced density matrix reads

$\left\langle
x_{2}\middle|\hat{\rho}_{2}\left(\lambda_{1},\lambda_{2}\right)\middle|x_{2}^{\prime}\right\rangle=\left(\frac{\gamma-\beta}{\pi}\right)^{1/2}\exp\left[-\frac{\gamma}{2}\left(y_{2}^{2}+y^{\prime
2}_{2}\right)+\beta
y_{2}y^{\prime}_{2}+i\delta\left(y_{2}-y^{\prime}_{2}\right)\right],$ (3.11)

where we used the arguments $\lambda_{1}$ and $\lambda_{2}$ to indicate that
this is the reduced density matrix that corresponds to the coherent state
(3.7) and

$\displaystyle y_{2}$ $\displaystyle=x_{2}-\mathrm{Re}\lambda_{2},\qquad
y^{\prime}_{2}=x^{\prime}_{2}-\mathrm{Re}\lambda_{2},$ (3.12)
$\displaystyle\beta$
$\displaystyle=\frac{\left(\omega_{+}-\omega_{-}\right)^{2}}{4\left(\omega_{+}+\omega_{-}\right)},$
(3.13) $\displaystyle\gamma$
$\displaystyle=\frac{\omega_{+}^{2}+6\omega_{+}\omega_{-}+\omega_{-}^{2}}{4\left(\omega_{+}+\omega_{-}\right)},$
(3.14) $\displaystyle\delta$
$\displaystyle=\frac{1}{2}\left[\left(\omega_{+}-\omega_{-}\right)\mathrm{Im}\lambda_{1}+\left(\omega_{+}+\omega_{-}\right)\mathrm{Im}\lambda_{2}\right].$
(3.15)

The spectrum of a density matrix $\hat{\rho}$ is in principle completely
determined by the sequence of the traces
$\mathrm{Tr}\left(\hat{\rho}^{n}\right)$ for all $n\in\mathbb{N}$. As the
parameters $\beta$ and $\gamma$ _do not depend_ on the values of values of
$\lambda_{1}$ and $\lambda_{2}$, it is evident that

$\mathrm{Tr}\left[\hat{\rho}_{2}^{n}\left(\lambda_{1},\lambda_{2}\right)\right]=\mathrm{Tr}\left[\hat{\rho}_{2}^{n}\left(0,0\right)\right],\quad\forall
n\in\mathbb{N}.$ (3.16)

Thus, the spectrum of the reduced density matrix coincides with the spectrum
of the reduced density matrix when the overall system lies at its ground
state. It follows that the entanglement entropy is also identical to the
entanglement entropy in the case of the ground state. Nevertheless, this
neither implies that the former reduced density matrix coincides with the
latter nor that their dynamics coincide.

We can verify the above and obtain the spectrum of the reduced density matrix
in a more direct way, which also provides its eigenfunctions. Mehler’s formula
for Hermite polynomials states that

$\sum_{n=0}^{\infty}\frac{1}{n!}\left(\frac{\xi}{2}\right)^{n}H_{n}\left(x\right)H_{n}\left(y\right)\exp\left[-\frac{1}{2}\left(x^{2}+y^{2}\right)\right]\\\
=\frac{1}{\sqrt{1-\xi^{2}}}\exp\left[-\frac{\left(1+\xi^{2}\right)\left(x^{2}+y^{2}\right)-4\xi
xy}{2\left(1-\xi^{2}\right)}\right].$ (3.17)

So identifying

$\xi=\frac{\beta}{\gamma+\alpha},\qquad x=\sqrt{\alpha}y_{2},\qquad
y=\sqrt{\alpha}y_{2}^{\prime},$ (3.18)

where we have defined

$\alpha:=\sqrt{\gamma^{2}-\beta^{2}}=\sqrt{\omega_{+}\omega_{-}},$ (3.19)

we may directly apply Mehler’s formula to the reduced density matrix (3.11) in
order to re-express it as

$\left\langle
x_{2}\middle|\hat{\rho}_{2}\middle|x_{2}^{\prime}\right\rangle=\sqrt{\frac{\alpha}{\pi}}\sum_{n=0}^{\infty}\frac{1-\xi}{n!}\left(\frac{\xi}{2}\right)^{n}H_{n}\left(\sqrt{\alpha}y_{2}\right)H_{n}\left(\sqrt{\alpha}y^{\prime}_{2}\right)\\\
\times\exp\left[-\frac{\alpha}{2}\left(y_{2}^{2}+y^{\prime
2}_{2}\right)+i\delta\left(y_{2}-y^{\prime}_{2}\right)\right].$ (3.20)

As a result

$\left\langle
x_{2}\middle|\hat{\rho}_{2}\middle|x_{2}^{\prime}\right\rangle=\left(1-\xi\right)\sum_{n=0}^{\infty}\xi^{n}\psi_{n}\left(x_{2}-\mathrm{Re}\lambda_{2}\right)\psi_{n}\left(x^{\prime}_{2}-\mathrm{Re}\lambda_{2}\right)\exp\left[i\delta\left(x_{2}-x^{\prime}_{2}\right)\right],$
(3.21)

where $\psi_{n}$ are the energy eigenstates of an effective simple harmonic
oscillator with eigenfrequency equal to $\omega_{\text{eff}}=\alpha$. We may
express the reduced density matrix as

$\left\langle
x_{2}\middle|\hat{\rho}_{2}\middle|x_{2}^{\prime}\right\rangle=\left(1-\xi\right)\sum_{n=0}^{\infty}\xi^{n}\Psi_{n}\left(x_{2}\right)\bar{\Psi}_{n}\left(x^{\prime}_{2}\right),$
(3.22)

where

$\Psi_{n}\left(x\right)=\psi_{n}\left(x-\mathrm{Re}\lambda_{2}\right)\exp\left[i\delta\left(x-\mathrm{Re}\lambda_{2}\right)-i\alpha
t\left(n+\frac{1}{2}\right)\right].$ (3.23)

The orthonormality of the eigenstates of the harmonic oscillator, combined
with the fact that the phase $\delta$ is _the same_ for all states $\Psi_{n}$
implies that the states $\Psi_{n}$ are also orthonormal. It follows that
equation (3.22) provides the expansion of the reduced density matrix into its
eigenfunctions.

It is obvious that the eigenvalues of the reduced density matrix

$p_{n}=\left(1-\xi\right)\xi^{n}$ (3.24)

are independent of $\lambda_{\pm}$, implying that they are time-independent
and equal to the eigenvalues of the reduced system when the composite system
lies in the vacuum state.

### 3.2 Time Evolution of the Reduced System

The eigenfunctions (3.23) of the reduced density matrix (3.11) belong to the
family of wave-functions (2.23) that solve the generalized oscillator system,
which we presented in Section 2. In particular they are identical under the
identification

$m\omega^{R}=\alpha,\quad\omega^{I}=0,\quad x_{0}=\mathrm{Re}\lambda_{2},\quad
p_{0}=\delta,$ (3.25)

where $\hbar$ is set to one. Therefore, all the eigenfunctions (3.23) are
solutions of the same time-dependent Schrödinger equation, where the
corresponding Hamiltonian is given by (2.12), i.e.

$\hat{H}_{2}=\frac{1}{2m}\left(\hat{p}-\delta+m\mathrm{Re}\dot{\lambda}_{2}\right)^{2}+\frac{1}{2m}\alpha^{2}\left(\hat{x}-\mathrm{Re}\lambda_{2}\right)^{2}\\\
-\dot{\delta}\left(\hat{x}-\mathrm{Re}\lambda_{2}\right)+\mathrm{Re}\dot{\lambda}_{2}\left(\delta-\frac{1}{2}m\mathrm{Re}\dot{\lambda}_{2}\right).$
(3.26)

It is interesting to find the time derivative of the function
$x_{0}\left(t\right)$. Bearing in mind that the coherent state parameters for
the normal modes have a very simple time dependence, namely
$\lambda_{\pm}\left(t\right)=\lambda_{\pm}\left(0\right)e^{-i\omega_{\pm}t}$,
it is not difficult to show that

$\begin{split}\dot{x}_{0}&=\mathrm{Re}\dot{\lambda}_{2}=\frac{1}{\sqrt{2}}\mathrm{Re}\left(\dot{\lambda}_{+}-\dot{\lambda}_{-}\right)\\\
&=\frac{1}{\sqrt{2}}\mathrm{Re}\left(-i\omega_{+}\lambda_{+}+i\omega_{-}\lambda_{-}\right)\\\
&=\frac{1}{2}\left[\omega_{+}\mathrm{Im}\left(\lambda_{1}+\lambda_{2}\right)-\omega_{-}\mathrm{Im}\left(\lambda_{1}-\lambda_{2}\right)\right]=\delta=p_{0},\end{split}$
(3.27)

This implies that the Hamiltonian (3.26) is simplified if we chose $m=1$.
Then, it reads

$\hat{H}_{2}=\frac{1}{2}\hat{p}_{2}^{2}+\frac{1}{2}\alpha^{2}\left(\hat{x}_{2}-\mathrm{Re}\lambda_{2}-\frac{\dot{\delta}}{\alpha^{2}}\right)^{2}+\frac{1}{2}\delta^{2}-\frac{1}{2}\frac{\dot{\delta}^{2}}{\alpha^{2}}.$
(3.28)

After tracing out the first harmonic oscillator, the second one comprises an
open quantum system. The time evolution of such a system is non-unitary in
general. Nevertheless, since _all_ the eigenstates of the reduced density
matrix obey the equation

$i\hbar\dot{\Psi}_{n}=\hat{H}_{2}\Psi_{n},$ (3.29)

where $\hat{H}_{2}$ is given by (3.28), equation (3.22) directly implies that

$i\hbar\dot{\hat{\rho}}_{2}=\left[\hat{H}_{2},\hat{\rho}_{2}\right],$ (3.30)

meaning that when the overall system lies in a coherent state, the time
evolution of the reduced system _is unitary_.

How can we conceive this effective Hamiltonian at a classical level? This
Hamiltonian describes a force from an effective spring with time-independent
spring constant equal to $a^{2}$, but time-dependent equilibrium position
which is equal to $x_{0}+\dot{p}_{0}/a^{2}$, i.e. the spring force reads

$F_{\mathrm{eff}}=-a^{2}\left(x_{2}-x_{0}\right)-\dot{p}_{0}.$ (3.31)

If we associate the coherent state to a classical motion, so that the
classical values of the coordinates are given by their mean values of their
quantum counterparts, which are described by the coherent state, (this is
well-known to be indeed a solution of the classical equations of motion),
then, this force will simulate exactly the force from both real springs that
exert forces on the second oscillator, i.e.

$F_{\mathrm{eff}}=-k_{0}x_{2}-k_{1}\left(x_{2}-x_{1}\right)=F_{2},$ (3.32)

including the one which connects it to the invisible first oscillator.

So classically this effective Hamiltonian is a nice one-particle model that
recovers the correct motion of the second observable oscillator as part of the
composite system, accumulating the effect from the first invisible oscillator
in the time-dependence of the equilibrium position of the single effective
spring. Although this effective Hamiltonian _does depend on the state of the
overall system_ , the interesting observation is that this dependence appears
solely in the effective equilibrium position of the effective spring and not
at all at its constant. No matter what is the coherent state that we prepare
the composite system, the constant of the effective spring is the same; only
the specific time-dependence of its equilibrium position is altered.

### 3.3 The Modular Hamiltonian

We recall that the modular Hamiltonian is defined as the logarithm of the
reduced density matrix

$\hat{\rho}=e^{-\hat{H}_{M}}.$ (3.33)

So employing equation (3.22), we obtain

$\left\langle
x^{\prime}_{2}\middle|H_{M}\middle|x_{2}\right\rangle=-\sum_{n=1}^{\infty}\left(n\ln\xi+\ln\left(1-\xi\right)\right)\Psi_{n}\left(x_{2}\right)\bar{\Psi}_{n}\left(x^{\prime}_{2}\right).$
(3.34)

As a direct consequence of the relations (2.31) and (2.32) the modular
Hamiltonian assumes the form

$\hat{H}_{M}=-\ln\xi\hat{A}^{\dagger}\hat{A}-\ln\left(1-\xi\right)\hat{I},$
(3.35)

where the parameters in the generalized creation and annihilation operators
are given by equation (3.25).

Employing equation (2.8), we may write the modular Hamiltonian in terms of the
position and momentum operators as,

$\hat{H}_{M}=-\ln\xi\left(\frac{1}{2m\alpha}\left(\hat{p}_{2}-p_{0}\right)^{2}+\frac{\alpha}{2m}\left(\hat{x}_{2}-x_{0}\right)^{2}\right)+\frac{1}{2}\ln\xi-\ln\left(1-\xi\right).$
(3.36)

This Hamiltonian, when seen as a classical Hamiltonian, reproduces the correct
time evolution of the classical system, like the effective Hamiltonian that we
presented in Section 3.2, only when $\ln\xi=\alpha$. However, we should not
expect this to happen. The physical meaning of the modular Hamiltonian is
rather thermodynamic and assumes form through the first law of entanglement
thermodynamics.

## 4 System of Harmonic Oscillators at Coherent State

In this section, we treat an arbitrary system of coupled oscillators. The
analysis follows closely the original presentation of [1].

### 4.1 The Reduced Density Matrix

Assume a system of $N$ coupled harmonic oscillators described by the quadratic
Hamiltonian

$H=\frac{1}{2}\mathbf{p}^{T}\mathbf{p}+\frac{1}{2}\mathbf{x}^{T}K\mathbf{x},$
(4.1)

where $\mathbf{x}$ and $\mathbf{p}$ stand for $N$-dimensional column matrices
containing the coordinates and momenta of the $N$ harmonic oscillators. The
matrix $K$ is symmetric and has positive eigenvalues, so that the system is
stable. As $K$ is positively defined, we can define the matrix $\Omega$ as the
square root of $K$, which has all its eigenvalues positive.

In the following, we consider as subsystem $A$ the $N-n$ oscillators, which
are described by the coordinates $x_{i}$ with $i>n$. The complementary
subsystem $A^{C}$ comprises of the $n$ oscillators, which are described by the
coordinates $x_{i}$ with $i\leq n$.

Following the usual normal mode analysis, there is an orthogonal
transformation $O$, relating the coordinates $\mathbf{x}$ to the normal
coordinates $\tilde{\mathbf{x}}$, which diagonalizes the matrix $K$, reducing
the system to a set of decoupled harmonic oscillators, one for each normal
mode. In other words,

$\tilde{\mathbf{x}}=O\mathbf{x},\quad\tilde{\mathbf{p}}=O\mathbf{p},$ (4.2)

such that

$H=\frac{1}{2}\tilde{\mathbf{p}}^{T}\tilde{\mathbf{p}}+\frac{1}{2}\tilde{\mathbf{x}}^{T}K^{D}\tilde{\mathbf{x}}.$
(4.3)

The diagonalized matrix $K^{D}$

${K^{D}_{ij}}\equiv{\omega_{i}}^{2}{\delta_{ij}},$ (4.4)

is related to the initial matrix $K$ as,

$K={O^{T}}{K^{D}}O.$ (4.5)

As in Section 3, we want to study states of the overall system, where _all_
normal modes lie in a coherent state. Such a state reads

$\left\langle\mathbf{x}\middle|\Psi\left(\tilde{\boldsymbol{\lambda}}\right)\right\rangle=\sqrt[4]{\frac{\det\Omega_{D}}{\pi^{N}}}\exp\left[-\frac{1}{2}\left(\tilde{\mathbf{x}}-\tilde{\boldsymbol{\lambda}}\right)^{T}\Omega^{D}\left(\tilde{\mathbf{x}}-\tilde{\boldsymbol{\lambda}}\right)\right.\\\
\left.+\frac{1}{4}\left(\tilde{\boldsymbol{\lambda}}^{T}\Omega^{D}\tilde{\boldsymbol{\lambda}}-\tilde{\boldsymbol{\lambda}}^{\dagger}\Omega^{D}\tilde{\boldsymbol{\lambda}}\right)-\frac{i}{2}\mathrm{Tr}\Omega^{D}t\right],$
(4.6)

where $\Omega^{D}_{ij}=\omega_{i}\delta_{ij}$, is the square root of $K^{D}$.
In an obvious manner, it is connected to the matrix $\Omega$ as

$\Omega=O^{T}\Omega^{D}O.$ (4.7)

The column matrix $\tilde{\boldsymbol{\lambda}}$ contains the complex
“canonical” moduli of the coherent states of the $N$ normal modes. It evolves
with time as

$\tilde{\boldsymbol{\lambda}}\left(t\right)=e^{-i\Omega^{D}t}\tilde{\boldsymbol{\lambda}}\left(0\right).$
(4.8)

We define the moduli $\boldsymbol{\lambda}$, so that they are connected to the
“canonical” moduli as

$\boldsymbol{\lambda}=O^{T}\tilde{\boldsymbol{\lambda}},$ (4.9)

i.e. in exactly the same way as the original coordinates are related to the
canonical ones. We rewrite the wavefunction in terms of the original
coordinates and the moduli $\boldsymbol{\lambda}$ to get

$\left\langle\mathbf{x}\middle|\Psi\left(\boldsymbol{\lambda}\right)\right\rangle=\sqrt[4]{\frac{\det\Omega}{\pi^{N}}}\exp\left[-\frac{1}{2}\left(\mathbf{x}-\boldsymbol{\lambda}\right)^{T}\Omega\left(\mathbf{x}-\boldsymbol{\lambda}\right)+\frac{1}{4}\left(\boldsymbol{\lambda}^{T}\Omega\boldsymbol{\lambda}-\boldsymbol{\lambda}^{\dagger}\Omega\boldsymbol{\lambda}\right)-\frac{i}{2}\mathrm{Tr}\Omega
t\right].$ (4.10)

The density matrix reads

$\left\langle\mathbf{x}\middle|\rho\left(\boldsymbol{\lambda}\right)\middle|\mathbf{x}^{\prime}\right\rangle=\sqrt{\frac{\det\Omega}{\pi^{N}}}\exp\left[-\frac{1}{2}\left(\mathbf{x}-\boldsymbol{\lambda}\right)^{T}\Omega\left(\mathbf{x}-\boldsymbol{\lambda}\right)-\frac{1}{2}\left(\mathbf{x}^{\prime}-\bar{\boldsymbol{\lambda}}\right)^{T}\Omega\left(\mathbf{x}^{\prime}-\bar{\boldsymbol{\lambda}}\right)\right]\\\
\times\exp\left[\frac{1}{4}\left(\boldsymbol{\lambda}-\bar{\boldsymbol{\lambda}}\right)^{T}\Omega\left(\boldsymbol{\lambda}-\bar{\boldsymbol{\lambda}}\right)\right].$
(4.11)

In the following, we use the block notation

$\mathbf{x}=\begin{pmatrix}\mathbf{x}_{C}\\\ \mathbf{x}_{A}\\\
\end{pmatrix},\quad\boldsymbol{\lambda}=\begin{pmatrix}\boldsymbol{\lambda}_{C}\\\
\boldsymbol{\lambda}_{A}\\\
\end{pmatrix},\quad\Omega=\begin{pmatrix}\Omega_{C}&\Omega_{B}\\\
\Omega_{B}^{T}&\Omega_{A}\\\ \end{pmatrix}.$ (4.12)

The matrix columns $\mathbf{x}_{C}$ and $\boldsymbol{\lambda}_{C}$ are
$n$-dimensional and contain the coordinates and corresponding coherent state
moduli of the degrees of freedom of subsystem $A^{C}$. Similarly, the matrix
columns $\mathbf{x}_{A}$ and $\boldsymbol{\lambda}_{A}$ are
$\left(N-n\right)$-dimensional and contain the coordinates and corresponding
coherent state moduli of the degrees of freedom of subsystem $A$. Finally, in
an obvious manner, the blocks of the matrix $\Omega$ have been defined so that
they have the appropriate dimensions: $\Omega_{C}$ is an $n\times n$ matrix,
$\Omega_{B}$ is $n\times\left(N-n\right)$ and $\Omega_{A}$ is
$\left(N-n\right)\times\left(N-n\right)$.

Once we have introduced this block formalism, we can proceed to find the
reduced density matrix of subsystem $A$. This reads

$\left\langle\mathbf{x}_{A}\middle|\rho_{A}\middle|\mathbf{x}^{\prime}_{A}\right\rangle=\int
dx_{C}\left\langle\mathbf{x}_{A},\mathbf{x}_{C}\middle|\rho\middle|\mathbf{x}^{\prime}_{A},\mathbf{x}_{C}\right\rangle.$
(4.13)

It is a matter of tedious but straightforward algebra to show that this yields

$\left\langle\mathbf{x}_{A}\middle|\rho_{A}\left(\boldsymbol{\lambda}\right)\middle|\mathbf{x}^{\prime}_{A}\right\rangle=\sqrt{\frac{\det\Omega}{\pi^{N}}}\exp\left[-\frac{1}{2}\left(\mathbf{x}_{A}-\mathbf{x}_{0A}\right)^{T}\gamma\left(\mathbf{x}_{A}-\mathbf{x}_{0A}\right)\right.\\\
\left.-\frac{1}{2}\left(\mathbf{x}^{\prime}_{A}-\mathbf{x}_{0A}\right)^{T}\gamma\left(\mathbf{x}^{\prime}_{A}-\mathbf{x}_{0A}\right)+\left(\mathbf{x}^{\prime}_{A}-\mathbf{x}_{0A}\right)^{T}\beta\left(\mathbf{x}_{A}-\mathbf{x}_{0A}\right)+i\left(\mathbf{x}_{A}-\mathbf{x}^{\prime}_{A}\right)^{T}\boldsymbol{\delta}\right]\\\
\times\int
d\mathbf{x}_{C}\exp\left[-\left(\mathbf{x}_{C}-\mathbf{x}_{0C}\right)^{T}\Omega_{C}\left(\mathbf{x}_{C}-\mathbf{x}_{0C}\right)\right],$
(4.14)

where

$\displaystyle\mathbf{x}_{0C}$
$\displaystyle=\mathrm{Re}\boldsymbol{\lambda}_{C}+\Omega_{C}^{-1}\Omega_{B}\left(\mathrm{Re}\boldsymbol{\lambda}_{A}-\frac{\mathbf{x}_{A}+\mathbf{x}^{\prime}_{A}}{2}\right),$
(4.15) $\displaystyle\mathbf{x}_{0A}$
$\displaystyle=\mathrm{Re}\boldsymbol{\lambda}_{A},$ (4.16)
$\displaystyle\beta$
$\displaystyle=\frac{1}{2}\Omega_{B}^{T}\Omega_{C}^{-1}\Omega_{B},$ (4.17)
$\displaystyle\gamma$ $\displaystyle=\Omega_{A}-\beta,$ (4.18)
$\displaystyle\boldsymbol{\delta}$
$\displaystyle=\Omega_{B}^{T}\mathrm{Im}\boldsymbol{\lambda}_{C}+\Omega_{A}\mathrm{Im}\boldsymbol{\lambda}_{A}.$
(4.19)

Notice that the matrices $\gamma$ and $\beta$ are independent of the moduli
$\boldsymbol{\lambda}$. Therefore, they do not depend on the coherent state of
the overall system; they are the same for all coherent states.

The Gaussian integral can trivially be performed to obtain

$\left\langle\mathbf{x}_{A}\middle|\rho_{A}\left(\boldsymbol{\lambda}\right)\middle|\mathbf{x}^{\prime}_{A}\right\rangle=\sqrt{\frac{\det\left(\gamma-\beta\right)}{\pi^{\left(N-n\right)}}}\exp\left[-\frac{1}{2}\left(\mathbf{x}_{A}-\mathbf{x}_{0A}\right)^{T}\gamma\left(\mathbf{x}_{A}-\mathbf{x}_{0A}\right)\right.\\\
\left.-\frac{1}{2}\left(\mathbf{x}^{\prime}_{A}-\mathbf{x}_{0A}\right)^{T}\gamma\left(\mathbf{x}^{\prime}_{A}-\mathbf{x}_{0A}\right)+\left(\mathbf{x}^{\prime}_{A}-\mathbf{x}_{0A}\right)^{T}\beta\left(\mathbf{x}_{A}-\mathbf{x}_{0A}\right)+i\left(\mathbf{x}_{A}-\mathbf{x}^{\prime}_{A}\right)^{T}\boldsymbol{\delta}\right].$
(4.20)

The matrices $\gamma$ and $\beta$ are by construction real symmetric
$\left(N-n\right)\times\left(N-n\right)$ matrices. Let $\mathbf{v}_{i}$,
$i=1,2,\ldots,N-n$, be a set of orthonormal eigenvectors of the symmetric
matrix

$B:=\gamma^{-\frac{1}{2}}\beta\gamma^{-\frac{1}{2}},$ (4.21)

i.e.

$I_{N-n}=\sum_{i=1}^{N-n}\mathbf{v}_{i}\mathbf{v}_{i}^{T},\quad
B=\sum_{i=1}^{N-n}\beta^{D}_{i}\mathbf{v}_{i}\mathbf{v}_{i}^{T},$ (4.22)

where $\beta^{D}_{i}$ is the eigenvalue corresponding to the eigenvector
$\mathbf{v}_{i}$.

We define

$\tilde{x}_{Ai}:=\mathbf{v}_{i}^{T}\gamma^{\frac{1}{2}}\mathbf{x}_{A},\quad\tilde{x}^{\prime}_{Ai}:=\mathbf{v}_{i}^{T}\gamma^{\frac{1}{2}}\mathbf{x}^{\prime}_{A},\quad\tilde{x}_{0Ai}:=\mathbf{v}_{i}^{T}\gamma^{\frac{1}{2}}\mathbf{x}_{0A},$
(4.23)
$\tilde{y}_{Ai}:=\tilde{x}_{Ai}-\tilde{x}_{0Ai},\quad\tilde{y}^{\prime}_{Ai}:=\tilde{x}^{\prime}_{Ai}-\tilde{x}_{0Ai}$
(4.24)

and finally

$w_{i}:=\mathbf{v}_{i}^{T}\gamma^{-\frac{1}{2}}\boldsymbol{\delta}.$ (4.25)

Using the above definitions, the reduced density matrix assumes the form

$\left\langle\mathbf{x}_{A}\middle|\rho_{A}\left(\boldsymbol{\lambda}\right)\middle|\mathbf{x}^{\prime}_{A}\right\rangle\\\
=\sqrt{\frac{\det\left(I-\beta^{D}\right)}{\pi^{\left(N-n\right)}}}\exp\left[\sum_{i=1}^{N-n}\left(-\frac{1}{2}\tilde{y}_{Ai}^{2}-\frac{1}{2}\tilde{y}^{\prime
2}_{Ai}+\beta^{D}_{i}\tilde{y}_{Ai}\tilde{y}^{\prime}_{Ai}+i\left(\tilde{y}_{Ai}-\tilde{y}^{\prime}_{Ai}\right)w_{i}\right)\right].$
(4.26)

We observe that the density matrix has assumed a form, where it can be written
as the tensor product of $N-n$ matrices; each of these matrices is a function
of only one of the coordinates $\tilde{x}_{Ai}$. For this reason we name the
coordinates $\tilde{x}_{Ai}$ as “the normal coordinates of the reduced
system”. Having written down the reduced density matrix in this form, we may
take advantage of Mehler’s formula (3.17) and write the reduced density matrix
as

$\left\langle\mathbf{x}_{A}\middle|\rho_{A}\left(\boldsymbol{\lambda}\right)\middle|\mathbf{x}^{\prime}_{A}\right\rangle=\prod_{i=1}^{N-n}\left(\frac{\alpha_{i}}{\pi}\right)^{1/2}\sum_{n_{i}=0}^{\infty}\frac{1-\xi_{i}}{n!}\left(\frac{\xi_{i}}{2}\right)^{n}H_{n}\left(\sqrt{\alpha_{i}}\tilde{y}_{Ai}\right)H_{n}\left(\sqrt{\alpha_{i}}\tilde{y}^{\prime}_{Ai}\right)\\\
\times\exp\left[-\frac{\alpha_{i}}{2}\left(\tilde{y}_{Ai}^{2}+\tilde{y}^{\prime
2}_{Ai}\right)+iw_{i}\left(\tilde{y}_{Ai}-\tilde{y}^{\prime}_{Ai}\right)\right].$
(4.27)

where we have defined

$\alpha_{i}:=\sqrt{1-\left(\beta^{D}_{i}\right)^{2}},\quad\xi_{i}=\frac{\beta^{D}_{i}}{1+\alpha_{i}}.$
(4.28)

In other words,

$\left\langle\mathbf{x}_{A}\middle|\rho_{A}\left(\boldsymbol{\lambda}\right)\middle|\mathbf{x}^{\prime}_{A}\right\rangle=\prod_{i=1}^{N-n}\left(1-\xi_{i}\right)\sum_{n_{i}=0}^{\infty}\xi_{i}^{n_{i}}\psi_{n_{i}}\left(\tilde{y}_{Ai}\right)\psi_{n_{i}}\left(\tilde{y}^{\prime}_{Ai}\right)\exp\left[iw_{i}\left(\tilde{y}_{Ai}-\tilde{y}^{\prime}_{Ai}\right)\right],$
(4.29)

where $\psi_{n_{i}}$ are the eigenfunctions of the Hamiltonian of an effective
simple harmonic oscillator with eigenfrequency equal to
$\omega_{\text{eff}}=\alpha_{i}$. We may further write the reduced density
matrix as

$\left\langle\mathbf{x}_{A}\middle|\rho_{A}\left(\boldsymbol{\lambda}\right)\middle|\mathbf{x}^{\prime}_{A}\right\rangle=\prod_{i=1}^{N-n}\left(1-\xi_{i}\right)\sum_{n_{i}=0}^{\infty}\xi_{i}^{n_{i}}\Psi_{n_{i}}\left(\tilde{y}_{Ai}\right)\bar{\Psi}_{n_{i}}\left(\tilde{y}^{\prime}_{Ai}\right),$
(4.30)

where

$\Psi_{n_{i}}\left(x\right)=\psi_{n_{i}}\left(x\right)\exp\left[iw_{i}x-i\alpha_{i}t\left(n+\frac{1}{2}\right)\right].$
(4.31)

The orthonormality of the eigenstates of the harmonic oscillator, combined
with the fact that the phase $w_{i}$ is _the same_ for all states
$\Psi_{n_{i}}$ implies that the states $\Psi_{n_{i}}$ are also orthonormal. It
follows that equation (4.30) provides the expansion of the reduced density
matrix into its eigenfunctions, which read

$\Psi_{\left\\{n_{i}\right\\}}=\prod_{i=1}^{N-n}\Psi_{n_{i}}.$ (4.32)

It is obvious that the corresponding eigenvalues

$p_{\left\\{n_{i}\right\\}}=\prod_{i=1}^{N-n}\left(1-\xi_{i}\right)\xi_{i}^{n_{i}}$
(4.33)

are independent of $\boldsymbol{\lambda}$, implying that they are time-
independent and equal to the eigenvalues of the reduced system when the
composite system lies in the vacuum state.

### 4.2 Time Evolution of the Reduced System

The eigenfunctions (3.23) of the reduced density matrix (4.29) belong to the
family of wave-functions (2.23) that solve the generalized oscillator system
that we presented in Section 2. In particular they coincide under the
identification

$m_{i}\omega_{i}^{R}=\alpha_{i},\quad\omega_{i}^{I}=0,\quad\tilde{x}_{0Ai}=\mathrm{Re}\left(\mathbf{v}_{i}^{T}\gamma^{\frac{1}{2}}\boldsymbol{\lambda}_{A}\right),\quad\tilde{p}_{0Ai}=w_{i}=\mathbf{v}_{i}^{T}\gamma^{-\frac{1}{2}}\boldsymbol{\delta},$
(4.34)

where $\hbar$ is set to one. It follows that the corresponding Hamiltonian is

$\hat{H}_{A}=\sum_{i=1}^{N-n}\frac{1}{2m_{i}}\left(\hat{\tilde{p}}_{Ai}-w_{i}+m_{i}\dot{\tilde{x}}_{0Ai}\right)^{2}+\frac{1}{2m_{i}}\alpha_{i}^{2}\left(\hat{\tilde{x}}_{Ai}-\tilde{x}_{0Ai}\right)^{2}\\\
-\dot{w}_{i}\left(\hat{\tilde{x}}_{Ai}-\tilde{x}_{0Ai}\right)+\dot{\tilde{x}}_{0Ai}\left(w_{i}-\frac{1}{2}m_{i}\dot{\tilde{x}}_{0Ai}\right).$
(4.35)

We have the freedom of determining freely the value of each $m_{i}$
independently.

As a direct consequence of (4.8), we get

$\mathrm{Re}\dot{\tilde{\boldsymbol{\lambda}}}=\Omega^{D}\mathrm{Im}\tilde{\boldsymbol{\lambda}}.$

This implies that

$\mathrm{Re}O^{T}\dot{\tilde{\boldsymbol{\lambda}}}=O^{T}\Omega^{D}O\mathrm{Im}O^{T}\tilde{\boldsymbol{\lambda}},$

i.e.

$\mathrm{Re}\dot{\boldsymbol{\lambda}}=\Omega\mathrm{Im}\boldsymbol{\lambda}.$

Employing the block form notation (4.12) yields

$\mathrm{Re}\dot{\boldsymbol{\lambda}}_{A}=\Omega_{B}^{T}\mathrm{Im}\boldsymbol{\lambda}_{C}+\Omega_{C}\mathrm{Im}\boldsymbol{\lambda}_{A}=\boldsymbol{\delta}.$

Finally, we need to revert back to the original coordinates $\mathbf{x}_{0A}$.
It is a matter of tedious algebra to show that

$\hat{H}_{A}=\frac{1}{2}\left(\mathbf{p}_{A}-\mathrm{Re}\dot{\boldsymbol{\lambda}}_{A}\right)^{T}\gamma^{-\frac{1}{2}}\left[\sum_{i=1}^{N-n}\mathbf{v}_{i}\frac{1}{m_{i}}\mathbf{v}_{i}^{T}\right]\gamma^{-\frac{1}{2}}\left(\mathbf{p}_{A}-\mathrm{Re}\dot{\boldsymbol{\lambda}}_{A}\right)\\\
+\frac{1}{2}\left(\mathbf{x}_{A}-\mathrm{Re}\boldsymbol{\lambda}_{A}\right)^{T}\gamma^{\frac{1}{2}}\left[\sum_{i=1}^{N-n}\mathbf{v}_{i}\frac{\alpha_{i}^{2}}{m_{i}}\mathbf{v}_{i}^{T}\right]\gamma^{\frac{1}{2}}\left(\mathbf{x}_{A}-\mathrm{Re}\boldsymbol{\lambda}_{A}\right)\\\
+\mathbf{p}_{A}^{T}\mathrm{Re}\dot{\boldsymbol{\lambda}}_{A}-\left(\mathbf{x}_{A}-\mathrm{Re}\boldsymbol{\lambda}_{A}\right)^{T}\mathrm{Re}\ddot{\boldsymbol{\lambda}}_{A}.$
(4.36)

Let us make two special choices for the parameters $m_{i}$. First if we choose
all $m_{i}$ so that

$m_{i}=1,$ (4.37)

then the effective Hamiltonian will assume the form

$\hat{H}_{A}=\frac{1}{2}\left(\mathbf{p}_{A}-\mathrm{Re}\dot{\boldsymbol{\lambda}}_{A}\right)^{T}\gamma^{-1}\left(\mathbf{p}_{A}-\mathrm{Re}\dot{\boldsymbol{\lambda}}_{A}\right)+\mathbf{p}_{A}^{T}\mathrm{Re}\dot{\boldsymbol{\lambda}}_{A}\\\
+\frac{1}{2}\left(\mathbf{x}_{A}-\mathrm{Re}\boldsymbol{\lambda}_{A}\right)^{T}\left(\gamma-\beta\gamma^{-1}\beta\right)\left(\mathbf{x}_{A}-\mathrm{Re}\boldsymbol{\lambda}_{A}\right)-\left(\mathbf{x}_{A}-\mathrm{Re}\boldsymbol{\lambda}_{A}\right)^{T}\mathrm{Re}\ddot{\boldsymbol{\lambda}}_{A}.$
(4.38)

If we choose (notice that this is a time-dependent mass, and, thus, the
equation (2.10) applies)

$m_{i}=\frac{\mathbf{v}_{i}^{T}\gamma^{-\frac{1}{2}}\mathrm{Re}\dot{\boldsymbol{\lambda}}_{A}}{\mathbf{v}_{i}^{T}\gamma^{\frac{1}{2}}\mathrm{Re}\dot{\boldsymbol{\lambda}}_{A}},$
(4.39)

then the linear terms in the momentum will vanish, and the effective
Hamiltonian will assume the form

$\hat{H}_{A}=\frac{1}{2}\mathbf{p}_{A}^{T}\gamma^{-\frac{1}{2}}\left[\sum_{i=1}^{N-n}\mathbf{v}_{i}\frac{\mathbf{v}_{i}^{T}\gamma^{\frac{1}{2}}\mathrm{Re}\dot{\boldsymbol{\lambda}}_{A}}{\mathbf{v}_{i}^{T}\gamma^{-\frac{1}{2}}\mathrm{Re}\dot{\boldsymbol{\lambda}}_{A}}\mathbf{v}_{i}^{T}\right]\gamma^{-\frac{1}{2}}\mathbf{p}_{A}+\frac{1}{2}\mathrm{Re}\dot{\boldsymbol{\lambda}}_{A}^{T}\mathrm{Re}\dot{\boldsymbol{\lambda}}_{A}\\\
+\frac{1}{2}\left(\mathbf{x}_{A}-\mathrm{Re}\boldsymbol{\lambda}_{A}\right)^{T}\gamma^{\frac{1}{2}}\left[\sum_{i=1}^{N-n}\mathbf{v}_{i}\frac{\alpha_{i}^{2}\mathbf{v}_{i}^{T}\gamma^{\frac{1}{2}}\mathrm{Re}\dot{\boldsymbol{\lambda}}_{A}}{\mathbf{v}_{i}^{T}\gamma^{-\frac{1}{2}}\mathrm{Re}\dot{\boldsymbol{\lambda}}_{A}}\mathbf{v}_{i}^{T}\right]\gamma^{\frac{1}{2}}\left(\mathbf{x}_{A}-\mathrm{Re}\boldsymbol{\lambda}_{A}\right)\\\
-\left(\mathbf{x}_{A}-\mathrm{Re}\boldsymbol{\lambda}_{A}\right)^{T}\mathrm{Re}\ddot{\boldsymbol{\lambda}}_{A}.$
(4.40)

Defining $\gamma^{D}$ as

$\gamma^{D}:=\sum_{i=1}^{N-n}\mathbf{v}_{i}\frac{\mathbf{v}_{i}^{T}\gamma^{\frac{1}{2}}\mathrm{Re}\dot{\boldsymbol{\lambda}}_{A}}{\mathbf{v}_{i}^{T}\gamma^{-\frac{1}{2}}\mathrm{Re}\dot{\boldsymbol{\lambda}}_{A}}\mathbf{v}_{i}^{T},$
(4.41)

Then, the latter assumes the form

$\hat{H}_{A}=\frac{1}{2}\mathbf{p}_{A}^{T}\gamma^{-\frac{1}{2}}\gamma^{D}\gamma^{-\frac{1}{2}}\mathbf{p}_{A}+\frac{1}{2}\mathrm{Re}\dot{\boldsymbol{\lambda}}_{A}^{T}\mathrm{Re}\dot{\boldsymbol{\lambda}}_{A}\\\
+\frac{1}{2}\left(\mathbf{x}_{A}-\mathrm{Re}\boldsymbol{\lambda}_{A}\right)^{T}\left(\gamma^{\frac{1}{2}}\gamma^{D}\gamma^{\frac{1}{2}}-\beta\gamma^{-\frac{1}{2}}\gamma^{D}\gamma^{-\frac{1}{2}}\beta\right)\left(\mathbf{x}_{A}-\mathrm{Re}\boldsymbol{\lambda}_{A}\right)-\left(\mathbf{x}_{A}-\mathrm{Re}\boldsymbol{\lambda}_{A}\right)^{T}\mathrm{Re}\ddot{\boldsymbol{\lambda}}_{A}.$
(4.42)

Having traced out the subsystem $A^{C}$, the subsystem $A$ constitutes an open
quantum system. The time evolution of such a system is non-unitary in general.
Nevertheless, since _all_ the eigenstates of the reduced density matrix obey
the equation

$i\hbar\dot{\Psi}_{\left\\{n_{i}\right\\}}=\hat{H}_{A}\Psi_{\left\\{n_{i}\right\\}},$
(4.43)

where $\hat{H}_{A}$ is given by (4.36), equation (3.22) directly implies that

$i\hbar\dot{\rho}_{A}=\left[\hat{H}_{A},\rho_{A}\right],$ (4.44)

meaning that when the overall system lies in a coherent state, the time
evolution of the reduced system _is unitary_ and it is driven by the
Hamiltonian (4.36).

### 4.3 The Modular Hamiltonian

In a similar manner to the specification of the effective Hamiltonian, which
evolves the reduced density matrix, we may express the modular Hamiltonian as

$\hat{H}_{M}=-\sum_{i=1}^{N-n}\ln\xi_{i}\left(\frac{1}{2m_{i}\alpha_{i}}\left(\hat{\tilde{p}}_{Ai}-w_{i}\right)^{2}+\frac{\alpha_{i}}{2m_{i}}\left(\hat{\tilde{x}}_{Ai}-\tilde{x}_{0Ai}\right)^{2}\right)+\frac{1}{2}\ln\xi_{i}-\ln\left(1-\xi_{i}\right).$
(4.45)

We need to revert back to the original coordinates $\mathbf{x}_{0A}$. It is a
matter of tedious algebra to show that

$\hat{H}_{A}=\frac{1}{2}\left(\mathbf{p}_{A}-\mathrm{Re}\dot{\boldsymbol{\lambda}}_{A}\right)^{T}\gamma^{-\frac{1}{2}}\left[\sum_{i=1}^{N-n}\mathbf{v}_{i}\frac{\ln\xi_{i}}{m_{i}\alpha_{i}}\mathbf{v}_{i}^{T}\right]\gamma^{-\frac{1}{2}}\left(\mathbf{p}_{A}-\mathrm{Re}\dot{\boldsymbol{\lambda}}_{A}\right)\\\
+\frac{1}{2}\left(\mathbf{x}_{A}-\mathrm{Re}\boldsymbol{\lambda}_{A}\right)^{T}\gamma^{\frac{1}{2}}\left[\sum_{i=1}^{N-n}\mathbf{v}_{i}\frac{\ln\xi_{i}\alpha_{i}}{m_{i}}\mathbf{v}_{i}^{T}\right]\gamma^{\frac{1}{2}}\left(\mathbf{x}_{A}-\mathrm{Re}\boldsymbol{\lambda}_{A}\right)\\\
-\sum_{i=1}^{N-n}\frac{1}{2}\ln\xi_{i}-\ln\left(1-\xi_{i}\right).$ (4.46)

If we choose all parameters $m_{i}$ to be equal to 1,

$m_{i}=1,$ (4.47)

then the modular Hamiltonian will assume the form

$\hat{H}_{A}=-\frac{1}{2}\left(\mathbf{p}_{A}-\mathrm{Re}\dot{\boldsymbol{\lambda}}_{A}\right)^{T}K_{\mathrm{mod}}\left(\mathbf{p}_{A}-\mathrm{Re}\dot{\boldsymbol{\lambda}}_{A}\right)\\\
-\frac{1}{2}\left(\mathbf{x}_{A}-\mathrm{Re}\boldsymbol{\lambda}_{A}\right)^{T}V_{\mathrm{mod}}\left(\mathbf{x}_{A}-\mathrm{Re}\boldsymbol{\lambda}_{A}\right)\\\
-\sum_{i=1}^{N-n}\frac{1}{2}\ln\xi_{i}-\ln\left(1-\xi_{i}\right),$ (4.48)

where

$\displaystyle K_{\mathrm{mod}}$ $\displaystyle=\gamma^{-\frac{1}{2}}\left[\ln
B-\ln\left(I+\left(I-B^{2}\right)^{\frac{1}{2}}\right)\right]\left(I-B^{2}\right)^{-\frac{1}{2}}\gamma^{-\frac{1}{2}},$
(4.49) $\displaystyle V_{\mathrm{mod}}$
$\displaystyle=\gamma^{\frac{1}{2}}\left[\ln
B-\ln\left(I+\left(I-B^{2}\right)^{\frac{1}{2}}\right)\right]\left(I-B^{2}\right)^{\frac{1}{2}}\gamma^{\frac{1}{2}}.$
(4.50)

## 5 Discussion

We have managed to calculate the reduced density matrix and the modular
Hamiltonian in an arbitrary harmonic system lying at an arbitrary coherent
state. Furthermore, we showed that for these specific states, the time
evolution of the reduced density matrix is unitary and we specified the
effective Hamiltonian, which generates this unitary evolution. For a general
coherent state, the reduced density matrix has non-trivial time evolution.
This is natural, since the state of the composite system also evolves with
time. In the special case that the composite system lies at its ground state,
the state of the composite system and the reduced density matrix are time-
independent. However, even in this case the system does not lie in a thermal
state. Each normal mode of the reduced system actually lies in a thermal
state, but the corresponding temperatures are different. These are determined
by the eigenvalues $\beta^{D}_{i}$ of the matrix $B$, namely,

$T_{i}=\frac{\alpha_{i}}{\ln\xi_{i}},$ (5.1)

where $\alpha_{i}$ and $\xi_{i}$ are given as functions of $\beta^{D}_{i}$ by
equation (4.28). This relation between thermodynamics and entanglement in
harmonic systems at their ground state has been exploited in [13]. This
temperature spectrum is not contradictory to the thermalization hypothesis
[14, 15]. The latter applies only to non-integrable systems.

The form of the temperatures above implies that the modular Hamiltonian,
expressed in terms of the normal modes of the reduced system (4.45) assumes
the form

$\hat{H}_{M}=-\sum_{i=1}^{N-n}\frac{1}{T_{i}}\left(\frac{1}{2m_{i}}\left(\hat{\tilde{p}}_{Ai}-w_{i}\right)^{2}+\frac{\alpha_{i}^{2}}{2m_{i}}\left(\hat{\tilde{x}}_{Ai}-\tilde{x}_{0Ai}\right)^{2}\right)+\frac{1}{2}\ln\xi_{i}-\ln\left(1-\xi_{i}\right).$
(5.2)

In other words, the modular Hamiltonian contains the Hamiltonian of each
normal mode of the reduced system, which appears in the effective Hamiltonian
(4.35), but weighted with its inverse temperature.

The above expression of the modular Hamiltonian is naively a non-local
function of the stress-energy tensor of the underlying harmonic system, since
it is expressed in terms of the normal modes of the reduced system. It is
well-known that in general the modular Hamiltonian is indeed a non-local
function of the stress-energy tensor, except for cases with very specific
characteristics. Spherical entangling surfaces in conformal field theory at
the ground state is such an exceptional case [17]; this derivation is based on
the long-known local expression of the modular Hamiltonian of the Rindler
wedge [18, 19], which holds in any quantum field theory. If any of these
conditions is violated, e.g. by introducing a mass term [20], non-local terms
emerge. The expression (5.2) can be useful in the investigation of whether and
under which conditions the modular Hamiltonian can be expressed as a local
function of the stress-energy tensor. Such investigations have advanced more
in strongly constrained systems, such as two-dimensional conformal field
theory [21].

Continuing the above discussion, the expression (5.2) provides an insight of
the interrelation between temperatures and the modular Hamiltonian at the
level of the normal modes of the reduced system, i.e. in a non-local way. It
is interesting to investigate whether this relation can be used to further
probe the relation between local temperature and the modular Hamiltonian, or
at least the local terms of the latter [22, 23].

We would like to make a comment about the special case of the ground state of
the composite system. In this case, the effective Hamiltonian (4.35) is simply
the sum of the quadratic simple harmonic oscillator Hamiltonians of each
normal mode of the reduced system. The modular Hamiltonian is the same sum,
but weighted with the appropriate temperatures. What the effective Hamiltonian
should obey? it should provide the time evolution of the reduced density
matrix, which in this case is trivial, since it is time-independent. Since the
normal coordinates of the reduced system are orthogonal, any linear
combination of the quadratic simple harmonic oscillator Hamiltonians of each
normal mode of the reduced system commutes with the modular Hamiltonian and,
thus, with the reduced density matrix. So how one chooses the effective
Hamiltonian in this case? The limit of (4.35) as the state of the composite
system tends to the ground state is the specific selection which is
continuously connected to the corresponding effective Hamiltonians that
generate the dynamics of the reduced system for all classicalmost states
including the ground state. So this is a natural choice for the ground state
effective Hamiltonian. However, even when making this special selection, the
freedom of arbitrary choice of the parameters $m_{i}$, results once again in
this Hamiltonian being an arbitrary linear combination of the aforementioned
normal modes. Notice that the freedom of these linear combination is not that
innocent. The multiplication of a harmonic Hamiltonian with a constant greatly
alters the dynamics it describes as it alters the eigenfrequency of the
corresponding oscillator.

There is also another property of the spectrum of the reduced density matrix,
which complicates the calculation of the modular Hamiltonian, via the formula
(4.45). In the case that the subsystem $A^{C}$ has fewer degrees of freedom
than subsystem $A$, the spectrum of the matrix $B$ contains at least as many
vanishing eigenvalues as the difference of the numbers of the degrees of
freedom of the two subsystems. It is simple to prove this statement.
Performing a Schmidt decomposition of the state of the composite system state
yields

$\left|\Psi\right>=\sum_{i=1}^{m}\sum_{j=1}^{m}c_{i}\left|\psi_{Ai}\right>\otimes\left|\psi_{A^{C}i}\right>,$
(5.3)

where the states $\left|\psi_{Ai}\right>$ and $\left|\psi_{A^{C}i}\right>$ are
sets of orthonormal states of the subsystems $A$ and $A^{C}$, which contain
equal number of states at most equal to the dimensionality of the smaller
subsystem, i.e.

$m\leq\min\left(\dim A,\dim A^{C}\right).$ (5.4)

This directly implies that

$\rho_{A}=\sum_{i=1}^{m}c_{i}c^{*}_{i}\left|\psi_{Ai}\right>\left<\psi_{Ai}\right|,\quad\rho_{A^{C}}=\sum_{i=1}^{m}c_{i}c^{*}_{i}\left|\psi_{A^{C}i}\right>\left<\psi_{A^{C}i}\right|.$
(5.5)

It directly follows that the two reduced density matrices have exactly the
same spectrum up to vanishing eigenvalues. These vanishing eigenvalues
correspond to the vectors of the Hilbert space of the larger subsystem, which
are orthogonal to the vectors $\left|\psi_{Ai}\right>$ or
$\left|\psi_{A^{C}i}\right>$, depending on which is the larger subsystem.

This means that the larger subsystem necessarily contains modes which are
disentangled to the smaller subsystem, and, thus, remain frozen upon tracing
out the smaller subsystem. Equation (5.2) directly implies that the modular
Hamiltonian of the larger subsystem diverges. This apparent problem can be
resolved via the introduction of an IR cutoff such as temperature. If the
composite system lies at a thermal state instead of its ground state, these
vanishing temperatures will tend to become equal to the temperature of the
composite system yielding the modular Hamiltonian finite.

As a final comment, in the special case that all temperatures $T_{i}$ become
equal to each other, the quadratic parts of the effective Hamiltonian and the
modular Hamiltonian become proportional to each other. In the special case of
the ground state there are no linear terms and the two Hamiltonians become
literally proportional to each other.

These results have a profound application in the spirit of [1]. In the seminal
work of Srednicki, free scalar field theory is treated as a harmonic system
upon discretization of the degrees of freedom. The determination of the
corresponding modular Hamiltonian in this case may provide more hints for the
similarity between entanglement and gravity beyond the similarity of
entanglement and black hole entropies. More specifically, we propose that the
gravitational dynamics may be encoded in the structure of the modular
Hamiltonian. The modular Hamiltonian obeys the first law of entanglement
thermodynamics with respect to entanglement entropy. The latter is dominated
by an area law as shown in the same work [1]. As a direct consequence of [16]
the modular Hamiltonian should somehow embody Einstein gravity dynamics.
Having written down the modular Hamiltonian in the form (4.48) allows to check
whether the above argument indeed implies that the modular Hamiltonian
resembles the discretized Hamiltonian for scalar degrees of freedom
propagating on a curved background, greatly enhancing the understanding of the
relation between gravity and entanglement.

## Acknowledgements

D. K. was supported by FAPESP Grant No. 2021/01819-0. The authors would like
to thank K. Papadodimas for useful discussions.

##  Appendix A Coherent and Squeezed States of the Simple Harmonic Oscillator

A simple application of the formulation of the generalized, time-dependent
oscillator that we developed in Section 2 is the construction of the Gaussian
states of the simple harmonic oscillator, namely the coherent and squeezed
states [12]. Although this material is well-known and easily accessible in the
literature, we present it for completeness and in order to demonstrate how
easily it emerges in this formalism.

In Section 2 we showed that the Gaussian state (2.16) solves the time-
dependent Schrödinger equation with Hamiltonian given by (2.12). This
Hamiltonian coincides with the Hamiltonian of the simple harmonic oscillator
with eigenfrequency equal to $\omega$, provided that

$\displaystyle\omega^{I}$
$\displaystyle=\frac{\dot{\omega}^{R}}{2\omega^{R}},$ (A.1)
$\displaystyle\omega^{2}$
$\displaystyle=\frac{d}{dt}\left(\frac{\dot{\omega}^{R}}{2\omega^{R}}\right)-\left(\frac{\dot{\omega}^{R}}{2\omega^{R}}\right)^{2}+\left(\omega^{R}\right)^{2},$
(A.2) $\displaystyle p_{0}$
$\displaystyle=m\dot{x}_{0},\quad\ddot{x}_{0}+\omega^{2}x_{0}=0.$ (A.3)

The equations (A.3) are solved by

$x_{0}=\alpha\cos\left(\omega t+\phi_{0}\right),\qquad
p_{0}=-m\omega\alpha\sin\left(\omega t+\phi_{0}\right).$ (A.4)

In order to solve equation (A.2), we set

$\omega^{R}=\frac{\omega}{f^{2}}.$ (A.5)

Then, $f$ obeys the equation

$\frac{1}{\omega^{2}}\ddot{f}+f-\frac{1}{f^{3}}=0.$ (A.6)

We can integrate this equation once to obtain

$\frac{1}{4\omega^{2}}\left[\frac{\partial}{\partial
t}\left(f^{2}\right)\right]^{2}=-f^{4}+2\cosh\gamma f^{2}-1,$ (A.7)

where we defined the integration constant as $2\cosh\gamma$. The integration
constant has been selected so that the right-hand-side of the above equation
can be positive, allowing the existence of a real solution. For this purpose,
the integration constant has to be larger than 2, and thus we defined it as
$2\cosh\gamma$. Solving this equation yields

$f^{2}=\cosh\gamma+\sin\left(2\omega\left(t-t_{0}\right)\right)\sinh\gamma,$
(A.8)

implying that $\omega^{R}$ is given by

$\omega^{R}=\frac{\omega}{\cosh\gamma+\sin\left(2\omega\left(t-t_{0}\right)\right)\sinh\gamma}.$
(A.9)

Equation (A.1) directly implies that

$\omega^{I}=-\frac{\omega\cos\left(2\omega\left(t-t_{0}\right)\right)\tanh\gamma}{1+\sin\left(2\omega\left(t-t_{0}\right)\right)\tanh\gamma}.$
(A.10)

Finally, the phase $\phi$, which appears in the wavefunctions (2.23), is
simply the time integral of $\omega^{R}$, as given by (2.17). It turns out
that it is equal to

$\phi=\arctan\left[\frac{\tan\left(\omega\left(t-t_{0}\right)\right)+\tanh\frac{\gamma}{2}}{1+\tan\left(\omega\left(t-t_{0}\right)\right)\tanh\frac{\gamma}{2}}\right].$
(A.11)

The potential of the generalized harmonic oscillator is not exactly the
potential of the simple harmonic oscillator; it contains a c-number term too.
This term is eliminated by

$\Psi_{n}\left(x,t\right)\rightarrow e^{-i\frac{m\omega
a^{2}}{4h}\sin\left(2\omega t+\phi\right)}\Psi_{n}\left(x,t\right).$ (A.12)

Thus, the time-dependent Schrödinger equation

$i\hbar\frac{\partial}{\partial
t}\Psi=-\frac{\hbar^{2}}{2m}\frac{\partial^{2}}{{\partial
x}^{2}}\Psi+\frac{1}{2}m\omega^{2}x^{2}\Psi,$ (A.13)

has the following family of solutions

$\Psi_{n}=\frac{1}{\sqrt{2^{n}n!}}\left(\frac{m\omega^{R}}{\pi\hbar}\right)^{\frac{1}{4}}H_{n}\left(\sqrt{\frac{m\omega^{R}}{\hbar}}\left(x-x_{0}\right)\right)\\\
\times\exp\left[-\frac{m\left(\omega^{R}+i\omega^{I}\right)}{2\hbar}\left(x-x_{0}\right)^{2}+i\frac{p_{0}}{\hbar}\left(x-x_{0}\right)-i\phi_{n}\right],$
(A.14)

where

$\phi_{n}=\left(n+\frac{1}{2}\right)\arctan\left[\frac{\tan\left(\omega\left(t-t_{0}\right)\right)+\tanh\frac{\gamma}{2}}{1+\tan\left(\omega\left(t-t_{0}\right)\right)\tanh\frac{\gamma}{2}}\right]+\frac{m\omega}{4h}\alpha^{2}\sin\left(2\left(\omega
t+\phi\right)\right).$ (A.15)

Notice, that the $n=0$ wavefunction is Gaussian and actually it describes a
squeezed state. In the special case $\gamma=0$, this reduces to a coherent
state.

Using the usual creation and annihilation operator techniques, it is
straightforward to calculate the following expectation values

$\displaystyle\langle x\rangle$ $\displaystyle=x_{0},$ (A.16)
$\displaystyle\langle x^{2}\rangle$
$\displaystyle=\left(n+\frac{1}{2}\right)\frac{\hbar}{m\omega^{R}}+x_{0}^{2},$
(A.17) $\displaystyle\langle p\rangle$ $\displaystyle=p_{0}$ (A.18)
$\displaystyle\langle p^{2}\rangle$
$\displaystyle=\left(n+\frac{1}{2}\right)\hbar
m\omega^{R}\left(1+\left(\frac{\omega^{I}}{\omega^{R}}\right)^{2}\right)+p_{0}^{2},$
(A.19)

when the system is described by the wave-functions (2.23). Thus, we obtain

$\Delta x\Delta
p=\hbar\left(n+\frac{1}{2}\right)\sqrt{1+\left(\frac{\omega^{I}}{\omega^{R}}\right)^{2}}.$
(A.20)

In the case of squeezed states ($n=0$), it follows that

$\Delta x\Delta
p=\frac{\hbar}{2}\sqrt{1+\cos^{2}\left(2\omega\left(t-t_{0}\right)\right)\sinh^{2}\gamma},$
(A.21)

which is the well-known uncertainty relation of the squeezed states.

Summing up, the algebraic construction of the generalized time-dependent
oscillator, which is presented in Section 2, trivializes the construction of
the Gaussian states of the simple harmonic oscillator.

## References

  * [1] M. Srednicki, “Entropy and area,” Phys. Rev. Lett. 71, 666 (1993) [hep-th/9303048].
  * [2] M. Van Raamsdonk, “Comments on quantum gravity and entanglement”, [arXiv:0907.2939 [hep-th]].
  * [3] M. Van Raamsdonk, “Building up spacetime with quantum entanglement”, Gen. Rel. Grav. 42, 2323-2329 (2010) [arXiv:1005.3035 [hep-th]].
  * [4] N. Lashkari, M. B. McDermott and M. Van Raamsdonk, “Gravitational dynamics from entanglement ’thermodynamics’ ”, JHEP 04, 195 (2014) [arXiv:1308.3716 [hep-th]].
  * [5] K. Husimi, “Miscellanea in Elementary Quantum Mechanics, II” Progr. Theoret. Phys., 9 (1953), p. 381
  * [6] X. Chapuisat, J. Manz, “A closed-form solution of the Schroedinger equation describing the harmonic oscillator with time dependent frequency and acted on by a time dependent perturbative force” Mol. Phys., 26 (1973), p. 1577
  * [7] H.-D. Meyer, “On the forced harmonic oscillator with time-dependent frequency”, Chemical Physics, Volume 61, Issue 3 (1981), 365-383,
  * [8] C. J. Ballhausen, “Time-dependent raising and lowering operators for the ageing, forced and damped oscillator in quantum mechanics”, Chem. Phys. Lett. 192,491992
  * [9] K. B. Moller and T. G. Jorgensen, “Wave packet dynamics and photofragmentation in time-dependent quadratic potentials”, J. Chem. Phys. 105, 5037 (1996).
  * [10] E. Benedict and S. Y. Pi, “Entanglement entropy of nontrivial states,” Annals Phys. 245, 209 (1996) [hep-th/9505121].
  * [11] C. G. Callan, Jr. and F. Wilczek, “On geometric entropy,” Phys. Lett. B 333, 55 (1994) [hep-th/9401072].
  * [12] D. Stoler, “Equivalence Classes of Minimum Uncertainty Packets”, Physical Review D, 1 (12), 3217–3219 (1970)
  * [13] V. Eisler, G. Di Giulio, E. Tonni and I. Peschel, “Entanglement Hamiltonians for non-critical quantum chains”, J. Stat. Mech. 2010 (2020), 103102 [arXiv:2007.01804 [cond-mat.stat-mech]].
  * [14] M. Srednicki, “Chaos and Quantum Thermalization”, Physical Review E. 50 (2), 888–901 (1994) [arXiv:cond-mat/9403051 [cond-mat]].
  * [15] J.M. Deutsch, “Quantum statistical mechanics in a closed system”, Physical Review A. 43 (4): 2046–2049 (1991)
  * [16] T. Jacobson, “Thermodynamics of space-time: The Einstein equation of state”, Phys. Rev. Lett. 75, 1260-1263 (1995) [arXiv:gr-qc/9504004 [gr-qc]].
  * [17] H. Casini, M. Huerta and R. C. Myers, “Towards a derivation of holographic entanglement entropy”, JHEP 05, 036 (2011) [arXiv:1102.0440 [hep-th]].
  * [18] J. J. Bisognano and E. H. Wichmann, “On the Duality Condition for Quantum Fields”, J. Math. Phys. 17, 303-321 (1976)
  * [19] J. J. Bisognano and E. H. Wichmann, “On the Duality Condition for a Hermitian Scalar Field” J. Math. Phys. 16, 985-1007 (1975)
  * [20] R. Longo and G. Morsella, “The massive modular Hamiltonian,” [arXiv:2012.00565 [math-ph]].
  * [21] J. Cardy and E. Tonni, “Entanglement hamiltonians in two-dimensional conformal field theory”, J. Stat. Mech. 1612 (2016) no.12, 123103 [arXiv:1608.01283 [cond-mat.stat-mech]].
  * [22] R. Arias, D. Blanco, H. Casini and M. Huerta, “Local temperatures and local terms in modular Hamiltonians”, Phys. Rev. D 95 (2017) no.6, 065005 doi:10.1103/PhysRevD.95.065005 [arXiv:1611.08517 [hep-th]].
  * [23] R. Arias, H. Casini, M. Huerta and D. Pontello, “Anisotropic Unruh temperatures”, Phys. Rev. D 96 (2017) no.10, 105019 doi:10.1103/PhysRevD.96.105019 [arXiv:1707.05375 [hep-th]].