Entanglement production in the dynamical Casimir effect at parametric resonance

TL;DR

This paper analytically studies entanglement growth in the dynamical Casimir effect at parametric resonance, revealing dimension-dependent entropy scaling and unique properties of bosonic field theories.

Contribution

It provides the first analytical computation of entanglement entropy evolution in the dynamical Casimir effect at resonance across arbitrary dimensions.

Findings

In (1+1)D, entanglement entropy grows logarithmically over time.

In higher dimensions, entropy growth is linear with time.

The logarithmic growth has a unique 1/2 pre-factor specific to bosonic fields.

Abstract

The particles produced from the vacuum in the dynamical Casimir effect are highly entangled. In order to quantify the correlations generated by the process of vacuum decay induced by moving mirrors, we study the entanglement evolution in the dynamical Casimir effect by computing the time-dependent R\'enyi and von Neumann entanglement entropy analytically in arbitrary dimensions. We consider the system at parametric resonance, where the effect is enhanced. We find that, in (1+1) dimensions, the entropies grow logarithmically for large times, $S_A(\tau)\sim\frac{1}{2}\log(\tau)$, while in higher dimensions (n+1) the growth is linear, $S_A(t)\sim \lambda\,\tau$ where $\lambda$ can be identified with the Lyapunov exponent of a classical instability in the system. In $(1+1)$ dimensions, strong interactions among field modes prevent the parametric resonance from manifesting as a Lyapunov instability, leading to a sublinear entropy growth associated with a constant rate of particle production in the resonant mode. Interestingly, the logarithmic growth comes with a pre-factor with $1/2$ which cannot occur in time-periodic systems with finitely many degrees of freedom and is thus a special property of bosonic field theories.

Equations

$$L(t)=L_1[1+ϵ\sin (2\pi t/T)].$$

$$L=-\frac{1}{2}\int d^{n}x({\partial }^{\mu }\varphi )({\partial }_{\mu }\varphi ),$$

$$\frac{{\partial }^2\varphi }{\partial t^2}-{\nabla }^2\varphi =0.$$

$$\pi (t,x)={\partial }_t\varphi (t,x),$$

$$H=\frac{1}{2}\int d^{n}x\left[{\pi }^2+(\nabla \varphi {)}^2\right].$$

$$\varphi (t,x){|}_{\partial V(t)}=0.$$

$$\varphi (t,x^1=0)=\varphi (t,x^1=L(t))=0.$$

$$L(t)=L_1,\text{ for }t\le 0\text{ and }t\ge T.$$

$$(\varphi ,\psi )=i{\int }_{V(t)}{d}^{n}x\left[{\varphi }^{\ast }({\partial }_t\psi )-({\partial }_t{\varphi }^{\ast })\psi \right],$$

$$[\hat{\varphi }(t,x),\hat{\pi }(t,x^{\prime })]$$

$$[\hat{\varphi }(t,x),\hat{\varphi }(t,x^{\prime })]$$

$$\hat{\varphi }(t,x)=\sum _n[u_n(t,x){\hat{a}}_n+u_n^{\ast }(t,x){\hat{a}}_n^{†}],$$

$$[{\hat{a}}_m,{\hat{a}}_n^{†}]={\delta }_{mn},[{\hat{a}}_m,{\hat{a}}_n]=[{\hat{a}}_m^{†},{\hat{a}}_n^{†}]=0.$$

$$u_n^{(in)}(t,x)={\psi }^n(t,x)\propto e^{-i{\omega }_{n}t},\text{for }t\le 0.$$

$$u_k^{(out)}(t,x)={\psi }_k(t,x)\propto e^{-i{\omega }_{k}t},\text{for }t\ge T.$$

$$\hat{\varphi }(t,x)$$

$$=\sum _n[{\psi }^n(t,x){\hat{b}}_n+{\psi }^{n\ast }(t,x){\hat{b}}_n^{†}],$$

$${\psi }^n=\sum _k({\alpha }_{nk}{\psi }_k+{\beta }_{nk}{\psi }_k^{\ast }),$$

$${\hat{a}}_k$$

$${\hat{a}}_k^{†}$$

$$M=\left(\begin{array}{cc}{\alpha }^T& {\beta }^{†}\\ {\beta }^T& {\alpha }^{†}\end{array}\right).$$

$${\hat{b}}_n|0,in⟩=0,\forall n.$$

$$⟨N_k⟩=⟨0,in|{\hat{a}}_k^{†}{\hat{a}}_k|0,in⟩=\sum _n|{\beta }_{nk}{|}^2.$$

$${\psi }^n=\left[\prod _{i=1}^d\frac{1}{\sqrt{\pi n_i}}\sin \left(\frac{n_i\pi }{L_i}{x}^i\right)\right]e^{-i{\omega }_{n}t},\text{for }t\le 0,$$

$${\omega }_n=\pi \sqrt{\sum _{i=1}^d{\left(\frac{n_i}{L_i}\right)}^2}.$$

$${\psi }^n=$$

$$\times \sqrt{\frac{L_1}{L(t)}}\frac{1}{\sqrt{\pi n_1}}\sin \left(\frac{n_1\pi }{L(t)}{x}^1\right)Q_k^n(t).$$

$$Q_k^n(0)={\delta }_{kn},{\dot{Q}}_k^n(0)=-i{\omega }_n{\delta }_{kn}.$$

$${\psi }_k=\left[\prod _{i=1}^d\frac{1}{\sqrt{\pi k_i}}\sin \left(\frac{k_i\pi }{L_i}{x}^i\right)\right]e^{-i{\omega }_{k}t},\text{for }t\ge T,$$

$$Q_k^n(t)={\alpha }_{nk}{e}^{-i{\omega }_{k}t}+{\beta }_{nk}{e}^{i{\omega }_{k}t},\text{for }t\ge T.$$

$${\psi }^n=\sum _k\sqrt{\frac{L_1}{L(t)}}\frac{1}{\sqrt{\pi n}}\sin \left(\frac{n\pi }{L(t)}x\right)Q_k^n(t).$$

$$\begin{array}{rl}{\ddot{Q}}_k^n+{\Omega }_k^2(t)Q_k^n& =2\lambda (t)\sum _j{g}_{kj}{\dot{Q}}_j^n+\dot{\lambda }(t)\sum _j{g}_{kj}{Q}_j^n\\ & +{\lambda }^2(t)\sum _{j,l}{g}_{jk}{g}_{jl}{Q}_l^n,\end{array}$$

Full text

I Introduction

The dynamical Casimir effect (DCE) describes the creation of particles from the vacuum by moving mirrors Moore (1970); Dodonov (2010). Many of its basic features, as particle production and the dynamics of quantum fields with time-dependent boundary conditions, closely resemble analogous effects in sophisticated phenomena as the Hawking effect, the Unruh effect, or particle production in an expanding Universe. As a result, the DCE has been extensively investigated as a simple model to explore physical aspects of these more complicated systems in a more manageable context Davies and Fulling (1977); Carlitz and Willey (1987); Brevik et al. (2000); Dodonov (2001); Lock and Fuentes (2017). In addition, the DCE has been observed experimentally in a superconducting circuit Wilson et al. (2011), opening an avenue for the investigation of processes of particle production from the vacuum in tabletop experiments111In the experiment Wilson et al. (2011), the idealized time-dependent boundary conditions imposed by the moving mirror are implemented in practice by changing the properties of the superconducting circuit by the application of a magnetic flux..

Because the DCE effects are usually small in laboratory conditions, it is interesting to analyze it in a condition of resonance Dodonov and Klimov (1992, 1996); Lambrecht et al. (1996); Plunien et al. (2000). In this work, we focus on the case where a mode of the field is at parametric resonance with the mirror oscillation, following Dodonov and Klimov (1996). The general setup of the DCE at parametric resonance is depicted in Fig. 1. A one dimensional cavity is bounded by two mirrors, one of which is static, while the other is allowed to move. If the moving mirror is set to oscillate harmonically with twice the frequency of some mode of the field, such a mode becomes resonant. This leads to a strong enhancement of the effect, with particles being continuously produced from the vacuum, providing a natural scenario for eventual further experimental explorations of the DCE Dodonov (2010).

