PHYSICAL REVIEW E 85, 036208 (2012)

Classical dynamics of quantum entanglement

Giulio Casati,1,2 Italo Guarneri,3,4 and Jose Reslen5

1Consorzio Nazionale Italiano di Struttura della Materia, Consiglio Nazionale delle Ricerche, Istituto Nazionale per la Fisica della Materia,and Center for Nonlinear and Complex Systems, Universita degli Studi dell’Insubria, Via Valleggio 11, I-22100 Como, Italy

2Istituto Nazionale di Fisica Nucleare, Sezione di Milano, via Celoria 16, I-20133 Milano, Italy3Center for Nonlinear and Complex Systems, Universita degli Studi dell’Insubria, Via Valleggio 11, I-22100 Como, Italy

4Istituto Nazionale di Fisica Nucleare, Sezione di Pavia, via Bassi 6, I-27100 Pavia, Italy5Department of Physics, National University of Singapore, Singapore 117543, Singapore

(Received 16 September 2011; published 19 March 2012)

We analyze numerically the dynamical generation of quantum entanglement in a system of two interactingparticles, started in a coherent separable state, for decreasing values of h. As h → 0 the entanglement entropy,computed at any finite time, converges to a finite nonzero value. The limit law that rules the time dependenceof entropy is well reproduced by purely classical computations. Its general features can be explained by simpleclassical arguments, which expose the different ways entanglement is generated in systems that are classicallychaotic or regular.

DOI: 10.1103/PhysRevE.85.036208 PACS number(s): 05.45.Mt, 03.67.Mn, 03.65.Ud

I. INTRODUCTION

The development of the study of quantum chaos has shownthat quantum dynamics is deeply influenced by classicalchaotic motion. It is then natural to inquire how the dynamicalgeneration of quantum entanglement is affected by the chaoticor regular nature of the underlying classical dynamics [1–8];an especially intriguing aspect of this question is that quantumentanglement does not have a classical counterpart. Since theearly paper by Furuya et al. [1], most investigations about thissubject have supported the view that regular classical dynamicsleads to slower dynamical generation of entanglement (DGE)[2]. The opposite point of view, that the DGE is scarcelyaffected by the nature of classical dynamics, has also beensurmised [3]; nevertheless, there are convincing reasons tobelieve that the way quantum entanglement is dynamicallygenerated does indeed have a classical counterpart [4,9,10].On the one hand, quasiclassical analytical methods have beenimplemented to investigate the decay of purity of reducedstates in integrable systems [5,6]. On the other hand, DGEmay be studied by resorting to a phase-space formulation ofquantum mechanics [4] and then the purity of reduced states,which is a measure of entanglement (or disentanglement),can be directly obtained from the phase-space functions thatrepresent such quantum states in the classical phase spaces.By letting such functions evolve according to the classicalLiouville dynamics, a classical counterpart to DGE is obtained[4]. Along this line of thought, it was then observed [9]that dynamically generated entanglement, as measured by thereduced linear entropy, may not tend to 0 in the limit h → 0and instead approaches a purely classical quantity in the caseof Gaussian initial states.

In this paper we study the h → 0 limit of DGE in a systemof two interacting particles, which has well-defined quantumand classical Hamiltonians. The classical dynamics can beturned from completely integrable to essentially chaotic bysuitably tuning parameters. The quantum system is initiallyprepared in a separable coherent state so in the limit h → 0,the evolving quantum state tends to a classical pure state(i.e., a Dirac δ function in phase space, supported in a

point of the classical trajectory). In spite of that, we observethat the quantum entanglement does not vanish in the limith → 0. Here we measure entanglement by the von Neumannentropy of reduced states and find that it is well approximatedby a coarse-grained classical entropy, computed by solvingHamilton’s equations. Thus our simple-minded approach toDGE in terms of classical coarse-grained entropy opens theway to an intuitive cell-counting analysis of how DGE isaffected by regular or chaotic classical motion and allowsfor predictions that would hardly be possible within a strictlyquantum mechanical formalism.

II. MODEL: QUANTUM AND CLASSICAL ENTROPIES

We consider a Toda model of two interacting particlesmoving in a line under the following Hamiltonian:

H = p21

2m1+ p2

2

2m2+ e−q1 + e−(q2−q1) + eq2 − 3. (1)

