# Evolution of capacity of entanglement and modular entropy in harmonic chains
and scalar fields

K. Andrzejewski111e-mail<EMAIL_ADDRESS>
Faculty of Physics and Applied Informatics
University of Lodz
Pomorska 149/153, 90-236, Lodz, Poland

###### Abstract

We examine the temporal evolution of the modular entropy and capacity (in
particular, the fluctuation of the entanglement entropy) for systems of time-
dependent oscillators coupled by a (time-dependent) parameter. Such models,
through the discretization procedure, fit into field theory problems arising
from quench phenomena or non-static spacetimes.

First, we compare the dynamics of the modular and Rényi entropies and derive
the form of the modular capacity for the single time-dependent oscillator as
well as chains with bipartite decompositions. In the latter case we analyse
distinguished periodicities during the evolution and the role of various
boundary conditions. Next, we focus on the dynamics of the capacity
(fluctuation) of entanglement. We compare the results obtained with the
predictions of quasiparticles models; in particular, we obtain a theoretical
value of the initial slope of the capacity for abrupt quenches. We study also
continuous protocols with the frequency that vanishes at plus (and minus)
infinity, including a model in which the frequency tends to the Dirac delta.
All the above issues are discussed with the emphasis on the analytical
methods.

##  1 Introduction

The notion of the entropy appears in many branches of physics and thus various
approaches and generalizations to it have been proposed. The von Neumann (vN)
entropy plays the prominent role due to its applications in quantum physics,
especially in the study of the entanglement which, in turn, is fundamental for
quantum information processing and technology. Obviously, the vN entropy
provides only an overall characterization of the state (in particular, the
entanglement of the state). In order to get more information some other
measures, such as the Rényi or Tsallis entropies, were proposed. The
motivation for various measures of the entanglement arises also from attempts
to understand this notion in quantum field theory; this is challenging due to
divergences (and thus the necessity of a regularization procedure, e.g. UV
cutoff). However, the progress has been made in the recent years, see e.g. [1,
2, 3] for review. One of the most intensively studied issues was the relation
between the entropy and area [4, 5] emerging from black holes physics (see [6,
7] and references therein) and its generalization in the context gauge/gravity
duality (holographic entanglement entropy [8]).

An interesting view on the notion of entropy can be obtained by considering a
(semipositive) operator defined as the minus logarithm of the density
operator. Then, the vN entropy becomes the expectation value of such an
operator, which is therefore called the entanglement (or modular) Hamiltonian,
see, among others, [9]-[13]. By having the modular Hamiltonian we can apply
thermodynamical approach and construct thermodynamical quantities, in
particular the entropy and capacity. Such an approach leads to new notions,
the so-called modular entropy and modular capacity [12, 14]. Obviously, these
notions will be useful if they allow us to understand better the physical
phenomena. For the modular entropy such a situation appears in the study of
the relation between the geometry of spacetime and strongly coupled field
theories via the AdS/CFT correspondence (in the Ryu-Takayanagi approach [8]).
It turns out that in some cases the modular entropy seems a natural
generalization of the vN entropy [14, 15, 16].

At the same time, the notion of the modular capacity gained also some
attention. Originally it appears in the condensed mater physics [12] and
recently in other contexts (in particular, quantum gravitational effects)
[17]-[24] (for more references see [24]). It turns out that, for some
specified value of parameter the modular capacity reduces to the variance of
entropy, thereby it measures the fluctuation of the entanglement entropy [25,
26]. In this way the capacity of the entanglement gives another qualitative
measure of the entanglement spectrum [27]-[30].

In this work we analyse and compare modular entropy and capacity with the
special emphasis put on their dynamics. We will focus on the harmonic chains
with time-dependent parameters as examples of continuous-variable quantum
systems. Such systems are interesting due to discretization procedure of
various time-dependent field models resulting from quench phenomena and non-
static spacetimes. We start by reviewing briefly the recent results on this
topic and outlining the plan of the paper. Namely, the dynamics of the vN
entanglement entropy of time-dependent harmonic chains was analysed in various
approaches, see [31]-[40]; despite this effort some misunderstanding appeared.
We clarify these issues and then extend them to the modular entropy and/or
continuous quenches, see Secs. 3 and 4.

On the other hand, the evolution of the capacity (fluctuation) of the
entanglement has been recently analysed for infinite harmonic chains (and the
abrupt protocol) in Ref. [24]. Such a model corresponds to 1+1 dimensional
(conformal) field theory on the infinite line. In particular, the problem of
different slopes for the linear growth of capacity and entropy has been
indicated. In Sec. 5 we consider these issues for the finite size systems by
considering various boundary conditions as well as the quasiparticles
approach. In Sec. 6 we study the above points for continuous quenches by means
of some analytical examples.

##  2 Preliminaries

The von Neumann (vN) entropy of the density matrix $\rho$

$S_{vN}=S=-\textrm{tr}(\rho\ln(\rho)),$ (2.1)

is one of the basic information-theoretic measures. However, the vN entropy
provides only a first characterization of $\rho$ (e.g. entanglement), in order
to get more information some other measures were proposed. One of them is the
Rényi entropy defined as

$R_{\alpha}=\frac{1}{(1-\alpha)}\ln(\textrm{tr}(\rho^{\alpha})),$ (2.2)

for which $\lim_{\alpha\rightarrow 1}R_{\alpha}=S$ holds. Another
generalization, motivated by the fact that trace of the power of density
matrix can be used to characterize its eigenvalues, is given by the Tsallis
entropy $T_{\alpha}=(\textrm{tr}(\rho^{\alpha})-1)/(1-\alpha)$.

In order to get more insight into the notion of entropy let us consider the
logarithm of the density matrix

$K=-\ln(\rho).$ (2.3)

Then, the vN entropy is the expectation value of $K$, $S=\textrm{tr}(\rho
K)=\langle K\rangle$. Since $K$ is semipositively defined, $K\geq 0$, it is
called the modular Hamiltonian (due to the relation $\rho=\exp(-K)$). For the
density matrix describing a part of a bipartite system (i.e. the reduced
density) the vN entropy becomes the basic measure of the entanglement; thus
sometimes the modular Hamiltonian is called the entanglement Hamiltonian or
entropy operator (such a nomenclature appears especially in the context of the
condensed matter physics).

Following the above observation we can consider analogues of thermodynamical
quantities with respect to $K$; in particular, the notion of entropy for $K$.
To this end we denote the inverse $\alpha=1/T$ of “temperature” $T$ (we use
$\alpha$ instead of thermodynamical $\beta$ index and the Boltzmann constant
is put one) and we define the partition function

$Z_{\alpha}=\textrm{tr}(e^{-\alpha K})=\textrm{tr}(\rho^{\alpha}),$ (2.4)

as well as the corresponding free energy

$F_{\frac{1}{T}}=F_{\alpha}=-\frac{1}{\alpha}\ln(Z_{\alpha}).$ (2.5)

Then, following the standard thermodynamical approach we construct the
corresponding entropy

$S_{\alpha}=-\partial_{T}F_{\frac{1}{T}}=\alpha^{2}\partial_{\alpha}F_{\alpha}.$
(2.6)

The new quantity $S_{\alpha}$ is called the modular entropy (due to its
relation to the modular Hamiltonian $K$) and can be easily compared with the
Rényi entropy

$S_{\alpha}=\alpha^{2}\partial_{\alpha}(\frac{\alpha-1}{\alpha}R_{\alpha})=R_{\alpha}+\alpha(\alpha-1)\partial_{\alpha}R_{\alpha}.$
(2.7)

Equivalently, it is the vN entropy of the modified density operator
$\rho_{\alpha}=\rho^{\alpha}/Z_{\alpha}$, i.e. $S_{\alpha}=\langle
K\rangle_{\alpha}$ where $\langle\cdot\rangle_{\alpha}$ stands for the
expectation value at $\rho_{\alpha}$. The last observation gives immediately
$S_{\alpha}\geq 0$ and $\lim_{\alpha\rightarrow 1}S_{\alpha}=S_{vN}$. Some
more properties of the modular entropy as well as its meaning in the context
of gravitational holography can be found in Refs. [14, 15, 16].

The above thermodynamical analogy can be continued by introducing the
(modular) capacity
$C_{\alpha}=-\partial^{2}_{T}F_{\frac{1}{T}}=-\alpha\partial_{\alpha}S_{\alpha}$.
Then one gets

$C_{\alpha}=\alpha^{2}\partial_{\alpha}^{2}((1-\alpha)R_{\alpha})=\alpha^{2}(\langle
K^{2}\rangle_{\alpha}-\langle K\rangle_{\alpha}^{2})\geq 0,$ (2.8)

i.e. $C_{\alpha}$ measures the variance of the modular Hamiltonian at
$\rho_{\alpha}$. Since for $\alpha=1$ the modular entropy coincides with the
vN one ($S_{1}=S_{vN}$), the capacity reduces to the variance of the modular
Hamiltonian $C\equiv C_{1}=\langle K^{2}\rangle-\langle K\rangle^{2}$, see
[25]-[30]. In consequence, it gives the quantum fluctuation with respect to
the original state. For the reduced density matrix $C$ measures entanglement
entropy fluctuations. The magnitude of the relative entanglement fluctuations
is defined as $\delta K=\sqrt{C}/S$ [25, 17]. As a border between strong and
weak fluctuations the condition $\delta K=1$ can be used, i.e. $C=S^{2}$.
Similarly, for an arbitrary $\alpha$ we have

$\langle K\rangle_{\alpha}=\partial_{\alpha}((\alpha-1)R_{\alpha}),$ (2.9)

and we can also easily express
$\delta_{\alpha}K=\sqrt{C_{\alpha}}/(\alpha\langle K\rangle_{\alpha})$ by
means of the Rényi entropy.