Several works deal with the most immediate aspects of the phenomenon, such as the time evolution of particle production and the Casimir force generated between the plates (for a review, see section II of Dodonov (2001)). However, a highly unexplored direction is the dynamics of the quantum correlations that emerge due to particle creation and the mixing of field modes induced by the moving mirror. In this work, we contribute to fill this gap by analitically determining the time evolution of entanglement measures for the DCE at parametric resonance.

We compute the Rényi and the entanglement entropy production from an initial vacuum state for the DCE at parametric resonance in arbitrary dimensions. We consider a massless scalar field in a hyperrectangular cavity where the field is required to vanish on the boundaries. For $t\leq 0$, the cavity is kept fixed. For $t>0$, we allow the length of the cavity to vary in one direction for a finite amount of time by keeping one mirror fixed and setting the other to oscillate according to

[TABLE]

After a finite number of oscillations, we then let the mirror return to its starting position, and the cavity is kept fixed again. The frequency of the mirror oscillation is set to twice that of some cavity mode, giving rise to parametric resonance. We select the lowest energy mode to be the resonant, since this is the most likely scenario to be reproduced in a laboratory. The time evolution of particle production in such conditions was fully analyzed in Dodonov and Klimov (1996), where, under reasonable approximations, the Bogoliubov transformation between the initial and final states was computed analytically for all times in both $(1+1)$ and $(3+1)$ dimensions.

In order to compute the entropies, we apply recently developed symplectic techniques for the description of the dynamics of entanglement for Gaussian states in linear bosonic systems Bianchi et al. (2015, 2018). In this framework, the Rényi entropy is first computed as a volume associated with the subsystem in phase space. Then, according to a theorem derived in Bianchi et al. (2018), the asymptotic behavior of the Rényi and entanglement entropies will converge. Moreover, in the presence of Lyapunov instabilities, the asymptotic production of entanglement entropy is completely determined by the Lyapunov exponents. We show that in $(1+1)$ dimensions the entanglement entropy grows logarithmically for large times. In higher dimensions, the resonance is associated with a Lyapunov instability, leading to an asymptotic exponential particle production and linear growth of the entanglement entropy, with a production rate equal to the Lyapunov exponent of the system.

This paper is structured as follows. In Section II, we review the Hamiltonian formalism of the dynamical Casimir effect and the analytical solution for the case of parametric resonance obtained in Dodonov and Klimov (1996). In Section III, we discuss the necessary background on symplectic techniques for the dynamics of entanglement of Gaussian states. In Section IV, we compute analytically the time evolution of the Rényi and entanglement entropies in the dynamical Casimir effect at parametric resonance. Finally, we summarize and discuss our results in Section V.

II Hamiltonian formalism for the dynamical Casimir effect

II.1 The model

Let $\phi(t,x)$ be a neutral massless scalar field in $(d+1)$-dimensions described by the Lagrangian

[TABLE]

where we adopt the mostly plus convention for the Minkowski metric, $\eta_{\mu\nu}=\operatorname{diag}(-,+,\dots,+)$. Requiring the action to be stationary, one finds that the field satisfies the Klein-Gordon equation

[TABLE]

The associated momentum is

[TABLE]

and the Hamiltonian is

[TABLE]

We are interested in the situation where the field is confined to a finite spatial volume $V(t)$. More specifically, we consider a field restricted to a time-dependent parallelepiped bounded in some direction $x^{1}$ to a finite interval $x^{1}\in[0,L(t)]$, where the first boundary is fixed and the second boundary is allowed to move arbitrarily. In the remaining directions, the boundaries are fixed, so that $x^{i}\in[0,L_{i}]$, $i=2,\dots,d$, for $x\in V(t)$. Dirichlet boundary conditions are imposed at $\partial V(t)$:

[TABLE]

In particular,

[TABLE]

If the field $\phi$ represents a component of the electromagnetic field, then these boundary conditions correspond to perfect mirrors, and the field is confined to a cavity bounded by perfectly reflecting plates. The boundary is assumed to move only for a finite period of time $T$, returning at the end to the starting position $L_{1}$:

[TABLE]

In the quantum theory, the motion of the mirrors will lead to particle production, and under this restriction the number of produced particles remains finite.

Let $\mathcal{S}$ be the space of solutions of the Klein-Gordon equation (3) with boundary conditions (6). This space is equipped with an invariant bilinear form

[TABLE]

which is not positive definite. The space of solutions can be decomposed into a direct sum of subspaces of positive and negative energy, $\mathcal{S}=\mathcal{S}^{+}\oplus\mathcal{S}^{-}$. The decomposition is required to be such that (i) positive-energy solutions have positive squared norm $(\phi,\phi)$ and (ii) the complex conjugate of a positive-energy solution is a negative-energy solution, $\mathcal{S}^{-}=\overline{\mathcal{S}^{+}}$. The decomposition is not unique, and distinct decompositions will correspond to distinct choices of vacuum in the quantum theory. The restriction of the invariant bilinear form to $\mathcal{S}^{+}$ defines an inner product in this subspace, turning it into a Hilbert space $\mathcal{H}$. We denote by $\{u_{n}(t,x)\}$ an orthonormal basis of $\mathcal{H}$. The solutions $u_{n}$ are positive-energy normal modes for the Klein-Gordon equation.

The canonical quantization of the scalar field $\phi(t,x)$ and its associated momentum $\pi(t,x)$ provides field operators satisfying the canonical comutation relations at equal times:

[TABLE]

and the equation of motion (3). For each choice of normal modes $u_{n}$, the field operator can be expanded as

[TABLE]

where the $\hat{a}_{n},\hat{a}_{n}^{\dagger}$ are creation and annihilation operators satisfying the usual commutation relations:

[TABLE]

The quantized canonical momentum operator is obtained by taking the time derivative of the field operator in the representation (12). Distinct choices of normal modes correspond to distinct representations of the quantum field. The vacuum state $\ket{0}$ is the state annihilated by all $\hat{a}_{n}$’s and depends on the choice of $\mathcal{H}=\mathcal{S}^{+}$ (although not on the specific basis $\{u_{n}\}$ chosen for $\mathcal{H}$).

We introduce two special representations, called the in- and out-representations. For $t<0$, the Hamiltonian is time-independent, and one can find an orthonormal basis of positive energy $in$-solutions for which the time-dependence is factored out in this region:

[TABLE]

Similarly, one can find an orthonormal basis of positive energy $out$-solutions satisfying

[TABLE]

Both sets $\mathcal{B}_{in}=\{\psi^{n},\psi^{n*}\}$ and $\mathcal{B}_{out}=\{\psi_{k},\psi_{k}^{*}\}$ are linear bases of $\mathcal{S}$, but define distinct decompositions into positive and negative energy subspaces. The field operator can be expanded in either set of normal modes:

[TABLE]

where the $\hat{b}_{n},\hat{b}_{n}^{\dagger}$ are annihilation and creation operators in the $in$-representation, and the $\hat{a}_{k},\hat{a}_{k}^{\dagger}$ are annihilation and creation operators in the $out$-representation. Expanding the $in$-modes in the basis of $out$-modes,

[TABLE]

and substituting into Eq. (17), we find

[TABLE]

This is a Bogoliubov transformation describing the relation between the two representations. The constants $\alpha_{nk},\beta_{nk}$ are the Bogoliubov coefficients of the transformation. The transformation is encoded in a symplectic matrix:

[TABLE]

The evolution of any property of the quantum field due to the motion of the boundaries for $t\in[0,T]$ is completely characterized by the Bogoliubov coefficients. In particular, the process of particle creation from the initial vacuum has a simple description. Let the initial state for $t\leq 0$ be the $in$-vacuum annihilated by all $\hat{b}_{n}$:

[TABLE]

Since we are in the Heisenberg representation, the state of the field remains unchanged. However, the number of particles $N_{k}$ in a mode $k$ observed after the motion of the boundary is described by the $out$-operators:

[TABLE]

