Classical evolution of quantum systems

For studying quantum phenomena in time, it is vital to have a profound understanding of the classical dynamics. For this reason, we derive equations of motions describing the classical propagation of a quantum system. A comparison of this classical evolution with the actual temporal behavior enables us to identify quantum effects of the evolution itself and distinguish them from static quantum features and quantum phenomena for a single point in time. As applications of our universal technique, we analyze nonlinear processes in quantum optics, semi-classical models, and the multipartite entanglement dynamics of macroscopic ensembles.

For studying quantum phenomena in time, it is vital to have a profound understanding of the classical dynamics. For this reason, we derive equations of motions describing the classical propagation of a quantum system. A comparison of this classical evolution with the actual temporal behavior enables us to identify quantum effects of the evolution itself and distinguish them from static quantum features and quantum phenomena for a single point in time. As applications of our universal technique, we analyze nonlinear processes in quantum optics, semi-classical models, and the multipartite entanglement dynamics of macroscopic ensembles.

I. INTRODUCTION
Quantum physics is one cornerstone of modern science. The laws governing the quantum domain led to new and surprising insights into the principles of nature. For example, interference patterns observed in a double-slit experiment are inconsistent with the propagation of classical particles. Yet, this and other phenomena are easily explained on the basis of the Schrödinger equation [1,2]. In contrast, other effects are only superficially related to quantum properties, such as interferences resulting from classical optics [3]. Therefore, to decide whether a process is genuinely quantum, the classical counterpart to a quantum evolution is required.
Moreover, quantum features are essential for the development of technologies [4,5]. New applications exploit the resources provided by quantum systems to perform tasks not achievable by classical means, e.g., quantum teleportation [6] and dense coding [7]. Other protocols rely on quantum properties to improve the communication security, such as quantum key distribution [8,9]. Thus, the classical limitations of the underlying processes have to be known in order to assess the benefit of utilizing quantum phenomena.
An early approach to certify that quantum evolutions are incompatible with classical dynamics is the seminal Leggett-Garg inequality [10]; see Ref. [11] for a recent review. Similarly to the concept of nonlocallity, the violation of this and related inequalities uncovers temporal quantum correlations [12,13]. The question of causality in quantum physics is a closely related concept, which addresses another aspect of the quantumness of a process [14][15][16][17]. Furthermore, the non-commutative features in the mathematical formulation of quantum physics lead to the quantum effect of time-ordering [18,19]. For example, the resulting nonclassical propagation of light yields multi-time quantum correlations for nonlinear optical processes [20][21][22].
Alongside with temporal correlations, time-independent quantum phenomena are frequently studied. For instance, the concept of quantum coherence characterizes the quantum nature in terms of quantum superpositions [23][24][25][26][27]; see Ref. [28] for a review. For different physical scenarios, different families of classical states can be considered [29][30][31], and quantum coherence is defined as the inability to describe the state of a * jan.sperling@physics.ox.ac.uk system in terms of such classical reference states and classical statistics. Moreover, this notion includes quantum entanglement between different degrees of freedom [32][33][34], which is arguably the most prominent and versatile quantum phenomena and resource for quantum protocols [35].
Furthermore, a process maps an initial state to a state at a later time. In this input-output formalism, the quantumness is characterized via the potential of the underlying map to create coherences from an initially classical state [36][37][38]. This renders it possible, for example, to characterize quantum correlations between a system and a bath [39] after a given time and to study the temporal propagation of entanglement [40][41][42][43]. Still, such methods are limited as they do not address the quantum nature of the dynamics itself. Rather, the generation and degradation of quantum coherence are quantified at each final time separately.
To overcome such a restriction, we recently introduced a method to analyze the entanglement dynamics [44]. This was achieved by constraining the quantum trajectories to nonentangled (i.e., separable) ones. Yet, this approach does not encompass other forms of temporal quantum phenomena based on general notions of quantum coherence.
In this contribution, we formulate a method to determine the classical evolution of a quantum system. Our technique is based on Hamilton's equations restricted to a domain of classical states under study. This results in equations of motion whose solution can be compared with the quantum propagation to directly uncover the quantum properties of the dynamics. The practicability of our approach is demonstrated with several applications assessing the time-dependent quantum behavior in different systems. Moreover, our technique bridges the gap between classical and quantum systems by consistently modeling semi-classical dynamics. Thus, we devise a versatile method which enables us to predict the classical evolution of quantum systems.

