Computing Quantum Mean Values in the Deep Chaotic Regime

Semiclassical approximations to quantum mechanics, which we distinguish from purely “classical approximations” in the sense that they should account for interference effects, have a long and venerable history since they have been introduced almost immediately after the invention of quantum mechanics [1] (and remarkably somewhat before [2]). The successes of these semiclassical techniques in providing both effective approximation tools and physical understanding of the mechanisms at play have been demonstrated over time in a large number of systems in the integrable or nearly integrable regime [3, *tannor2000semiclassical, *brodier2001resonance, *tomsovic2018post, *ray2016dynamics, *lando2019quantum, *schlagheck2022resurgent] as well as in the deep chaotic regime [10, *tomsovic1993long, *sieber2001correlations, *heusler2007periodic].

Surprisingly, the semiclassical computation of some very simple and common quantities poses both conceptual and practical problems that have not been solved until today. The conceptual problem has to do with the internal coherence of semiclassical approximations, namely the fact that different semiclassical approaches, or the use of different coordinate systems, must lead to the same result up to negligibly small terms. The derivation of this equivalence almost always implies the use of the stationary-phase (or steepest-descent) approximation (SPA), which appears therefore as one of the keystones on which we strongly rely to ensure the soundness of semiclassical approaches.

The problem is the following. For some physical quantities, among which the simplest is the time evolution of a smooth operator’s mean value, a blind application of the SPA leads to the incorrect result that only the classical contribution remains, and that all interference effects are washed out. This failure of the SPA, which goes beyond the need for uniform approximations [14], can actually be seen as one further example of the lack of commutation between the semiclassical ($\hbar\to 0$) and long-time ($t\to\infty$) limits [15]. The practical consequence of this conceptual difficulty is that the computation of several quantities always relies on some form of initial or final value representation [16, 17, *kay2005semiclassical, *ozorio2013initial], the most famous one being Herman and Kluk’s (HK) [20], which amounts essentially to evaluating all integrals using a numerical Monte Carlo technique. This family of approaches, and in particular HK, are easy to employ in practice and usually provide good approximations to correlations and spectra in systems whose classical dynamics is relatively close to integrability. However, the computation of mean values using such methods converges slowly and has larger errors when compared to correlation functions [21]. This is further worsened by the fact that, due to an inherent prefactor inversion, initial value representations fail as one enters the deep chaotic regime [16, 22]. The main consequence of such problems is that there is nowadays no semiclassical tool capable of computing something as simple as the time evolution of a smooth operator’s mean value in the deep chaotic regime – even in the near-integrable regime, no efficient methods exist. Note that a pure quantum numerical approach is also usually computationally expensive in the deep semiclassical regime (for which semiclassical expressions are most accurate), due to the need of very fine grids [23].

The conceptual difficulty mentioned above has been addressed in a recent paper by some of us [24]. The source of failure of the SPA was identified, and a canonically invariant expression for the mean value of smooth operators was derived in which only the integrals for which the SPA cannot be used remain to be performed numerically. Moreover, these integrals should be sampled only on a classical scale, and are therefore much easier to evaluate than the highly oscillatory (on a quantum scale) integrals in HK [25]. However, [24] did not assess the numerical efficiency or semiclassical accuracy of the approach, nor did it compare it with other widely used methods such as HK. The goal of this Letter is to illustrate the effectiveness of this approach on a simple system, while simultaneously providing a rather general conceptual background on why it works. In the regime where one can apply HK, we demonstrate that our approach does at least as well and at a significantly lower numerical cost. More importantly, our approach retains its remarkable accuracy as one goes deeper in the chaotic regime, while the HK results become essentially numerical noise.

As typical when employing semiclassical methods, our conclusions come accompanied by a wealth of geometrical insight into the underpinnings of mean-value calculations, tracing the source of quantum coherence to the self-interference of an evolving curve associated with the system’s initial state. Here, we assume evolution to be generated by the kicked rotor system (KRS),