We see that the vacuum is unstable if the coefficients $\beta_{nk}$ do not vanish. In the next section we will discuss formulas for the Rényi and entanglement entropies in the $out$ region in terms of the coefficients of the Bogoliubov transformation. Before that, we review the calculation of the Bogoliubov coefficients for the dynamical Casimir effect at parametric resonance performed in Dodonov and Klimov (1996).

II.2 Dynamical Casimir effect at resonance

Let $n$ be a multi-index $n=(n_{1},\dots,n_{d})$. The $in$-modes $\psi^{n}$ are solutions of the Klein-Gordon equation (3) that reduce to standing waves in the $in$-region:

[TABLE]

where

[TABLE]

While the mirror is moving, the modes evolve in a nontrivial way. At each time $t$, we can expand them in a basis of instantaneous Fourier modes that satisfy the imposed boundary conditions,

[TABLE]

Substituting this expansion into the Klein-Gordon equation (3), one finds equations of motion for the Fourier amplitudes $Q^{n}_{k}(t)$, which we will discuss in a moment. The Fourier amplitudes satisfy the initial conditions

[TABLE]

In the $out$-region, the evolution is again trivial. The positive energy $out$-modes are defined as

[TABLE]

The $in$-mode $\psi^{n}$ evolves into a superposition of $out$-modes of positive and negative energy, as written in Eq. (18), with constant Bogoliubov coefficients $\alpha,\beta$ determined by:

[TABLE]

Up to this point, we considered a generic motion $L(t)$ of the boundary. Now let the moving boundary oscillate harmonically around its initial position. When the frequency of the boundary oscillation is twice that of some normal mode, the system is said to be at parametric resonance. Explicit formulas for the Bogoliubov coefficients can then be obtained under an approximation of slow variation of the coefficients, as first shown in Dodonov and Klimov (1996), and we review the relevant results in the following subsections.

II.2.1 (1+1) dimensions

In one spatial dimension, all multi-indices reduce to integer numbers, and the expansion (II.2) of the $in$-modes in the instantaneous basis becomes

[TABLE]

The evolution equation for the Fourier amplitudes reads

[TABLE]

where

[TABLE]

and the coefficients $g$ form an antisymmetric matrix with components

[TABLE]

Consider the case where the first normal mode is resonant, that is, let the moving boundary oscillate harmonically around its initial position

[TABLE]

with twice the frequency of the first normal mode,

[TABLE]

The functions $\Omega_{k}(t)$ and $\lambda(t)$ can then be determined from Eqs. (34) and (32). The oscillation amplitude is assumed to be small compared to the initial width $L_{1}$ of the boundary, $\epsilon\ll 1$. To first order in $\epsilon$, the parameters $\lambda$ and $\dot{\lambda}$ in the evolution equation oscillate harmonically with frequency $2\omega_{1}$, and $\lambda^{2}\sim\epsilon^{2}$ can be neglected.

One can look for solutions of the evolution equation (31) of the form

[TABLE]

where the Bogoliubov coefficients are now allowed to be time-dependent, reaching their asymptotic values $\alpha_{nk},\beta_{nk}$ at $t=T$. The time-dependent coefficients are assumed to vary slowly in time: after substituting (36) into Eq. (31), terms proportional to $\ddot{\alpha}$ and $\ddot{\beta}$ are neglected, and the Bogoliubov coefficients are considered approximately constant during one period of oscillation of any normal mode. Multiplying the resulting equation of motion by $e^{i\omega_{k}t}$ or $e^{-i\omega_{k}t}$ and averaging over intervals $T_{k}=2\pi/\omega_{k}$, one obtains a set of coupled differential equations for $\alpha$ and $\beta$, respectively. These are solved in Dodonov and Klimov (1996). Let us gather the relevant results for our purposes.

The complete elliptic integrals of the first and second kind are defined, respectively, as

[TABLE]

It is convenient to introduce the new time variable

[TABLE]

and define the quantities

[TABLE]

The lowest Bogoliubov coefficients are then given by

[TABLE]

Coefficients with one larger odd index are obtained from the recurrence relations

[TABLE]

where the dots represent derivatives with respect to $\tau$, and the relations

[TABLE]

All coefficients with some even index vanish. In order to compute the evolved state of the subsystem formed only by the resonant mode, it is enough to know the Bogoliubov coefficients with some index equal to $1$, which are all determined by the relations above.

In our analysis of entanglement production, the following formulas for infinite sums of products of Bogoliubov coefficients will play an important role:

[TABLE]

II.2.2 (d+1) dimensions

In higher dimensions, the expansion of the $in$-modes in the instantaneous basis has the general form (II.2). Upon substitution of this expansion into the Klein-Gordon equation, one obtains equations of motion for the Fourier coefficients $Q^{n}_{k}$ similar to (31). The difference is that $j,k,n$ are now multi-indices, the frequencies of the normal modes are

[TABLE]

and the coefficients $g_{jk}$ have a more complicated form, which is not relevant for our purposes, except for the fact that they are still constant.

A mode $r$ is set at resonance by letting the boundary oscillate with twice the time-independent frequency $\omega_{r}$ of the mode. One may take, for instance,

[TABLE]

with

[TABLE]

as done in Dodonov and Klimov (1996). Note that, in contrast to the one-dimensional case, the frequencies are not equidistant for $d>1$. To first order in $\epsilon$, the parameters $\lambda$ and $\dot{\lambda}$ in the evolution equation oscillate harmonically with frequency $2\omega_{r}$, and $\lambda^{2}$ can be neglected, as before.

One can look again for solutions of the form (36) under a slow-variation approximation. The resulting equations for the $Q$’s are quite different, however, as compared to the one-dimensional case. The reason for that is the following. When multiplied by $e^{\pm i\omega_{k}t}$, to first order in $\epsilon$, the right-hand-side of Eq. (31) becomes a sum of terms proportional to $\exp[i(\pm\omega_{j}\pm\omega_{k}\pm 2\omega_{r})]$. The integral over a period $T_{k}=2\pi/\omega_{k}$ of each such term vanishes unless the sum of frequencies in the exponential vanishes. But since they are not equally spaced, this never happens, and the equations for the distinct modes all decouple. One is left with an infinite set of equations for independent oscillators with time-dependent frequencies

[TABLE]

with

[TABLE]

The calculation of the Bogoliubov coefficients can be done independently for each mode. Put $Q_{k}=Q^{k}_{k}$. Fourier amplitudes $G^{n}_{k}$ with $k\neq n$ vanish, since the modes are decoupled. We need to solve

[TABLE]

the equation of a time-dependent oscillator with a harmonically oscillating frequency. One can look again for solutions of the form (36), which we now write as

[TABLE]

with $\alpha_{k}=\alpha_{kk}$ and $\beta_{k}=\beta_{kk}$. Under a slow variation approximation, terms proportional to $\ddot{\alpha}_{k}$ and $\ddot{\beta}_{k}$ are neglected. After averaging over the fast oscillations, one finds for the resonant mode

[TABLE]

while the evolution of the non-resonant modes is trivial in this approximation: $\alpha_{k},\beta_{k}=\textrm{constant}$, for $k\neq r$. Initial conditions corresponding to the $in$-vacuum are given by

[TABLE]

Integrating Eq. (53) with these initial conditions, we obtain the desired Bogoliubov coefficients

[TABLE]

The time-dependent symplectic transformation associated with such Bogoliubov coefficients, which describes the evolution of the resonant mode, is given by

[TABLE]

and has the simple form

[TABLE]

with a time-independent symplectic generator $K_{r}$.

III Entanglement dynamics of Gaussian states

In this section, we review the basic properties of Gaussian states. We focus on systems with a finite number of bosonic degrees of freedom, which we can later take to infinity to recover a bosonic field theory. In particular, we review compact expressions for the entropy and the von Neumann entanglement entropy in terms of the covariance matrix of the Gaussian state. Our conventions closely follow Bianchi et al. (2015, 2018); Hackl (2018), while other reviews include Eisert and Plenio (2003); Weedbrook et al. (2012).

III.1 Gaussian states

Before considering the field theory case, we focus on a system with $N$ degrees of freedom. Here, the classical phase space $V\simeq\mathbb{R}^{2N}$ is equipped with a symplectic form $\Omega^{ab}$, i.e. an antisymmetric and non-degenerate bilinear form. For the quantum theory, we choose the basis

