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

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 Renyi 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.


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 the 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.

Preliminaries
The von Neumann (vN) entropy of the density matrix ρ S vN = S = −tr(ρ ln(ρ)), (2.1) is one of the basic information-theoretic measures.However, the vN entropy provides only a first characterization of ρ (e.g.entanglement), in order to get more information some other measures were proposed.One of them is the Rényi entropy defined as for which lim α→1 R α = 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 In order to get more insight into the notion of the entropy let us consider the logarithm of the density matrix Then, the vN entropy is the expectation value of K, S = tr(ρK) = ⟨K⟩.Since K is semipositively defined, K ≥ 0, it is called the modular Hamiltonian (due to the relation ρ = 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 α = 1/T of "temperature" T (we use α instead of thermodynamical β index and the Boltzmann constant is put one) and we define the partition function as well as the corresponding free energy Then, following the standard thermodynamical approach we construct the corresponding entropy The new quantity S α is called the modular entropy (due to its relation to the modular Hamiltonian K) and can be easily compared with the Rényi entropy Equivalently, it is the vN entropy of the modified density operator ρ α = ρ α /Z α , i.e. S α = ⟨K⟩ α where ⟨•⟩ α stands for the expectation value at ρ α .The last observation gives immediately S α ≥ 0 and lim α→1 S α = 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 i.e.C α measures the variance of the modular Hamiltonian at ρ α .Since for α = 1 the modular entropy coincides with the vN one (S 1 = S vN ), the capacity reduces to the variance of the modular Hamiltonian C ≡ C 1 = ⟨K 2 ⟩ − ⟨K⟩ 2 , see [25]- [30].In consequence, it gives the quantum fluctuation with respect to the original state.For the reduced density matrix C measures the entanglement entropy fluctuations.The magnitude of the relative entanglement fluctuations is defined as δK = √ C/S [25,17].As a border between strong and weak fluctuations the condition δK = 1 can be used, i.e.C = S 2 .Similarly, for an arbitrary α we have and we can also easily express δ α K = √ C α /(α⟨K⟩ α ) 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 α , S α (for any α) coincide and C α = δ α 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.

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 ω(t).It is described by the Hamiltonian for which the equation of motion reads ..
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] ..
where c is a constant (we assume c ̸ = 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 where H n for n = 0, 1 . . .denote the Hermite polynomials and τ (t) = 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 and will be assumed in what follows.
Then the density function of the state ψ n (x, t) reads while the density of the Fourier transform of ψ n (x, t) is given by the formula We start with the coordinate Rényi entropy as well as with their momentum counterpart R α,p n (t).For the TDHO they read respectively where e −αy 2 dy is the so-called entropic moment of the Hermite polynomials [44].For lower states the latter quantities can be given explicitly (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 and analogously for the momentum case.From eq. (3.11) we see that the time dependence of the modular entropies S α,n is the same as for the Rényi ones; in particular, the increase of the modular entropy for an arbitrary state ψ n does not depend on n, e.g.
,n (t 0 ) (obviously it vanishes for the ordinary harmonic oscillator, b(t) = 1).In summary, for the states ψ n (t) both kinds of entropies differ by a constant (depending on α and n) only.Now, let us pass to the notion of the modular capacity.By virtue of eq.(2.8) we have thus the capacity is time-independent (and bounded for fixed α and n).For the lowest (initial) state we have for any α while for the excited state the α-dependence emerges where Υ 1 is trigamma function (the first derivative of the digamma function).For the most interesting case α = 1 (corresponding to the fluctuation of the vN entropy) we have For n ≥ 2 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. ..
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.

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.

Two coupled oscillators
Let us consider the Hamiltonian of the form By means of the transformation x = Ry where R is the 2-dimensional rotation matrix with the rotation angle π/4, we transform the Hamiltonian (4.1) into the following one where now p's denote the canonical momenta associated with y's and ω 2 1 (t) = ω 2 (t) + 2k(t), ω 2 2 (t) = ω 2 (t).The frequencies ω 1,2 (t) determine the parameters of the initial Hamiltonian (4.1) as follows 3) The evolution ψ 0 (x, t) of the ground state ψ 0 (x, t 0 ) of the Hamiltonian operator H(t 0 ) as well as the reduced density matrix )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] where and ζ(t) > χ(t) ≥ 0 are given by , (4.6) while the functions b 1,2 (t) satisfy the EMP equation (3.3) with the frequencies ω 1,2 (t) and the constants c 1 = ω 2 (t 0 ) + 2k(t 0 ) and c 2 = ω(t 0 ), respectively.In consequence, the Rényi entropy and the vN entropy of the reduced density ρ red 0 can be easily found, see Ref. [34] where . Now, we are in the position to find the modular entropy and capacity.First, the straightforward computations give i.e. S α is obtained by the replacement ξ(t) → ξ α (t) in the vN entropy, cf.eq.(4.7) (this can be also seen from the fact the spectrum of ρ red 0 consists of the powers of ξ).Moreover, we have the following relation 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 in consequence, we have 0 ≤ C α (t) ≤ 1.In particular, the capacity (fluctuation) of the entanglement is given by the formula From the above we see that the discussed quantities are determined by the solutions of the two EMP equations with the frequencies ω 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 ω(t) = ω in the Hamiltonian (4.1).Then b 1 (t) = b 2 (t) = 1 and ξ(t) = ξ is a constant which can be easily expressed as a function of the ratio k/ω 2 (the same concerns S α , C α ); in particular, ξ → 1 for k → ∞.In this case, the reduced matrix ρ red 0 is equivalent to the thermal density operator for a single harmonic oscillator specified by the frequency ζ 2 − χ 2 and temperature T = − ζ 2 − χ 2 / ln(ξ).The modular entropy is also the vN entropy of the harmonic oscillator provided we define the temperature T = − ζ 2 − χ 2 /(α ln(ξ)).Moreover, the capacity becomes the ordinary heat capacity of a quantum harmonic oscillator.Thus when ξ tends to one (k to infinity) the entanglement fluctuations are close 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 such systems in the more general time-dependent case.Then all entropies are time-dependent and measure the evolution of entanglement (non-local correlations).This is interesting from the thermodynamical point of view since such models can be considered as many-body non-equilibrium systems.Although the unqualified definition of the entropy for non-equilibrium systems is still problematic [45] such information can give some insight into how thermodynamics arises in isolated quantum systems as well as probe thermalization processes, see e.g.[46,47,48] and [49] for experimental analysis.Moreover, such considerations are relevant for quantum information processing and technology and can lead to some universal scaling properties as for the correlations functions in the Kibble-Zurek scenario [50].
A typical example of such a situation are many-body quantum quench phenomena characterized by sudden (fast) changes of global parameters of the Hamiltonian.Then it turns out that the entanglement entropy starts to grow linearly in time.Next it is saturated to a volume-law behavior suggesting in this way the thermalization of the subsystem, see [10,51] and references therein.A deeper analysis of such a scenario (see e.g.[52,53,54]) implies that the resulting reduced density matrix for a subsystem relaxes to the Gibbs ensemble or in the case of integrable systems, a generalized Gibbs ensemble.Moreover, the other approaches to this problem are also related to the entanglement entropy.For example, holographic considerations presented in Ref. [55] suggest that entanglement entropy can be treated also as a kind of coarse-grained entropy for time-dependent system.Similar conclusions, but in a different approach [56], appear for the cosmological spacetimes and non-equilibrium processes.

