Particle Trajectories for Quantum Maps

We study the trajectories of a semiclassical quantum particle under repeated indirect measurement by Kraus operators, in the setting of the quantized torus. In between measurements, the system evolves via either Hamiltonian propagators or metaplectic operators. We show in both cases the convergence in total variation of the quantum trajectory to its corresponding classical trajectory, as defined by propagation of a semiclassical defect measure. This convergence holds up to the Ehrenfest time of the classical system, which is larger when the system is less chaotic. In addition, we present numerical simulations of these effects. In proving this result, we provide a characterization of a type of semi-classical defect measure we call uniform defect measures. We also prove derivative estimates of a function composed with a flow on the torus.

Each grey curve represents the trajectory of a quantum particle (under the evolution procedure discussed in this paper) on a torus with identical initial conditions and a fixed Hamiltonian.A single trajectory is plotted in red and the corresponding classical trajectory is plotted as a blue dotted line.Up until time 0.4, most of the positions of the quantum particles are within 0.1 of the classical trajectory, but shortly after become decoherent.More numerics are provided in §4.This simulation is discussed at the end of §1.1 and §4.

Introduction
The framework of "quantum trajectories" has gotten increasing attention in recent years.Instead of treating measurement as a one-time event, researchers consider many measurements obtained over time and study their effects on the particle being measured, described mathematically via quantum instruments.The sequence of measured values is called the quantum trajectory, and is a random variable worthy of study.
In this paper, we take quantum trajectories into the setting of the semiclassical quantized torus, as studied in Bouzouina-De Biévre [BD96], Schenck [Sch09], Dyatlov-Jezequel [DJ21], and others.The spaces involved are all finite-dimensional, which makes numerical simulations easier than in PDE-based models.Our result is inspired by the recent work of Benoist, Fraas, and Fröhlich [BFF22], who prove in the PDE setting that the trajectory of a repeatedly observed quantum particle in the semiclassical limit approaches the natural notion of a classical trajectory.See Figure 1 for a simulation of this quantum-classical correspondence of particle trajectories on the quantized torus and §4 for more numerics.We strengthen the result in [BFF22] in our setting by allowing the time to depend on the semiclassical parameter, and we relate the rate of convergence to the "chaotic behavior" of the underlying classical system.This matches the intuition that a less chaotic transformation should be well-approximated for longer by its classical counterpart.Additionally, we illustrate the results numerically to demonstrate in the chaotic case that our time restriction is essentially optimal.
1.1.Model.In our setting, a quantum particle is modeled by a density operator (see Definition 2.13) ρ N on the quantized 2d−dimensional torus H N (see (2.8)) where N ∈ N is proportional to the reciprocal of the semi-classical parameter h.This particle is "indirectly measured" by conjugating by Kraus operators (see Definition 2.14) Op N (f q ) where {f q } q∈Ω are fixed elements of C ∞ (T 2d ) with uniformly bounded derivatives (such that Op N (f q ) provide a resolution of the identity), Ω is a compact metric space with finite Borel measure ν, and Op N is the quantization of functions on the torus as described in §2.4.The measurement procedure is further described in §2.5.
After observing the particle, we evolve it for one unit of time by conjugating ρ N by a unitary operator U , which is either a quantization of a fixed symplectic matrix M (defined in §2.2) or Schrödinger evolution e − i h P for a fixed Hamiltonian P .If we repeat this process n times, we get a resulting measure on the space of trajectories.Suppose F ⊂ Ω n+1 , then the probability of obtaining the set of trajectories F is given by: where Φ * N,q i (ρ) := U Op N (f q i )ρ N Op N (f q i ) * U * .Therefore, we have a "quantum probability measure" given by We now describe the trajectory of a corresponding classical particle.Let µ be the defect measure (defined in Definition 2.5) of ρ N .We interpret µ as the probability distribution of the initial position and momentum of a classical particle.We apply an "approximate measurement," which computes an observable quantity q to be in E 0 ⊆ Ω with probability E 0 |f q 0 (ζ)| 2 dν(q 0 ).The particle is then allowed to classically evolve for one unit of time by a flow ϕ t on T 2d which is either multiplication by the symplectic matrix M or else the nonlinear Hamiltonian evolution exp(H p ).We again repeat the process of alternating measurement and evolution, though we must remark that unlike in the quantum case, the measurement does not affect the location of the particle.Conditioned on initially having position/momentum ζ, we have the probability of classically measuring (q 0 , . . .q n ) in a measurable set F ⊂ Ω n+1 to be One can easily check that (q n ) n is a process of independent random variables (although not identically distributed).As the initial value of ζ was taken to be randomly chosen from the distribution µ, we take the classical probability of measuring (q 0 , . . ., q n ) ∈ F to be We therefore have a "classical probability measure" given by (1.2) Note that the functions f q are allowed to overlap so that P (n) µ is a probability of trajectories in Ω n+1 while ϕ t is the unique trajectory on the torus.
The goal of this paper is to describe in what sense, and at what quantitative rate, the measure P (n) N,ρ approaches P (n) as N grows large.We can state our main result.
In proving this, we show uniform convergence of P (n) N,ρ to P (n) µ and hence µ in total variation.For a stronger version of Theorem 1.1, see Theorem 1, in which we provide a quantitative rate of convergence in the case where the number of steps n approaches the Ehrenfest time (defined in 2.4).
Figure 1 numerically samples 60 quantum trajectories from such an evolution procedure (see the end of §4 for more details) for 100 time steps of 0.01 units of time.Observe that up until t = 0.4, most trajectories are within 0.1 of the classical trajectory 1 , while shortly after this time, they disagree.This paper's main result estimates how long they agree, and to what distribution they agree with.1.2.Previous Work.The understanding that a measurement affects a quantum state goes back to the early days of quantum mechanics.In his paper on the uncertainty principle, Heisenberg [Hei27] (see [Hei83] for an English translation) discussed the socalled "collapse" of the wave function, which was later rephrased in mathematical terms by von Neumann [Von13].When a system is only partially or indirectly measured, the necessary framework is that of positive operator-valued measures, which were first studied by Naimark [Neu43].The analog of wave-function collapse was introduced by Davies and Lewis [DL70] with the development of quantum instruments.
Further research focused on multiple measurements as a probabilistic process, which could be modeled in either discrete or continuous time.The first analysis of continuoustime measurements was due to Davies [Dav69], and required extensive machinery from stochastic calculus (See [BG09] or [Hol03] for a brief introduction).Meanwhile, the discrete-time model consisting of repeated applications of a quantum instrument, while requiring substantially fewer technicalities, has nonetheless demonstrated a rich variety of behavior and is a subject of current study.Kümmerer and Maassen [KM03] demonstrated ergodicity of a repeated measurement process and showed that the states approach a pure state provided the measurement operators do not have a "dark" subspace [KM06].Ballesteros et al. [Bal+18] show that when the dynamic is trivial, the state of the system localizes in space.
If the state is allowed to evolve in-between measurements, the situation becomes more complex.The mathematical analysis of such a quantum trajectory, which approximates a classical one, was done by Ballesteros et al. [Bal+21] following extensive physical evidence of the phenomenon, see the exposition by Figari [FT13].That paper studied particles, initially in the set of normal states, evolving under a quadratic 1 Our result shows that before the Ehrenfest time, at each fixed time, the distribution of quantum positions is a Gaussian centered at the classical trajectory position with variance given by the precision of the indirect measurement.
Hamiltonian.The motivation for our work came from the study of the semiclassical case under a more general Hamiltonian by Benoist, Fraas, and Fröhlich in [BFF22].We adapt their framework to the setting of quantum maps.This has two advantages: we can improve on the number of measurements exhibiting classical/quantum correspondence and we can provide (more easily than in the PDE setting) numerical simulations.In particular, the numerics indicate the accuracy of the dynamical bound on the number of measurements.