[TABLE]

of creation and annihilation operators. This basis also fixes the form of the symplectic form to get the correct commutation relations

[TABLE]

The basis is not Hermitian, which means that we can define the $2N\times 2N$ matrix $C$, such that

[TABLE]

Given a linear observable $\mathcal{O}=f_{a}\hat{\xi}^{a}$, we need to require $f_{a}^{*}C^{a}{}_{b}=f_{b}$ for $\mathcal{O}$ to be a Hermitian operator.

We can now introduce a special class of states $\ket{\psi}$ that are fully characterized by their displacement vector $z^{a}$ and their covariance matrix $G^{ab}$, defined as

[TABLE]

We refer to a state $\ket{\psi}$ as Gaussian state and label it as $\ket{G,z}$ if its linear complex structure $J$ satisfies the condition

[TABLE]

where we introduced the inverse symplectic form $\Omega^{-1}_{ab}$ with $\Omega^{ac}\Omega^{-1}_{cb}=\delta^{a}{}_{b}$. The linear map $J:V\to V$ is called linear complex structure because it represents the imaginary unit on the classical phase space, i.e. it squares to minus identity. We can use directly $J$ to encode the state $\ket{G,z}$ as solution to the equation

[TABLE]

While complex structures have been used to describe vacua of quantum fields in curved spacetime Ashtekar and Magnon (1975); Wald (1994) for many years, only recently they became a useful tool in quantum information to parametrize both bosonic Bianchi et al. (2015, 2018); Hackl (2018); Hackl and Jonsson (2019) and fermionic Gaussian states Vidmar et al. (2017, 2018); Hackl and Myers (2018); Hackl et al. (2019) in a unified manner. As an example, we can choose $z^{a}=0$ and the covariance matrix to be given by

[TABLE]

which gives rise to the operator-valued vector

[TABLE]

that annihilates the state $\ket{G_{0},0}$, i.e. we just confirmed that the Gaussian state $\ket{G_{0},0}$ is nothing else than the vacuum associated to the annihilation operators $\hat{b}_{i}$. For other choices of $G^{ab}$ and $z^{a}$, we will find that (64) gives rise to a Bogoliubov transformation of the form

[TABLE]

where $\hat{a}_{i}$ represent a new set of annihilation operators annihilating the state $\ket{G,z}$ under consideration.

A defining property of bosonic Gaussian states is the well-known Wick’s theorem, i.e. that higher order $n$-point correlation functions can be efficiently computed from $1$- and $2$-point correlation functions. To state Wick’s theorem, it is useful to define the connected $n$-point correlation function of the state $\ket{G,z}$ as

[TABLE]

We can then use the connected $2$-point function given by

[TABLE]

to state Wick’s theorem as the equality

[TABLE]

i.e. the connected $n$-point correlation function for odd $n$ vanishes, while for even $n$, it is given by a sum over products of $C_{2}^{ab}$.

For the rest of this paper, we will focus on Gaussian states $\ket{G,z}$ with $z^{a}=0$, i.e. they are centered at the origin of phase space. We will write $\ket{G}=\ket{G,0}$ to simplify notation.

III.2 Quadratic Hamiltonians

The most general class of Hamiltonians that preserve the family of Gaussian states, i.e. for which $e^{t\hat{H}}\ket{G,z}$ is again Gaussian, takes the form

[TABLE]

If we only look at Gaussian states $\ket{G,z}$ with $z^{a}=0$, we need to choose $f_{a}=0$ in our Hamiltonian to ensure that $z^{a}$ remains equal to zero under time evolution. We will therefore restrict ourselves to the class of quadratic Hamiltonians

[TABLE]

where we can derive a condition on $h_{ab}$ to ensure that $\hat{H}$ is Hermitian. For $\hat{\xi}^{\dagger a}=C^{a}{}_{b}\hat{\xi}^{b}$, we find the condition

[TABLE]

Furthermore, we can require $h_{ab}$ to be symmetric, i.e. $h_{ab}=h_{ba}$, because the anti-symmetric part will only contribute a constant to the energy. Defining $U(t)=e^{-\mathrm{i}t\hat{H}}$, we can use the well-known Baker-Campbell-Hausdorff identity to compute

[TABLE]

where the symplectic transformation is given by

[TABLE]

For the dynamical Casimir effects, we will need to consider quadratic Hamiltonians with explicit time dependence

[TABLE]

which leads to a time evolution determined by the time-ordered exponentials

[TABLE]

In practice, we will not compute $M(t)$ as a time-ordered exponential, but rather by solving the underlying equations of motion explicitly. Using (75), which also continues to be valid for Hamiltonian with explicit time dependence, we can compute the time evolution of the covariance matrix to be given by

[TABLE]

For our class of Gaussian states, we therefore find

[TABLE]

A special class of time-dependent Hamiltonians are Floquet Hamiltonians, i.e. time-periodic Hamiltonians with period $T$ satisfying

[TABLE]

Such Hamiltonians can be described stroboscopically by only looking at the evolved states at $t=n\,\tau$ with $n\in\mathbb{Z}$. For this, it is sufficient to compute $U(\tau)$ which allows us to define a time-independent effective Hamiltonian $\hat{H}_{\text{eff}}=\frac{1}{\tau}\log U(\tau)$. For quadratic Hamiltonians $\hat{H}(t)$, we find $\hat{H}_{\text{eff}}=\frac{1}{2}h^{\text{eff}}_{ab}\hat{\xi}^{a}\hat{\xi}^{b}$ to be also quadratic with

[TABLE]

where $M(\tau)$ can be evaluated as the time-ordered exponential (79) or by integrating the equations of motion.

III.3 Entanglement measures

Given a pure Gaussian state $\ket{G}$ and a choice of subsystem $A\subset V$ (inducing a tensor product decomposition $\mathcal{H}=\mathcal{H}_{A}\otimes\mathcal{H}_{B}$), we can quantify the amount of quantum correlations between $A$ and its complement $B$ using different entanglement measures. We are particularly interested in the von Neumann entropy $S_{A}(\ket{G})$ and the Rényi entropy $R_{A}(\ket{G})$ of the reduced state

[TABLE]

where we trace out the degrees of freedom in $B$.

One can show that the mixed state $\rho_{A}$ is also Gaussian, i.e. it continues to satisfy Wick’s theorem and can be efficiently parametrized by the restriction of $G^{ab}$ to the subsystem $A$. Mathematically speaking, we have a decomposition of the classical phase space $V$ as a direct sum $V=A\oplus B$, as well as of its dual $V^{*}=A^{*}\oplus B^{*}$, such that the bilinear form $G:V^{*}\times V^{*}\to\mathbb{R}$ can be meaningfully restricted to

[TABLE]

In practice, we use the matrix representation of $G$ with respect to a basis $(\hat{\xi}^{a}_{A},\hat{\xi}^{a}_{B})$

[TABLE]

with $\hat{\xi}^{a}=(\hat{\xi}^{a}_{A},\hat{\xi}^{a}_{B})$ and $N=N_{A}+N_{B}$. Note that this basis is complex, which implies that also $G$ correspond to the analytically continued bilinear form on the complexified phase space $V_{\mathbb{C}}$ and its dual $V^{*}_{\mathbb{C}}$.

The key result is that $\rho_{A}$ is fully characterized by $G_{A}$, which implies that any function of $\rho_{A}$ can be written purely in terms of $G_{A}$. The von Neumann entropy $S_{A}(\ket{G})$ and the Rényi entropy $R^{(n)}_{A}(\ket{G})$ of order $n$ are

[TABLE]

and we define $R_{A}(\ket{G}):=R_{A}^{(2)}(\ket{G})$ as the Rényi entropy, dropping the reference to order 2 in this case. To express these functions in terms of $G_{A}$, it is useful to define the restricted linear complex structure

[TABLE]

where $\Omega_{A}$ is the restriction of $\Omega$ to the subsystem $A$. An important difference between the full linear complex structure $J$ (describing a pure Gaussian state) and the restricted linear complex structure $J_{A}$ (describing in general a mixed Gaussian state) lies in the fact that $-J_{A}^{2}$ is not necessarily the identity—instead, we have the inequality