Concluding, let us note that for a finite-dimensional system and the density
matrix proportional to the identity (e.g. for the maximally entanglement
states) all the entropies $S,R_{\alpha},S_{\alpha}$ (for any $\alpha$)
coincide and $C_{\alpha}=\delta_{\alpha}K=0$. On the other hand, the maximal
value of the capacity appears for partially mixed states [25, 17]. For the
continuous variables systems the situation is more complicated and will be
discussed below.

##  3 Warm-up: time-dependent oscillator

In this section we examine the temporal evolution of the above quantities in a
very simple case. Namely, let us take the harmonic oscillator with the time-
dependent frequency $\omega(t)$. It is described by the Hamiltonian

$H(t)=\frac{p^{2}}{2}+\frac{1}{2}\omega^{2}(t)x^{2},$ (3.1)

for which the equation of motion reads

$\overset{..}{x}(t)=-\omega^{2}(t)x(t).$ (3.2)

The time-dependent harmonic oscillator (TDHO) appears in many physical models
and has been studied in various contexts. It turns out that eq. (3.2) is
equivalent to the so-called Ermakov-Milne-Pinney (EMP) equation [41, 42, 43]

$\overset{..}{b}(t)+\omega^{2}(t)b(t)=\frac{c^{2}}{b^{3}(t)},$ (3.3)

where $c$ is a constant (we assume $c\neq 0$ to ensure the non-vanishing of
the function $b(t)$). At the quantum level, the general solution to the
Schrödinger equation of the TDHO is a superposition of the following wave
functions

$\psi_{n}(x,t)=\frac{1}{\sqrt{2^{n}n!b(t)}}\sqrt[4]{\frac{c}{\pi}}e^{-ic(n+1/2)\tau(t)}H_{n}\left(\frac{\sqrt{c}x}{b(t)}\right)e^{-\frac{cx^{2}}{2b^{2}(t)}+\frac{i\dot{b}(t)x^{2}}{2b(t)}},$
(3.4)

where $H_{n}$ for $n=0,1\ldots$ denote the Hermite polynomials and
$\tau(t)=\int b^{-2}(t)$.

In physical considerations we usually assume that the states (3.4) at initial
time $t=t_{0}$ are eigenstates of the instantaneous Hamiltonian $H(t_{0})$.
This holds for the following initial conditions

$c=\omega(t_{0})\equiv\omega_{0},\quad b(t_{0})=1,\quad\dot{b}(t_{0})=0,$
(3.5)

and will be assumed in what follows.

Then the density function of the state $\psi_{n}(x,t)$ reads

$\rho_{n}(x,t)=\frac{1}{2^{n}n!b(t)}\sqrt{\frac{\omega_{0}}{\pi}}H_{n}^{2}\left(\frac{x\sqrt{\omega_{0}}}{b(t)}\right)e^{-\frac{\omega_{0}x^{2}}{b^{2}(t)}},$
(3.6)

while the density of the Fourier transform of $\psi_{n}(x,t)$ is given by the
formula

$\rho_{n}(p,t)=\frac{b(t)}{2^{n}n!\sqrt{\omega_{0}^{2}+b^{2}(t)\dot{b}^{2}(t)}}\sqrt{\frac{\omega_{0}}{\pi}}e^{-\frac{\omega_{0}b^{2}(t)p^{2}}{\omega_{0}^{2}+b^{2}(t)\dot{b}^{2}(t)}}H_{n}^{2}\left(\frac{\sqrt{\omega_{0}}b(t)p}{\sqrt{\omega_{0}^{2}+b^{2}(t)\dot{b}^{2}(t)}}\right).$
(3.7)

We start with the coordinate Rényi entropy

$R^{\alpha,x}_{n}(t)=\frac{1}{1-\alpha}\ln\left(\int\rho^{\alpha}_{n}(x,t)dx\right),$
(3.8)

as well as with their momentum counterpart $R^{\alpha,p}_{n}(t)$. For the TDHO
they read respectively

$\begin{split}R_{\alpha,n}^{x}(t)&=\ln\left(\frac{b(t)}{\sqrt{\omega_{0}}}\right)+\frac{1}{1-\alpha}\ln(W_{\alpha,n}),\\\
R_{\alpha,n}^{p}(t)&=\frac{1}{2}\ln\left(\frac{\omega_{0}^{2}+b^{2}(t)\dot{b}^{2}(t)}{\omega_{0}^{2}b^{2}(t)}\right)+\frac{1}{1-\alpha}\ln(W_{\alpha,n}),\end{split}$
(3.9)

where
$W_{\alpha,n}=\int_{-\infty}^{\infty}\left(\frac{H_{n}^{2}(y)}{\sqrt{\pi}2^{n}n!}\right)^{\alpha}e^{-\alpha
y^{2}}dy$ is the so-called entropic moment of the Hermite polynomials [44].
For lower states these quantities can be given explicitly

$W_{\alpha,0}=\sqrt{\frac{\pi^{1-\alpha}}{\alpha}},\quad
W_{\alpha,1}=\frac{2^{\alpha}}{\pi^{\alpha/2}\alpha^{\alpha+1/2}}\Gamma(\alpha+1/2);$
(3.10)

for higher $n$ they are a more complicated combination of the Bell polynomials
[44]. By virtue of eq. (2.7) we conclude that the coordinate modular entropy
is given by

$S_{\alpha,n}^{x}(t)=R_{\alpha,n}^{x}(t)+\frac{\alpha}{\alpha-1}\ln(W_{\alpha,n})-\alpha\partial_{\alpha}\ln(W_{\alpha,n}),$
(3.11)

and analogously for the momentum case. From eq. (3.11) we see that the time
dependence of the modular entropies $S_{\alpha,n}$ is the same as for the
Rényi ones; in particular, the increase of the modular entropy for an
arbitrary state $\psi_{n}$ does not depend on $n$, e.g.
$S_{\alpha,n}^{x}(t)-S_{\alpha,n}^{x}(t_{0})=\ln(b(t))=R_{\alpha,n}^{x}(t)-R_{\alpha,n}^{x}(t_{0})$
(obviously it vanishes for ordinary harmonic oscillator, $b(t)=1$). In
summary, for the states $\psi_{n}(t)$ both kinds of entropies differ by a
constant (depending on $\alpha$ and $n$) only.

Now, let us pass to the notion of the modular capacity. By virtue of eq. (2.8)
we have

$C_{\alpha,n}(t)=\alpha^{2}\partial^{2}_{\alpha}\ln(W_{\alpha,n});$ (3.12)

thus it is time-independent (and bounded for fixed $\alpha$ and $n$). For the
lowest (initial) state we have for any $\alpha$

$C_{\alpha,0}(t)=\frac{1}{2},$ (3.13)

while for the excited state the $\alpha$-dependence emerges

$C_{\alpha,1}(t)=\frac{1}{2}-\alpha+\alpha^{2}\Upsilon_{1}(\alpha+\frac{1}{2}),$
(3.14)

where $\Upsilon_{1}$ is trigamma function (the first derivative of the digamma
function). For the most interesting case $\alpha=1$ (corresponding to the
fluctuation of the vN entropy) we have

$C_{1,1}(t)=\frac{1}{2}(\pi^{2}-9).$ (3.15)

For higher $n$ only numerical results can be obtained. Finally, let us note
that the above considerations can be extended to the TDHO driven by an
external time-dependent force, i.e.

$\overset{..}{x}(t)=-\omega^{2}(t)x(t)+f(t).$ (3.16)

Namely, after direct calculations, we find that all the entropies reduce to
the ones for the force-free case.

In summary, in this section we showed that the dynamics of the modular entropy
of the basic solutions of the TDHO is the same as for the Rényi entropy. This
supports the point of view that the modular entropy can be considered as
another candidate for a generalization of the vN entropy. Moreover, we showed
that the capacity (fluctuation) of the vN entropy for such solutions of the
TDHO is time-independent and we found its form. In the next section we will
see that the situation changes when we consider more complicated systems
consisting of more oscillators and reduced densities.

##  4 Time-dependent coupled oscillators

In this section we investigate the notion of the modular entropy and the
capacity in the context of entanglement of the systems. To this end we
consider the system of TDHO coupled by a time-dependent parameter and next we
analyse various bipartite decompositions and corresponding reduced density
operators.

###  4.1 Two coupled oscillators

Let us consider the Hamiltonian of the form

$H(t)=\frac{1}{2}(p_{1}^{2}+p_{2}^{2})+\frac{1}{2}\omega^{2}(t)\left((x^{1})^{2}+(x^{2})^{2}\right)+\frac{1}{2}k(t)\left(x^{1}-x^{2}\right)^{2}.$
(4.1)

By means of the transformation ${x}=R{y}$ where $R$ is the 2-dimensional
rotation matrix with the rotation angle $\pi/4$, we transform the Hamiltonian
(4.1) into the following one

$H_{y}(t)=\frac{1}{2}(p_{1}^{2}+p_{2}^{2})+\frac{1}{2}\left(\omega^{2}_{1}(t)(y^{1})^{2}+\omega^{2}_{2}(t)(y^{2})^{2}\right),$
(4.2)

where now $p$’s denote the canonical momenta associated with $y$’s and
$\omega_{1}^{2}(t)=\omega^{2}(t)+2k(t),\quad\omega_{2}^{2}(t)=\omega^{2}(t)$.
The frequencies $\omega_{1,2}(t)$ determine the parameters of the initial
Hamiltonian (4.1) as follows

$\omega(t)=\omega_{2}(t),\quad
k(t)=\frac{1}{2}(\omega_{1}^{2}(t)-\omega_{2}^{2}(t)).$ (4.3)

The evolution $\psi_{0}(x,t)$ of the ground state $\psi_{0}(x,t_{0})$ of the
Hamiltonian operator $H(t_{0})$ as well as the reduced density matrix
$\rho_{0}^{red}(x^{1},\tilde{x}^{1},t)=\int\psi_{0}(x^{1},x^{2},t)\psi^{*}_{0}(\tilde{x}^{1},x^{2},t)dx_{2}$
can be easily computed when we take into account the form of the Hamiltonian
(4.2) and then return to the $x=(x^{1},x^{2})$ variable. The final result
reads [34]