II. EQUATIONS OF MOTION
In mathematical physics, a general technique to formulate equations of motion is based on Hamilton's equations, Here, H is the Hamiltonian, q is a generalized position, and p is its conjugate momentum. For the purpose of this work, it arXiv:1802.02807v1 [quant-ph] 8 Feb 2018 is also useful to define a complex degree of freedom, with a positive constant γ.
we can now combine the equations for x and p. In general, for a system consisting of N degrees of freedom, ξ = (ξ 1 , . . . , ξ N ) T , we readily obtain Moreover, Hamilton's equation apply to scenarios beyond classical physics. For example, the Hamiltonian H relates to the Hamilton operatorĤ as We emphasize that the distinction between Hamiltonian H and the Hamilton operatorĤ is important for our work. Further, we identify ξ = |ψ and γ =h −1 . Then, the Schrödinger equation is directly retrieved from Eq. (3), using ∂ ψ|Ĥ|ψ /∂ ψ| =Ĥ|ψ [1].
For analyzing the evolution with respect to the theory of quantum coherence, we may consider a family of classical states, |ψ ζ , parametrized by an, in general, complex vector ζ . Different types of states can be selected, connecting to diverse classical scenarios [28]; specific examples are studied later. Further, the concept of a classical evolution requires that a classical initial state remains classical for all times. The actual quantum evolution in terms of Eq. (5) can achieve this in some approximation only [45]. Thus, we have to require that the state |ψ ζ (0) propagates as |ψ ζ (t) for all times t ≥ 0. In comparison, one typically assumes and input-output formalism in which an initially classical state is mapped onto a classical state after a fixed time, independent of the possibility to have nonclassical states for intermediate times. By contrast, we demand classicality over the entire duration of the process to define a classical evolution. Now, we formulate the necessary equation of motion for a dynamics constrained to the given family of classical states. Applying the concepts discussed earlier [see Eqs. (3) and (4)], we can directly describe the desired evolution via This simple, yet useful approach is capable of determining the classical evolution of any quantum system. On the one hand, the quantum nature of the underlying system is represented by the Hamiltonian defined in terms of the Hamilton operator, On the other hand, the classicality of the dynamics is ensured by the fact that that the solutions ζ (t) of Eq. (6) yield a classical state |ψ ζ (t) for any t. In Appendix 1, we provide additional information and relations, including the treatment of mixed states.