[TABLE]

i.e. $J_{A}^{2}$ has negative eigenvalues whose magnitude is at least $1$, but possibly larger. The magnitude of these eigenvalues or, put differently, the failure of $J_{A}$ to square to minus identity provides a measure of entanglement between $A$ and $B$ in the state $\ket{G}$. The formulas for $S_{A}(\ket{G})$ and $R_{A}(\ket{G})$ are then given by Sorkin (1983); Bianchi et al. (2015); Hackl et al. (2018)

[TABLE]

where the absolute value $|\cdot|$ in the second equation is meant in terms of eigenvalues. Note that in the basis from (58) leads to the standard form of $\Omega$ (and $\Omega_{A}$) from (59), which leads to $\det\Omega_{A}=(-1)^{N_{A}}$.

III.4 Theorems on entanglement asymptotics

For classically unstable quadratic systems and Gaussian initial states, there exist several theorems that describe the asymptotic behavior of the entanglement entropy and the Rényi entropy for large times. For this, let us introduce the concept of Lyapunov exponents and unstable quadratic systems.

Definition 1** (Lyapunov exponents).**

Given a quadratic Hamiltonian system $H(t)=\frac{1}{2}h(t)_{ab}\xi^{a}\xi^{b}$ with linear Hamiltonian flow

[TABLE]

(where $K(t)^{a}{}_{b}=\Omega^{ac}h(t)_{cb}$) on the classical phase space $V$, we define the Lyapunov exponent $\lambda_{\ell}$ of a linear observable $\ell\in V^{*}$ as the limit

[TABLE]

where the norm $\lVert\cdot\rVert$ can be induced by any positive definite inner product on the dual phase space. Note that the action on the dual phase space is given by the transpose $M^{\intercal}(t)_{a}{}^{b}=M(t)^{b}{}_{a}$.

Definition 2** (Unstable quadratic system).**

A classical quadratic Hamiltonian $H(t)=\frac{1}{2}h(t)_{ab}\xi^{a}\xi^{b}$ is called unstable if there exists at least one linear observable $\ell\in V^{*}$ with positive Lyapunov exponent $\lambda_{\ell}>0$. There are at most $2N$ distinct Lyapunov exponents and we can always find a symplectic basis $\mathcal{D}_{L}=(\ell_{1},\cdots,\ell_{2N})$ with associated Lyapunov exponents $\lambda_{i}:=\lambda_{\ell_{i}}$ satisfying

[TABLE]

and $\lambda_{2N+1-i}=-\lambda_{i}$, i.e. all Lyapunov exponents come in conjugate pairs.

Proof.

The proof can be found in definition 2 of appendix A.2 in Bianchi et al. (2018). ∎

With these definitions at hand, we can now review the key theorems that we will use in subsequent sections to study dynamics of Rényi entropy (of order 2) and von Neumann entanglement entropy.

Theorem 1** (Rényi entropy as phase space volume).**

The Rényi entropy (of order $2$) $R_{A}(\ket{G})=R_{A}^{(2)}(\ket{G})$ is given by

[TABLE]

where $\mathcal{V}_{A}\subset A^{*}\subset V^{*}$ represents an arbitrary parallelepiped of symplectic volume $1$.

Proof.

The proof can be found in Sec. 6.2 of Bianchi et al. (2018). ∎

Theorem 2** (Asymptotic entanglement production).**

We consider an unstable and regular222The word “regular” is used to exclude some pathological cases discussed in Appendix A.3 of Bianchi et al. (2018). Hamiltonian system with $N$ degrees of freedom and quadratic Hamiltonian $\hat{H}(t)=\frac{1}{2}h(t)_{ab}\hat{\xi}^{a}\hat{\xi}^{b}$, whose Lyapunov exponents are

[TABLE]

An initial Gaussian state $\ket{G_{0}}$ will evolve into $\ket{G_{t}}$, where the time-dependent covariance matrix is given by

[TABLE]

A system decomposition $V=A\oplus B$ of the classical phasespace $V$ (with $\dim A=2N_{A}$ and $\dim B=2N_{B}$) induces a tensor product decomposition $\mathcal{H}=\mathcal{H}_{A}\otimes\mathcal{H}_{B}$ that we can use to compute entanglement entropy $S_{A}(\ket{G_{t}})$. The asymptotics of the Rényi entropy $R_{A}(\ket{G_{t}})$ and the entanglement entropy $S_{A}(\ket{G_{t}})$ is given by

[TABLE]

where we assume the subsystem to be generic, i.e. above statement applies to all subsystems except a set of measure zero.

Proof.

A detailed proof can be found in Sec. 2.1 and 2.2 of Bianchi et al. (2018). ∎

Next, we will compare the rate of entanglement entropy production in the dynamical Casimir effect with these predictions.

IV Entanglement production at parametric resonance

In this section, we compute the time-dependent Rényi and entanglement entropies at resonance. In $(1+1)$-D, the normal modes are strongly coupled by the evolution and become highly entangled. We compute the entropies explicitly for the subsystem associated with the resonant mode. In higher dimensions, the modes decouple under the assumed approximations, and the normal modes are not entangled. In this case, we consider a generic subsystem that intersects nontrivially the resonant mode.

IV.1 (1+1) dimensions

Let the scalar field be initially in its vacuum state $\ket{0,in}$. We work in the Heisenberg representation, so that the state of the system is constant and the time evolution is encoded in the observables. Annihilation and creation operators $\hat{a}_{k}(t),\hat{a}_{k}^{\dagger}(t)$ for a mode $k$ at each time $t$ are defined in terms of the time-dependent Bogoliubov coefficients as in Eq. (19) :

[TABLE]

Such operators are associated with the instantaneous basis at $t$. It is convenient to use the variable $\tau=\epsilon\omega t/2$ to compute the evolution of the entropies.

The covariance matrix at time $\tau$ is given by

[TABLE]

where $\hat{\xi}^{a}=(\hat{a}_{k}(\tau),\hat{a}_{k}^{\dagger}(\tau))$ in the complex representation. The restriction of the covariance matrix to a subsystem $A$ defines the reduced covariance matrix $G_{A}$ associated with it. The Rényi entropy $R_{A}(\ket{G})$ can be computed from $G_{A}$ as

[TABLE]

where we used Eq. (92) and the fact that $\det(\Omega_{A})=-1$ for a single degree of freedom ($N_{A}=1$). Similarly, the entanglement entropy is given by

[TABLE]

With $J_{A}=-G_{A}\,\omega^{-1}_{A}$ as explained in Eq. (93). We will now compute these quantities.

Expressing $\hat{\xi}^{a}$ in terms of creation and annihilation $in$-operators $\hat{b}_{m}^{\dagger}$ and $\hat{b}_{m}$ using Eq. (101), and representing the covariance matrix as

[TABLE]

we find for the matrix components

[TABLE]

From now on, we will omit the arguments of the time-dependent Bogoliubov coefficients for conciseness.

Since $\alpha_{nm}$ and $\beta_{nm}$ are real for any $m$ and $n$, we have

[TABLE]

Differentiating with respect to $\tau$, these sums can be evaluated using Eqs. (43)–(45) and then integrated analytically.

Let us start with the term $G_{A}^{11}$. Taking its derivative with respect to $\tau$ yields

[TABLE]

Then, according to Eq. (44), we have

[TABLE]

Similarly, taking the derivative of the term $G^{12}_{A}$, we find:

[TABLE]

From Eq. (43), we conclude that

[TABLE]

Since $G^{11}_{A}=G^{22}_{A}$ and $G^{12}_{A}=G^{21}_{A}$, the only remaining step is to integrate Eqs. (110) and (112).

The explicit form of the Bogoliubov coefficients $\alpha_{11}$ and $\beta_{11}$ is presented in Eq. (41). Inserting these formulas in Eq. (110) we find

[TABLE]

where the complete elliptic integrals of the first and second kind, $K(\kappa)$ and $E(\kappa)$, are defined in Eqs. (38) and (37), and $\kappa$ and $\tilde{\kappa}$ are given by Eq. (40).