Many coupled oscillators
The system of N coupled time-dependent oscillators (with the nearest neighbour interaction) is described by the Hamiltonian Imposing boundary conditions the Hamiltonian can be rewritten in the form where x = (x 1 , . . ., x N ) and Λ(t) is a symmetric N × N matrix with the eigenvalues λ j (t), for j = 1, . . ., 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 Λ(t) are of the form λ j (t) = ω 2 (t) + 4k(t) sin 2 ( jπ N ), j = 1, . . ., N .
• Finally, we can impose the Dirichlet boundary conditions (DBC), x 0 = x N +1 = 0, then λ j (t) = ω 2 (t) + 4k(t) cos 2 ( jπ 2(N +1) ), j = 1, . . ., 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 inconsistency 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,57] where x = (x 1 , . . ., x N ) and Ω(t respectively, while b j (t) are the solutions of the EMP equations with the frequencies λ j (t) and, finally, U is a time-independent matrix diagonalizing Λ(t): Next, we split the whole system into two parts: the first one A consisting of n first oscillators (x 1 , . . ., x n ) and the second one B consisting of the remaining N − n ones described by x = (x n+1 , . . ., x N ) and rewrite Ω and B in the form 1 where Ω 1 , B 1 are n × n matrices.Tracing (integrating) over the subsystem A we obtain the reduced density of the subsystem B. Namely, after straightforward computations we arrive at the formula where Z, Υ, ∆ are (N − n) × (N − n) matrices given by 1 For simplicity of notation we omit the time parameter t in the matrices when it does not cause ambiguity.
Such a form of the reduced density coincides with the one obtained in Ref. [57].The matrices Ω 1 , B 1 are real and symmetric thus Θ is skewsymmetric.In consequence, for N = 2 (and n = 1) Θ vanishes and we obtain the results from previous section for two oscillators.However, for higher N (and n) ∆ is a complex (but Hermitian) matrix.In consequence, we cannot directly apply methods from Ref. [5] to obtain the spectrum of ρ B (i.e.simultaneously diagonalize both Υ and ∆ by means of an orthogonal matrix).The term with Θ 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 A and S B were not equal (this is also in contrast to results of Ref. [38]).For similar reasons an approach in Ref. [31] seems also to miss something; namely, in our notation the matrix E therein takes the form Υ −1 ∆ and thus has a real (not complex as in [31]) spectrum (in fact the same as the Hermitian matrix ∆ below, eq.(4.24)); moreover, the matrices E and its complex conjugate are, in general, not simultaneously diagonalizable and thus an extension of the reasoning from Ref. [58] is more involved.The problem of the spectrum of the density operator in the case of Hermitian, but not real, ∆ has been recently discussed, in the context of squeezed state, in Ref. [59].Here we adopt essential results and for more details we refer to [59].
Namely, it turns out that the spectrum of the reduced density matrix with Hermitian ∆ is of the form where ξ's are the inverse of the eigenvalues (larger than one) of the following matrix where while Π is an orthogonal matrix diagonalizing Υ, i.e.ΠΥΠ T = Υ; let us note that ∆ is a Hermitian matrix.Obviously, such an approach involves the invertibility of the matrix ∆ thus we can consider the subsystem B provided N − n ≤ n, in the other case we should take the subsystem A and use the fact that for the pure state the spectrum of ρ B is equivalent (up to irrelevant zeros) to the spectrum of ρ A , see also the discussion in [59]; later on, using a different approach, we will avoid this technical problem.Now, from (4.22) we can immediately find the Rényi entropies where R α [ξ j (t)] has the form as in (4.7).
In consequence, the modular entropy as well as the capacity read where S α [ξ j (t)], C α [ξ j (t)] are of the form (4.8) and (4.10) (i.e. with the replacement ξ(t) → ξ j (t)).In particular, we immediately obtain the formula for entanglement fluctuation C(t) ≡ 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 the entanglement entropy and abrupt quenches, see e.g.[32,33,38].First, we define the time-dependent correlations (covariance) matrix where Using the results of Sec. 3 we obtain, after straightforward computations, that where U is defined as in (4.16) and Q, P , R are diagonal matrices with the following diagonal elements for k = 1, . . ., N .Since the matrix Θ is skew-symmetric the covariance matrix of the reduced density of the system B is given by the restriction of the initial one to the subsystem B where (Q B ) kl = Q kl , (R B ) kl = R kl and (P B ) kl = P kl for k, l = n+1, . . ., N .Next, following the procedure based on the symplectic transformation, see e.g.Refs.[33,38,60], we construct the matrix ΓB = iJΓ B , where J is given by and I is (N − n) × (N − n) identity matrix.Now, the spectrum of the matrix ΓB consists of elements ±γ k , k = n + 1, . . ., N and the Rényi entropies take the form By means of (4.32) we can easily find the modular entropy and capacity; for example, the entanglement fluctuation reads 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 Finally, let us note that 0 ≤ C(t) ≤ N − n.Moreover, the state is pure thus 0 ≤ C(t) ≤ n.
In consequence, for the initial ground state the capacity remains bounded 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 the results for the finite-dimensional spaces and fit into the field theory problems.

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 EMP 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 ω(t), k(t) change, at time t 0 = 0, from constant values (ω(i), k(i))) to another constant values (ω(f ), k(f )).For the abrupt quench the solutions of the EMP equations (4.15) with the initial conditions (3.5) read b j (t) = r j cos(2t λ j (f )) + s j , (4.36) where λ j (i), λ j (f ) are the eigenvalues of Λ before and after quench and r j = (λ j (f ) − λ j (i))/(2λ j (f )), s j = (λ j (f ) + λ j (i))/(2λ 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 λ N (f ) = ω 2 (f ) is distinguished.Indeed, it tends to zero as ω(f ) is small, in contrast to the previous frequency λ N −1 (f ) = ω 2 (f )+4k(f ) sin 2 ((N −1)π/N ) when k(f ) ≫ ω 2 (f ).Under these assumptions r N , s N tend to infinity as ω(f ) is close to zero; thus due to eq. ( 4.36) we expect the distinguished periodicity T = π/ω(f ) (in particular, it does not depend on k(f )).This situation is presented in Fig. 1, where N = 4 and ω(f ) = 0.01 (i.e.T ≃ 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 → ∞ the situation changes: λ N is not distinguished, and the above reasoning breaks down; we return to this issue in the next sections.For the DBC the situation is quite different.Namely, for ω(f ) = 0 the λ's remain nonzero.However, for large N the frequency λ N tends to zero implying that the distinguished periodicity should be of the form 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 T ≃ 100 while for k(f ) = 9 we obtain T ≃ 33; see Fig. 4 and 5 in the Sec.5.2.1.In general, we have T ≃ N √ k for large N , thus for k = 1 (lattice approximation) we have T ≃ N , see the next section.

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 ϕ.Namely, we discretize the field Hamiltonian, in general with time-dependent mass m(t), into a chain of harmonic oscillators by imposing a UV cutoff ϵ and IR cutoff N ϵ.Then we obtain a model described by the Hamiltonian (4.12) with k = 1 and ω 2 (t) = m 2 (t)ϵ 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.

Constant frequency
For start we consider the constant mass, which can be realized as a special case of the abrupt quench with ω(i) = ω(f ) = ω = const.We will focus on more subtle case when ω → 0 which is related 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 ω equals zero (the divergence is in the limit N → ∞ only).For the entanglement capacity the situation is quite different.Namely, as we noted above the entanglement capacity is bounded thus for ω = 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 ϵ the capacity and entropy are equal to each other.Namely, for the PBC (circle of the length L) we have In particular, for the infinite length L → ∞ we obtain S = 1 3 ln(l/ϵ) + 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 ϵ (and l = nϵ).For the DBC instead of eq. ( 5.2) we expect the relation S = 1  6 ln L πϵ sin πl L + O(1) = C.In the former case we take ω = 0.001 and next analyse the limit as ω → 0, in the Dirichlet case we can put directly ω = 0 from very beginning.
The results for the PBC are presented in the left panel of Fig. 2 (for N = 100).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 ω 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 ω = 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).

Quenches
Now, let us pass to the case of oscillatory chains with time-dependent ω(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.Due to zero modes we will consider the DBC and PBC separately; in the latter case, the final frequency will be small but not equal to zero.

Abrupt quench
First, we consider the abrupt quenches; namely, we start with some ω(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] 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 references 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).
Let us compare these results with the ones for the DBC.In this case we can put the final frequency zero ω(f ) = 0 (the entropy is finite).Then, in agreement with the discussion presented in Sec.4.3, we observe the distinguished periodicity T ≃ N/ √ k for large N ; namely, for k = 1 (N = 100) we have T ≃ 100 , see Fig. 5, while for k = 9 it gives T ≃ 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 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.
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 more 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(πn/N ) appears (see green data points in Fig. 6).Moreover, the symmetry n → N −n is preserved for entropy and capacity.In consequence, the entropy and capacity of the subsystems A and B coincide (the state is pure).Finally in order to compare better the dynamics of the entropy and capacity we plot the magnitude of relative quantum entanglement fluctuations δK = √ 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).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.ω(i) ≥ 0.1 for the DBC and PBC (in the latter case ω(f ) = 0.01 is fixed).The results are presented in Fig. 8. First, we see that the slope coefficient behaves linearly for 0.1 ≤ ω(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 ≤ ω(i) ≤ 0.5.More preciously, the dashed lines in Fig. 8 for the PBC are given by 0.575ω(i) + 0.003 and 0.657ω(i) + 0.0004 for the entropy and capacity, respectively; while for the DBC we have 0.287ω(i) + 0.001 and 0.329ω(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 t S = (π − 2)/2 ≃ 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. the left and right panels in Fig. 8. Next, as ω(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 ω(i) the entropy slope is greater than for the capacity.Let us return to the slope of the capacity.Using the considerations from Ref. [33] (basing on the quasiparticles approach [61]) and formula (2.8) we find the initial slope for the capacity and the PBC: we present this situation using a t S and a t C 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 ω(i) = 1 for the entropy and ω(i) = 0.5 for the capacity, see Fig. 9 for ω(i) = 0.5; such a restriction is related to the validity of lattice approximation (ω(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.

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.The first one with 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 (the massless field).For better transparency we describe the DBC and we mention only the necessary changes for the PBC or NBC.
To this end we put k = 1 and take in the Hamiltonian (4.12).Then ω(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 In order to obtain a continuous (differentiable) quench protocol we take Then the initial frequency ω(i) is quenched to zero at infinity.To analyse such a model we need the solutions of the EMP equations (4.15) with λ j (t) = ω2 (t) + a j , (6.4) where, for the DBC, a j = 4k cos 2 ( jπ 2(N +1) ) > 0, j = 1, . . ., N .The corresponding functions bj (t) can be given explicitly.Indeed, by means of the results of Ref. [62], we obtain bj where b j (t), for j = 1, . . ., N , are of the form Note that for the PBC we have a N = 0 and thus we have to modify b N as follows 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 ε = √ 2/3 which corresponds to the initial frequency ω(i) = 3.Let us compare it with the abrupt case.First, we should note that ω(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 ω(i) i.e. large ε).For larger times, we observe the oscillations about a saturation value similarly to the abrupt case.The second interesting case is given by the initial conditions (3.5) imposed at t 0 = −∞, 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/ε 2 at time t = 0.In this case the functions b j (t) are given by b for j = 1, . . ., N .A remarkable property of the functions (6.8) is that they satisfy the condition (3.5) also at t = ∞; there is no oscillatory behavior for a j > 0 at plus infinity (independently of the value of parameter ε).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.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 ε → 0 the frequency tends to the Dirac delta.Then the numerical results for ε = 0.0001 and t ̸ = 0 give the entropy and capacity as in Sec.5.1, while for t = 0 we obtain the similar behavior; 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.

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 the 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 ≪ N and the abrupt quenches they coincide with the observations made in Ref. [24].Obviously, we observed (for n ∼ N/2) some differences due to finite-size assumption (e.g.oscillations, in general non-linear behavior w.r.t. the subsystem size, and local area laws), see Figs. 3-7.We confirmed, in agreement with general theory, that the capacity (as the 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 the 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 issues have been discussed in the analytical manner, due to the explicit forms (given by elementary functions) of the solutions of the EMP equations (numerical computations were used, for higher N , to the standard matrix algebra only).
Turning to possible further developments, it would be interesting to extend the above considerations 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, [59,63].Additionally, some extensions to larger dimensions are also interesting (see [33,64] and references therein) or higher derivative theories [65].In this context, particularly interesting are applications in field theory in expanding spacetimes [66,67].Finally, basing on the methods presented in Refs.[40,68], a deeper analysis of the capacity (fluctuations of the entropy) would be worthwhile.
) 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