$\rho_{0}^{red}(x^{1},\tilde{x}^{1},t)=\frac{1}{\sqrt{\pi(\zeta(t)-\chi(t))}}e^{\chi(t)x^{1}\tilde{x}^{1}+i((x^{1})^{2}-(\tilde{x}^{1})^{2})\varphi(t)-\frac{\zeta(t)}{2}((x^{1})^{2}+(\tilde{x}^{1})^{2})},$
(4.4)

where

$\varphi(t)=\frac{\dot{b}_{1}(t)}{4b_{1}(t)}+\frac{\dot{b}_{2}(t)}{4b_{2}(t)}-\frac{\frac{c_{1}}{b_{1}^{2}(t)}-\frac{c_{2}}{b_{2}^{2}(t)}}{\frac{c_{1}}{b_{1}^{2}(t)}+\frac{c_{2}}{b_{2}^{2}(t)}}\left(\frac{\dot{b}_{1}(t)}{4b_{1}(t)}-\frac{\dot{b}_{2}(t)}{4b_{2}(t)}\right),$
(4.5)

and $\zeta(t)>\chi(t)\geq 0$ are given by

$\zeta(t)=\frac{\left(\frac{c_{1}}{b_{1}^{2}(t)}+\frac{c_{2}}{b_{2}^{2}(t)}\right)^{2}+4\frac{c_{1}c_{2}}{b_{1}^{2}(t)b_{2}^{2}(t)}+\left(\frac{\dot{b}_{1}(t)}{b_{1}^{2}(t)}-\frac{\dot{b}_{2}(t)}{b_{2}^{2}(t)}\right)^{2}}{4\left(\frac{c_{1}}{b_{1}^{2}(t)}+\frac{c_{2}}{b_{2}^{2}(t)}\right)},\quad\chi(t)=\frac{\left(\frac{c_{1}}{b_{1}^{2}(t)}-\frac{c_{2}}{b_{2}^{2}(t)}\right)^{2}+\left(\frac{\dot{b}_{1}(t)}{b_{1}^{2}(t)}-\frac{\dot{b}_{2}(t)}{b_{2}^{2}(t)}\right)^{2}}{4\left(\frac{c_{1}}{b_{1}^{2}(t)}+\frac{c_{2}}{b_{2}^{2}(t)}\right)},$
(4.6)

while the functions $b_{1,2}(t)$ satisfy the EMP equation (3.3) with the
frequencies $\omega_{1,2}(t)$ and the constants
$c_{1}=\sqrt{\omega^{2}(t_{0})+2k(t_{0})}$ and $c_{2}=\omega(t_{0})$,
respectively. In consequence, the Rényi entropy and the vN entropy of the
reduced density $\rho_{0}^{red}$ can be easily found, see Ref. [34]

$R_{\alpha}(t)=\frac{1}{1-\alpha}\ln\frac{(1-\xi(t))^{\alpha}}{1-\xi^{\alpha}(t)},\quad
S(t)=-\ln(1-\xi(t))-\frac{\xi(t)}{1-\xi(t)}\ln\xi(t),$ (4.7)

where $\xi(t)=\frac{\chi(t)}{\zeta(t)+\sqrt{\zeta^{2}(t)-\chi^{2}(t)}}$.

Now, we are in the position to find the modular entropy and capacity. First,
the straightforward computations give

$S_{\alpha}(t)=-\ln(1-\xi^{\alpha}(t))-\frac{\xi^{\alpha}(t)}{1-\xi^{\alpha}(t)}\ln(\xi^{\alpha}(t)),$
(4.8)

i.e. $S_{\alpha}$ is obtained by the replacement
$\xi(t)\rightarrow\xi^{\alpha}(t)$ in the vN entropy, cf. eq. (4.7) (this can
be also seen from the fact the spectrum of $\rho_{0}^{red}$ consists of the
powers of $\xi$). Moreover, we have the following relation

$S_{\alpha}(t)=R_{\alpha}(t)+\frac{\alpha}{\alpha-1}\left(\ln\frac{1-\xi(t)}{1-\xi^{\alpha}(t)}+\frac{(1-\alpha)\xi^{\alpha}(t)\ln(\xi(t))}{1-\xi^{\alpha}(t)}\right);$
(4.9)

thus, in contrast to the single TDHO, the dynamics of the modular entropy for
the coupled oscillators is different from the Rényi one.

Similarly, for the capacity we get

$C_{\alpha}(t)=\frac{\xi^{\alpha}(t)\ln^{2}(\xi^{\alpha}(t))}{(1-\xi^{\alpha}(t))^{2}};$
(4.10)

in consequence, we have $0\leq C_{\alpha}(t)\leq 1$. In particular, the
capacity (fluctuation) of the entanglement is given by the formula

$C(t)=C_{1}(t)=\frac{\xi(t)\ln^{2}(\xi(t))}{(1-\xi(t))^{2}}.$ (4.11)

From the above we see that the discussed quantities are determined by the
solutions of the two EMP equations with the frequencies $\omega_{1,2}(t)$,
respectively.

In order to analyse the above results, let us start with the simplest case,
namely constant both the coupling $k(t)=k$ and the frequency
$\omega(t)=\omega$ in the Hamiltonian (4.1). Then $b_{1}(t)=b_{2}(t)=1$ and
$\xi(t)=\xi$ is a constant which can be easily expressed as a function of the
ratio $k/\omega^{2}$ (the same concerns $S_{\alpha},C_{\alpha}$); in
particular, $\xi\rightarrow 1$ for $k\rightarrow\infty$. In this case, the
reduced matrix $\rho_{0}^{red}$ is equivalent to the thermal density operator
for a single harmonic oscillator specified by the frequency
$\sqrt{\zeta^{2}-\chi^{2}}$ and temperature
$T=-\sqrt{\zeta^{2}-\chi^{2}}/\ln(\xi)$. The modular entropy is also the vN
entropy of the harmonic oscillator provided we define the temperature
$T=-\sqrt{\zeta^{2}-\chi^{2}}/(\alpha\ln(\xi))$. Moreover, the capacity
becomes the ordinary heat capacity of a quantum harmonic oscillator. Thus when
$\xi$ tends to one ($k$ to infinity) the entanglement fluctuations tend to one
while the entanglement entropy tends to infinity; this is in contrast with the
finite-dimensional case where the capacity vanishes for the maximally
entanglement states. For more coupled oscillators the reduced matrix is not
(in general) equivalent to a thermal density matrix for a system of
oscillators [5]. In what follows we will consider more general, time-
dependent, case.

###  4.2 Many coupled oscillators

The system of $N$ coupled time-dependent oscillators (with the nearest
neighbour interaction) is described by the Hamiltonian

$H(t)=\frac{1}{2}\sum_{l=1}^{N}(p_{l}^{2}+\omega^{2}(t)(x^{l})^{2})+\frac{1}{2}k(t)\sum_{l=0}^{N}(x^{l}-x^{l+1})^{2}.$
(4.12)

Imposing boundary conditions the Hamiltonian can be rewritten in the form

$H(t)=\frac{1}{2}\sum_{j=1}^{N}p_{j}^{2}+\frac{1}{2}x^{T}\Lambda(t)x,$ (4.13)

where $x=(x^{1},\ldots,x^{N})$ and $\Lambda(t)$ is a symmetric $N\times N$
matrix with the eigenvalues $\lambda_{j}(t)$, for $j=1,\ldots,N$. The most
popular boundary conditions are the following three:

  * •

The periodic boundary conditions (PBC) defined by $x^{0}=x^{1}$ and
$x^{N+1}=x^{0}$. Then the eigenvalues of $\Lambda(t)$ are of the form
$\lambda_{j}(t)=\omega^{2}(t)+4k(t)\sin^{2}(\frac{j\pi}{N})$, $j=1,\ldots,N$.

  * •

The Neaumann boundary conditions (NBC) given by $x^{0}=x^{1}$ and
$x^{N+1}=x^{N}$ (corresponding to the system of coupled pendulums in small-
angle approximation); then the eigenvalues read
$\lambda_{j}(t)=\omega^{2}(t)+4k(t)\sin^{2}(\frac{j\pi}{2N})$,
$j=0,\ldots,N-1$ (for $N=2$ we obtain the previous case, see Sec. 4.1).

  * •

Finally, we can impose the Dirichlet boundary conditions (DBC),
$x^{0}=x^{N+1}=0$, then
$\lambda_{j}(t)=\omega^{2}(t)+4k(t)\cos^{2}(\frac{j\pi}{2(N+1)})$,
$j=1,\ldots,N$.

As in the case of two oscillators, we can easily find the evolution of the
density matrix of the ground state of the instantaneous Hamiltonian
$H(t_{0})$. However, in contrast to the case of two oscillators, the reduced
density for more oscillators is substantial more complicated and some
misunderstanding appeared in the literature. In view of this, first, we
discuss this problem.

We start with the (time-dependent) density matrix of the whole system, see
e.g. [34, 45]

$\rho(x,x^{\prime},t)=\sqrt{\det(\Omega(t)/\pi)}\exp(ix^{T}B(t)x-i{x^{\prime}}^{T}B(t)x^{\prime}-\frac{1}{2}x^{T}\Omega(t)x-\frac{1}{2}{x^{\prime}}^{T}\Omega(t)x^{\prime}),$
(4.14)

where $x=(x^{1},\ldots,x^{N})$ and
$\Omega(t)=U^{T}\sqrt{\tilde{\Lambda}(t)}U$, $B=U^{T}\tilde{B}(t)U$ where
$\tilde{B}(t),\tilde{\Lambda}(t)$ are diagonal matrices with elements
$(\tilde{\Lambda}(t))_{ij}=\lambda_{i}(t_{0})/b_{i}^{4}(t)\delta_{ij}$ and
$(\tilde{B}(t))_{ij}=\dot{b}_{i}(t)/(2b_{i}(t))\delta_{ij}$, respectively,
while $b_{j}(t)$ are the solutions of the Ermakov-Milne-Pinney (EPM) equations
with the frequencies $\lambda_{j}(t)$