Relevant parameters are the energy E and the mass ratiom1/m2. For m1 = m2 the model is completely integrable nomatter how large the energy is [Fig. 1(a)], while for anyother mass ratio (except limit values 0 and ∞) the dynamicsexhibits a transition from integrability to chaos as the energyE is increased [11]. In particular, for m2/m1 = 0.54 motion ispredominantly chaotic already at energy E = 7 [Fig. 1(b)].

To compute the quantum evolution we first calculatedanalytically matrix elements of the quantum Hamiltonian[which is obtained from Eq. (1) by canonical quantization] overthe energy eigenbasis of two uncoupled harmonic oscillators ofunit mass and then diagonalized the matrix thus obtained. Withthe chosen matrix size, the computed eigenbasis was verifiedto support the chosen initial states within negligible errors.Such initial states were separable coherent states, centered atq1 = q2 = 0 and p1 = p2 = √

E in the regular case (m1 =m2 = 1) and q1 = q2 = 0, p1 = √

E, and p2 = −√m2E in

the chaotic case (m1 = 1 and m2 = 0.54). The latter initialcondition is inside the chaotic sea. As time goes on theparticles become entangled due to interactions. A measure

036208-11539-3755/2012/85(3)/036208(5) ©2012 American Physical Society

GIULIO CASATI, ITALO GUARNERI, AND JOSE RESLEN PHYSICAL REVIEW E 85, 036208 (2012)

-4

-2

0

2

4

-2.5 -1.5 -0.5 0.5 1.5 2.5

p1

q1

(a)

-4

-2

0

2

4

-2.5 -1.5 -0.5 0.5 1.5 2.5

p1

q1

(b)

FIG. 1. Surfaces of section for the Toda model at energy E = 7:(a) integrable case with masses m1 = m2 = 1 and (b) mixed phasespace for m1 = 1 and m2 = 0.54.

of their entanglement at time t is given by the von Neumannentropy S(t) = −Tr(ρ1lnρ1), where ρ1 is the reduced densitymatrix of particle 1. The thus-computed time dependence ofentanglement is shown in Fig. 2 for three different valuesof h. One observation is that integrable and chaotic casescan hardly be told from each other by just comparing thedegrees of entanglement to which they give rise. In this respectit is worth noting that the dependence of entanglement oninitial conditions, which was observed in several papers (e.g.,Ref. [7]) as the initial state was moved from the chaotic sea toan integrable island, is not a signature of chaos. Instead, it ismerely due to a decrease of the number of accessible quantumstates as one approaches the periodic orbit in the center of theisland. Indeed, a similar, significant decrease of entanglementis recorded even in the completely integrable case of our model(not shown), starting with initial conditions that approach thecentral periodic orbit of the stable island, e.g., along a verticalline (q1 = 1.2) in Fig. 1(a).

However, the following are the most interesting observa-tions: First, as h decreases, the curves S(t,h), which describethe entanglement of the quantum Toda system as a functionof time t at a given value of h, come closer and closer to anontrivial limit curve (Fig. 2, full curve). Second, this curveis well approximated by a purely classical construction, to bedescribed next.

From the classical Gaussian ensemble in the four-dimensional phase space that corresponds to the initial quan-tum coherent state we have generated a sample of M initialpoints. To compute the evolution of this initial ensemble we

0

1

2

3

4

5

6

0 10 20 30 40 50

Time

(a)

δ=0.02-h=0.0375-h=0.1-h=0.2

0

1

2

3

4

5

6

7

0 10 20 30 40 50

Time

(b)

δ=0.02-h=0.0375-h=0.1-h=0.2

FIG. 2. (Color online) Entanglement entropy vs time for E = 7and several values of h in (a) a regular and (b) a chaotic case. The fullcurve is classical entropy.

have integrated the classical Hamilton equation, using a fourth-order Runge-Kutta routine. The four-dimensional ensemblesthus computed at times t > 0 were projected onto reducedensembles in the (q1,p1) plane. To obtain a classical analogfor the von Neumann entropy of quantum reduced states, theBoltzmann-Shannon entropies Scl(δ,t) were computed using afixed partition of the two-dimensional phase space in squarecells of side

√δ:

Scl(δ,t) = −∑

i

wi(t)

Mln

wi(t)

M, (2)