III. NONLINEAR QUANTUM OPTICS
Let us apply our approach to show its general function. A coherent state |α behaves similarly to a classical harmonic oscillator [46]. A monochromatic radiation mode, represented by the annihilation operatorâ, is an example of such a system, where the coherent amplitude α relates to the classical electromagnetic field [47][48][49]. Here, we focus on a single mode as the generalization to multimode light is straightforward.
To demonstrate this, let us study the propagation of light in a Kerr medium, given by the Hamilton operator [51,52] H =hωâ †â +h κ 2â †2â2 , where κ is the coupling constant. Consequently, the classical Hamiltonian reads Applying our method and the Heisenberg equation, we get respectively, where χ = κ/ω quantifies the intensity dependent contribution to the refractive index. Moreover, the solution of the classical dynamics reads Similarly, the quantum field propagates according tô wheren =â †â is the photon-number operator.
FIG. 1. Phase-space representation at t = π/κ, rotated by e iωt to counter the free evolution. The red (blue) parts depict the positive (negative) contributions in the Wigner function. The first and second letter of a panel's label indicates the type of initial state and dynamics, respectively; both can be either classical (C) or quantum (Q). The curves depict the trajectory of the mean value, â e iωt . Bullet points mark the times t = 0 (right) and t = π/(2κ) (left).
The classical position in phase space is given by Eq. (13). By contrast, for the initial state |α(0) , the quantum solution (14) yields using e xn |α = |e x α and α|α = e −(|α| 2 +|α | 2 )/2 e α * α . For small times, implying e −iκt − 1 ≈ −iκt, the classical position in phase space coincides with the quantum one. However, this is not true for all time. Moreover, let us mention that we can apply the classical solution even to an initially nonclassical state. This enables us to distinguish quantum effects resulting from a nonclassical dynamics from nonclassicality which was already present at the initial time, t = 0. In Fig. 1, we compare the classical (top row) and quantum (bottom row) dynamics for states which can be either classical (left column) or nonclassical (right column) for t = 0. The trajectories (solid curves) significantly differ from one another, cf. Eqs. (13) and (15). Thus, the quantum evolution exhibits characteristics (here, additional loops) not present in their classical counterparts.
In addition, the Wigner function of the evolved states for tκ = π/2 is depicted in Fig. 1. Note that negativities of the Wigner function (red areas) are a sufficient indication of the nonclassicality in terms of the Glauber-Sudarshan function [53,54], which is the actual figure of merit but highly singular in the considered scenario [38]. The classical evolution preserves-up to a π rotation-the features of the initial state, cf. panels CC for |α(0) = 3 and QC for the superposition |ψ(0) ∝ |3e −iπ/4 − |3e iπ/4 . In case of the quantum evolution, nonclassicality is generated for the initially classical state (panel CQ) and additional, nonclassical interference patterns occur for the initially nonclassical state (panel QQ). Thus, we can discriminate quantum effects originating from the dynamics from the initial nonclassicality.

IV. SEMI-CLASSICAL MODELS
Another application of our technique are semi-classical models. Such hybrid systems are important for coupling a classical degree of freedom to a quantum-mechanical one, e.g., to describe measurements with macroscopic devices. However, the quantum evolution typically yields unphysical results as quantum properties can "leak" into the classical domain due to the interaction; see Ref. [55] for a recent study. Our approach can overcome such a deficiency.
As a proof of concept, we study a light field interacting with a two-level atom using the Jaynes-Cummings model [56], were |e and |g are the excited and ground state of the atom and κ is a positive coupling constant. Here, the resonant scenario is considered in which the energy difference between the atomic states is proportional to the frequency of the field,hω. To treat the system in a semi-classical framework, we say that a classical coherent state represents the optical mode. The state of the atom, |φ , is treated in the quantum domain. Thus, we have to restrict to states |ψ ζ = |α ⊗ |φ together with a semi-classical parametrization ζ = (α, |φ ). Now, the Hamiltonian [Eq. (7)] reads Applying the equation of motion (6) to this Hamiltonian for each component of the pair ζ = (α, |φ ), we find On the one hand, Eq. (18) is a classical wave equation which includes a source term proportional to the scalar, timedependent function φ (t)|g e|φ (t) . On the other hand, the dynamics of the quantum atom is described in terms of the Schrödinger equation (19), including an effective interaction proportional to the operator iα(t)|e g| − iα(t) * |g e|. Let us stress that the classical part interacts with the quantum part. Still, the solutions α(t) are classical field amplitudes for all times t.
where |n denotes the n-photon Fock state, which is a nonclassical for n > 0 [50]. For sin(κt) = 0, this solution is additionally entangled. Thus, the joint quantum system evolves into a state in which the optical mode cannot be considered in a classical frame, despite being a classical for κt = 0 mod π. By contrast, our semi-classical equations of motion, (18) and (19), yield a proper semi-classical state |α(t) ⊗ |φ (t) . Therein, the propagated classical field amplitude is and the quantum state of the atom reads Here, ϑ (x) denotes a Jacobi elliptic function defined by dϑ /dx = (1 + sin 2 ϑ ) 1/2 and ϑ (0) = 0, and it can be approximated with ϑ (x) ≈ µx (µ = 1.198 140 . . .). Thus, our method enables us to formulate consistent semi-classical models for studying the joint evolution of a classical and quantum system and their mutual influence for all times.

