Entanglement dynamics of coupled oscillators from Gaussian states

Entanglement is one of the most distinctive features of the quantum mechanics. The counter-intuitive behavior of quantum entanglement, which produces various paradoxical physical phenomena, has confounded the physics community for a very long time. Yet, quantum entanglement has gained a revived interest with the developments in quantum information processing, as it plays a indispensable role, being the source of quantum computation, quantum communication adesso2014continuous ; braunstein2005quantum ; horodecki2009quantum ; nielsen2002quantum ; weedbrook2012gaussian , quantum cryptography ekert1991quantum and metrology nichols2018multiparameter .

Regarding the entanglement between subsystems constituting a composite system, over the last couple of decades, there have been several attempts to quantify the degree of entanglement in a given bipartite quantum state. One of the measures proposed is the entropy of entanglement or entanglement entropy, which is a natural measure of the degree of quantum correlations for pure bipartite states in composite systems. The first motivations for studying the features of entanglement entropy came from the physics of black holes. In the seminal work of Sorkin et al. Bombelli:1986rw , the analytical expression for the ground state entanglement entropy of a system of coupled harmonic oscillators has been derived. Following this milestone, the analysis of Bombelli:1986rw has been extended to the bosonic scalar fields in reference Srednicki:1993im . In his pioneering work pointing out a possible connection between black hole thermodynamics and entanglement entropy, Srednicki has derived an area law, which scales as the boundary surface dividing the two subsytems and shows a high similarity to the intrinsic entropy formula of a black hole.

On the other hand, since its widespread adoption, there has been a robust interest from the quantum information community to gain enhanced understanding about entanglement entropy due to its capability of providing accurate descriptions of various physical phenomena related to quantum many-body physics osborne2002entanglement ; osterloh2002scaling ; vidal2003entanglement . Traditionally, the study of entanglement in the scope of quantum information theory has primarily concerned itself with the research in discrete physical systems with quantum states defined in finite dimensional Hilbert spaces. Nevertheless, although a daunting complexity arises due to the structure of entangled states in infinite dimensional Hilbert spaces, there has been a growing curiosity for investigating entanglement in continuous variable systems. In this regard, it is essential to remark that a particular family of continuous states called the Gaussian states, which are quantum states that are represented by Gaussian-Wigner quasiprobability functions, have provided the means to explore the dynamics of entanglement in continuous variable systems Schumaker1986317 . These states naturally arise as the thermal or ground states of any bosonic physical system described by at most second order quadrature Hamiltonians in canonical coordinates. Despite the enormous complexity of the underlying infinite dimensional Hilbert space, they can be completely characterized by finite number of degrees of freedom. Besides, Gaussian states play an important role in quantum optics, as they can be easily created and manipulated in a vast number of optical setups. Furthermore, they naturally appear in areas of physics like ion traps, gases of cold atoms and optomechanics adesso2014continuous ; olivares2012quantum .

In this paper, our main interest is to analyze the dynamics of entanglement for an isolated quantum system consisting of two coupled harmonic oscillators after a sudden quench. The departure point of our analysis is the use of a numerical method (the so-called Gaussian state approximation) that relies on the estimation of the system’s Wigner representation by a Gaussian function that is characterized by time-dependent wave packet centers and dispersion Buividovich:2017kfk ; Buividovich:2018scl . The paper is structured as follows. Section 2 starts out with brief introductions to the system of two coupled oscillators and the formalism of covariance matrices, which is followed by the derivation of the initial covariance matrix that is employed in the simulations. In section 3, we first describe how the covariance matrix representations of Gaussian states can be utilized to numerically compute entanglement entropy. Subsequently, the derivation of the analytical expression for the time-dependent von Neumann entropy, which is originated from the exact solutions of the time-dependent Schrödinger equation, is reviewed. In section 4, we investigate the time evolution of the bipartite entanglement entropy within the Gaussian state approximation by utilizing the covariance matrix representations of Gaussian states and compare the numerically obtained results with the findings from the analytical entropy expression. Lastly, section 5 is devoted to conclusions and outlook.