$\overset{..}{b_{j}}(t)+\lambda_{j}(t)b_{j}(t)=\frac{\lambda_{j}(t_{0})}{b^{3}_{j}(t)},\quad
j=1,\ldots,N;$ (4.15)

and, finally, $U$ is a time-independent matrix diagonalizing $\Lambda(t)$:

$U\Lambda(t)U^{T}=Diag(\lambda_{1}(t),\ldots,\lambda_{N}(t)).$ (4.16)

Next, we split the whole system into two parts: the first one $\mathcal{A}$
consisting of $n$ first oscillators $(x^{1},\ldots,x^{n})$ and the second one
$\mathcal{B}$ consisting of the remaining $N-n$ ones described by
${\bf{x}}=(x^{n+1},\ldots,x^{N})$ and rewrite $\Omega$ and $B$ in the form222
For simplicity of notation we omit the time parameter $t$ in the matrices when
it does not cause ambiguity.

$\Omega=\left(\begin{array}[]{cc}\Omega_{1}&\Omega_{2}\\\
\Omega_{2}^{T}&\Omega_{3}\end{array}\right),\quad
B=\left(\begin{array}[]{cc}B_{1}&B_{2}\\\ B_{2}^{T}&B_{3}\end{array}\right),$
(4.17)

where $\Omega_{1},B_{1}$ are $n\times n$ matrices.

Tracing (integrating) over the subsystem $\mathcal{A}$ we obtain the reduced
density of the subsystem ${\mathcal{B}}$. Namely, after straightforward
computations we arrive at the formula

$\rho_{\mathcal{B}}({\bf{x}},{\bf{x}}^{\prime},t)=A\exp(i{\bf{x}}^{T}Z{\bf{x}}-i{{\bf{x}}^{\prime}}^{T}Z{\bf{x}}^{\prime}-\frac{1}{2}{\bf{x}}^{T}\Upsilon{\bf{x}}-\frac{1}{2}{{\bf{x}}^{\prime}}^{T}\Upsilon{\bf{x}}^{\prime}+{\bf{x}}^{T}\Delta{\bf{x}}^{\prime}),$
(4.18)

where $Z,\Upsilon,\Delta$ are $(N-n)\times(N-n)$ matrices given by

$\displaystyle Z$ $\displaystyle=B_{3}-B_{2}^{T}\Omega_{1}^{-1}\Omega_{2},$
(4.19) $\displaystyle\Upsilon$
$\displaystyle=\Omega_{3}-\frac{1}{2}\Omega_{2}^{T}\Omega_{1}^{-1}\Omega_{2}+2B_{2}^{T}\Omega_{1}^{-1}B_{2},$
(4.20) $\displaystyle\Delta$
$\displaystyle=\frac{1}{2}\Omega_{2}^{T}\Omega_{1}^{-1}\Omega_{2}+2B_{2}^{T}\Omega_{1}^{-1}B_{2}+i\Theta,$
(4.21)

with
$\Theta=\Omega_{2}^{T}\Omega_{1}^{-1}B_{2}-B_{2}^{T}\Omega_{1}^{-1}\Omega_{2}$.
Such a form of the reduced density coincides with the one obtained in Ref.
[45]. The matrices $\Omega_{1},B_{1}$ are real and symmetric thus $\Theta$ is
skew-symmetric. In consequence, for $N=2$ (and $n=1$) $\Theta$ vanishes and we
obtain the results from previous section for two oscillators. However, for
higher $N$ (and $n$ ) $\Delta$ is a complex (but hermitian) matrix. In
consequence, we cannot directly apply methods from Ref. [5] to obtain the
spectrum of $\rho_{{\mathcal{B}}}$ (i.e. simultaneously diagonalize both
$\Upsilon$ and $\Delta$ by means of an orthogonal matrix). The term with
$\Theta$ has been missed in Refs. [34, 39] leading to some problems, e.g.,
despite the fact the initial state is pure the corresponding entropies
$S_{\mathcal{A}}$ and $S_{\mathcal{B}}$ were not equal to each other. For
similar reasons an approach presented in Ref. [31] seems also to miss
something; namely, in our notation the matrix $\mathcal{E}$ therein takes the
form $\Upsilon^{-1}\Delta$ and thus has a real (not complex as in [31])
spectrum (in fact the same as the hermitian matrix $\tilde{\Delta}$ below, eq.
(4.24)); moreover, the matrices $\mathcal{E}$ and its complex conjugate are,
in general, not simultaneously diagonalizable and thus an extension of the
reasoning from Ref. [46] is more involved. The problem of the spectrum of the
density operator in the case of hermitian, but not real, $\Delta$ has been
recently discussed, in the context of squeezed state, in Ref. [47]. Here we
adopt essential results and for more details we refer to [47].

Namely, it turns out that the spectrum of the reduced density matrix with
hermitian $\Delta$ is of the form

$(1-\xi_{1})(1-\xi_{2})\ldots(1-\xi_{N-n})\xi_{1}^{m_{1}}\xi_{2}^{m_{2}}\ldots\xi_{N-n}^{m_{N-n}},$
(4.22)

where $\xi$’s are the inverse of the eigenvalues (larger than one) of the
following matrix

$\left(\begin{array}[]{cc}2\tilde{\Delta}^{-1}&-\tilde{\Delta}^{-1}\tilde{\Delta}^{T}\\\
I&0\end{array}\right),$ (4.23)

where

$\tilde{\Delta}=(\tilde{\Upsilon})^{-1/2}\Pi\Delta\Pi^{T}(\tilde{\Upsilon})^{-1/2},$
(4.24)

while $\Pi$ is an orthogonal matrix diagonalizing $\Upsilon$, i.e.
$\Pi\Upsilon\Pi^{T}=\tilde{\Upsilon}$; let us note that $\tilde{\Delta}$ is a
hermitian matrix. Obviously, such an approach involves the invertibility of
the matrix $\Delta$ thus we can consider the subsystem ${\mathcal{B}}$
provided $N-n\leq n$, in the other case we should take the subsystem
${\mathcal{A}}$ and use the fact that for the pure state the spectrum of
$\rho_{\mathcal{B}}$ is equivalent (up to irrelevant zeros) to the spectrum
$\rho_{\mathcal{A}}$, see also the discussion in [47]; using a different
approach we will avoid this technical problem later. Now, from (4.22) we can
immediately find the Rényi entropies

$R_{\alpha}(t)=\sum_{j=n+1}^{N}R_{\alpha}[\xi_{j}(t)],$ (4.25)

where $R_{\alpha}[\xi_{j}(t)]$ has the form as in (4.7).

In consequence, the modular entropy as well as the capacity read

$S_{\alpha}(t)=\sum_{j=n+1}^{N}S_{\alpha}[\xi_{j}(t)],\quad
C_{\alpha}(t)=\sum_{j=n+1}^{N}C_{\alpha}[\xi_{j}(t)],$ (4.26)

where $S_{\alpha}[\xi_{j}(t)],C_{\alpha}[\xi_{j}(t)]$ are of the form (4.8)
and (4.10) (i.e. with the replacement $\xi(t)\rightarrow\xi_{j}(t)$). In
particular, we immediately obtain the formula for entanglement fluctuation
$C(t)\equiv C_{1}(t)$ of the system.

In what follows we will use another approach based on the correlation matrix
and symplectic spectrum; it avoids a direct computation of the spectrum and
has been successfully applied in the study of entanglement entropy and abrupt
quenches, see e.g. [32, 33, 38]. First, we define the time-dependent
correlations (covariance) matrix

$\Gamma=\left(\begin{array}[]{cc}Q&R\\\ R^{T}&P\end{array}\right),$ (4.27)

where $Q_{ij}=\langle x^{i}x^{j}\rangle$, $P_{ij}=\langle p_{i}p_{j}\rangle$,
$R_{ij}=1/2\langle\\{x^{i},p_{j}\\}\rangle$ are real and symmetric $N\times N$
matrices. Using the results of Sec. 3 we obtain, after straightforward
computations, that

$Q=U^{T}\tilde{Q}U,\quad P=U^{T}\tilde{P}U,\quad R=U^{T}\tilde{R}U;$ (4.28)

where $U$ is defined as in (4.16) and $\tilde{Q},\tilde{P},\tilde{R}$ are
diagonal matrices with the following diagonal elements

$\tilde{Q}_{kk}=\frac{b_{k}^{2}}{2\sqrt{\lambda_{k}(t_{0})}},\quad\tilde{P}_{kk}=\frac{1}{2}\left(\frac{\sqrt{\lambda_{k}(t_{0})}}{b_{k}^{2}}+\frac{\dot{b}_{k}^{2}}{\sqrt{\lambda_{k}(t_{0})}}\right),\quad\tilde{R}_{kk}=\frac{b_{k}\dot{b}_{k}}{2\sqrt{\lambda_{k}(t_{0})}},$
(4.29)

for $k=1,\ldots,N$. Since the matrix $\Theta$ is skew-symmetric the covariance
matrix of the reduced density operator of the system $\mathcal{B}$ is given by
the restriction of the initial one to the subsystem $\mathcal{B}$ i.e.

$\Gamma_{\mathcal{B}}=\left(\begin{array}[]{cc}Q_{\mathcal{B}}&R_{\mathcal{B}}\\\
R^{T}_{\mathcal{B}}&P_{\mathcal{B}}\end{array}\right),$ (4.30)

where $(Q_{\mathcal{B}})_{kl}=Q_{kl}$, $(R_{\mathcal{B}})_{kl}=R_{kl}$ and
$(P_{\mathcal{B}})_{kl}=P_{kl}$ for $k,l=n+1,\ldots,N$. Next, following the
procedure based on the symplectic transformation, see e.g. Refs. [33, 38, 48],
we construct the matrix $\tilde{\Gamma}_{\mathcal{B}}=iJ\Gamma_{\mathcal{B}}$,
where $J$ is given by