where wi(t) is the number of points in the reduced ensemble,which were found in the ith cell at time t . We used enoughcells to guarantee hosting all points at all times. This classicalcoarse-grained entropy is expected to somehow approximatethe quantum entropy provided the classical ensemble populatesa large number of states.

We numerically observe that, as δ → 0, the curve Scl(δ,t)approaches a limit curve S∞(t). This seemingly crude classicalconstruction works surprisingly well; indeed, the classicallimit curve is approached by the curves of quantum entan-glement as h → 0. At nonzero h and δ the quantum-classicalagreement is worse at the minima of the curves and alsoin the early stage of evolution. A possible explanation isthat, as already noted, the classical coarse-grained entropy isexpected to well approximate the quantum entropy only whenthe ensemble populates a large number of cells. Regardless,the overall approximation is so good that even details of

036208-2

CLASSICAL DYNAMICS OF QUANTUM ENTANGLEMENT PHYSICAL REVIEW E 85, 036208 (2012)

the quantum curves, such as superimposed, almost regularoscillations (Fig. 2), can be explained in terms of classicalphase-space structures, as we shall presently explain.

First, however, it is important to stress that our findingsdo not involve any contradiction. Measures of quantumentanglement need not tend to 0 in the h → 0 limit, as shownby simple examples; e.g., a system of two identical harmonicoscillators in the pure state 1√

2H0 ⊗ H0 + 1√

2H1 ⊗ H1, where

H0 and H1 are the ground state and the first excited stateof a harmonic oscillator, has entanglement entropy equal toln(2), independently of the value of h > 0, even though in thelimit h → 0 its Wigner function tends to a Dirac δ functionsupported in the origin of the four-dimensional phase space.In contrast, the meaning of the above-described classicalconstruction is an interesting question in purely classicalterms, which may be related to a recently proposed classicalseparability entropy [12].

III. THEORETICAL INTERPRETATION

In Figs. 3 and 4 we show the numerically computed classicaland quantum entropies, respectively. All the curves shownshare some general qualitative features: Notably, at long timesthey approach saturation values via an average monotonic

0

1

2

3

4

5

6

7

8

1 10 100 1000

Scl

Time

δ=0.02

δ=0.1

δ=0.5

δ=1

Scl ~ln(t2) (a)

16

32

64

128

256

1 10 100

4

5

6

7

8

tD-Scl

1/δ

0

1

2

3

4

5

6

7

8

0 50 100 150 200

Scl

Time

δ=0.02

δ=0.1

δ=0.5

δ=1

(b) Scl ~t/10

20

40

60

80

0 1 2 4 8 16 32 64

4

8

12

16

tD -Scl

1/δ

FIG. 3. (Color online) Classical entropy Scl(t) vs time for differ-ent values of the cell area δ for E = 7: (a) an integrable case, withthe dashed line representing const + 2 ln(t), and (b) a chaotic case,with the dashed line representing const + t/10. The insets show thesaturation time tD (circles, left axis) and saturation value Scl (squares,right axis) versus 1/δ. In the regular case, tD ≈ const + δ−0.43 andScl ≈ const − 0.86 ln(δ). In the chaotic case, tD ≈ const − 9.3 ln(δ)and Scl ≈ const − 0.93 ln(δ).

-1

0

1

2

3

4

5

6

7

0 50 100 150 200

Time

Scl ~ln(t2)

(a)

-h=0.2

-h=0.1

-h=0.0375

δ=0.02

8

16

32

4 8 16 32

2

3

4

5

tD-S

1/-h

-1

0

1

2

3

4

5

6

7

8

0 50 100 150 200

Time

(b) Scl ~t/10

-h=0.2

-h=0.1

-h=0.0375

δ=0.02

10

20

30

40

4 8 16 32

2

3

4

5

tD -S

1/-h

FIG. 4. (Color online) Entanglement entropy S(t) vs time fordifferent values of h, in an integrable case (above) and in a chaoticcase (below). The case δ = 0.02 corresponds to the classical entropyScl shown also in fig 3. Insets: saturation time tD (circles - left axis)and saturation value S (squares - right axis) of entanglement, versus1/h. In the regular case: tD ≈const+h−0.43, S ≈ const−0.86 ln(h). Inthe chaotic case: tD ≈const−8.6 ln(h); S ≈const−0.86 ln(h)