Background
In this section, we provide a background on the measurement procedure as well as the semiclassical analysis tools required to prove our main result.Throughout this paper, we consider h = (2πN ) −1 > 0 as the semiclassical parameter where N ∈ N.
2.1.Semiclassical Analysis on the Real Line.Physically, we are often concerned with how our quantum mechanical world approximately produces classical mechanics.This limit is achieved by taking a semiclassical parameter h ∈ (0, 1) to be very small, and the math involved is called semiclassical analysis.
Definition 2.1 (symbol class).For each δ ∈ [0, 1/2), the symbol class S δ is defined as Of most importance is the case where δ = 0, though technical steps of the proof of Theorem 1 will require us to consider general δ.Nevertheless, if δ is omitted it will be taken to be 0. Elements of S δ (called symbols) may depend on h, though we require that the constants C α,β in the definition be uniform in h.Real-valued symbols correspond to classical observables, which can be "quantized" to get quantum observables as follows.
Definition 2.2 (Weyl quantization).For δ ∈ [0, 1/2), the Weyl quantization of a symbol a ∈ S δ is an operator Op h a : S(R d ) → S ′ (R d ) (sometimes written a(x, hD)) given by A theorem of Calderón and Vaillancourt states that if δ ∈ [0, 1/2) and a ∈ S δ , then Op h (a) is a bounded operator on L 2 (R d ), with bound where C > 0 and K ∈ N depends only on d.(In fact, we remark that (2.1) even holds for the threshold scaling δ = 1 2 , though in this case, all derivatives of a in the sum will contribute equal orders of magnitude.) where in these formulas σ(x, ξ, y, η) := ⟨ξ, y⟩ − ⟨x, η⟩, and a where r has S δ seminorms at most O(h 1−2δ ), with constants depending only on finitely many derivatives of a and b.More generally, given a 1 , . . ., a n ∈ S δ , we have by an easy induction together with (2.1) that (2.3) Finally, we will need a theorem of Egorov, which gives a sense in which quantum operators evolve in time analogously to their classical counterparts.We provide the quantitative version due to Bouzouina and Robert [BR02], see also Zworski [Zwo12,Chapter 11.4].We first require the following definition.
Definition 2.3 (Lyapunov exponent).For any dynamical system given by a flow ϕ t on T * R d or T 2d (smooth in time), define the Lyapunov exponent Γ to be and analogously for iterative maps.
Given p ∈ S, let ϕ t be evolution by the classical Hamiltonian flow of p, induced by the Hamiltonian vector field X p = ∂ ξ p∂ x − ∂ x p∂ ξ .On the quantum side, a particle evolves by the Schrödinger equation hD t ψ = P ψ, so ψ(t) = e − i h tP ψ(0) and observables A evolve in the Heisenberg picture by A(t) = e i h tP Ae − i h tP .Egorov's theorem states that this evolved operator is equal, up to order h, to the quantization of the classically observed symbol of A, for t less than Ehrenfest time t E ∼ log 1 h .Specifically, we have that if Γ is the Lyapunov exponent of ϕ t and for some δ < 1 2 and ε > 0, then for a ∈ S, and r ∈ S δ has all seminorms of order O(h 2−3δ ).In cases where Γ = 0, such as a completely integrable classical system, some stronger results are available, see [BR02].
2.2.Fourier Integral Operators.As a generalization of the operators e − i h tP in Egorov's theorem, we briefly discuss Fourier integral operators.One may view the operator e − i h P as quantizing a transformation (or equivalently, a "change of variables") on phase space; namely the symplectomorphism on T * R d given by (x, ξ) → ϕ t (x, ξ).It is then natural to generalize the quantization to arbitrary symplectomorphisms.This abstraction will be rewarded when we work on the quantized torus in §2.4, on which the most natural classical transformations do not arise from a flow at all.Let ϕ : T * R d → T * R d be any symplectomorphism.We consider only unitary Fourier integral operators here, and only consider their action in the case δ = 0.For the purposes of this paper, define a Fourier integral operator quantizing ϕ in the following way.
Definition 2.4 (quantized flow).For ϕ : T * R d → T * R d , the Fourier integral operator quantizing ϕ is the operator holds for any a ∈ S with b = ϕ * a + O S (h).
In particular (using these definitions), Egorov's theorem is precisely the statement that e i h tP is a one-parameter family of Fourier integral operators, which quantize the classical Hamiltonian flow ϕ t .
Perhaps the most important Fourier integral operators are the metaplectic operators [Zwo12, Chapter 11.3], which is the special case when ϕ is a linear symplectomorphism.In this case for M ∈ Sp(2d, R), we write the corresponding metaplectic operator as M .The metaplectic operators are unitary and benefit from an exact version of (2.4), ie they satisfy for all a ∈ S, with no error.We may remark that while the metaplectic operators define a projective unitary representation of Sp(2d, R), they are in fact a unitary representation of its double cover, which is known as the metaplectic group.This is analogous to projective representations of rotation operators on spinor spaces and leads similarly to a non-canonical choice of phase for a metaplectic operator.As we will only encounter metaplectic operators when we are conjugating by them, this phase vanishes and we need not concern ourselves with it.
2.3.Semiclassical Defect Measures.We now discuss semiclassical defect measures, which give a quantitative answer to where in phase space a quantum particle "is" in the semiclassical limit as h → 0.
Definition 2.5 (defect measure).We say an h-dependent set of density operators we say it has defect measure µ if the density operator |ψ h ⟩ ⟨ψ h | has defect measure µ.
Having a semiclassical defect measure is a very special property of a sequence, corresponding to being "essentially classical" as h → 0. Nevertheless, such a limit necessarily exists along a subsequence of h's, as the following proposition states.
Proposition 2.6.Let ρ h be an h-dependent family of density operators on L 2 (R d ).
Then there is a nonnegative Radon measure µ on T * R d and a sequence h j → 0 such that The proof of Proposition 2.6 is essentially the same as [Zwo12, Theorem 5.2], which proves it in the special case of pure states.
In many natural cases, we may take the limit in (2.6) to be uniform in the symbol a, in which case we refer to µ as a uniform defect measure.For a precise characterization of aforesaid measures, see Proposition A.1.In particular, in this paper, we will work on the quantized torus H N (see §2.4) where all defect measures are uniform.
For semiclassical pseudodifferential operators on the real line, there are standard classes of states with known defect measures [BFF22, Proposition 3.3] which we restate in Proposition 2.7.For a full proof (in French) see [LP93, §3] and a proof for less general coherent states see [Zwo12,Chapter 5 Then ρ h has defect measure µ, given as follows: Here ĝ is Fourier transform of g.
In the last case of Proposition 2.7, the states ψ h are called generalized coherent states.
In fact, as we will be working with S δ classes and must control the error of our estimates, we will need a slightly stronger definition to accommodate h-dependent symbols.
We easily see that if µ is a (δ, θ)-defect measure, it is necessarily a (δ ′ , θ ′ )-defect measure for all δ ′ ⩽ δ and θ ′ ⩽ θ.To address a potential point of confusion, note that in (2.7) we are enlarging the space of symbols, a, used in (2.6).These symbols are allowed to depend on h with derivatives growing in h, provided they do not grow too fast, and their support may not be uniformly bounded in h.
We can now state a refinement of Proposition 2.7 using this new terminology.
Proposition 2.9.Let ρ h be as in Proposition 2.7.Then the defect measures µ given in Proposition 2.7 are uniform (δ, θ)-defect measures with δ, θ given as follows: The proof of Proposition 2.9 is straightforward and given in Appendix B.1.
2.4.The Quantized Torus.We now give an introduction to the quantized torus, and the pseudodifferential calculus on the Hilbert space of corresponding quantum states.
For each N ∈ N, define We may remark that H N consists precisely of the elements of S ′ (R d ) which are periodic in physical and Fourier space under the semiclassical Fourier transform F h , given by where h = (2πN ) −1 .For simplicity define and define a Hilbert space structure on H N such that (Q n ) is an orthonormal basis of H N .Though it is not necessary here, we may note that this is an elementary example of the general geometric Toeplitz quantization studied in [Del19].
We now give an overview of the pseudodifferential calculus on H N , which also may be found in [CZ10], [DJ21], [Sch09].Let a ∈ C ∞ (T 2d ).Intuitively, a corresponds to a classical observable, with T 2d playing the role of phase space.We see that a lifts to a doubly-periodic function ã on where â is the (non-semiclassical) Fourier series in both variables.Furthermore, we may define an analogue of S δ symbols as which the obvious generalization of (2.10) can be defined.
Most results from pseudodifferential calculus carry over to H N , see [CZ10; DJ21] for details.In particular, the formulas for compositions and Egorov's theorem are the same and follow immediately, and we have a version of the Calderón-Vaillancourt theorem [CZ10, Proposition 2.7] that says where ∥•∥ H N is the norm given by the Hilbert space structure.In fact, one can directly bound Op N (a) by the corresponding map on L 2 (R d ): (2.11) see [DJ21, §2.2.3] for an explanation using the direct integral decomposition of L 2 (R d ).
A particular class of metaplectic operators can also be defined on H N .Let M ∈ Sp(2d, Z) be an integer symplectic matrix, and let N be even.Define M N = M | H N , which can be shown [DJ21, §2.2.4] to be a unitary map on H N .We then have an analogous exact Egorov's theorem: (2.12) The definition of a defect measure on H N is also identical.Given an N -dependent sequence of density operators on H N , we say they have defect measure µ if for every The proof of Proposition A.1 shows that every defect measure on the torus is uniform.
In particular, we have the following result to be an immediate consequence.
Proposition 2.10.If ρ N are density operators on H N with defect measure µ, then µ is a probability measure, and for every ε > 0 there is an h 0 such that if h < h 0 and ∥a (α) ∥ ⩽ 1 for all |α| < K(d), then The notion of a (δ, θ)-defect measure (2.7) also carries over to the quantized torus case, where the symbols in this case are in S δ (T 2d ).Proposition 2.7, which gives the defect measures of generalized coherent states, has a natural analog on the quantized torus.We do not claim these values of δ, θ to be sharp.
To prepare for the next proposition, we define the map Π (2.13) Proposition 2.12.Given β ∈ (0, 1) and g ∈ S(R d ), define Next, define ψ N := Π N g N , and Remark 1.We can easily see that . Indeed, by Riemann sums: ).In particular, this shows that the constants C N , C ′ N in Propositions 2.11 and 2.12 can actually be taken to be 1 without affecting the value of the defect measure, or even the possible values of (δ, θ).They are only included so the density operators match our physical motivation of states with norm exactly 1.
The proofs of Propositions 2.11 and 2.12 are given in Appendix B.2. 2.5.Measurement Procedure.We provide a brief background on the quantum mechanics used in this paper.For a more comprehensive overview, see [BG09] [Hol03].To model particles in quantum mechanics, we begin by fixing a Hilbert space H, known as the state space.Particles in our model are described by density operators.Definition 2.13 (density operator).A density operator on H is a symmetric, positive semidefinite, trace-class operator on H with trace 1.Quantum particles in this paper are measured to have some value in Ω (a locally compact metric space) by a type of quantum instrument called Kraus operators.
Definition 2.14 (Kraus operators).Suppose Ω is a locally compact metric space with Borel measure ν.And suppose for q ∈ Ω, A q is a bounded operator on H, q → A q is measurable, and Ω A * q A q dν(q) = I.
Then we call A q 's Kraus operators.
Given a quantum particle described by a density operator ρ and Kraus operators A q on Ω, the probability of measuring ρ in F ⊂ Ω (a measurable set) is: After measurement within F , the state changes to: For each measurable set F ⊂ Ω, we can define I(F ) := ρ → F A q ρA * q dν(q).In this case, I is an example of a quantum instrument.This paper considers a specific class of Kraus operators.The Hilbert space we work on is the quantized torus H N (recall §2.4).We fix a compact metric space Ω with finite Borel measure ν(q), and let where Op N is defined in (2.10).Additionally, assume that f q are in S (recall Definition 2.1) uniformly in q, meaning ∥∂ α f q ∥ L ∞ ⩽ C α with constant independent of q.These Op N (f q ) are the Kraus operators we will consider and we will refer to f q as the symbols of the Kraus operators.
As a natural example of the above construction, the reader may consider the example where Ω = T d , dν(q) = dq, and f q (x, ξ) = f (x − q), where T d |f (x)| 2 dx = 1.Then conjugation by Op N (f q ) corresponds intuitively to an "approximate position measurement" (see also [Bal+18]), where the approximation approaches a (nonexistent) "exact" position measurement as f approaches a delta function.
2.6.Evolution Procedure.This paper models quantum particles which are repeatedly measured and evolved.The previous section discussed the measurement procedure, here we discuss the evolution procedure.
We consider two different evolution procedures.Particles will either be evolved by (1) a metaplectic operator or (2) exponentiation of a quantization of a Hamiltonian on the torus.In the first case, we fix M ∈ Sp(2d, Z) and let M N be the corresponding metaplectic operator acting on H N (recall §2.2).In the second case, we fix p ∈ C ∞ (T 2d ), and let P = Op N (p) (recall (2.10)).
Define an operator Φ * N,q on density operators given by Φ * N,q (ρ) : where U is a unitary operator given either by U = M N or U = e −2πiN P .We remark that Φ * N,q (ρ) is U ρ(q)U −1 where ρ(q) = Op N (f q )ρ Op N (f q ) * is an un-normalized a posteriori state of the instrument determined by Kraus operators Op N (f q ).We stress again that Φ * N,q (ρ) is not a density operator, as it is not normalized; the trace will be taken at the end and interpreted as the probability density of the trajectory (q 0 , . . ., q n ).
The evolution of the quantum system by U corresponds to a classical dynamical system ϕ t : T 2d → T 2d .In the case of evolving by the metaplectic operator M N , ϕ t := M t (where t takes integer values).In the case of evolving by e −2πiN t Op N (p) , the corresponding classical dynamical system is ϕ t := exp(tH p ) where H p is the Hamiltonian vector field generated by p.
The quantitative convergence rates of the measures on the quantum and classical trajectories depend on the Lyapunov exponents (recall Definition 2.3) of these corresponding classical dynamical systems.A system with a larger Lyapunov exponent is more exponentially sensitive to initial conditions, and hence seen as "more chaotic."In particular, completely integrable systems have Lyapunov exponents equal to zero.
Lyapunov exponents appear crucially in the following bound on symbol seminorms of a time-evolved symbol.
Proposition 2.15.Suppose ϕ t is a smooth dynamical system on T 2d where t either takes values in Z ⩾0 or R ⩾0 (in which case we assume ϕ t is smooth in time), with Lyapunov exponent Γ. Suppose a ∈ C ∞ (T 2d ), and define In particular if a ∈ S and t ⩽ δ Γ+ε log 1 h , then a t ∈ S δ with seminorms independent of t.
A version of Proposition 2.15 on manifolds is discussed in [AN07, §5.2], [DG14, Appendix C], and [DJN22, §5.2].We present a proof of Proposition 2.15 in Appendix C, adapting the arguments of the mentioned references to the simpler case of the torus.