$J=\left(\begin{array}[]{cc}0&I\\\ -I&0\end{array}\right),$ (4.31)

and $I$ is $(N-n)\times(N-n)$ identity matrix. Now, the spectrum of the matrix
$\tilde{\Gamma}_{\mathcal{B}}$ consists of elements $\pm\gamma_{k}$,
$k=n+1,\ldots,N$ and the Rényi entropies take the form

$R_{\alpha}=\frac{1}{\alpha-1}\sum_{k=n+1}^{N}\ln\left((\gamma_{k}+1/2)^{\alpha}-(\gamma_{k}-1/2)^{\alpha}\right).$
(4.32)

By means of (4.32) we can easily find the modular entropy and capacity; for
example, the entanglement fluctuation reads

$C_{1}=C=\sum_{k=n+1}^{N}(\gamma_{k}^{2}-\frac{1}{4})(\ln(\gamma_{k}+1/2)-\ln(\gamma_{k}-1/2))^{2}.$
(4.33)

In view of this and eqs. (4.28), (4.29) the temporal evolution of the modular
entropy and capacity is algebraically determined by the solutions (and their
derivatives) of the EMP equations. Moreover, the quantities obtained in this
way coincide with the formulae (4.25) and (4.26) after the identification

$2\gamma_{k}=\frac{1+\xi_{k}}{1-\xi_{k}}.$ (4.34)

Finally, let us note that $0\leq C(t)\leq N-n$. Moreover, the state is pure
thus $0\leq C(t)\leq n$. In consequence, for the initial ground state the
capacity remains bounded

$0\leq C(t)\leq\min(n,N-n).$ (4.35)

Concluding, eqs. (4.8) for $N=2$ and (4.26) (equivalently (4.32) and (4.33))
for higher $N$ enable us to analyse the evolution more explicitly provided we
have solutions $b_{j}(t)$ to the EMP equation. Such a situation holds for some
special forms of frequencies. In what follows we will discus both abrupt and
continuous examples and concentrate on the entanglement fluctuation and its
relation to the vN entropy; this completes considerations for finite-
dimensional spaces and fit into the field theory considerations.

###  4.3  Frequency and coupling jumps