trend with strong oscillations superimposed. Such oscillationslook surprisingly regular even in chaotic cases and graduallydecrease in time; in the classical cases, at saturation they are onthe order of fluctuations due to finite statistics. As the cell areaδ is decreased, for times less than a time scale tD(δ), a limitcurve is approached. This curve is different in the regular andthe chaotic cases and the saturation time tD may be thought ofas the time when a finite δ curve departs from the limit behavior.For times t > tD(δ), entropy curves approach a saturation valueand this value scales proportionally to ln(1/δ) (see the insets inFig. 3). To infer estimates of saturation times from numericaldata, tD(δ) is defined as the time when the average growth(represented by dashed lines) of the limit curve reaches thesaturation value at the given δ.

Such numerical observations are easily explained by in-terpreting entropy as the logarithm of an effective numberof populated cells. To this end, one may first replace theinitial Gaussian ensemble by a single four-dimensional cellof size

√δ. Saturation is then an obvious consequence of

the fact that only a finite number of cells is available dueto energy conservation and compactness of the energy shell.The saturation value of entropy is expected to scale with δ likethe logarithm of that number, i.e., in the leading order, likeC ln(1/δ) with C = 1 in both regular and chaotic cases. This

036208-3

GIULIO CASATI, ITALO GUARNERI, AND JOSE RESLEN PHYSICAL REVIEW E 85, 036208 (2012)

is consistent with numerically observed values C � 0.86 andC = 0.93 for the regular and chaotic cases, respectively.

The average growth of entropy at large times yet still farfrom saturation may be understood similarly. The initial cellintersects different three-dimensional manifolds of constantenergy. In the chaotic case, evolution stretches the three-dimensional section at energy E along the unstable directionand contracts it along a stable one, at an exponential rate givenby the maximal Lyapunov exponent λE . Hence the numberof cells that are populated by the projection of this sectiongrows roughly like exp(λEt). The number of cells populatedby the projections of all sections grows like a continuoussuperposition of exponentials with rates λE varying in asmall interval of size proportional

√δ, so in the leading order

entropy grows proportionally to t times the average Lyapunovexponent, independently of δ. Equating this to the saturationvalue the order of magnitude estimate tD ∼ λ−1ln(1/δ) isobtained, consistently with numerical results.1

In the regular case, the initial four-dimensional cell inter-sects several invariant two-tori. Each two-dimensional sectionglides uniformly in time along the torus where it belongs,carried by the linear flow that is associated with the frequencyvector of that torus; so its projection onto the single-particlephase space is just deformed quasiperiodically in time, with noaverage growth of the number of populated two-dimensionalcells. The number of two-dimensional cells populated byall such projections nevertheless increases because differentsections move with different frequency vectors due to nonlin-earity; so, although all sections were originally piled up overthe same two-dimensional cell, their projections slide apartlinearly in time, proportionally t

√δ (because the mismatch of

frequency vectors is on the order of the linear size of the cell).As sections are a two-parameter family, the two-dimensionalarea covered by such projections increases quadratically intime: ∼ constδt2 in the leading order. Hence the number ofoccupied cells of area δ increases ∼ constt2 and entropy grows,in the leading order, like 2 ln(t), independently of δ. At timetD this should be on the order of magnitude of the saturationvalue ln(1/δ), whence the observed scaling of the saturationtime tD ∼ 1/

√δ follows.

As shown in Fig. 4, the behavior of quantum entanglementis quite similar. The same scaling with h of deviation timesand saturation values emerges as the classical limiting curvesare approached, in both the integrable and chaotic cases [8].

The oscillatory pattern that is seen in all figures, superim-posed on the just-discussed average behavior, is qualitativelyexplained as follows. In the regular case, the projection ofeach invariant torus on the (q1,p1) plane has caustics at itsboundaries because there the tangent plane of the torus isvertical. As sections of the initial ensemble slide around theirtori, the areas of their projections change quasiperiodically

1In some chaotic systems the exponential growth of the numberof cells may be followed by a diffusive growth proportional to

√t .

Diffusion may continue forever if, unlike the present case, motionis not confined to compact energy surfaces and no saturation occurs.Also some conservative systems may display a diffusive window afterthe Lyapunov stage and prior to saturation. However, in the presentcase no such evidence could be detected.

-4

-3

-2

-1

0

1

2

3

4

-3 -2 -1 0 1 2 3 4

p1

q1

FIG. 5. Projected (p1,q1) ensemble for the chaotic case with E =7 and m2/m1 = 0.54, at time t = 200.