Main Result and Proof
We are now able to provide a quantitative proof of our result.
Theorem 1. Suppose H N is the quantized torus of dimension 2d, ρ = ρ N : H N → H N is a set of density operators depending on N with semiclassical defect measure µ on T 2d , Ω is a compact metric space with finite Borel measure ν, f q ∈ C ∞ (T 2d ) satisfy (2.14), P (n) µ is defined by (1.2), and P (n) N,ρ is defined by (1.1), with Φ * N,q defined by (2.15).Then for every n ∈ N we have as N → ∞ the convergence in total variation Furthermore, we have the quantitative estimate that if µ is a (δ, θ)-defect measure as defined in (2.7) and ε > 0, then there are C 0 , C 1 independent of n such that if then where Γ is the Lyapunov exponent of M (metaplectic case) or ϕ t (Hamiltonian evolution case).
Remark 2. The condition (3.1) can be recognized as the Ehrenfest time given in (2.4).
Proof of Theorem 1.This proof goes along the lines of [BFF22, Theorem 3.1], but gives quantitative bounds.The version for U = e −2πiN P requires the quantitative version of Egorov's theorem given in (2.5).
We first consider when U = M N .We have dP (3.3) The first three equalities follow from the definition of Φ * N,q j , commutativity of the trace, the exact Egorov property (2.12) of M N , and the unitarity of M * N .We then used Proposition 2.15 to see that for each i, f so that we could use the composition rule for pseudodifferential operators (2.3) together with (2.11) to get the fourth equality.We lastly used that ρ has (δ, θ)-defect measure µ to get the last equality.Note that we may remove the poly log(N ) term by adjusting ε and replacing δ with a slightly smaller δ ′ .
On the other hand, We may now integrate in q and apply the dominated convergence theorem (using that the ∥f q ∥ are uniformly bounded in q) to get the theorem for the case that U = M N .
When U = e −2πiN P , the proof is the same but uses the quantitative version of Egorov's theorem (2.5) on the second line of (3.3), to give The O(N −(2−3δ) ) error is dominated by the O(N −(1−2δ) ) error later on in (3.3), and the remainder of the proof is the same.■