Although the general dynamics of the entropy and capacity in the quench
phenomena is complicated, some qualitative description can be obtained if we
know the explicit form of the solutions to the EPM equations. To see this let
us consider the abrupt change. Namely, let us analyse this evolution for the
model (4.12) with a sudden quench, i.e. where $\omega(t),k(t)$ change, at time
$t_{0}=0$, from constant values $(\omega(i),k(i)))$ to another constant values
$(\omega(f),k(f)$. For the abrupt quench the solutions of the EPM equations
(4.15) with the initial conditions (3.5) read

$b_{j}(t)=\sqrt{r_{j}\cos(2t\sqrt{\lambda_{j}(f)})+s_{j}},$ (4.36)

where $\lambda_{j}(i),\lambda_{j}(f)$ are the eigenvalues of $\Lambda$ before
and after quench and
$r_{j}=(\lambda_{j}(f)-\lambda_{j}(i))/(2\lambda_{j}(f))$,
$s_{j}=(\lambda_{j}(f)+\lambda_{j}(i))/(2\lambda_{j}(f))$.

Both the entanglement and capacity are functions of $b(t)$’s and their
derivatives. Thus, due to the formula (4.36), we expect some distinguished
periodicities in the temporal evolution.

For the PBC the last frequency $\lambda_{N}(f)=\omega^{2}(f)$ is
distinguished. Indeed, it tends to zero as $\omega(f)$ is small, in contrast
to the previous frequency
$\lambda_{N-1}(f)=\omega^{2}(f)+4k(f)\sin^{2}((N-1)\pi/N)$ when
$k(f)\gg\omega^{2}(f)$. Under these assumptions $r_{N},s_{N}$ tend to infinity
as $\omega(f)$ is close to zero; thus due to eq. (4.36) we expect a
distinguished periodicity ${\mathbb{T}}=\pi/\omega(f)$ (in particular, it does
not depend on $k(f)$). This situation is presented in Fig. 1, where $N=4$ and
$\omega(f)=0.01$ (i.e. ${\mathbb{T}}\simeq 314$). Let us note that for the
capacity (which is bounded) this is less evident, especially, as we pointed
out above, for small $k(f)$, see the right panel in Fig. 1. Finally, it is
worth to notice that for $N\rightarrow\infty$ the situation changes:
$\lambda_{N}$ is not distinguished, and the above reasoning breaks down; we
return to this issue in the next sections.

Figure 1: Abrupt quench - temporal evolution: $N=4,n=1$ and PBC
($\omega(i)=3$, $\omega(f)=0.01$ ). Lines: blue - entropy, red - capacity. The
left panel $k(f)=10$. The right panel $k(f)=0.5$

For the DBC the situation is quite different. Namely, for $\omega(f)=0$ the
$\lambda$’s remain non-zero. However, for large $N$ the frequency
$\lambda_{N}$ tends to zero implying that the distinguished periodicity should
be of the form

${\mathbb{T}}=\frac{\pi}{2\sqrt{k(f)}\cos(\pi N/(2(N+1)))},$ (4.37)

thus it depends on $k(f)$ (and does not depend on $k(i)$; in particular, we
can put $k(i)=k(f)=k$). For example, for $N=100$ and $k(f)=1$ we have
${\mathbb{T}}\simeq 100$ while for $k(f)=9$ we obtain ${\mathbb{T}}\simeq 33$;
see Fig. 4 and 5 in the Sec. 5.2.1. In general, we have an approximation
${\mathbb{T}}\simeq\frac{N}{\sqrt{k}}$ for large $N$, thus for $k=1$ (lattice
approximation) we have ${\mathbb{T}}\simeq N$, see the considerations in the
next section.

##  5 Fluctuation of entanglement and field theory

In this section we consider larger values of $N$. This is interesting due to
the lattice approximation to quantum fields, e.g. a time-dependent bosonic
scalar field $\phi$. Namely, we discretize the field Hamiltonian, in general
with time-dependent mass $m(t)$,

$H=\frac{1}{2}\int dx(\pi^{2}+(\nabla\phi)^{2}+m^{2}(t)\phi^{2}),$ (5.1)

into a chain of harmonic oscillators by imposing a UV cutoff $\epsilon$ and IR
cutoff $N\epsilon$. Then we obtain a model described by the Hamiltonian (4.12)
with $k=1$ and $\omega^{2}(t)=m^{2}(t)\epsilon^{2}$. To complete the
discretization procedure we need boundary conditions; we take the PBC for the
spacetime cylinder; the other choices are the DBC or NBC.

###  5.1 Constant frequency

For start we consider the constant mass, which can be realized as a special
case of the abrupt quench with $\omega(i)=\omega(f)=\omega=const$. We will
focus on more subtle case when $\omega\rightarrow 0$ which corresponds to the
2-dimensional CFT. In this case the PBC (or NBC) leads to zero modes which
cause the divergence of the entanglement entropy. For the DBC there are no
zero modes and no infrared divergences arise when setting $\omega$ equals zero
(the divergence is in the limit $N\rightarrow\infty$ only). For the
entanglement capacity the situation is quite different. Namely, as we noted
above the entanglement capacity is bounded thus for $\omega=0$ we expect a
finite value.

Using methods of the 2-dimensional CFT it has been shown [24] that at the
leading order of the UV cutoff parameter $\epsilon$ the capacity and entropy
are equal to each other. Namely, for the PBC (circle of the length $L$) we
have

$S=\frac{1}{3}\ln\left(\frac{L}{\pi\epsilon}\sin\left(\frac{\pi
l}{L}\right)\right)+O(1)=C.$ (5.2)

In particular, for the infinite length $L\rightarrow\infty$ we obtain
$S=\frac{1}{3}\ln(l/\epsilon)+O(1)=C$; such a behavior has been confirmed in
Ref. [24] by considering the infinite harmonic chain. Here, we analyse this
problem for $L$ finite: $L=N\epsilon$ (and $l=n\epsilon$). For the DBC instead
of eq. (5.2) we expect the relation
$S=\frac{1}{6}\ln\left(\frac{L}{\pi\epsilon}\sin\left(\frac{\pi
l}{L}\right)\right)+O(1)=C$. In the former case we take $\omega=0.001$ and
next analyse the limit as $\omega\rightarrow 0$, in the Dirichlet case we can
put directly $\omega=0$ from very beginning.

The results for the PBC are presented in the left panel of Fig. 2 (for
$N=100$).

Figure 2: $N=100$ and PBC ($\omega=0.001$ and $k=1$). The left panel: entropy
- blue data points (green line corresponds to (5.2)), capacity - red data
points (yellow line corresponds to (5.2)). The right panel: omega dependence,
entropy - blue date points, capacity - red date points.

We see that the entropy and capacity coincide with eqs. (5.2) (represented
here by green and yellow lines, respectively). For the small ratio $n/N=l/L$
this formula exhibits a logarithmic scaling, i.e. it diverges, in contrast to
the area law. Moreover, for the PBC the entropy increases (for fixed $n$ and
$N$) as $\omega$ tends to zero as is presented in the right panel of Fig. 2,
in contrast to the capacity which is bounded (cf. eq. (4.35)).

For the DBC we can put directly $\omega=0$, then both entropy and capacity are
finite and their behavior agrees with the above theoretical prediction (i.e.
eq. (5.2) with $1/3$ replaced by $1/6$).

###  5.2 Quenches

Now, let us pass to the case of oscillatory chains with time-dependent
$\omega(t)$. Such models can approximate quench protocols in the field theory
where the mass varies in time or in expanding spacetimes. We will focus on the
more interesting case where the final theory is massless i.e. the two-
dimensional conformal field theory. Due to zero modes we will consider DBC and
PBC separately; in the latter case, the final frequency will be small but not
equal to zero.

####  5.2.1 Abrupt quench

First, we consider the abrupt quenches; namely, we start with some $\omega(i)$
and at time $t_{0}=0$ then there is a sudden change of the frequency to a
small or zero value. Such a situation is motivated by the CFT methods
developed in the study of the temporal evolution of the Rényi entropies. For
the global quench and the subsystem which is an interval on the infinite line
the field approach yields [24]

$S=\textrm{const}+\left\\{\begin{array}[]{cc}a_{S}t\quad t<t^{*},\\\
b_{S}n\quad t>t^{*};\end{array}\right.\quad
C=\textrm{const}+\left\\{\begin{array}[]{cc}a_{C}t\quad t<t^{*},\\\
b_{C}n\quad t>t^{*};\end{array}\right.$ (5.3)

with some constants $a_{S},b_{S},a_{C},b_{C}$ and $t^{*}=n/2$; moreover, the
relation $a_{S}=a_{C}$ holds. However, the analysis performed in [24] suggests
that the situation is more subtle and the slopes $a_{C}$ and $a_{S}$ are not
equal (and, consequently, the saturation values too). Here, using the results
presented in Sec. 4, we analyse this problem in the finite case and various
boundary conditions; in particular, we obtain some theoretical predictions for
these constants.

We start with temporal evolution of the entropy and capacity for various $n$
and PBC. In Fig. 3 (top panels), we see that there is indeed the linear growth
of the both quantities (more precisely, just after the quench there is
parabolic-like growth); however, as in Ref. [24], we observe that the slopes
$a_{S}$ and $a_{C}$ are, in general, not equal; it depends on the initial
frequency. Approximately, after time $t>n/2$ and for lower $n$ there is a
logarithmic-type growth instead of the constant value (both for the entropy
and capacity); this can be related to zero modes appearing for the PBC, see
e.g. the discussion in [33] and referees therein. Next, both quantities
oscillate; after longer time these oscillations are around saturation values
(bottom panels in Fig. 3). Note that the values of both the entropy and
capacity increase with $n$ ($l$); namely, for $n<N/2$ they are less than for
$n=N/2$ (for $n>N/2$ they coincide with the ones for $N-n$). Finally, let us
note that the dynamics of capacity is bounded as in eq. (4.35).

Figure 3: Abrupt quench - temporal evolution: $N=100$ and PBC ($\omega(i)=3$,
$\omega(f)=0.01$ and $k=1$). Lines: black $n=6$ , blue $n=10$, green $n=20$,
yellow $n=36$, red $n=50$. The left panels entropy. The right panels capacity.

Let us compare these results with the ones for the DBC. In this case we can
put the final frequency zero $\omega(f)=0$ (the entropy is finite). Then, in
agreement with the discussion presented in Sec. 4.3, we observe the
distinguished periodicity ${\mathbb{T}}\simeq N/\sqrt{k}$ for large $N$;
namely, for $k=1$ ($N=100$) we have ${\mathbb{T}}\simeq 100$ , see Fig. 5,
while for for $k=9$ it gives ${\mathbb{T}}\simeq 33$, see Fig. 4. This special
periodicity is related to local minima, such that the quaintness coincide for
all $n$ (especially when the initial frequency is not too large, see see Fig.
4), thus at these points we can observe the area law ($n^{0}$). In contrast to
this, at local maxima or for further times the area law is broken. Moreover,
the initial linear growth is up to the first maximum which, see Sec. 4.3,
depends on $k$ (in contrast to the PBC case), it is around the half of the
period ${\mathbb{T}}$, i.e. $t=N/2$ (in agreement with the maximal linear
growth resulting from the field approach, i.e. $t=n=N/2$). For small $n$ the
dynamics exhibits plateaus (instead of maxima) and finally oscillates around
saturation values.

Figure 4: Abrupt quench - temporal evolution: $N=100$ and DBC ($\omega(i)=3$,
$\omega(f)=0$ and $k=9$). Lines: black $n=6$ , blue $n=10$, green $n=20$,
yellow $n=36$, red $n=50$. The left panels entropy. The right panels capacity.
Figure 5: Abrupt quench - temporal evolution: $N=100$ and DBC ($\omega(i)=3$,
$\omega(f)=0$ and $k=1$). Lines: black $n=6$ , blue $n=10$, green $n=20$,
yellow $n=36$, red $n=50$. The left panels entropy. The right panels capacity.

These results can be confirmed by considering time-fixed slices, see Fig. 6.
Then just after the quench the behavior of the entropy and capacity resembles
the one for the finite mass, i.e. it is constant w.r.t. $n$; the area law
holds. For larger time the entropy (capacity) behaves linearly for larger $n$;
however, around $n=N/2$ this behavior changes. This is related to the finite
size of the whole system; instead of the linear growth rather $\sin(\pi n/N)$
appears (see green data points in Fig. 6). Moreover, the symmetry
$n\rightarrow N-n$ is preserved for entropy and capacity. In consequence, the
entropy and capacity of the subsystems ${\mathcal{A}}$ and ${\mathcal{B}}$
coincide (this agrees with the results obtained for entropy in [38] and the
fact that the state is pure).

Figure 6: Abrupt quench - time slices, $N=100$ and DBC ($k=1$, $\omega(f)=0$)
and $t=10150$. Date points: blue - entropy, red - capacity. Left panel
$\omega(i)=1$ and green date points: $10\sin(\pi n/100)+0.1$. Right panel
$\omega(i)=10$.

Finally in order to compare better the dynamics of the entropy and capacity we
plot the magnitude of relative quantum entanglement fluctuations $\delta
K=\sqrt{C}/S$, see Sec. 2. In Fig. 7 we see that after some initial time the
relative entanglement fluctuations are small and exhibit some revival time;
moreover, they are larger for lower $n$ (for the DBC the maxima are related to
the distinguished periodicities discussed above).

Figure 7: Abrupt quench - temporal evolution of relative fluctuations, $N=100$
($k=1$). Lines: black $n=6$ , blue $n=10$, green $n=20$, yellow $n=36$, red
$n=50$. The left panel: PBC ($\omega(i)=3$, $\omega(f)=0.01$). The right
panel: DBC ($\omega(i)=3$, $\omega(f)=0$).

Now, let us return to the field theory picture. Based on the quasiparticles
considerations we expect the relation (5.3) to hold (at least for small $n$).
To this end let us compute the initial slope of $S(t)-S(0)$ as well as
$C(t)-C(0)$ w.r.t. $\omega(i)\geq 0.1$ for the DBC and PBC (in the latter case
$\omega(f)=0.01$ is fixed). The results are presented in Fig. 8. First, we see
that the slope coefficient behaves linearly for $0.1\leq\omega(i)<1$ and it is
smaller for the entropy than for the capacity. However, for the capacity which
is slowly increasing, the linearity is broken earlier; thus to analyse it we
restrict ourselves to the interval $0.1\leq\omega(i)\leq 0.5$. More
preciously, the dashed lines in Fig. 8 for the PBC are given by
$0.575\omega(i)+0.003$ and $0.657\omega(i)+0.0004$ for the entropy and
capacity, respectively; while for the DBC we have $0.287\omega(i)+0.001$ and
$0.329\omega(i)+0.0004$, respectively. In view of this the slope $a_{S}=0.575$
for the entropy and PBC is two times greater than for the DBC ($2a_{S}=0.574$)
and agrees with theoretical prediction $a_{S}^{t}=(\pi-2)/2\simeq 0.571$
presented in [33]; for the capacity we have also the same relation ($0.657$
and $2a_{C}=0.658$). Even more in the non-linear region the slope for the PBC
is twice the one for the DBC (both for the entropy and capacity), cf. left and
right panels in Fig. 8. Next, as $\omega(i)$ becomes greater than one the
linearity is more and more broken and both slopes meet; then the entropy
increases while capacity stabilizes. In consequence, for large $\omega(i)$ the
entropy slope is greater than for the capacity.

Figure 8: The initial slope for the abrupt quench, $N=100$. Lines: blue -
entropy, red - capacity. Left panel PBC ($\omega(f)=0.01$); right panel DBC
($\omega(f)=0$).

Let us return to the slope of the capacity. Using the considerations from Ref.
[33] (basing on the quasiparticles approach [49]) and formula (2.8) we find
the initial slope for the capacity and the PBC:

$a_{C}^{t}=\partial_{\alpha}^{2}(\alpha-\cot(\pi/4\alpha))(1)=\pi(1-\frac{\pi}{4})\simeq
0.674,$ (5.4)

thus it quite well agrees with the above numerical result $a_{C}=0.657$
obtained by the lattice approximation. It remains to analyse the behavior for
further time; then for the PBC $b_{S},b_{C}$ we should have $b_{S}=a_{S}/2$
and $b_{C}=a_{C}/2$ while for the DBC $b_{S}=a_{S}$ and $b_{C}=a_{C}$. In Fig.
9 we present this situation using $a_{S}^{t}$ and $a_{C}^{t}$ coefficients
(for the PBC and DBC, respectively). Let us stress that these theoretical
predictions match quite well only for some range of the initial frequencies,
i.e. around $\omega(i)=1$ for the entropy and $\omega(i)=0.5$ for the
capacity, see Fig. 9 for $\omega(i)=0.5$; such a restriction is related to the
validity of lattice approximation ($\omega(i)$ should be not too small or
large, see e.g. the discussion in Refs. [32, 33]) and the fact the
quasiparticles model gives a qualitative description. A further investigation
of the field theoretical considerations presented recently in Refs. [24, 40]
can give more insight into this problem.

Summarizing, for the abrupt quench and small $n$ the relation (5.3)
approximates the dynamics of the entropy and capacity. More precisely, we have
a linear initial growth and for further times the oscillations are about
saturation values (up to logarithmic period for the PBC); however, the initial
slope coefficients are different for the entropy and capacity. For some range
of initial frequencies the capacity slope can be also obtained by means of the
quasiparticles picture.

Figure 9: Abrupt quench - temporal evolution: $N=100$ ($\omega(i)=0.5$ and
$k=1$). The top panels PBC ($\omega(f)=0.01$). The bottom panels DBC
($\omega(f)=0$). The left panels entropy. The right panels capacity. Lines:
black $n=4$ , blue $n=8$, green $n=12$, yellow $n=16$. Dashed lines correspond
to eq. (5.3) with theoretical parameters $a^{t}_{S},a^{t}_{C}$ .

##  6 Continuous protocols

In this section we analyse the aforementioned issues in the case of continuous
protocols. To this end we consider two examples which enable us more analytic
considerations. First, there is a continuous frequency change at time
$t_{0}=0$ such that it asymptotically tends to zero. In the second case we
consider the situation when at minus and plus infinities the frequency is zero
(i.e. approximation of CFT ). For better transparency we describe the DBC and
we mention only the necessary changes for the PBC.

To this end we put $k=1$ and take

$\omega^{2}(t)=\frac{2}{\varepsilon^{2}\cosh^{2}(t/\varepsilon)},$ (6.1)

in the Hamiltonian (4.12). Then $\omega(t)$ is a bell shaped function with the
maximum at $t=0$ and tending to zero at infinities. In particular, we obtain
the Dirac delta limit

$\lim_{\varepsilon\rightarrow 0}\omega(t)=\pi\sqrt{2}\delta(t).$ (6.2)

In order to obtain a continuous (differentiable) protocol we take

$\tilde{\omega}^{2}(t)=\left\\{\begin{array}[]{cc}\omega^{2}(i)\equiv
2/\varepsilon^{2}&\quad\textrm{for}\quad t\leq 0,\\\
\omega^{2}(t)&\quad\textrm{for}\quad 0<t.\end{array}\right.$ (6.3)

Then the initial frequency $\omega(i)$ is quenched to zero at infinity. To
analyse such a model we need the solutions of the EPM equations (4.15) with

$\lambda_{j}(t)=\tilde{\omega}^{2}(t)+a_{j},$ (6.4)

where, for the DBC, $a_{j}=4k\cos^{2}(\frac{j\pi}{2(N+1)})>0$, $j=1,\ldots,N$.
The corresponding functions $\tilde{b}_{j}(t)$ can be given explicitly.
Indeed, by means of the results of Ref. [50], we obtain

$\tilde{b}_{j}(t)=\left\\{\begin{array}[]{cc}1&\quad\textrm{for}\quad t\leq
0,\\\ b_{j}(t)&\quad\textrm{for}\quad 0<t;\end{array}\right.$ (6.5)

where $b_{j}(t)$, for $j=1,\ldots,N$, are of the form

$b_{j}^{2}(t)=\left(1+\frac{\tanh^{2}(t/\varepsilon)}{a_{j}^{2}\varepsilon^{2}}\right)\left(1-\frac{\sin^{2}\left(a_{j}t+\tan^{-1}(\frac{\tanh(t/\varepsilon)}{a_{j}\varepsilon})\right)}{(1+a_{j}^{2}\varepsilon^{2})^{2}}\right).$
(6.6)

Note that for the PBC we have $a_{N}=0$ and thus we have to modify $b_{N}$ as
follows

$b^{2}_{N}(t)=\left(1-\frac{t}{\varepsilon}\tanh(t/\varepsilon)\right)^{2}+2\tanh^{2}(t/\varepsilon).$
(6.7)

Now, using the results of Sec. 4 we will analyse the entropy and capacity of
the entanglement. To make contact with our previous considerations we take
$\varepsilon=\sqrt{2}/3$ which corresponds to the initial frequency
$\omega(i)=3$. Let us compare it with the abrupt case. First, we should note
that $\omega(t)$ tends (is not equal) to zero at infinity, so we cannot use
the argument presented in Sec. 4.3. In consequence, it is more difficult to
identify a distinguished periodicity and the temporal evolution is more
complicated as shown in Figs. 10 and 11. Despite the identical initial
conditions the values of the capacity are smaller than those for the abrupt
quench; moreover, the transition from the linear to the oscillatory regime is
much smoother and occurs earlier compared to the abrupt quench. However, there
are minima where the area law approximately holds (especially for smaller
$\omega(i)$ i.e. large $\varepsilon$). For larger times, we observe the
oscillations about a saturation value similarly to the abrupt case.

Figure 10: (6.3) quench - temporal evolution: $N=100$ and DBC ($\omega(i)=3$
and $k=1$). Lines: black $n=6$ , blue $n=10$, green $n=20$, yellow $n=36$, red
$n=50$. The left panel entropy. The right panel capacity.  Figure 11: (6.3)
quench - temporal evolution: $N=100$ and DBC ($\omega(i)=3$ and $k=1$). Lines:
black $n=6$ , blue $n=10$, green $n=20$, yellow $n=36$, red $n=50$. The left
panels entropy. The right panels capacity.

The second interesting case is given by the initial conditions (3.5) imposed
at $t_{0}=-\infty$, i.e. when the frequency is given directly by eq. (6.1).
Then, we start and end with massless (in the field theory picture) case with
the peak $2/\varepsilon^{2}$ at time $t=0$. In this case the functions
$b_{j}(t)$ are given by

$b_{j}^{2}(t)=\frac{a_{j}^{2}\varepsilon^{2}+\tanh^{2}(t/\varepsilon)}{1+a_{j}^{2}\varepsilon^{2}},$
(6.8)

for $j=1,\ldots,N$. A remarkable property of the functions (6.8) is that they
satisfy the condition (3.5) also at $t=\infty$; there is no oscillatory
behavior for $a_{j}>0$ at plus infinity (independently of the value of
parameter $\varepsilon$). Such a property implies that we observe only a peak
in both the entropy and capacity, see Fig. 12. More insight can be obtained
when we consider the time-fixed slices. For sufficiently large (positive and
negative) time the behavior of both the entropy and capacity coincide with the
ones presented in Sec. 5.1, which corresponds to the CFT picture. However,
near $t=0$ this behavior becomes disordered, as depicted in the left panel
Fig. 13. At $t=0$ it resembles the massless case again, see the right panel in
Fig. 13.

Finally, as $\varepsilon\rightarrow 0$ the frequency tends to the Dirac delta.
Then the numerical results for $\varepsilon=0.0001$ and $t\neq 0$ give the
entropy and capacity as in Sec. 5.1, while for $t=0$ we obtain behavior
similar to the conformal case; however, with different parameters. Namely,
instead of the factor $1/6$ (see Sec. 5.1) we have $N$-dependent parameter;
e.g. for $N=100$ we have $0.77$ for entropy and $0.2$ for capacity,
respectively; the suitable curves are depicted in Fig. 14.

Figure 12:  (6.1) quench with $\varepsilon=5$ \- temporal evolution: $N=100$
and DBC ($k=1$). Lines: black $n=6$ , blue $n=10$, green $n=20$, yellow
$n=36$, red $n=50$. The left panel entropy. The right panel capacity.  Figure
13: (6.1) quench with $\varepsilon=5$ \- time slices, $N=100$ and DBC ($k=1$).
Lines: blue - entropy, red - capacity. Left panel $t=1$. Right panel $t=0$.
Figure 14: (6.1) quench with $\varepsilon=0.0001$ (Dirac delta limit) - time
slices $t=0$, $N=100$ and DBC ($k=1$). Date points: blue - entropy, red -
capacity. Green line: $0.77\ln((N/\pi)\sin(\pi n/100))+1.3$, yellow line:
$0.2\ln((N/\pi)\sin(\pi n/100))+0.95$.

##  7 Summary and final discussion

In this work we have studied the notion of the modular entropy and capacity
(in particular the fluctuations of the entanglement entropy) in more detail.
This is motivated by the recent investigations of both the quantities in
various physical contexts; e.g. quantum gravitational effects. To this end we
considered systems of oscillators with the time-dependent frequency coupled by
a parameter. Such models, due to the discretization procedure, can be used to
analyse field theory problems resulting from quench phenomena or non-static
gravitational metrics. First, we showed that the modular and Rényi entropies
have the same dynamics for the basic solutions of the single TDHO. Moreover,
in this case the capacity is time-independent and we found its form (eqs.
(3.12)-(3.15)).

Next, we studied the aforementioned notions for a finite number of the TDHO
and various bipartite decompositions. As a result we obtained analytic
formulae for them; in contrast to the single oscillator case the dynamics of
modular entropy is different from the Rényi one and the capacity becomes time-
dependent (eqs. (4.8) and (4.10)). Next, we focused on the capacity
(fluctuation) of the entanglement, which has recently gained increasing
attention. In addition to the work [24] we considered various boundary
conditions related to the approximation of the finite size system (in field
theory picture). Taking $n\ll N$ and abrupt quenches they coincide with the
observations made in Ref. [24]. Obviously, we observed (for $n\sim N/2$) some
differences due to finite-size assumption (e.g. oscillations, in general non-
linear behavior w.r.t. subsystem size, and local area laws), see Figs. 3-7. We
confirmed, in agreement with general theory, that the capacity (as entropy) is
symmetric with respect to both subsystems; moreover, for fixed time they
increase with $n$ (for $n<N/2$). We analysed also the relative fluctuation of
entanglement and, in Sec. 4.3, some conditions which give rise to
distinguished periodicities. Special attention has been paid to the
theoretical prediction resulting from the quasiparticles model; the initial
slope coefficients, obtained in this way, are different for the entropy and
capacity and agree with the ones from lattice approximation (at least for some
range of parameters).

Finally, we studied some of the above issues in the case of continuous
protocols. To this end, in Sec. 6, we considered frequencies which vanish at
plus (and minus) infinity. In particular, we examined the model in which the
frequency tends to the Dirac delta. In such a case the behavior of entropy and
capacity, except one point (time zero), coincide with the CFT results (see
Sec. 5.1); fairly surprisingly, at time zero we obtain similar behavior,
however, with $N$-dependent coefficient, see Fig. 14.

All the above models have been discussed in the analytical manner, due to the
explicit forms (given by elementary functions) of solutions of the EMP
equations (for higher $N$ numerical computations were used to the standard
matrix algebra only).

Turning to possible further developments, it would be interesting to extend
the considerations presented to the transition from the DBC to NBC (such a a
scenario can be used to simulate the dynamical Casimir effect [39]) and/or
multiply quenched harmonic chains [36]. The other states can be also
considered, e.g. the coherent or squeezed ones, [47, 51]. Additionally, some
extensions to larger dimensions are also interesting (see [33, 52] and
references therein) or higher derivative theories [53]. In this context,
particularly interesting are applications in field theory in expanding
spacetimes [54, 55]. Finally, basing on the methods presented in Refs. [40,
56], a deeper analysis of the capacity (fluctuations of the entropy) would be
worthwhile.

Acknowledgment

The author would like to thank Piotr Kosiński for useful comments.

## References

  * [1] P. Calabrese, J. Cardy, “Entanglement Entropy and Quantum Field Theory” J. Stat. Mech. 0406 (2004) P06002
  * [2] P. Calabrese, J. Cardy, “Entanglement entropy and conformal field theory” J. Phys. A 42 (2009) 504005
  * [3] H. Casini, M. Huerta, “Entanglement entropy in free quantum field theory” J. Phys. A 42 (2009) 504007
  * [4] L. Bombelli, R. Koul, J. Lee, R. Sorkin, “Quantum source of entropy for black holes” Phys. Rev. D 34 (1986) 373
  * [5] M. Srednicki, “Entropy and Area” Phys. Rev. Lett. 71 (1993) 666
  * [6] S. Solodukhin, “Entanglement entropy of black holes” Living Rev. Relat. 14 (2011) 8
  * [7] J. Eisert, M. Cramer, M. Plenio, “Colloquium: Area laws for the entanglement entropy” Rev. Mod. Phys. 82 (2010) 277
  * [8] S. Ryu, T. Takayanagi, “Holographic Derivation of Entanglement Entropy from AdS/CFT” Phys. Rev. Lett. 96 (2006) 18160
  * [9] H. Li, F. Haldane, “Entanglement Spectrum as a Generalization of Entanglement Entropy: Identification of Topological Order in Non-Abelian Fractional Quantum Hall Effect States” Phys. Rev. Lett. 101 (2008) 010504
  * [10] P. Calabrese, A. Lefevre, “Entanglement spectrum in one-dimensional systems” Phys. Rev. A 78 (2008) 032329
  * [11] G. Di Giulio, E. Tonni, “On entanglement hamiltonians of an interval in massless harmonic chains” J. Stat. Mech. (2020) 033102
  * [12] H. Yao, X.-L. Qi, “Entanglement entropy and entanglement spectrum of the Kitaev model” Phys. Rev. Lett. 105 (2010) 080501
  * [13] D. Blanco, H. Casini, L.-Y. Hung, R. Myers, “Relative Entropy and Holography” JHEP 2013 (2013) 60
  * [14] X. Dong, “The gravity dual of Rényi entropy” Nature Commun. 7 (2016) 12472
  * [15] Y. Nakaguchi, T. Nishioka, “A holographic proof of Rényi entropic inequalities” JHEP 12 (2016) 129
  * [16] K. Kawabata, T. Nishioka, Y. Okuyamaa, K. Watanabec, “Replica wormholes and capacity of entanglement” JHEP 10 (2021) 227
  * [17] J. Boer, J. Järvelä, E. Keski-Vakkuri, “Aspects of capacity of entanglement” Phys. Rev. D 99 (2019) 066012
  * [18] K. Kawabata, T. Nishioka, Y. Okuyama, K. Watanabe, “Probing Hawking radiation through capacity of entanglement” JHEP 05 (2021) 062
  * [19] K. Okuyama, “Capacity of entanglement in random pure state” Phys. Lett. B 820 (2021) 136600
  * [20] B. Bhattacharjee, P. Nandy, T. Pathak, “Eigenstate capacity and Page curve in fermionic Gaussian states” Phys. Rev. B 104 (2021) 214306
  * [21] P. Nandy, “Capacity of entanglement in local operators” JHEP 07 (2021) 019
  * [22] D. Shrimali, S. Bhowmick, V. Pandey, A. Pati, “Capacity of Entanglement for Non-local Hamiltonian” Phys. Rev. A 106 (2022) 042419
  * [23] L. Wei, “Average capacity of quantum entanglement” J. Phys. A 56 (2023) 015302
  * [24] R. Arias, G. Di Giulio, E. Keski-Vakkuri, E. Tonni, “Probing RG flows, symmetry resolution and quench dynamics through the capacity of entanglement” JHEP 03 (2023) 175
  * [25] E. Fel’dman, M. Yurishchev, “Fluctuations of Quantum Entanglement” JETP Letters 90 (2009) 70
  * [26] M. Yurishchev, “Entanglement entropy fluctuations in quantum Ising chains” J. Exp. Theor. Phys. 111 (2010) 525
  * [27] C.-H. Wen, J.-C. Wang, “Fluctuations of Entanglement Entropy in Noninertial Frames” Commun. Theor. Phys. 62 (2014) 327
  * [28] S.-X. Zhang , T.-H. Liu, S. Cao, Y.-T. Liu, Y.-J. Lia “Quantum fluctuation of entanglement for accelerated two-level detectors” Chinese Phys. B 29 (2020) 050402
  * [29] E. Verlinde, K. Zurek, “Modular Fluctuations from Shockwave Geometries” Phys. Rev. D 106 (2022) 106011
  * [30] E. Verlinde, K. Zurek, “Spacetime Fluctuations in AdS/CFT” JHEP 04 (2020) 209
  * [31] M. Nezhadhaghighi, M. Rajabpour, “Entanglement dynamics in short and long-range harmonic oscillators” Phys. Rev. B 90 (2014) 205438
  * [32] A. Coser, E. Tonni, P. Calabrese, “Entanglement negativity after a global quantum quench” J. Stat. Mech. (2014) P12017
  * [33] J. Cotler, M. Hertzberg, M. Mezei, M. Mueller “Entanglement growth after a global quench in free scalar field theory” JHEP 11 (2016) 166
  * [34] S. Ghosh, K. Gupta, S. Srivastava, “Entanglement dynamics following a sudden quench: an exact solution” EPL 120  (2017) 50005
  * [35] S. Santhosh Kumar, S. Shankaranarayanan, “Evidence of quantum phase transition in real-space vacuum entanglement of higher derivative scalar quantum field theories” Scientific Rep. 7 (2017) 15774
  * [36] S. Ghosh, K. Gupta, S. Srivastava, “Exact relaxation dynamics and quantum information scrambling in multiply quenched harmonic chains” Phys. Rev. E 100 (2019) 012215
  * [37] G. Di Giulio, R. Arias, E. Tonni, “Entanglement hamiltonians in 1D free lattice models after a global quantum quench” J. Stat. Mech. (2019) 123103
  * [38] G. Di Giulio, E. Tonni, “Subsystem complexity after a global quantum quench” JHEP 05 (2021) 022
  * [39] S. Chandran, S. Shankaranarayanan, “Dynamical scaling symmetry and asymptotic quantum correlations for time-dependent scalar fields” Phys. Rev. D 107 (2023) 025003
  * [40] B. Estienne, Y. Ikhlef, A. Rotaru, E. Tonni, “Entanglement entropies of an interval for the massless scalar field in the presence of a boundary” arXiv:2308.00614 (2023)
  * [41] V. Ermakov, “Second order differential equations. Conditions of complete integrability” Univ. Izv. Kiev, Series III 9 (1880) 1 (English translation: A. Harin, under redaction by P. Leach, Appl. Anal. Discrete Math. 2 (2008) 123)
  * [42] E. Milne, “The Numerical Determination of Characteristic Numbers” Phys. Rev. 35 (1930) 863
  * [43] E. Pinney, “The nonlinear differential equation $y^{\prime\prime}+p(x)y+cy^{-3}=0$” Proc. Am. Math. Soc. 1 (1950) 681
  * [44] P. Sánchez-Moreno, J. Dehesa, D. Manzano, R. Yánez, “Spreading lengths of Hermite polynomials” J. Comput. Appl. Math. 233 (2010) 2136
  * [45] R. Unanyan, M. Fleischhauer, “Entanglement dynamics in harmonic-oscillator chains” Phys. Rev. A 89 (2014) 062330
  * [46] C. Callan, F. Wilczek, “On Geometric Entropy” Phys. Lett. B 333 (1994) 55
  * [47] D. Katsinis, G. Pastras, N. Tetradis, “Entanglement of Harmonic Systems in Squeezed States” arXiv:2304.04241 (2023)
  * [48] V. Beaud, J. Sieber, S. Warzel,“Bounds on the bipartite entanglement entropy for oscillator systems with or without disorder” J. Phys. A 52 (2019) 235202
  * [49] H. Casini, H. Liu, M. Mezei, “Spread of entanglement and causality” JHEP 07 (2016) 077
  * [50] K. Andrzejewski, “Dynamics of entropy and information of time-dependent quantum systems: exact results” Quantum Inf. Process. 21 (2022) 117
  * [51] D. Katsinis, G. Pastras, “Entanglement in harmonic systems at coherent states” arXiv:2206.05781 (2022)
  * [52] S. Braunstein, S, Das, S. Shankaranarayanan, “Entanglement entropy in all dimensions” JHEP 1307 (2013) 130
  * [53] S. Santhosh Kumar, S. Shankaranarayanan, “Role of spatial higher order derivatives in momentum space entanglement” Phys. Rev. D 95 (2017) 065023
  * [54] K. Boutivas, G. Pastras, N. Tetradis, “Entanglement and expansion” JHEP 2023 (2023) 199
  * [55] S. M. Chandran, K. Rajeev, S. Shankaranarayanan, “Real-space quantum-to-classical transition of time dependent background fluctuations” arXiv:2307.13611 (2023)
  * [56] D. Katsinis, G. Pastras, “An inverse mass expansion for entanglement entropy in free massive scalar field theory” Eur. Phys. J. C 78 (2018) 282