V. MULTIPARTITE ENTANGLEMENT
As a final application, let us study quantum correlations in multipartite systems. An N-partite quantum system can be decomposed into K parties, each including N j individual subsystems ( j = 1, . . . , K and N 1 + · · · + N K = N). For instance, the composite system {A, B,C, D} can be decomposed into a bipartition {A, B}:{C, D}, i.e., N = 4, K = 2, and N 1 = N 2 = 2. With increasing number of parties, the number of such possible decompositions increases exponentially, which already makes the time-independent characterization of multipartite correlations a sophisticated problem [35].
Here, as an example, we study the case in which the initial state is the thermal state ofĤ without interaction, cf. Eq. (26) for κ = 0. This yields a mixed N-fold tensor-product state. Further, the observable under study is the mean momentum, For t ≥ 0 and a non-vanishing coupling, κ > 0, we obtain the expected separable evolution, which results in the following propagation of the variance ofP: for a balanced partitioning (i.e., N 1 = · · · = N K = N/K), the temperature T , the Boltzmann constant k B , and the parameter In Fig. 2, we show the evolution of V P (t)/V P (0) for different K partitions of a macroscopic ensemble consisting of N = 10! ≈ 3.6 × 10 6 oscillators; see Refs. [59,60] for a time-independent analysis. The depicted curves apply to all initial temperatures. The variances are not below the initial one, V P (t)/V P (0) ≥ 1, which eliminates the possibility of squeezing as a single-time quantum effect. Here, the actual entanglement-generating evolution (red curve), resembling the trivial partition with K = 1, has the smallest initial increment of the variance. This is clearly different when compared to the bipartition, K = 2, separating the ensemble into two parts with about 1.8 × 10 6 oscillators each. In fact, with increasing separation K (increasingly lighter blue curves), the initial increment of V P (t)/V P (0) increases as well, which enables the discrimination of those cases. More generally, the difference between the K-separable dynamics and the actual propagation becomes increasingly pronounced in both the amplitude and period of the oscillation. Therefore, our general method enables us to study the classically correlated dynamics (for arbitrary degrees of separation, K) to identify measurable characteristics genuine to the quantum-correlated dynamics.

VI. CONCLUSION
In summary, we devised a technique to describe the classical evolution of quantum systems. This allows us to compare the quantum dynamics with the classically predicted behavior for identifying nonclassical properties in time. Based on Hamilton's equations, we derived the equations of motions for the parameters which characterize the set of classical states under study. Consequently, the solutions are confined to those classical states for all times. Whenever the actual quantum dynamics is incompatible with the classical behavior, the quantum nature of the evolution is verified. To demonstrate the broad applicability, different applications were studied.
Firstly, we analyzed the notion of nonclassicality in terms of coherent states to study processes in nonlinear quantum optics. We demonstrated that our classical equations of motion and the quantum-physical Heisenberg equations formally share the same structure. However, their solutions exhibit distinguishing features. This was shown by analyzing the propagation in a Kerr medium. For the classical and quantum dynamics, we compared states which are initially classical or nonclassical to separate initial nonclassicalities from quantum features due to the quantum evolution.
Secondly, we considered semi-classical models in which one part is a quantum system and the other one is classical. This is important, for instance, for studying quantum measurements processes with a macroscopic device in a consistent manner. Our method renders it possible to perform such a task as demonstrated for the example of a classical light field coupled to a two-level quantum atom. We could show that, despite of interactions with quantum subsystem, the field remains classical for all times, which cannot be achieved when applying standard approaches.
Finally, we studied the entanglement dynamics of a system consisting of a macroscopic number of quantum-mechanical oscillators. In this scenario, the classical states are identified with separable ones. The classical dynamics resulting from our method is described in terms of Schrödinger-type equations for each subsystem separately. This is consistent with our recent approach [44] and underlines the the general potential of our method introduced here. As we have shown for different forms of separability (or, conversely, partial entanglement), the entanglement-generating evolution can be clearly distinguished from the classical one, even for initially fully separable mixed states.
In conclusion, we developed a general and easily accessible approach to certify quantum properties in time. Our technique is not restricted to a specific type of classical states and applies to arbitrary time scales and interaction regimes. Thus, we believe that our method provides a powerful tool and versatile starting point for studies of temporal quantum phenomena.