Numerical Illustrations
Here we present numerical simulations to illustrate the results of this paper.In these numerics, we simulate the evolution of a state on the quantization of the 2-dimensional torus.The measurement procedure will use Kraus operators f q (x, ξ) = f (x − q), with where σ > 0 is a parameter to be chosen 2 and c > 0 is such that T |f | 2 = 1.At each time step, after observation, we evolve the particle with the metaplectic operator M N (as defined in §2.4).For now, we let M be Arnold's cat map, 2 σ can be viewed as the measurement precision.The smaller the value of σ, the more precise.
Let Q i be the observed location of the state at time i.These Q i 's are random variables, and their joint distribution is described by (1.1).To numerically simulate this distribution, we note that formally P ( n i=0 {Q i = q i }) can be written Therefore to sample Q i from this distribution, we first sample Q 0 with law P(Q 0 = q 0 ) dν(q 0 ), and get a value, which we call q0 .Then we sample Q 1 with law We continue this process to get values of qi .
For our choice of initial state and Kraus operators, the probability density functions of these distributions greatly simplify.Working in the basis of H N described in §2.4,we write |ψ N ⟩ as a vector v(0) ∈ C N (the argument of v will denote the time at which the state is to be observed).It is a straightforward computation that, with respect to this basis, Op N (f q ) = diag(f (x 1 − q), f (x 2 − q), . . ., f (x N − q)) where x i = i/N .To ease notation, let M N := U and Op N (f q ) := A q .Then, by (1.1), Observe that this is the convolution of the probability density functions f (x) 2 and k 1 x k v(0) 2 k (call random variables with these distributions F and V (0) respectively).Here 1 x k is 1 at x k and zero everywhere else.Therefore to sample Q 0 , we sample from the random variable F + V (0) to get q0 .
To sample Q i , for i > 1, we proceed in a similar manner.Assuming we have already computed that Q j = qj for j < i, we must sample a random variable with distribution Like before, let V (i) be a random variable with probability density function k v(i) 2 k 1 x k .Therefore, to sample Q i , sample from V (i) + F .This proves the following proposition.
Proposition 4.1.Suppose ρ = |ψ N ⟩⟨ψ N | is an initial state (as defined by (4.3)), f (x − q) are the symbols of Kraus operators (given by (4.1)), and n > 0. Then a random trajectory of a quantum particle with distribution P (n) N,ρ , described by (1.1), can be sampled by applying Algorithm 1 by storing the values q 0 , q 1 , . . ., q n (representing the observed locations of the particle at each instance of time).
Algorithm 1 Quantum Particle Evolution g ← sampled random number from Gaussian distribution of mean 0, variance σ 2 5: In this algorithm, we use the fact that in the fixed basis of H N , the components of M N are given by as computed in [DJ21].
Theorem 1 states that this distribution, P N,ρ , will converge (as N → ∞) to P (n) µ .By Proposition 2.12, the semi-classical defect measure µ for our chosen ρ is δ x=x 0 ,ξ=p 0 , so that Therefore the distribution of q k according to P (n) µ is a Gaussian centered at the position component of M k (x 0 , p 0 ) t , with variance σ 2 .