From the definitions of $\kappa$ and $\tilde{\kappa}$, we can easily check that

[TABLE]

Substituting this relation in (113), we can integrate it analytically (see Dodonov and Klimov (1996)), obtaining:

[TABLE]

where $C_{1}$ is a constant of integration to be determined from the initial conditions.

In a similar way, substituting (41) in Eq. (112), we find

[TABLE]

This equation can again be integrated Dodonov and Klimov (1996), leading to

[TABLE]

where $C_{2}$ is a constant to be determined from the initial conditions.

To fix the integration constants, we need to evaluate the expressions in Eqs. (115) and (115) at $\tau=0$ and compare them with the components of the initial covariance matrix $G_{A}(0)$. But at $\tau=0$, the Bogoliubov transformation is trivial, $\alpha_{mn}(0)=\delta_{mn},\beta_{mn}(0)=0$. Hence, the initial covariance matrix is given by

[TABLE]

Now, in the small-time limit $\tau\ll 1$, the elliptic integrals can be expanded as Gradshteyn and Ryzhik (2014)

[TABLE]

In addition, in this limit we see that $\kappa\to 0$ and $\tilde{\kappa}\to 1$. Keeping only terms up to first order in $\tau$, we find the matrix elements

[TABLE]

Hence, the integration constants are $C_{1}=C_{2}=0$, and the matrix $G_{A}(\tau)$ assumes the form

[TABLE]

Furthermore, we can analyze the long-time behavior of the Rényi entropy. According to Gradshteyn and Ryzhik (2014), when $\tau\gg 1$, the leading terms of the asymptotic expansion of the elliptic integrals are

[TABLE]

Therefore, the covariance matrix becomes

[TABLE]

where $\zeta\equiv\log 4+4\tau$. We find that, in the asymptotic limit, the Rényi entropy is

[TABLE]

The Rényi entropy converges fast to its asymptotic behavior, as shown in Fig. 2. In particular, we see in Fig. 3 that the first asymptotic expansion from Eq. (127) approaches the exact solution exponentially fast, while the leading behavior from Eq. (128) keeps a finite offset.

For a single degree of freedom, the relation between the Rényi entropy $R_{A}$ and the entanglement entropy $S_{A}$ is given by

[TABLE]

with

[TABLE]

This function behaves as

[TABLE]

Hence, for large $x$, we have

[TABLE]

Combining Eqs. (128) and (132), we obtain the asymptotic behavior of the entanglement entropy given by

[TABLE]

where we find the constant to be $1+\frac{1}{2}\log{32/\pi^{2}}\approx 0.44$.

In short, we found that in one spatial dimension, the asymptotic growth of the Rényi and entanglement entropies is logarithmic. As discussed in Section III.2, to any system with a periodic quadratic Hamiltonian there is an associated effective time-independent Hamiltonian that describes its stroboscopic evolution. Moreover, if such an effective Hamiltonian displays a Lyapunov instability, then the asymptotic growth of the entropy is linear, according to the general results summarized in Section III.4. Therefore, we conclude that for $d=1$ the DCE displays no Lyapunov instability, despite the presence of parametric resonance. As stressed in Dodonov and Klimov (1996), the case of $d=1$ is special in the sense that, after the averaging of the fast oscillations, the field modes remain strongly coupled. The intermode interactions compete with the effect of parametric resonance, keeping the production of particles in the resonant mode linear, in opposition to what happens in higher dimensions, where it grows exponentially. Our result show that this suppression of particle creation is reflected in the absence of a Lyapunov instability for the stroboscopic evolution.

The possible asymptotic behaviors of the entanglement entropy for systems with time-independent quadratic Hamiltonians and a finite number of degrees of freedom was fully analyzed in Hackl et al. (2018). In the presence of instabilities, the asymptotic growth is linear, in agreement with the previous results obtained in Bianchi et al. (2018), discussed in Section III.4. In the absence of instabilities, the entanglement entropy either grows logarithmically, for so-called metastable systems, or is bounded. Hence, the DCE in $d=1$ behaves qualitatively as a metastable system, but there is an important difference. In Hackl et al. (2018), it was proved that the prefactor multiplying the logarithmic function must be an integer, i.e. $S_{A}(t)\sim n\log t$, $n\in\mathbb{N}$. In contrast, we found $S_{A}(t)\sim 1/2\log t$, according to Eq. (133). This provides evidence that the extension of the results of Hackl et al. (2018) to systems with an infinite number of degrees of freedom is nontrivial, with new classes of asymptotic behavior arising in this limit. In fact, in preliminary numerical investigations, we found that, by truncating the number of degrees of freedom in the one dimensional DCE, the entanglement entropy becomes bounded, with an initial regime of logarithmic growth. As the size of the truncated system is increased, the regime of logarithmic growth becomes longer, consistently with the fact that for the infinite system such a regime has an infinite duration. This mechanism illustrates how a new class of logarithmic asymptotic growth can arise in the nontrivial limit of infinite number of degrees of freedom and may be investigated elsewhere.

IV.2 (2+1) dimensions

We now proceed to the case of $2$ dimensions. The frequencies of the normal modes are determined by Eq. (46), where $k$ is now a set of two integers. We let the moving mirror oscillate with twice the frequency $\omega_{r}=\omega_{11}$ of the first mode $r=(1,1)$,

[TABLE]

so that the resonant mode is $r$. The length of the cavity along the $x^{2}$-direction is kept fixed. The resonant frequency is

[TABLE]

The Rényi entropy can be computed from Eq. (92) as before. The main difference is that, as discussed in section II.2.2, the evolution of the system is now much simpler. Under the slow-variation approximation, all normal modes decouple, and the evolution is nontrivial only for the resonant mode. The problem is thus reduced to that of a one dimensional parametric oscillator at resonance. The nontrivial Bogoliubov coefficients are given by Eq. (II.2.2)

[TABLE]

where $\gamma$, defined in (50), now reads

[TABLE]

We start again from the covariance matrix in order to compute the Rényi entropy. The sums in the formulas (IV.1) for the matrix components for each mode are now easily computed, since only the first term of each sum is nonzero. For the resonant mode $r=(1,1)$, we obtain

[TABLE]

with

[TABLE]

For the remaining modes, the covariance matrix keeps its initial form throughout the evolution,

[TABLE]

In order to compute the entanglement entropy via Eq. (93), we need to project the covariance matrix to a subspace of interest. Here, a special care must be taken so that a nontrivial evolution is obtained. Since distinct modes evolve independently, if we choose the subsystem to be a single mode, then the entropy will be zero. In order for the entropy to be nonzero, we should consider subspaces that intersect both the resonant mode and its complement.

As a simple example, we study a mixture of two modes. In order to define the subsystem $A$ of interest, let us first introduce an extended subsystem $E$ including the resonant mode and some other (arbitrary) mode $s$. The only mode with a nontrivial evolution is the resonant mode; for all others, the covariance matrix keeps its initial form. As a result, regardless of which mode $s$ one chooses to build the subsystem $E$, the projection of the covariance matrix to $E$ will always have the form

[TABLE]

We now define a new set of annihilation operators

[TABLE]

and their conjugates $\hat{\tilde{a}}^{\dagger}_{1}$ and $\hat{\tilde{a}}^{\dagger}_{s}$. This symplectic transformation defines a new basis on the subsystem $E$. The transformation matrix $B$ that connects the new and the old basis is

[TABLE]

In order to mix the modes and find a nonzero entropy, we apply this transformation to $G_{E}$

[TABLE]

and then restrict to the subspace spanned by the first transformed mode, which gives the $2\times 2$ restricted covariance matrix in the transformed basis, which we write as

[TABLE]

It is now simple to compute the Rényi entropy as

[TABLE]

In the limit of $t\gg 1$, the Rényi entropy becomes

[TABLE]

From Eq. (132), the entanglement entropy has the same asymptotic behavior, differing only by an offset

[TABLE]