ACKNOWLEDGMENTS
The project leading to this application has received funding from the European Union's Horizon 2020 research and innovation program under Grant Agreement No. 665148 (QCUM-bER).

Appendix
Here, we provide additional details. A brief review about relevant concepts in theoretical mechanics is given in Sec. 1. The exact solutions for the entanglement example are computed in Sec. 2.

Lagrangian mechanics
The classical and quantum dynamics of particles and fields can be formulated in terms of Lagrangian mechanics. Let us recall some specific properties of this formalism because of their importance for our work. More detailed and general introductions to the field quantization and calculus of variations may be found elsewhere [61,62]. For a closed system, the Lagrangian L = L (q,q) for N degrees of freedom is a function of generalized coordinates q = (q 1 , . . . , q N ) and their time derivatives. We use the notatioṅ x = dx/dt for a function x = x(t). The action, S = dt L , determines the evolution of the system. Namely, the least action principle-a vanishing variation of the action, δ S = 0yields the Euler-Lagrange equations, d dt The conjugate momenta, p = (p 1 , . . . , p N ) = ∂ L /∂q, allow one the define the Hamiltonian, H = p ·q − L . This yields the equivalent Hamilton's equations,q = ∂ H /∂ q anḋ p = −∂ H /∂ p. One consequence of Lagrangian mechanics is that two variations can be related, δ dtq · p = −δ dt q ·ṗ. This can be used to formulate an equivalent relation between the Hamiltonian and the Lagrangian without changing the equations of motion of the system, e.g., the symmetric form L = (q · p − q ·ṗ)/2 − H . The latter is useful when considering complex representation combining position and momentum, where γ is a positive constant. Namely, the Lagrangian reads L = (ξ † ξ − ξ †ξ )/(2iγ) − H . The complex description is specifically useful for field theories, such as electromagnetism where electric and magnetic field components can be combined into one complex-valued electric field. Further, the equations of motions for q and p can be combined, One application can be found in the derivation of the Schrödinger equation (see Ref. [1]), where ξ = |ψ , γ = 1/h, and where H = ψ|Ĥ|ψ andĤ being the Hamilton operator. Then, the complex form of the Euler-Lagrange equation yields the desired equation of motion, ih|ψ =Ĥ|ψ . It is also worth pointing out that a contribution proportional to1 in the Hamilton operator leads to a global phase which can be ignored. Statistical ensembles in a closed system can be studied with the Lagrange formalism in a similar manner. For instance, the density operatorρ = ∑ n p n |ψ n ψ n | can be represented by the partial trace tr aux |ψ ψ| =ρ, using the pure state |ψ = ∑ n √ p n |ψ n ⊗ |n aux and an auxiliary system expanded in the orthonormal basis {|n aux } n . As the initial system (i.e, without the auxiliary part) is closed, the Hamiltonian takes the form H = ψ|Ĥ ⊗1 aux |ψ . Eventually, one observes the separation of the components for each n in the Lagrangian, L = ∑ n p n ih 2 ψ n |ψ n − ih 2 ψ n |ψ n − ψ n |Ĥ|ψ n . (A. 5) This means that for each component |ψ n , the equations of motions can be solved separately, and the ensemble average can be performed using the time-independent probabilities {p n } n . Note that this form can be directly extended to continuous distributions, ∑ n p n |ψ n ψ n | → dP(ψ)|ψ ψ|, and quasiprobility distributions, where P 0.