Now we present several numerical simulations.
First, Figure 2 demonstrates the evolution procedure.An initial state is given by (4.3) with N = 100, β = 0.5, σ ′ = 0.1, x 0 = 0.5, and p 0 = 0.The symbols of the Kraus operators are given by f (x − q) where f is given by (4.1), with σ = 0.1.The absolute value squared of the components of ψ (with respect to the basis of H N ) is plotted in blue.We then apply Algorithm 1 to get a value q 1 which is plotted as an orange dot.The state is then changed by observation (by step 6 of Algorithm 1), which is plotted in the next plot.Then the particle is evolved by applying M , where M is the cat map (4.2) to get a new state.The absolute value squared of the components of the new state are plotted in the next plot, and the process is repeated for 3 time steps.for n = 4.For this simulation, the initial state is given by (4.3) with N = 2000, β = 0.5, σ ′ = 0.1, x 0 = 0.1, and p 0 = 0.1.The symbols of the Kraus operators are given by f (x − q) where f is given by (4.1), with σ = 0.1.We then run the evolution procedure 2000 times (we call this the number of trials), saving the observed positions for each time step.The relative frequencies of the recorded positions are plotted as a histogram in blue.As the number of trials goes to infinity, these histograms should as an orange curve.This paper's main result states that the two distributions should agree until the Ehrenfest time.For the cat map, the Lyapunov exponent is Γ ≈ .96,so the Ehrenfest time is approximately log(N )/(2Γ) ≈ 3.95.