We initiate the developments with a consideration of the traditional system of two bosonic oscillators that are coupled by a quadratic Hamiltonian of the form²²2In this paper, we work in natural units by setting $\hbar=1$.

where $\kappa_{0}(t)$ and $\kappa_{1}(t)$ are arbitrary time-dependent parameters. Let us collect the quadrature pairs $\big{(}x_{\alpha},p_{\beta}\big{)}$ $(\alpha,\beta=1,2)$ into the vector operator $\bm{v}\coloneqq\big{(}x_{1},x_{2},p_{1},p_{2}\big{)}^{T}$, which enables us to express the canonical commutation relations satisfied by these operators, namely $[x_{\alpha},p_{\beta}]=i\delta_{\alpha\beta}$, in a more compact form as follows

where $\mu,\nu=1,2,3,4$ and the blocks $0_{2}$ and $\mathds{1}_{2}$ are the two-dimensional zero and identity matrices. $\Lambda_{\mu\nu}$ are the elements of the antisymmetric matrix $\Lambda$ that is known as the symplectic form. Let us also note that the matrix transformation $\Lambda$ is orthogonal, i.e. $\Lambda^{T}\Lambda=\mathds{1}_{4}$ Serafini .

Having now listed some of the essential features of (1), we move on to introduce the Gaussian states. A Gaussian state is defined as any quantum state whose Wigner representation³³3The Wigner representation (or function) is the quantum analog of a Liouville probability density. It is a quasiprobability distribution in the sense that it can take on negative values peres2006quantum . is a Gaussian function. Only two parameters, the first and second statistical moments, are required to completely determine a Wigner function. While the first moments are given by the expectation values $\bar{v}_{\mu}\coloneqq\expectationvalue{v_{\mu}}$, the second moments are defined as weedbrook2012gaussian ; olivares2012quantum ; adesso2014continuous

where $\delta v_{\mu}\coloneqq v_{\mu}-\bar{v}_{\mu}$ are the fluctuation operators and $C_{\mu\nu}$ are the elements of the real, symmetric matrix $C$, which is referred to as the covariance matrix. The Wigner function can be written by

where $v\in\mathbb{R}^{4}$.

To proceed to further, it is convenient to apply the Heisenberg equation of motion to the system of two coupled harmonic oscillators, which results in

The time derivatives of the second moments can be revealed by combining the averages of (5) with the averages of the operator products of the form $v_{\mu}v_{\nu}$ that are both taken over the Gaussian state characterized by (4) Buividovich:2017kfk . To illustrate, one may consider the elements of the first row of the covariance matrix $C$. After a straightforward calculation, the explicit time derivatives of $C_{1\nu}$ elements are evaluated and shown below

Besides, the remaining six equations describing the time derivatives of the covariance matrix are listed in Appendix 48.

In order to study the time evolution of the covariance matrix numerically, apart from determining how $C$ varies in time, we also need to specify an initial configuration. With the aim of preparing such a configuration, we start by defining $S\in Sp(4,\mathbb{R})$

such that $S^{T}\!\Lambda\hskip 0.56905ptS=S\hskip 0.56905pt\Lambda\hskip 0.56905ptS^{T}=\Lambda$. Equation (1) can be first cast into the matrix form, i.e.

and then by diagonalizing $H_{m}$, (8) can be rewritten as follows

where

Now, let us express $H_{d}$ as a direct sum of diagonal matrices

Since $H$ is a quadratic Hamiltonian, it is possible to deduce the ground state covariance matrix from the block diagonal form of $H_{d}$ as illustrated below schuch2006quantum

where

In the next section, a method for computing entanglement entropy from the covariance matrix representation of a Gaussian state is presented. As it will be discussed shortly, the significance of $C_{i}$ relies on the fact that it corresponds to a pure Gaussian state.