where $k$ counts the number of iterations, or kicks, an initial point $(m_{0},\theta_{0})$ in phase space is subjected to [26]. The kicking strength $\alpha$ is responsible for bringing the system from a near-integrable regime to a mixed one, and later to the deep chaotic regime in which the measure of regular islands tends to zero as phase space becomes almost fully chaotic. These regimes are reached for $\alpha\ll 1$, $\alpha\approx 1$ and $\alpha\gg 1$, respectively. Since the angular momentum $m$ is unbounded, a curve evolved according to (1) can diffuse in the deeply chaotic regimes, while an initial state $|\psi_{0}\rangle$ propagated using its exact quantum equivalent,

will remain bounded ¹¹1Note that while the angle operator is not well defined, its exponential (and, therefore, its sine and cosine) are perfectly rigorous, which is why we place the hat over the whole sine instead of only the angle. See [49] for a nice exposition of the construction of angle and angular momentum operators.. As is well-known, this renders the KRS a prototype for theoretical and experimental studies on dynamical and Anderson localizations [28, *floss2013anderson, *manai2015experimental].

We will focus on the time-dependent mean value $\langle\widehat{O}\rangle(t=k)=\langle\psi_{k}|\widehat{O}|\psi_{k}\rangle$ with the initial state given by a plane wave of initial angular momentum $M_{0}$,

with $M_{0}=n\,\hbar$, $n\in\mathbb{N}$. The operator $\widehat{O}$ will be assumed to act as a “classical detector” that measures the probability of finding the system in the neighborhood of some angular momentum $m_{c}$ and some angle $\theta_{c}$, within some classical scale $\sigma$. In practice $\widehat{O}$ is most easily defined through its Weyl symbol

which is just a symmetric, normalized phase-space Gaussian centered at $(m_{c},\theta_{c})$ and with standard deviation $\sigma$. The larger we pick $\sigma$ the more classical the detector becomes, as can be seen by measuring the effective area enclosed by (4) in units of $\hbar$. We shall fix a small value for $\sigma$, namely $\sigma=0.088$, in order to clearly observe a gradual improvement in the approximations as one moves deeper into the semiclassical regime.

Using (4) the quantum mean value of $\widehat{O}$ can then be calculated according to the phase-space average

with

the Wigner function of $\psi_{k}$ (for the KRS the angular momentum is quantized, so that the integrals over $m$ should be understood as discrete sums [31]). The semiclassical HK expression, $\langle\widehat{O}\rangle_{\text{HK}}$, is calculated analogously to $\langle\widehat{O}\rangle_{\text{QM}}$ but with the wavefunction $\psi_{k}$ replaced by its HK equivalent [31].

The WKB quantization of the initial Lagrangian manifold $\mathcal{L}_{0}=\{(M_{0},\theta_{0}),\,\theta_{0}\in[0,2\pi)\}$ results exactly in expression (3). By iterating $\mathcal{L}_{0}$ classically using (1), we produce an evolved Lagrangian manifold $\mathcal{L}_{k}=\{(m_{k},\theta_{k})\}$ that is the exact classical counterpart of (2) with initial state (3). As an example, Fig. 1 displays a snapshot of the classical manifold $\mathcal{L}_{k}$ (top) and the quantum Wigner function $\mathcal{W}[\psi_{k}]$ (bottom) after $k=200$ iterations for the KRS in the near-integrable regime. The classical mean value, which we call $\langle\widehat{O}\rangle_{\text{CL}}$, is then given by substituting $\mathcal{W}[\psi_{k}]$ by a delta function with support on $\mathcal{L}_{k}$ in (5) and is known in physics literature as the truncated Wigner approximation (e.g. [32, *polkovnikov2003quantum, *schlagheck2019enhancement]).

Following Maslov [35], our semiclassical mean value (SMV) approach is exclusively based on information contained in $\mathcal{L}_{k}$ and is given by