Hence, entropy is produced at a constant rate for large times, with a rate determined by the frequency $\omega_{r}$ of the resonant mode, the amplitude of oscillation $\epsilon$ of the mirror, and the dimensions $L_{1},L_{2}$ of the cavity, via the parameter $\gamma$. We see in Fig. 5 that, for large $t$, Eq. (149) gives the same asymptotic behavior of the exact entanglement entropy computed from Eqs. (129) and (IV.2).

This result can be interpreted as signaling the presence of an instability in the system, associated with the fast amplification of excitations in the resonant mode. As discussed in Bianchi et al. (2018), in any system displaying Lyapunov instabilities, the asymptotic production of entropy for a generic subsystem takes place at a constant rate $\Lambda_{A}$, which is determined by the Lyapunov exponents of the system. We will now show that $\Lambda_{A}=\omega_{r}\gamma$ in the present case, in agreement with Eq. (149).

The asymptotic rate of entanglement entropy production $\Lambda_{A}$ for a generic subsystem with $N_{A}$ degrees of freedom is given by the sum of the largest $2N_{A}$ Lyapunov exponents Bianchi et al. (2018). Since we are dealing with a single degree of freedom, $N_{A}=1$ here, and $\Lambda_{A}$ is the sum of the two largest Lyapunov exponents. We will later give a precise characterization of such generic subsystems.

The Lyapunov exponents can be computed from the Hamiltonian flow $M(t)$ of the system. In the special case where $M(t)=\exp(Kt)$, they are simply the eigenvalues of the time-independent symplectic generator $K$. In our case, the generator decomposes into a series of $2\times 2$ blocks $K_{k}$ associated with the independent modes $k$, as discussed in Section II.2.2. The evolution of the nonresonant modes is trivial, $K_{k}=0$, for $k\neq r$, and the associated Lyapunov exponents vanish. The symplectic generator of the resonant mode is given in Eq. (57). The Lyapunov exponents are the eigenvalues

[TABLE]

The two largest Lyapunov exponents of the system are then $\lambda^{(r)}_{1}=\omega_{r}\gamma$ and [math], leading to

[TABLE]

Comparing Eqs. (149), (150) and (151), we see that the asymptotic rate of growth of the entanglement entropy is precisely the value of the positive Lyapunov exponent, which is also equal to $\Lambda_{A}$.

The coefficent $\Lambda_{A}$ can also be computed directly from the equation of motion of the normal modes in the case of parametric resonance. A simple model of parametric resonance is provided by the Mathieu equation

[TABLE]

where

[TABLE]

is the time-dependent natural frequency, which oscillates around $\omega_{0}$ with frequency $\Omega=2\omega_{0}$ and amplitude $\alpha$. It is a well known fact that the Floquet exponents $\mu$ can be computed directly from this equation as

[TABLE]

In addition, the Floquet exponent is the positive Lyapunov exponent for the stroboscopic evolution in this case. Therefore, in the dynamical Casimir effect, we can write the time-dependent frequency $\Omega_{r}(t)$ of the resonant mode, given by Eq. (46), in the form (153) and compute the Lyapunov exponents directly from Eq. (154). With $r=(1,1)$, we have

[TABLE]

To first order in $\epsilon$,

[TABLE]

Comparing with the Mathieu equation, we see that $\alpha=2\epsilon\pi^{2}/L_{1}^{2}$. Substituting into Eq. (154), setting $\omega_{0}=\omega_{r}$, and using Eq. (137), we find

[TABLE]

which agrees with (150).

In the more explicit derivation of the entropy production rate leading to Eq. (149), we considered a special choice of subsystem $A$. The result remains valid in a much more general context, however. In fact, the arguments leading to the same result in Eq. (151) apply to a large class of generic finite-dimensional subsystems, which can be fully characterized from Theorem 2. In the next section, we provide a detailed characterization of such generic subsystems that applies to any spatial dimension $d\geq 2$.

IV.3 Higher dimensions

For any spatial dimension $d\geq 2$, the dynamical Casimir effect at resonance for the lowest energy mode has the same overall features in the slow-variation approximation, as discussed in Section II.2.2. The time evolution is trivial for all modes, except for the resonant one, which behaves as a single oscillator at parametric resonance. The evolution of the resonant mode is described by the symplectic transformation (56), with the natural frequency $\omega_{r}$ and the coefficient $\gamma$ given by Eqs. (48) and (50), respectively. A Lyapunov instability is associated with the resonance. For a generic subsystem, entropy is then produced asymptotically at a constant rate $\Lambda_{A}=\omega_{r}\gamma$, which is equal to the positive Lyapunov exponent, $\Lambda_{A}=\lambda_{1}^{(r)}$. We now present a simple characterization of such generic subsystems and prove that they satisfy $\Lambda_{A}=\omega_{r}\gamma$.

Let $\ell_{1},\ell_{2}\in V^{*}$ be the eigenvectors of the generator $K^{\intercal}$ whose eigenvalues $\lambda_{1},\lambda_{2}$ are the nonzero Lyapunov exponents $\lambda_{1}^{(r)}=-\lambda_{2}^{(r)}$,

[TABLE]

They form a canonical pair of $V$, $\Omega^{ab}(\ell_{1})_{a}(\ell_{2})_{b}=1$. The remaining eigenvectors of $K^{\intercal}$ have vanishing eigenvalues. We can construct a linear basis $\mathcal{D}_{L}=\{\ell_{i}\}$ of $V$ composed of canonically conjugate eigenvectors of $K$ that includes $\ell_{1}$ and $\ell_{2}$. The evolution of the basis vectors under the Hamiltonian flow has a simple form:

[TABLE]

This basis is adequate for analyzing the expansion of subsystems with time, which determines the asymptotic growth of the Rényi and entanglement entropies, as discussed in Section III.4.

Let $A$ be a proper symplectic subspace of $V$ representing a subsystem of interest. Its symplectic complement $B$ is defined as

[TABLE]

We have $V=A\oplus B$, and any vector can be uniquely decomposed as $v=v_{A}+v_{B}$, with $v_{A}\in A$ and $v_{B}\in B$. The symplectic projection of $v$ to the subsystem $A$ is the component $v_{A}$ in this decomposition. The decomposition $V=A\oplus B$ induces a unique decomposition $V^{*}=A^{*}\oplus B^{*}$ of the dual phase space based on the isomorphism induced by $\Omega_{ab}^{-1}$, i.e. we can define $A^{*}=\{\Omega^{-1}_{ab}v^{b}|v^{a}\in A\}$ and $B^{*}=\{\Omega^{-1}_{ab}v^{b}|v^{a}\in B\}$. We call $A$ (or, equivalently, $A^{*}$) a generic subsystem if:

• (i)

$\ell_{2}$ has a nonzero symplectic projection on $A^{*}$;

• (ii)

$\ell_{2}$ has a nonzero symplectic projection on $B^{*}$.

The condition (i) implies that there exists a vector $\theta_{1}\in A^{*}$ such that

[TABLE]

Expanding $\theta_{1}$ in the basis $\mathcal{D}_{L}$, $\theta_{1}=\sum_{k^{\prime}}\theta_{1}^{k^{\prime}}\ell_{k^{\prime}}$, and exploring the fact that $\mathcal{D}_{L}$ is a Darboux basis, we find that $\theta_{1}^{1^{\prime}}\neq 0$, that is, $\theta_{1}$ has a nonzero expansion coefficient in the first basis vector of $\mathcal{D}_{L}$.

Now let $\mathcal{D}_{A}=\{\theta_{r}\mid r=1,\dots,2N_{A}\}$ be a Darboux basis of the subsystem $A$ that includes $\theta_{1}$. Expanding its elements in the basis $\mathcal{D}_{L}$, we can write $\theta_{r}=\theta_{r}^{k^{\prime}}\ell_{k^{\prime}}$. By a simple procedure of Gaussian elimination, we can construct a new linear basis $\tilde{\mathcal{D}}_{A}=\{\tilde{\theta}_{r}\}$ of $A$ such that $\tilde{\theta}_{1}=\theta_{1}$, and $\tilde{\theta}_{r}^{1}=0$, for $r\neq 1$. Moreover, it follows from condition (ii) that, for all $r$, $\tilde{\theta}_{r}^{k^{\prime}}\neq 0$ for some $k^{\prime}\neq 2$, otherwise $\ell_{2}$ would be an element of $A$ and have a vanishing symplectic projection on $B$. In short, the only vector with a contribution from the expanding unstable mode is $\theta_{1}$, and all other modes have a contribution from some stable mode. As a result, for large times, we have