The entropy of entanglement is the natural measure of the degree of quantum correlations for pure bipartite states in composite systems. Let $\rho$ denote the density matrix of the quantum system $\mathcal{S}$ with Hilbert space $\mathcal{H}=\mathcal{H}_{a}\otimes\mathcal{H}_{b}$. The entanglement entropy of $\mathcal{S}$ can be determined by computing the von Neumann entropy of the reduced density matrix $\rho_{a}$, which can be obtained by taking the partial trace of $\rho$ over the subspace $\mathcal{H}_{b}$.

This section is devoted to an investigation of the dynamics of entanglement observed in the system of two coupled harmonic oscillators. In order to provide a comprehensive analysis, we carry out numerical simulations of the time evolution of (19) for various distinct frequency settings. After discretizing the equations of motion given in (6) and (48), an iterative algorithm is utilized to solve the discretized equations numerically. As it will be detailed shortly, we start the simulations with initial conditions that are formed from the elements of the matrix $C_{i}$, which is listed in (2). For comparison, aside from the numerical findings obtained from the evaluation of $S_{c}$, we also consider the bipartite entanglement entropy expression of (38) and calculate $S_{a}$ for the frequency combinations that are employed in the simulations of (19). From the numerical computations of $S_{c}$ and the comparisons of $S_{c}$ with $S_{a}$, we hope to gain valuable insight both into the dynamics of entanglement and the effectiveness of the Gaussian state approximation.

In the computations, a simulation code that is implemented in Matlab is used. The code is executed with a constant time step of $\Delta t=0.001$. Due to truncation of digits, errors are inevitable in numerical calculations. In this regard, although (as it is characterizing a pure Gaussian state) $C_{i}$ is expected to keep representing pure Gaussian states due to the quadratic structure of $(\ref{HamOsc})$ Buividovich:2017kfk ; Buividovich:2018scl , the cumulative effect of rounding errors could cause a violation of this preserving map and create mixed states. However, by constantly monitoring the symplectic eigenvalues during the the trial runs of the simulation, we made sure that no such effect is present.

Having now introduced the basic features of numerical computations, we move on to the details of starting configurations. The simulations are initiated from the previously determined covariance matrix $C_{i}$ defined by $(\ref{Ci})$. Namely, at the initial time $t=0$, we have

while the rest of the initial conditions are zero. Before moving on to the presentation of numerical findings, we need to point out that the simulations are started with frequencies $\omega_{p0}$ and $\omega_{m0}\,$. Upon the completion of the first time increment $\Delta t$, these frequencies are updated to $\omega_{p1}$ and $\omega_{m1}\hskip 0.56905pt$, which remain unchanged for the second time increment and the rest of the computations. This is how the quenching mechanism of $(\ref{QnchMec})$ is included in the implementation of the algorithm.

The comparison between two distinct measures of von Neumann entropy and the long-time behavior analysis of $S_{c}\hskip 0.56905pt$, both of which would shed light on the effectiveness of the Gaussian state approximation, are made through the use of figures and tables. In Figure 1, we present sample plots for the time series of $S_{a}$ and $S_{c}$ after a sudden quantum quench. While the initial frequencies are held fixed at $\omega_{p0}=0.9$ and $\omega_{m0}=4.9$, the final frequencies, i.e. $\omega_{p1}$ and $\omega_{m1}$, are varied. The first thing we can immediately observe from the plots is that the oscillation frequencies of both $S_{a}$ and $S_{c}$ decrease with reduced final frequencies. In Figure 1(a),