Separable dynamics
The Hamilton operator of the N-partite system under study readŝ where [x k ,p l ] = ihδ k,l with δ k,l denoting the Kronecker symbol. In particular, we can identify the specific form of canonical operators in position representation as Following the approach in Ref. [60], it is convenient to consider the dynamics in "natural" units of the system, leading to the following quantities: a rescaled time, τ = ωt; a coupling constant, R = κ/(mω 2 ); transformed canonical coordinates, ξ = (mω/h) 1/2 x and π = −i∂ /∂ ξ ; as well as a rescaled Hamilton operator, where [ξ k ,π l ] = iδ k,l and n = (1, . . . , 1) T . For R = 0, we obtain the non-interacting case. We aim to study the time-dependent behavior of secondorder correlations. Because of the non-commuting property of position and momentum, let us briefly mention the noncommuting components of the covariance matrix C, ξ kπl + π lξk /2 − ξ k π l . We haveξ kπl = (iδ k,l +ξ kπl +π lξk )/2 andπ lξk = (−iδ k,l +ξ kπl +π lξk )/2, which yields where E = diag(1, . . . , 1) is the N ×N identity. Here it is worth emphasizing that ( q T r T ) 0 E −E 0 q r = 0. Consequently, the dynamics is described by the following equations of motion: dξ k /dτ = −i[ξ k ,η] and dπ k /dτ = −i[π k ,η]. In addition, we have Because it is also useful for the separable case, let us assume for the time being that the Hamilton operator takes the more general formη with an N-dimensional vector v and a real-valued, symmetric, and positive definite N × N matrix G. Then we get with the symplectic matrix J = 0 E −E 0 = −J T . Furthermore, using the definition ∆ŷ =ŷ − ŷ for anyŷ, we get d dτ which is independent of v. For the example v = 0, the evolution is described in terms of the symplectic transformation S(τ) = cos(G 1/2 τ) G −1/2 sin(G 1/2 τ) −G 1/2 sin(G 1/2 τ) cos(G 1/2 τ) , (A.12) which satisfies S(0) = E 0 0 E and dS(τ)/dτ = J G 0 0 E S(τ) and applies to an initial covariance matrix C(0) as C(τ) = S(τ)C(0)S(τ) T . As a result of the spectral decomposition of the specific case G = µE + ν n n T , we get 13) for the N-partite system under study. For studying separable dynamics (see also Ref. [44] and its Supplemental Material for an detailed introduction), let us decompose the N-partite system into K subsystems, each consisting of N j individual oscillators, for j = 1, . . . , K. In this scenario it is useful to decompose, e.g., the generalized position as ξ = ( ξ T 1 , . . . , ξ T K ) T and the separable state as |ψ = |ψ 1 ⊗ · · · ⊗ |ψ K . The partially reduced Hamilton operator for the jth subsystem takes the form η ψ 1 ,...,ψ j−1 ,ψ j+1 ,...,ψ K = 1 2 π T j π j + where contributions proportional to1 j have been ignored. Applying our previous considerations, we find d dτ which is identical to the entangled case. However, the covariance of the jth component evolves in the separable scenario differently-namely, according to 16) with G j = (1 + RN)E j − R n j n T j . Let us emphasize that we have ∆ŷ ⊗ ∆ẑ ψ⊗φ = ∆ŷ ψ ∆ẑ φ = 0 · 0 for allŷ andẑ and all times. Thus, the evolution of the covariance in the separable case is fully described by the evolution of the jth block of the corresponding subsystem separately, which is given by Finally, the specific initial condition is the thermal state of the non-interacting system, R = 0. In this case, the covariance decomposes in one-partite blocks with the Hamiltonianη j = (π 2 j +ξ 2 j )/2. The partition function reads Z j = tr(e −βη j ) = e −β /2 /(1 − e −β ), where β =hω/(k B T ) in our natural units. Further, this yields the desired moments, ξ j = π j = ξ jπ j +π jξ j /2 = 0 and ξ 2 j = π 2 j = 1/2 + e −β /(1 − e −β ). Thus, the jth block of the propagated covariance matrix, C j (τ), reads C j (τ) = 1 2 + e −β 1 − e −β cos 2 (G 1/2 j τ) For the variance of the momentum operator parallel to n, i.e., Π = n T π/N, we find using N 1 = · · · = N K = N/K and r K = RN(1 − 1/K).