in time and in particular become quite small anytime theycome close to caustics. So long as different sections of theinitial ensemble have not yet strongly dephased from oneanother, this produces a regular sequence of minima in thenumber of occupied cells. In the chaotic case, observation ofthe projected ensembles in the (q1,p1) plane, at times thatcorrespond to minima in the regular sequence of oscillations,reveals that phase-space points are then concentrated nearcausticlike lines (Fig. 5). In our interpretation, these betraythe presence of remnants of a broken torus. The projectedensemble appears to quasiperiodically concentrate near theselines until the four-dimensional ensemble eventually leaksthrough the broken torus to invade most of the energy shell.

A final remark is in order about the classical and quantumreduced entropies. In the quantum case, the von Neumannentropies of the two reduced states that are obtained from apure state of a bipartite system, by tracing out either subsystem,have the same value. In contrast, it is fairly obvious that,in the classical case, the coarse-grained entropies (2) of thetwo ensembles that are obtained from a given ensemble byprojection onto the phase spaces of the subsystems need notbe equal in general. However, they are found to approach eachother in the limit δ → 0.

IV. CONCLUSION

In conclusion, we have shown that there is a close similaritybetween the laws that rule the dynamical generation ofquantum and classical reduced entropies. The classical lawsadmit of simple heuristic explanations, which provide cluesto understanding the quantum ones. Such arguments predictdifferences between the ways entanglement is generated inclassically regular and chaotic regimes, which should beobservable for generic quantum multipartite systems in thequasiclassical regime.

ACKNOWLEDGMENTS

G.C. and I.G. acknowledge support by Ministerodell’Istruzione dell’Universita e della Ricerca, PRIN 2008and by Regione Lombardia. J.R. acknowledges support fromGrant Projects No. R-144-000-276-112 and No. R-710-000-016-271.

036208-4

CLASSICAL DYNAMICS OF QUANTUM ENTANGLEMENT PHYSICAL REVIEW E 85, 036208 (2012)

[1] K. Furuya, M. C. Nemes, and G. Q. Pellegrino, Phys. Rev. Lett.80, 5524 (1998).

[2] Shi-hui Zhang and Quan-lin Jie, Phys. Rev. A 77, 012312 (2008);Y. S. Weinstein and L. Viola, Europhys. Lett. 76, 746 (2006);S. Ghose and B. C. Sanders, Phys. Rev. A 70, 062315 (2004);J. N. Bandyopadhyay and A. Lakshminarayan, Phys. Rev. E60, 016201 (2004); N. N. Chung and L. Y. Chew, ibid. 80,016204 (2009); Xiaoguang Wang, S. Ghose, B. C. Sanders, andBambi Hu, ibid. 70, 016217 (2004); R. M. Angelo, S. A. Vitiello,M. A. M. de Aguiar, and K. Furuya, Physica A 338, 458 (2004);S. Ghose, R. Stock, P. Jessen, R. Lal, and A. Silberfarb, Phys.Rev. A 78, 042318 (2008).

[3] M. Novaes, Ann. Phys. (N.Y.) 318, 308 (2005); M. Lombardiand A. Matzkin, Phys. Rev. E 83, 016207 (2011).

[4] J. Gong and P. Brumer, Phys. Rev. A 68, 022101 (2003).[5] P. Jacquod, Phys. Rev. Lett. 92, 150403 (2004).[6] M. Znidaric and T. Prosen, Phys. Rev. A 71, 032103 (2005).[7] S. Chaudhury, A. Smith, B. E. Anderson, S. Ghose, and P. S.

Jessen, Nature (London) 461, 769 (2009).[8] Exponential growth of entanglement at a rate given by the sum

of positive Lyapunov exponents was pointed out by P. A. Millerand S. Sarkar, Phys. Rev. E 60, 1542 (1999); see, however, acritical discussion by H. Fusijaki, T. Miyadera, and A. Tanaka,ibid. 67, 066201 (2003).

[9] R. M. Angelo and K. Furuya, Phys. Rev. A 71, 042321 (2005).[10] A. Matzkin, Phys. Rev. A 84, 022111 (2011).[11] G. Casati and J. Ford, Phys. Rev. A 12, 1702 (1975).[12] T. Prosen, Phys. Rev. E 83, 031124 (2011).

036208-5