the ratio of the oscillation periods of two entropy measures is approximately equal to $2.3\hskip 0.56905pt$. When both $\omega_{p1}$ and $\omega_{m1}$ are reduced by $0.05$, this ratio shows the tendency to converge to $1$ as it can be observed from Figure 1(b). Furthermore, in the same graph, the difference between $S_{a}$ and the numerically computed entanglement entropy $S_{c}$ is the smallest among all the plots. In particular, for $t\lessapprox 83\hskip 0.56905pt$, the Gaussian state approximation is highly effective for the frequency combination depicted in Figure 1(b). Similarly, in Figures 1(c) and 1(d), despite an apparent increase in the difference of the oscillation amplitudes, $S_{c}$ certainly shows the ability to follow the oscillatory motion of $S_{a}$.

Regarding the possible sources of error in the computations of the numerically calculated entropy, it is essential to note that despite its capability in conserving the von Neumann entropy by evolving pure states into pure states, the Gaussian state approximation describes a non-unitary evolution of the density matrices Buividovich:2018scl . This is certainly not true for the case of equation (38), which is emanated from the exact solutions of the time-dependent Schrödinger equation and describes a unitary time evolution.

On the other hand, in order to take the effects of the starting frequencies into consideration, we illustrate in Figure 2 the evolutions of the entanglement entropies with time at four different initial frequency combinations. Despite the presence of periodic ripples in Figure 2(d), as it can be clearly seen from all four plots, $S_{c}$ is capable of following $S_{a}\hskip 0.56905pt$, which shows an upward trend with time.

Lastly, we would like to report on the results of long-time behavior analysis of $S_{c}\hskip 0.56905pt$. Figure 3 shows the plots of $S_{c}$ versus time at the frequency combinations used for the preparation of Figure 2. This time we run the code for a sufficient amount of time

to clearly observe the values that the entanglement entropies converge to. The best-fitting model for the numerical data displayed in Figure is found to be a logarithmic function in the form

The fitting parameters of the best fit equations (47) are listed in Table 1. The adjusted R-squared values of the fitting curves depicted in Figures 3(a) - 3(d) are given by $0.9581$, $0.9394$, $0.9134$ and $0.884$ respectively. Although, due to the apparent increase in the variance of the data illustrated in Figures 3(c)-3(d), $\Phi_{3}$ and $\Phi_{4}$ fits are not as good in comparison to the

In this paper, we have considered the dynamics of entanglement of a system of two coupled harmonic oscillators following an instant quantum quench in the inherent parameters of the system. We have performed a detailed numerical analysis of the time evolution of the bipartite entanglement entropy within the Gaussian state approximation, whose accuracy has been verified qualitatively through the comparison with the analytical expression for the time-dependent entanglement entropy. Besides, from the results concerning the change of von Neumann entropy with time, we were able to show through a proper fitting function that entanglement entropy vary logarithmically in time under the same approximation.

Let us also mention a recent development regarding the possible applications of the Gaussian state approximation. Although calculating entanglement entropy in ordinary field theories is a rather difficult task, the numerical computations of the entanglement entropy in the Banks-Fischler-Shenker-Susskind (BFSS) matrix model were recently performed in reference Buividovich:2018scl by using the covariance matrix representations of Gaussian states. In this regard, a valuable direction of research would be to investigate the time-dependence of the entanglement entropy in a Yang-Mills matrix model with two distinct mass deformation terms, whose emerging chaotic motions and dynamics of thermalization have been recently studied in a series of papers Baskan:2019qsb ; Oktay:2020tzi . Aside from the mass deformation terms keeping the gauge invariance intact, this model contains the same matrix content as the bosonic part of the BFSS matrix model and provides an ideal system for the use of Gaussian state approximation. Moreover, with the help of the results obtained in Oktay:2020tzi , it would be interesting to investigate the possible use of entanglement entropy as a probe of thermalization. Another challenging direction of development is to generalize the approach and methods developed in this work for the analysis of the variation of entanglement entropy of two time-dependent, coupled harmonic oscillators to the case of $N$ time-dependent, coupled harmonic oscillators. We hope to report on progress in these directions elsewhere.