[TABLE]

We are interested in computing how the volume of the unit cube with sides $\theta_{r}$ evolves for large times, $\operatorname{Vol}(M^{\intercal}(t)\mathcal{D}_{A})$. Since the basis $\mathcal{D}_{A}$ is related to $\tilde{\mathcal{D}}_{A}$ by a fixed transformation, their volumes are related by a time-independent factor, $\operatorname{Vol}(M(t)\mathcal{D}_{A})\propto\operatorname{Vol}(M(t)\tilde{\mathcal{D}}_{A})$. But then, from Eq. (IV.3),

[TABLE]

and we find for the subsystem exponent:

[TABLE]

This completes the proof that $\Lambda_{A}=\omega_{r}\gamma$ for a generic subsystem satisfying conditions (i) and (ii).

Note that in $(2+1)$ dimensions and higher, the Lyapunov instability is restricted to the resonant mode, so that the condition (ii) is violated for the subsystem composed of the resonant mode only. This is the reason why we had to consider a subsystem mixing more than one mode to obtain a nontrivial behavior in Section IV.2.

V Discussion

We analyzed the production of Rényi and entanglement entropies for the dynamical Casimir effect at resonance in arbitrary dimensions. In our settings, a $d$-dimensional parallelepided has one of its perfectly reflecting boundaries oscillating harmonically, while all others are fixed. By setting the frequency of oscillation equal to twice that of some normal mode $r$, the system is set at parametric resonance and particles are continuously produced from the vacuum. We considered the regime where the oscillations of the mirror are small, with a relative amplitude $\epsilon=\Delta L/L\ll 1$. We combined analytical techniques introduced in Dodonov and Klimov (1996) for the study of pair creation under these circumstances with general symplectic techniques applicable to the analysis of entanglement production in general time-dependent linear systems Bianchi et al. (2015, 2018) in order to provide a thorough description of the build up of correlations and entanglement at resonance.

The system has two characteristic time scales: the period of the resonant mode, $\omega_{r}^{-1}$, and the time scale $(\omega_{r}\epsilon)^{-1}$ for which the system starts to depart considerably from the initial configuration at the ground state, which depends on the amplitude of oscillation of the mirror. Since $\epsilon\ll 1$, these scales are well separated. By averaging over the fast oscillations of normal modes, the equations of motion for the mode amplitudes are simplified and can be solved exaclty, as shown in Dodonov and Klimov (1996). The cases of $d=1$ and $3$ are treated in Dodonov and Klimov (1996), where a Heisenberg representation is used for $d=1$ and a Schrödinger representation for $d=3$. We extended the treatment for arbitrary dimensions, and showed that it can be implemented in the Heisenberg representation in any dimension.

The moving boundary induces a nontrivial evolution of the cavity modes, which in general includes two kinds of interactions: particle creation from the vacuum and mode coupling. The vacuum becomes unstable due to the time-dependent boundary conditions, and the created excitations can be transferred from one mode to another.

For any spatial dimension $d\geq 2$, the intermode interactions are strongly suppressed. This is due to the fact that the energies required to excite nonresonant modes do not match differences between energy levels of the resonant mode, making such transitions unlikely. The amplitude of the resonant mode then increases exponentially, being continuously pumped by the parametric resonance, providing an example of a Floquet instability. In the approximation where the fast oscillations are integrated out, the mode amplitude grows monotonically, and the original Floquet instability then corresponds to a Lyapunov instability. We found that, for any subsystem that intersects both the resonant mode and its complement, Rényi and entanglement entropies are produced at a constant rate for large times, $S_{A}\sim R_{A}\sim\Lambda_{A}t$, where the production rate is equal to the positive Lyapunov exponent of the system, as expected from general results of Bianchi et al. (2018). It can be explicitly written as $\Lambda_{A}=\lambda_{1}^{(r)}=\omega_{r}\gamma$, with $\omega_{r}$ and $\gamma$ given by Eqs. (48) and (50).

An alternate technique for studying a Floquet instability in a system with periodicity $\tau$ consists of looking at its stroboscopic evolution. That is, if $M(t)$ is its Hamiltonian flow, one focus on the discrete series of instants $M(n\tau)$, $n\in\mathbb{Z}$. An effective, time-independent Hamiltonian can then be defined in terms of $\tau^{-1}\log M(\tau)$ such that at full periods its evolution agrees with that of the original system. This effective Hamiltonian displays a Lyapunov instability associated with the original Floquet instability. Production of entropy in systems with Floquet instabilities was studied in detail in Bianchi et al. (2018) within the stroboscopic approach. Applying it to the dynamical Casimir effect at resonance for $d\geq 2$, we found the same results for the entanglement entropy production as in the approximation introduced in Dodonov and Klimov (1996), where the fast oscillations are integrated out. Our result suggests that, for systems satisfying a slow-variation condition, $M(\tau)\sim\mathbb{1}$, the technique of averaging out the fast oscillations provides a good approximation for computing the Lyapunov exponents of the effective Hamiltonian associated with the stroboscopic dynamics of the system.

The case of one spatial dimension displays special features that lead to a different asymptotic behavior for the Rényi and entanglement entropies. In constrast to what happens in higher dimensions, for $d=1$, the energy required to excite any normal mode is a multiple of the difference $\hbar\omega_{r}$ between energy levels of the lowest energy mode, which we set at resonance. All modes are then strongly coupled, and excitations produced at the resonant mode can transition to other modes. The energy injected in the field by the moving mirror thereby spreads throughout an infinite number of modes. As a result, two new effects take place. First, the number of particles in the resonant mode, instead of growing exponentially, increases only linearly with time, as first recognized in Dodonov and Klimov (1996). In addition, the resonant mode is strongly entangled with the nonresonant modes. We found that the associated Rényi and entanglement entropies of the resonant mode display a logarithmic growth, $S_{A}\sim 1/2\log\tau$, where $\tau=\epsilon\omega t/2$ is a dimensionless time variable. It is interesting that, even in the presence of parametric resonance, the system does not display features expected for unstable systems, i.e., a linear production of entropy. Instead, its behavior is qualitatively similar to that of metastable Hamiltonians discussed in Hackl et al. (2018), except for the fact that the logarithmic function is here multiplied by a factor of $1/2$, while for metastable systems the prefactor is always an integer.

The techniques here applied for the analysis of entropy production in the dynamical Casimir effect at resonance can also be applied to other mirror trajectories. For finite cavities, the adaptation of the procedure followed in this work should be straightforward. A natural question for further developments is how to implement our techniques for the dynamical Casimir effect with a single mirror. As is well known, one of the main motivations for the study of the dynamical Casimir effect is that it can mimic the Hawking effect Davies and Fulling (1976). The application of symplectic techniques for the study of entropy production in general field theories is discussed in Bianchi et al. (2018), and we hope our approach can be extended along such lines to the analysis of the time evolution of the entanglement entropy in models simulating aspects of particle creation in the Hawking effect Carlitz and Willey (1987); Good et al. (2013).

Acknowledgements.

We thank Eugenio Bianchi and Rodolfo R. Soldati for inspiring discussions during the development of this project. We thank Eugenio Bianchi for sharing his notes on Floquet instabilities and for many invaluable discussions during the conception and initial stages of this work. IR acknowledges support from CAPES and CNPq. LH thanks the members of the physics department at the Universidade Federal de Minas Gerais for their hospitality during his visit in the summer of 2018. LH acknowledges support through the Brazil-U.S. Exchange Program of the American Physical Society which made this visit possible. LH is funded by the the Max Planck Harvard Research Center for Quantum Optics and supported by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy – EXC-2111 – 39081486. NY acknowledges support from CNPq, Brazil, under the grant PQ 306744/2018-0, and from the Programa Institucional de Auxílio à Pesquisa de Docentes Recém-Contratados, PRPq/UFMG.