It is apparent that the convergence of P
as N → ∞ is relatively slow.As seen in Figure 3, after only 4 time steps, the distribution of observed positions starts to become uniform.This should be expected, as Theorem 1 only guarantees (3.2) for N ≳ exp((Γ + ε)n/δ).For the cat map, Γ = log((3 + √ 5)/2), and in this case δ < 1/4.Then, (3.2) requires (roughly) N ≳ 50 n .These estimates are far from sharp, but give some sense of the exponentially slow rate of convergence as seen in Figure 3.
In Figure 4, we run the same simulation as in Figure 3 but vary N from 100 to 4000.For each value of N , we compute the total variation of the (approximate) marginal distribution of P (time) N,ρ against P (time) µ for times between 1 and 6.Observe that increasing N leads to smaller total variation at an exponentially slow rate.µ .In Figure 5, the same evolution procedure is performed for 4 time steps but with 3 different quantized symplectic matrices.In each case, we plot the relative frequencies of the observed positions (approximating the marginal distributions of P µ .In all cases, we chose N = 3000, β = 0.5, σ ′ = σ = 0.1, x 0 = 0.3, p 0 = 0.9, and simulated 500 trajectories.The symplectic matrices chosen were: N,ρ and P (n) µ disagree.Because the Ehrenfest time grows logarithmically with N , and the computations grow like N 2 , it is not feasible to simulate numerics for which the two measures agree for many more time steps.
In Figure 6 we replace the evolution of the particle by M N at each time step, with evolution by exp(−2πitN Op N (a)), where a = cos(2πx) + cos(2πξ), and t = 0.05which should be thought of as the interval of time between observations.In this case, P (n) µ is a Gaussian whose mean is computed by evolving (x 0 , p 0 ) along the Hamiltonian flow generated by a(x, ξ) for time tn.Note that the convergence of P  more rapid in terms of N (as long as t is small).We chose the following parameters: N = 1000, β = 0.5, σ ′ = σ = 0.1, x 0 = 0.4, p 0 = 0.7, with 1000 trials.
Lastly, in Figure 7, we present numerics for the evolution of quantum particles, where the cat map is replaced by the matrix: which is an example of an elliptic matrix, having Lyapunov exponent equal to zero (in fact, E 3 = I).The parameters chosen here are N = 1000, β = 0.5, σ = σ ′ = 0.1, x 0 = 0.1, p 0 = 0.9, with 1000 trials.In this case P (n) N,ρ and P (n) µ agree for more time steps than the evolution of the other symplectic maps.The Ehrenfest time will subtly depend on the ε in (3.1).
The figure on the first page (Figure 1) is simulated in the following way.We chose an initial state ρ given in (4.3) with N = 500, σ ′ = 0.1, β = 0.1, x 0 = 0.4, p 0 = 0.7.We indirectly measured the particle location using Kraus operators with symbol given by (4.1) with σ = 0.1.Between each measurement, the particle is evolved by exp(−2πitN P ) with P = Op N (cos(2πx) + cos(2πξ)) and t = 0.01.This procedure is done for 100 time steps to sample a single trajectory with law P (100) N,ρ .We then repeat this 60 times to get 60 different trajectories, which are plotted on top of each other.We plot a single quantum trajectory in red.The corresponding classical probability  measure is a product of Gaussian distributions with variance σ 2 and at time t is centered at ϕ t ((x 0 , p 0 )).The flow ϕ t , the flow generated by the Hamiltonian vector field with Hamiltonian cos(2πx) + cos(2πξ), is numerically computed by solving with initial conditions x(0) = x 0 and ξ(0) = p 0 .We plot x(t) as a blue dotted line.
The reason for the name is given in the following proposition.We use the standard multi-index notation.For a ∈ C ∞ (T * R d ) and α ∈ N 2d , we let Proposition A.1.Let ρ h be a set of density operators with defect measure µ.The following are equivalent.
To show (2) =⇒ (3), let a ∈ S, and let χ k ∈ C ∞ 0 (T * R d ) be smooth partition of unity of T * R d , chosen as in [Zwo12, Theorem 4.23], with the first K(d) derivatives chosen small enough.By the Cotlar-Stein-Knapp lemma, we have the limit in the strong operator topology: for every ψ ∈ L 2 (R d ).We need the following brief lemma.
Lemma A.2. Let A K be a sequence of uniformly bounded operators converging to A in the strong topology, and let ρ be self-adjoint and trace class.Then Proof of Lemma A.2.By subtracting we may assume A = 0, and without loss of generality let sup K ∥A K ∥ = 1.Let ψ j be an orthonormal basis of eigenvectors of the compact operator ρ, so ρ = j λ j |ψ j ⟩ ⟨ψ j | with λ j → 0. as j → ∞.Let ε > 0, and choose J such that λ j < ε for j > J. Then ■ By Lemma A.2, we see that for fixed h, and similarly as µ is a finite measure lim Then letting ε > 0 and picking the h 0 (ε) from part (2) we get We also note that χ may be chosen so its first K(d) derivatives are small (possibly making R larger).Then we have that for C depending on the first K(d) derivatives of a. Next we apply (A.3) with the symbol 1 − χ (which is independent of a) to get that as h → 0. So that there exists h 0 > 0 such that for 0 < h < h 0 we have C tr(Op h (1 − χ)ρ h ) < ε/4.With this, the triangle inequality, and symbol calculus, we have for where for the last inequality we have chosen possibly smaller h 0 (ε) and used the fact that a has bounded derivatives.
We now prove (3) using a compactness argument.Observe that by Calderón-Vaillancourt, if (A.3) holds for a ∈ S, it must also hold for any a ∈ C K 0 (d) (T * R d ) for some sufficiently large K 0 .Given h j , let X j be the set of for all h < h j .Then for K 0 (d) large enough, condition (A.5) makes X j to be an open set in C K 0 (d) (B R (0)), so letting h j → 0, gives (X j ) j∈N to be an open cover of ), so there is a finite subset of the X ′ j s that covers Y , and hence an h 0 such that for h < h 0 ) completes the proof of (3).
To show (3) =⇒ (4), we simply see that Here we use that the function f (x) = 1 is in the symbol class S.
Finally, we show (4 as h → 0, so in particular (A.3) holds in the special case of the noncompact symbol 1 − χ.But that was the only symbol used in the proof of (1) =⇒ (3), so in particular shows (4) =⇒ (3) as well.This shows all equivalences and hence completes the proof.■

Appendix B. Coherent States
In this appendix, we prove Propositions 2.9, 2.11, and 2.12, which give examples of (δ, θ)-defect measures.B.1.Coherent States on R d .We prove Proposition 2.9.The proof is essentially the same as [Zwo12, §5.1 Examples 1 and 2], but is quantitative and keeps track of errors.
Proof of Proposition 2.9.We first consider the case when β = 0.In this case, we can without loss of generality assume x 0 = 0, by absorbing it into the definition of g.Then We evaluate the inner double integral using the explicit form of stationary phase given in [Zwo12, Theorem 3.17].
with ∼ referring to asymptotic summation and the last equality coming from a being an element of S δ .Then which gives the result for β = 0.The case of β = 1 follows from taking the semiclassical Fourier transform.
We now treat the case when 0 < β < 1.Here we have We again have that tr(ρ h Op h (a)) = ⟨Op h (a)ψ h , ψ h ⟩.We now approximate Op h (a) by its left quantization defined as: Because ∥g∥ L 2 (R d ) = 1, we have by the Fourier inversion formula that: By the mean value theorem and the fact that a ∈ S δ ∩C ∞ 0 (T * R d ), we have for x, ξ ∈ R d , So, using that g is Schwartz, Noting that 1 − 2δ = β − δ + 1 − β − δ shows the O(h 1−2δ ) term to be superfluous and gives the desired result.■ B.2. Coherent States on the Quantized Torus.We provide full proofs of Propositions 2.11 and 2.12, though we do not claim our parameters to be optimal.
Proof of Proposition 2.11.For this proof, we work in the orthonormal basis (Q n ).Then for a ∈ S δ (T 2d ), Op N (a) has a representation as a matrix A with Note that there exists a large constant M 0 > 0 (independent of N ) such that Therefore: Note that there exists c d > 0 such that for all r ̸ = 0, by the same reasoning as in (B.2).By Fourier inversion in the first variable, we have in that case where F 2 is the Fourier coefficient in the second component (so F 2 a acts on Then we compute But when |j − m| ⩽ M , we have and Therefore for all γ > δ. ■ Proof of Proposition 2.12.Let a ∈ S δ (T 2d ) and let A = Op N (a) given in the basis (Q n ).Once again choose M = ⌊N γ ⌋ for γ ∈ (δ, 1).We have by (B.2) and (B.3) that (writing ψ := ψ N for convenience and letting subscripts refer to the indices) But as g ∈ S(R d ), it has super-polynomial decay, and we get with the last line using that δ ⩽ β.We substitute and use |g| 2 = 1 to obtain This gives which completes the proof by taking ε and γ − δ small enough.■

Appendix C. Estimating Derivatives of a Flow
Here we present a proof of Proposition 2.15 which we restate here for convenience.
Proposition C.1.Suppose ϕ t : T 2d → T 2d is a smooth dynamical system on T 2d with Lyapunov exponent Λ where t either takes integer values, or values in R ⩾0 .Let a ∈ C ∞ (T 2d ).Then for all ε > 0 and α ∈ N 2d , there exists C = C(ε, α) such that: Here C k is the norm defined as This proof of Proposition C.1 essentially follows [DJN22, §5.2] but in the simpler setting of the torus.To our knowledge, such a proof in this setting is not present in the current literature.In an effort to improve readability, we include examples throughout the proof, at the cost of being rather voluminous.
Proof.We prove Proposition C.1 by strong induction on the order of α.We trivially have the claim when |α| = 0.If |α| = 1, let i be such that ∂ α = ∂ i .Then we have: where (ϕ t (x)) j is the jth component of ϕ t (x).By the definition of the Lyapunov exponent, for any ε > 0, there exists C = C(ε) > 0 such that: To ease notation in the remained of the proof, define λ := Λ + ε.
Before proving the inductive step, we show how to go from α of order 1 to order 2. We will use the following notation for Jacobians and tensors.If f : R → R 2d , then Df will be the vector with ith component ∂ i f .Note that Df will be a (0, 1) tensor.We define D 2 f as the (0, 2) tensor with i, j entry (D 2 f ) i,j = ∂ i ∂ j f .If ϕ : R 2d → R 2d , then we let Dϕ be the (1, 1) tensor with i, j entry (Dϕ) i j = ∂ j ϕ i .We let D 2 ϕ denote the (1, 2) tensor with i, j, k entry (D where the product of two tensors is the contraction of the upper and lower components.That is, (D(a • ϕ n−1 ) • ϕ 1 ) is a (0, 1) tensor which can be written A µ , and Dϕ 1 is a (1, 1) tensor which can be written B ν γ .In this case AB is the (0, 1) tensor contracting the lower component of A with the upper component of B: , so that By the upper triangular structure of M , it is easy to check that where: with (M 2 ) j (x) defined as where a ⊗2 := a ⊗ a.
We then have, by (C.1) and (C.2), where the absolute value of a tensor denotes the maximum of the absolute value of all entries.
But now observe that by the inductive hypothesis and setting a to be coordinate functions.And we also have, by (C.3), that where we use that ϕ 1 is a smooth function on a compact set so that ∥D 2 ϕ 1 ∥ L ∞ < ∞.
But we stress that this norm only comes up once in each (M 2 ) j term so that the constant at the end is independent of n.We therefore have that: We can absorb the (1 − e −λ ) −1 term into the constant C to get, using (C.4), (C.5), and (C.6), that Now if the flow is continuous in time and t ∈ R ⩾0 , we let n ∈ N be such that n ⩽ t ⩽ n + 1 and observe that: Note that s := t − n ∈ [0, 1] which is a compact set, so that: We can then absorb max s∈[0,1] ∥ϕ s ∥ C 2 into the constant C to get the result for continuous time.
Now we present the general inductive step, but the reader should keep in mind that morally the idea is the same as above, but with more notation.Suppose we have for all k ′ ⩽ k − 1 and |α| = k ′ : We then fix n ∈ N, and aim to bound D k (a • ϕ n )(x).Define A n (x) as the 1 × k vector of tensors whose ith entry is D i (a • ϕ n )(x).By writing and repeatedly applying the product rule, we can find a k × k matrix of tensors, M (x) such that A n (x) = A n−1 (ϕ 1 (x))M (x).Note that M (x) will be upper triangular.Indeed, for any i ⩽ k, D i (a • ϕ n ) is a linear combination of D j (a • ϕ n−1 ) • ϕ 1 , for j ⩽ i.Also, for each i, j = 1, . . ., k, the i, j entry of M (x) must be a (i, j) tensor because it will be paired with a (0, i) tensor to get a (0, j) tensor.Each component of (M (x)) i,j will contain j derivatives of the components of ϕ 1 (x) with a coefficient depending on k.Lastly note that for each i = 1, . . ., k, (M (x)) i,i = (Dϕ 1 (x)) ⊗i as this is the term in (C.8) when all derivatives fall on D(a • ϕ n−1 • ϕ 1 ).
We then recursively use the relation of A n to A n−1 to get that: M (ϕ n−j (x)) .
For i = 1, . . ., k, let M i be the (i, k) entry of n j=1 M (ϕ n−j ), so that: (D j (a) • ϕ n (x))M j (x).(C.9) It now suffices to provide appropriate estimates of M j (x).For i = 1, . . ., n let B i := M (ϕ n−i )(x).Then for each fixed j = 1, . . ., k: Because B i are upper triangular, this sum is nonzero only if ℓ m ⩽ ℓ m+1 for each m = 1, . . ., n − 1, and j ⩽ ℓ 1 .Each term in the sum can be expressed as a path on a grid of k × (n + 1) nodes.Labeling the bottom left node (1, 1), every path must begin at (j, 1) and end at (k, n + 1).The path can only go to the right or up.That is (m, n) can only go to (m ′ , n + 1) for m ′ ⩾ m (provided that the path terminates at (k, n + 1)).The edges of the path give the index of B i in the sum.That is, if the nodes (m, i), (n, i + 1) are on the path, then we include the term B i m,n in the sum.For example, if k = 4 and n = 5, then the following path represents the term B 1 1,2 B 2 2,2 B 3 2,2 B 4 2,4 B 5 4,4 .
Each path can be represented as a vector whose ith component is m where (m, i) is on the path, e.g. the example path is written as (1, 2, 2, 2, 4, 4).
Each path will consist of flat sections (where the i and i − 1 components of the vector representation are the same) and jumps (where the component of the vector representation is not flat).In the above example, the jumps are at indices 2 and 5.The corresponding entry of B i for the jumps will be individually norm-bounded.Recall that B i j,ℓ will be a (j, ℓ) tensor where all terms are derivatives of ϕ 1 .As long as (j, ℓ) ̸ = (1, k), at most k − 1 derivatives are taken of ϕ 1 .Note that B 1,k can only appear once in any term in the sum (C.10), and this is bounded by the C k norm of ϕ 1 .We therefore can apply the inductive hypothesis to get that |B i j,ℓ (x)| ⩽ Ce λℓ .Next, we bound the flat sequences.Suppose such a sequence is j 2 i=j 1 B i ℓ,ℓ .Recall that B i ℓ,ℓ = (Dϕ 1 ) ⊗ℓ • (ϕ n−i (x)), so that This follows by the inductive hypothesis by setting a to be coordinate functions.Now suppose we have a term of the sum (C.10) whose path has a vector representation (a 1 , a 2 , . . ., a n+1 ) with jumps at indices ℓ 1 , ℓ 2 , . . ., ℓ s .We assume that there is a jump at index 2 and the path ends in a flat section (as in the example), but all cases can be treated similarly.We can then bound this term of the sum (using (C.11)) by: a ℓs ,a ℓs ⩽ Ce λ(a ℓ 1 +•••+a ℓs ) e λa ℓ 1 (ℓ 2 −ℓ 1 −1) • • • e λa ℓs (n+1−ℓs−1) = Ce λ(a ℓ 1 (ℓ 2 −ℓ 1 )+•••+a ℓ s−1 (ℓs−ℓ s−1 )+a ℓs (n+1−ℓs)) = Ce λ(a 2 +•••+a n+1 ) .
The last equality follows by noting that ℓ i+1 − ℓ i is the length of each flat section and is where the sequence has values a ℓ i .Therefore: (an)∈A j,k e λ(a 2 +•••+a n+1 ) (C.12) where A j,k is the set of all sequences (a n ) of length n + 1 such that a i ∈ Z >0 , a 1 = j, a n+1 = k, and a i ⩽ a i+1 .Let V be the k × k upper triangular matrix such that V ℓ,m = e mλ for ℓ ⩽ m and zero otherwise.Then by the exact same argument for computing products of upper triangular matrices, we see that e λ(a 2 +•••+a n+1 ) .(C.13) Because the eigenvalues of V are distinct, V is diagonalizable.We can then write the entries of V n as linear combinations (with coefficients independent of n) of the eigenvalues of V raised to the nth power.That is there exist c j,ℓ,k ∈ R independent of n such that e λℓn c j,ℓ,k .
But then we immediately see that there exists a C = C(k) > 0 such that |(V n ) j,k | ⩽ Ce λkn for j = 1, . . ., k.We therefore have that |M j (x)| ⩽ Ce nλk , so that (C.9) becomes: Applying the exact same argument at the end of the case k = 2 to go from discrete time to continuous time completes the proof. ■

Figure 1 .
Figure 1.A numerical simulation demonstrating this paper's main result.Each grey curve represents the trajectory of a quantum particle (under the evolution procedure discussed in this paper) on a torus with identical initial conditions and a fixed Hamiltonian.A single trajectory is plotted in red and the corresponding classical trajectory is plotted as a blue dotted line.Up until time 0.4, most of the positions of the quantum particles are within 0.1 of the classical trajectory, but shortly after become decoherent.More numerics are provided in §4.This simulation is discussed at the end of §1.1 and §4.

Figure 2 .
Figure 2. Evolution Procedure Here we simulate the trajectory of a quantum particle for three time steps under the discussed evolution procedure.The specifics of this simulation are discussed in §4.

Figure 3
Figure 3 compares an approximation of the marginal distributions of P (n) N,ρ against P (n) µ

Figure 3 .
Figure 3. Marginal distributions of P (n) N,ρ vs P (n) µ : We perform Algorithm 1 (observation then evolution by the quantum cat map).We plot the relative frequencies of the observed positions (approximating the marginal distribution of P (n) N,ρ against the marginal distributions of P (n) µ ).See §4 for more details.

Figure 4 .
Figure 4. Convergence of marginal distributions of P (n) N,ρ → P (n) µ : We run the same simulation as in Figure 3, but vary N from 100 to 4000 and compare P (n) N,ρ to P (n) µ .See §4 for more details.

Figure 5 .
Figure 5. Changing the Lyapunov exponent: We present the same numerics as in Figure 3 but evolve by quantizations of three different matrices given by (4.6) with 3 different Lyapunov exponents.See §4 for more details.

Figure 6 .
Figure 6.Evolution of cos(2πx) + cos(2πξ): We present the same numerics as in Figure 3 but instead of evolving the quantum state by applying the quantum cat map, we multiply our state by exp(−2πitN Op N (cos(2πx) + cos(2πξ))).See §4 for more details.

Figure 7 .
Figure 7. Evolution by Elliptic Matrix: We present the same numerics as in Figure 3 but instead of evolving by the quantum cat map, we evolve by the quantization of the elliptic matrix E which has Lyapunov exponent zero.See §4 for more details.