as adapted from [24]. Here, each $\beta$ entering the sum represents a pair of sections of $\mathcal{L}_{k}$ that fall inside the detector, which we denote by $\mathcal{L}_{k}^{+}$ and $\mathcal{L}_{k}^{-}$. We refer to each of these sections as a filament. For each phase-space point $x_{\pm}=(m_{\pm},\theta_{\pm})\in\mathcal{L}_{k}^{\pm}$, the coordinates $\eta$ are formed as the centers $\eta=(x_{+}+x_{-})/2$. The center action $S^{\beta}$ is given by the symplectic area of the region enclosed by the path between $x_{+}$ and $x_{-}$ on $\mathcal{L}_{k}$ and the line segment joining $x_{+}$ and $x_{-}$ [36, 24]. Along the path on $\mathcal{L}_{k}$, one is expected to cross several points at which the tangent space to $\mathcal{L}_{k}$ is parallel to the $m$ axis, which semiclassically gives rise to divergences in the $\theta$ representation, known as caustics [37, *littlejohn1992van, *OzorioBook]. While (7) is not impacted by them, it is fundamental to count how many were traversed in the form of a Morse index $\mu^{\beta}$ [35, 40]. Lastly, the coefficients $a_{\pm}$ measure how much the manifold was stretched during evolution, i.e. they can be thought of as the ratio between the length of an initial line segment in $\mathcal{L}_{0}$ and of its image in $\mathcal{L}_{k}$. For a detailed discussion of these concepts, see [31].

The classical structure of Fig. 1(a) is standard for near-integrable systems: The initial manifold $\mathcal{L}_{0}$ is stretched due to the shearing action of (1) when far from fixed points, but also captured near them, spiraling into whorls [41, 42]. The semiclassical interpretation of the oscillations seen in Fig. 1(b) is that all the branches in $\mathcal{L}_{k}$ interfere with each other, forming the complicated patterns seen in the Wigner function. Before the so called Ehrenfest time [43], classical and quantum evolutions produce the same expectation values (that is, quantum interference cancels out when computing mean values). Since we pick a large number of kicks, we are far beyond the Ehrenfest time [43], and the classical features of Fig. 1(a) are blurred by quantum interferences in Fig. 1(b). By selecting a set of interfering branches in $\mathcal{L}_{k}$, detector (4) allows us to define a very objective notion of “quantumness”: the mean value is all the more quantum as more interfering segments of $\mathcal{L}_{k}$ fall inside the detector. The number of such filaments depends mostly on the curve defining the initial state and on the dynamics of the system, but for a high degree of chaoticity it will grow exponentially with $k$. In Fig. 1(a), for instance, there are 24 filaments within the yellow region defining 4 standard deviations of the detector.

The success of a semiclassical method will depend on how accurately it is able to reproduce pairwise interferences between filaments. This task becomes harder if the number of filaments increases, which can happen by increasing either the number of kicks, $k$, or the kicking strength, $\alpha$. From the detector’s point of view, the long-time (large $k$) semiclassical dynamics and the deeply chaotic regime (large $\alpha$) are very similar, and the only important measure is the number of filaments falling inside its domain. Both limits, however, are well-known to be problematic for semiclassical mechanics [15], and producing reasonable quantitative results for a large number of filaments is a formidable task. To have a clearer idea on how intricate such a task is, in Fig. 2 we zoom in on the detector in order to have a better view of the classical and quantum structures inside it, both in the near-integrable regime of Fig. 1 and in the deep chaotic regime. The rather structured Wigner function Fig. 2(b) in the near-integrable regime, semiclassically built on top of the simple filaments of Fig. 2(a), finds an intimidating analogue in Fig. 2(d). By emerging as a result of almost one thousand interfering filaments in Fig. 2(c), the corresponding Wigner function looks ergodic and structureless.

We expect a good match of both HK and SMV with the exact quantum mean values for the structure shown in Fig. 2(a). This is indeed confirmed in Fig. 3(a) and Fig. 3(c), where we plot mean values as a function of $\hbar$. Here, since the accuracy of semiclassical approximations improves as $\exp(i\,M_{0}\theta/\hbar)$ becomes highly oscillatory, instead of entering the semiclassical regime by increasing $M_{0}=n\,\hbar$ we define a tunable $\hbar=M_{0}/n$, tracking our results as a function of the integer $n$. The advantage is that by fixing $M_{0}$ the classical dynamics is frozen. We then see that in the near-integrable regime both HK and SMV become increasingly accurate as we decrease $\hbar$, with the former converging earlier but being slightly less accurate for smaller $\hbar$ values. On the other hand, in the deep-chaotic regime in Fig. 3(d) the HK calculation is essentially numerical noise, while the SMV in Fig. 3(b) displays a remarkable quantitative match with the quantum result. The result obtained from the HK method is also improved by the fact that we compute wavefunctions and only then the mean values. This allows us to renormalize the wavefunctions according to $\psi_{\text{HK}}\mapsto\psi_{\text{HK}}/\sqrt{\lVert\psi_{\text{HK}}\rVert}$ before plugging them into (5), correcting for normalization loss due to the fact that the HK propagator is not unitary. This procedure allows us to bypass one of the main flaws in the HK method, but can only be efficiently carried out for systems with a single degree of freedom. It is also important to be sure that the HK wavefunctions are properly converged before taking the mean values, and we provide extensive convergence testing in [31].

It can be argued that HK results could be improved by filtering trajectories or employing one of the many prefactor renormalization techniques developed throughout the years in computational chemistry [44, 22]. These techniques, however, are only successful when there exist regular trajectories that survive the filtering process [22]. In the deep chaotic regime phase space is completely dominated by the chaotic sea, and it is no surprise that even prefactor renormalization was not capable of improving the quality of HK calculations for the KRS – even when used in conjunction with wavefunction renormalization [45]. The fact that the HK wavefunctions are converged, yet highly inaccurate, confirms that the complexity of the classical structure in Fig. 2(c) is to blame: The HK method is simply incapable of accounting for all interfering filaments correctly. The large number of filaments, however, also leads to a combinatorial problem in a naive implementation of the SMV: Since all pairs enter (7) and the number of pairs is given by $N(N-1)/2$, the computational cost increases quadratically with the number of filaments, $N$. This difficulty can be soothed by keeping in mind that not all pairs need to be included in (7), only the ones for which the action is not highly oscillatory – in practice, of $\mathcal{O}(\hbar)$. This renders the computational cost linear in the number of filaments. Our implementation also relies on filaments not developing too much structure (such as coming out from the same direction as they came in, forming a finger). This is by no means a shortcoming of the method itself, and our algorithm can be improved and adapted to particular applications as our source code is made accessible to the public [31]. For dealing with the previously mentioned fingers, for instance, a uniformization procedure in terms of Airy functions can be applied to (7) [14].

In principle, our SMV method is ready to be generalized to the propagation of coherent states, where now complex Lagrangian manifolds need to be employed [46, *tomsovic2016generalized]. An implementation for higher-dimensional systems is also possible, and will mainly imply that the intersection of the Lagrangian manifold defining our propagated state with the detector will be of dimension $d\geq 2$. Although these manifolds can be quite complicated, in the deep chaotic regime they are locally $d$-dimensional planes, just as the filaments of Fig. 2 are essentially lines. This can extend the applicability of our method to complex chaotic systems (e.g. two- or more electron atoms), which lie beyond the applicability horizons of initial value representations such as Herman and Kluk’s. Moreover our approach can be adapted to address the dynamics of hybrid quantum-classical systems [48], or to compute other physical quantities with well-defined classical limits, such as the conductance matrix of quantum dots or the purity (linear entropy) of a bipartite quantum system.