Semi-Classical Discretization and Long-Time Evolution of Variable Spin Systems

Abstract

We apply the semi-classical limit of the generalized $SO\left(3\right)$ map for representation of variable-spin systems in a four-dimensional symplectic manifold and approximate their evolution terms of effective classical dynamics on $T^{*}{\mathcal{S}}_2$. Using the asymptotic form of the star-product, we manage to “quantize” one of the classical dynamic variables and introduce a discretized version of the Truncated Wigner Approximation (TWA). Two emblematic examples of quantum dynamics (rotor in an external field and two coupled spins) are analyzed, and the results of exact, continuous, and discretized versions of TWA are compared.

1. Introduction

Phase-space methods provide a very convenient framework for the analysis of large quantum systems (QS) [1,2,3,4]. In cases when the states of QS are elements of a Hilbert space $\mathbb{H}$ which carries a unitary irreducible representation of a Lie group G, trace-like maps from operators $\widehat{f}$ acting in $\mathbb{H}$ into their (Weyl) symbols $W_f\left(\mathsf{\Omega }\right)$ can be established. The symbols $W_f\left(\mathsf{\Omega }\right)$ are functions on the corresponding classical phase-space $\mathcal{M}$, $\mathsf{\Omega }\in \mathcal{M}$ being the phase-space coordinates. The most suitable in applications is frequently the Wigner (self-dual) map, allowing the average value of an observable to be computed as a convolution of its symbol with the symbol of the density matrix $W_{\rho }\left(\mathsf{\Omega }\right)$ (the Wigner function). The properties of the Wigner function are extremely useful for studying the quantum-classical correspondence. In particular, the quantum evolution is described by a partial differential equation (the Moyal equation [5]) for the Wigner function with a well-defined classical limit. In the non-harmonic case, the Moyal equation contains high-order derivatives, which turns the finding of its exact solution into quite a difficult task.

The advantage of the phase-space approach to quantum dynamics consists of the possibility of expanding the Moyal equation in powers of small parameters in the semi-classical limit. Such semi-classical parameters are related both to the symmetry of the Hamiltonian and of the map.

In general, different mappings can be performed for composite quantum systems. The type of map fixes the structure of the phase-space manifold, and thus the allowed set of shape-preserving transformations. In addition, the phase-space symmetry determines the physical semiclassical parameter $\epsilon \ll 1$.

The standard phase-space approach [6,7,8,9,10] is not applicable for the construction of covariant (under group transformations) invertible maps if the density matrix cannot be decomposed into a direct sum of components each acting in an irreducible representation of the dynamical group G. Physically, this happens when the Hamiltonian of the system induces transitions between irreducible subspaces i.e., the total angular momentum is changed in the course of evolution. This happens, for instance, in quantum systems with non-fixed (variable) spin, such as large interacting spins, a rigid rotor in an external field, coupled and externally pumped boson modes with and without decay, etc.

A suitable $SU\left(2\right)$ covariant map, establishing a one-to-one relation between operators and set of functions (discretely labelled symbols) in a three-dimensional space, can be found for variable spin systems [11]. In addition, there exists a continuous limit of such symbols for large values of the mean spin [12]. This allows the quantum operators to be put in correspondence with smooth functions on the four-dimensional cotangent bundle $T^{*}{\mathcal{S}}_2$, equipped with a symplectic structure. In particular, the semi-classical dynamics of the Wigner function of variable spin systems can be described in terms of effective “classical” trajectories ${\mathsf{\Omega }}^{cl}\left(t\right)$ in the phase-space $T^{*}{\mathcal{S}}_2$. As a rough approach, the evolution of average values can then be estimated within the framework of the so-called Truncated Wigner Approximation (TWA), which consists of propagating points of the initial distribution along the classical trajectories. Such an approximation has been successfully applied for studying a short-time evolution of QS with different dynamic symmetry groups [13,14,15,16,17,18,19,20,21,22,23,24,25,26,27]. Expectably, the TWA fails to describe the non-harmonic evolution beyond the Ehrenfest (or semiclassical) time ${\tau }_{sem}$, [28,29,30,31,32], where intrinsic quantum correlations effects start to emerge. The semi-classical time heavily depends both on the Hamiltonian and on the initial state (in particular on the stability of the classical motion). For instance, for spin S systems, the semiclassical time usually scales as the inverse power of the spin length, $gt\lesssim S^{-\alpha }$, $\alpha >0$, where g is a constant characterizing the non-harmonic dynamics. The semiclassical time for variable spin systems behaves in a similar way, where the effective spin size is proportional to the average total angular momentum.

Multiple attempts to improve the TWA [19,33,34] suggest that more sophisticated methods [35,36,37,38,39,40,41,42,43] should be applied in order to describe the quantum evolution in terms of continuously distributed classical phase-space trajectories. Alternatively, different types of discrete phase-space sampling were proposed [44,45,46,47,48,49,50] in order to emulate evolution of average values using the main idea of the TWA. In general, phase-space discretization is a tricky question, which has been addressed from different perspectives mainly focusing on the flat and torus phase-space manifolds [51,52].

It is worth noting that the phase-space analysis of some variable spin systems, such as, for example, two large interacting spins, can be formally performed using the Schwinger representation. Then, faithful mapping onto a flat $C^2\otimes C^2$ phase-space is carried out by applying the standard Heisenberg-Weyl $H\left(1\right)\times H\left(1\right)$ map [53,54,55,56]. However, the natural $SU\left(2\right)$ symmetry is largely lost in such an approach. This is reflected, in the fact that (a) the corresponding distributions are not covariant under rotations; and (b) the inverse excitation numbers in each of the boson modes play the role of dynamical semi-classical parameters. The explicit time-dependence of such semi-classical parameters may restrict the validity of the formal division of the Moyal equation to the classical part, containing only the Poisson brackets and the so-called “quantum corrections”. All of this leads to inefficiency of the standard semiclassical methods [19]. The map of the same systems onto the ${\mathcal{S}}_2\otimes {\mathcal{S}}_2$-spheres in the Stratonovich–Weyl framework [57,58,59] reveals only the local $SU\left(2\right)$ symmetry. The semiclassical parameters are the inverse spin lengths (constant in time). Thus, the $SU\left(2\right)$ TWA leads, in principle, to better results than its flat counterpart in $C^2\otimes C^2$. However, the standard discretization of the ${\mathcal{S}}_2$ sphere [60] does not lead to a significant improvement of the TWA, since the location of the initial distribution is not taken into account.

The situation is even more intricate when the number of involved invariant subspaces becomes formally infinite (or physically very large), as, for example, in the case of a highly-excited rigid rotor interacting with external fields. The use of the Schwinger representation leads to substantial complications in both the analytical and numerical calculations and the standard $SU\left(2\right)$ map is simply non-applicable. Thus, the analysis of the semi-classical limit becomes very challenging [61,62].

In the present paper, we show that there is a natural discretization of $T^{*}{\mathcal{S}}_2$ in the vicinity of the initial distribution, which allows the time-scale of validity of the TWA to be significantly increased, including the so-called revival times. Such a discretization is based on the asymptotic form of the star-product for spin-variable systems [12] and is directly applicable to calculations of the evolution of mean values of physical observables. Using our approach, we will be able to describe the long-time dynamics of molecule in an external field (modelled by a physical rotor) and a coupled two-spin system in the semi-classical limit. We restrict our study to quantum systems with $SO\left(3\right)$ symmetry, corresponding to integer spins.

The paper is organized as follows: In Section 2, we recall the basic results on the $SO\left(3\right)$ covariant mapping for variable spin systems. In Section 3, we discuss the asymptotic form of the star-product in the semi-classical limit. In Section 4, we develop a discretization scheme on the classical manifold of variable spin systems and apply it to computation of mean values in a “quantized” version of the TWA. Two applications of the proposed method with the corresponding numerical solutions are discussed in Section 5. A summary and conclusions are given in Section 6.

2. Variable Spin Wigner Function

Let us consider a QS whose states are elements of a Hilbert space $\mathbb{H}$ containing multiple $SO\left(3\right)$ irreps, so that the $SO\left(3\right)$ group, in general, does not act irreducibly on the density matrix of the system $\widehat{\rho }$, i.e.,

$$\widehat{\rho }=\sum _{L,L^{\prime }=0,1,..}\sum _{m,m^{\prime }}{c}_{m{m}^{\prime }}^{L{L}^{\prime }}|L,m\rangle \langle L^{\prime },m^{\prime }|.$$

The generalized Wigner-like map [11] from operators acting in $\mathbb{H}$ to a discrete set of functions, later called j-symbols,

$$\widehat{f}\iff \{W_f^j\left(\mathsf{\Theta }\right);j=0,1,\cdots \},$$

(1)

where

$$\mathsf{\Theta }=(\varphi ,\theta ,\psi )\in {\mathcal{S}}_3,0\le \varphi <2\pi ,0\le \theta <\pi ,0\le \psi <2\pi$$

(2)

are the Euler angles, is defined through a trace operation

$$W_f^j\left(\mathsf{\Theta }\right)=\mathrm{Tr}\left(\widehat{f}{\widehat{\omega }}_j\left(\mathsf{\Theta }\right)\right),$$

(3)

where the Hermitian $SO\left(3\right)$ covariant mapping kernels have the form

$${\widehat{\omega }}_j\left(\mathsf{\Theta }\right)=\sum _{K=0}^j\sum _{Q,Q^{\prime }=-K}^K\sqrt{\frac{2K+1}{j+1}}{D}_{Q{Q}^{\prime }}^K\left(\mathsf{\Theta }\right){\widehat{T}}_{KQ}^{\frac{j+Q^{\prime }}{2}\frac{j-Q^{\prime }}{2}},$$

(4)

where $D_{Q{Q}^{\prime }}^K\left(\mathsf{\Theta }\right)$ is the Wigner D-function, $D_{Q{Q}^{\prime }}^K\left(\mathsf{\Theta }\right)=\langle K,Q|e^{-i\varphi {\widehat{l}}_z}{e}^{-i\theta {\widehat{l}}_y}{e}^{-i\psi {\widehat{l}}_z}|K,Q^{\prime }\rangle$, here ${\widehat{l}}_{x,y,z}$ are generators of $SO\left(3\right)$ group, $[{\widehat{l}}_k,{\widehat{l}}_m]=i{\epsilon }_{kmn}{\widehat{l}}_n$,

$${\widehat{T}}_{KQ}^{J{J}^{\prime }}=\sum _{M,M^{\prime }}\sqrt{\frac{2K+1}{2J+1}}{C}_{J^{\prime }{M}^{\prime },KQ}^{JM}\left|J,M\right.〉\left.〈J^{\prime },M^{\prime }\right|,$$

(5)

are tensor operators [63] and $C_{a\alpha ,b\beta }^{c\gamma }$ are the Clebsch–Gordan coefficients. The map (3) is explicitly invertible

$$\begin{array}{ccc}\hfill \widehat{f}& =& \sum _{j=0,1\cdots }^{\infty }{\widehat{f}}_j,\hfill \end{array}$$

(6)

$$\begin{array}{ccc}\hfill {\widehat{f}}_j& =& \frac{j+1}{8{\pi }^2}\int d\mathsf{\Theta }{W}_f^j\left(\mathsf{\Theta }\right){\widehat{\omega }}_j\left(\mathsf{\Theta }\right),\hfill \end{array}$$

(7)

where $d\mathsf{\Theta }=sin\theta d\varphi d\theta d\psi$ is a volume element of $SO\left(3\right)$, leading to the overlap relation

$$\mathrm{Tr}\left(\widehat{f}\widehat{g}\right)=\sum _{j=0,1,2,\cdots }^{\infty }\frac{j+1}{8{\pi }^2}\int d\mathsf{\Theta }{W}_f^j\left(\mathsf{\Theta }\right)W_g^j\left(\mathsf{\Theta }\right).$$

(8)

It is worth noting that the operators ${\widehat{f}}_j$ correspond to the expansion of $\widehat{f}$ on the tensor operators ${\widehat{T}}_{KQ}^{J{J}^{\prime }}$ (5) in the sectors with fixed values of $j=J^{\prime }+J$:

$$\begin{array}{ccc}\hfill {\widehat{f}}_j& =& \sum _{K=0}^j\sum _{Q,Q^{\prime }=-K}^K{\widehat{T}}_{KQ}^{\frac{j+Q^{\prime }}{2}\frac{j-Q^{\prime }}{2}}{f}_{KQ}^{\frac{j+Q^{\prime }}{2}\frac{j-Q^{\prime }}{2}},\hfill \end{array}$$

(9)

$$\begin{array}{ccc}\hfill f_{KQ}^{\frac{j+Q^{\prime }}{2}\frac{j-Q^{\prime }}{2}}& =& Tr\left(\widehat{f}{\widehat{T}}_{KQ}^{\frac{j+Q^{\prime }}{2}\frac{j-Q^{\prime }}{2}†}\right).\hfill \end{array}$$

(10)

It should be noted that $Q^{\prime }$ in (4), running over even or odd integers depending on the parity of the index j, such that the restriction $Q^{\prime }\pm j$ is an even number, is fulfilled. This leads to the following symmetry properties of the kernel:

$$\begin{array}{ccc}\hfill {\widehat{\omega }}_j(\varphi ,\theta ,\psi )& =& {\widehat{\omega }}_j(\varphi ,\theta ,\psi +\pi ),\mathrm{even}j,\hfill \end{array}$$

(11)

$$\begin{array}{ccc}\hfill {\widehat{\omega }}_j(\varphi ,\theta ,\psi )& =& -{\widehat{\omega }}_j(\varphi ,\theta ,\psi +\pi ),\mathrm{odd}j.\hfill \end{array}$$

(12)

The advantage of the generalized map (1)–(3) consists of the possibility of a “classical” representation of the whole operator acting in $\mathbb{H}$ and not only its projections on the $SO\left(3\right)$ irreducible subspaces. For instance, for the orientation operators, $\widehat{\mathbf{r}}$,

$$\begin{array}{ccc}\hfill \widehat{\mathbf{r}}& =& \int d\phi sin\vartheta d\vartheta \mathbf{n}(\phi ,\vartheta )|\phi ,\vartheta \rangle \langle \phi ,\vartheta |,{\widehat{\mathbf{r}}}^2=\widehat{I},\hfill \end{array}$$

(13)

$$\begin{array}{ccc}\hfill |\phi ,\vartheta \rangle & =& \sum _{J=0,1,\cdots }\sum _{M=-J}^J{Y}_{JM}^{*}(\phi ,\vartheta )|J,M\rangle ,\hfill \end{array}$$

(14)

where $\phi$, $\vartheta$ are angles in the configuration space, $\mathbf{n}(\phi ,\vartheta )=\left(cos\phi sin\vartheta ,sin\phi sin\vartheta ,cos\vartheta \right)$, one obtains

$$W_{r_k}^j\left(\mathsf{\Theta }\right)=r_k\sum _{n\in {\mathbb{Z}}^{+}}{\delta }_{j,2n+1},$$

(15)

where $\mathbf{r} = \left(sin\varphi sin\psi -cos\varphi cos\theta cos\psi ,-cos\varphi sin\psi -sin\varphi cos\theta cos\psi ,sin\theta cos\psi \right)$, ${\mathbf{r}}^2=1$. The dependence of the symbol on the angle $\psi$ indicates that the corresponding operator mixes $SO\left(3\right)$ invariant subspaces. Vice-versa, symbols of the operators that preserve each $SO\left(3\right)$ irreducible subspace are independent of $\psi$, as, for instance, the angular momentum operators, $\widehat{\mathbf{l}}=({\widehat{l}}_x,{\widehat{l}}_y,{\widehat{l}}_z)$,

$$W_{l_k}^j\left(\mathsf{\Theta }\right)=\sqrt{\frac{j}{2}\left(\frac{j}{2}+1\right)}{n}_k(\varphi ,\theta )\sum _{n\in {\mathbb{Z}}^{+}}{\delta }_{j,2n},$$

(16)

where $\mathbf{n}(\varphi ,\theta )=\left(cos\varphi sin\theta ,sin\varphi sin\theta ,cos\theta \right)$ is a unitary vector in the parameter space (2), and, as a consequence,

$$W_{{\mathbf{l}}^2}^j\left(\mathsf{\Theta }\right)=\frac{j}{2}\left(\frac{j}{2}+1\right)\sum _{n\in {\mathbb{Z}}^{+}}{\delta }_{j,2n},$$

(17)

where ${\widehat{\mathbf{l}}}^2$ is the square angular momentum operator.

The standard Stratonovich–Weyl kernel ${\widehat{w}}_L(\varphi ,\theta )$ [57,58,59], used for mapping operators acting in a single $SO\left(3\right)$ subspace of dimension $2L+1=j+1$ (L is an integer),

$$\widehat{f}\iff W_f\left(\varphi ,\theta \right)=Tr\left(\widehat{f}{\widehat{w}}_{L=j/2}(\varphi ,\theta )\right),$$

(18)

is recovered from the generalized kernel ${\widehat{\omega }}_j\left(\mathsf{\Theta }\right)$ (4) by integrating over the angle $\psi$ (for even values of j):

$${\widehat{w}}_{L=j/2}(\varphi ,\theta )={\int }_0^{2\pi }\frac{d\psi }{2\pi }{\widehat{\omega }}_j\left(\mathsf{\Theta }\right)=\sqrt{\frac{4\pi }{2L+1}}\sum _{K=0}^{2L}\sum _{Q=-K}^K{Y}_{KQ}^{*}(\varphi ,\theta ){\widehat{T}}_{KQ}^L$$

(19)

where $Y_{KQ}(\varphi ,\theta )$ are spherical harmonics and ${\widehat{T}}_{KQ}^L$ are the standard (diagonal) tensor operators [64,65].

It is important to stress that the kernel (4) is not reduced to the direct product of the standard $SO\left(3\right)$ kernels (19). Therefore, the map (1), possessing the underlying global $SO\left(3\right)$ symmetry, allows us to faithfully represent operators in the form of c-functions that:

(a) act in two independent $SO\left(3\right)$ irreps, as, for example, a direct product of angular momentum operators ${\widehat{l}}_k^{\left(1\right)}\otimes {\widehat{l}}_m^{\left(2\right)}$. It should be observed that an alternative mapping can also be achieved with the kernel ${\widehat{w}}_{L_1}(\varphi ,\theta )\otimes {\widehat{w}}_{L_2}(\varphi ,\theta )$. However, in the latter case, the underlying symmetry group is $SO\left(3\right)\times SO\left(3\right)$. The advantage of one of the map over another is not obvious. It will be shown below that the map (1) admits a natural discretization in the semiclassical limit that significantly improves the range of applicability of the Truncated Wigner Approximation;

(b) mixes all $SO\left(3\right)$ irreps, as, for example, the orientation operator (13). This type of operators cannot be mapped into their classical counterparts in the framework of the standard Stratonovich–Weyl approach (18).

3. Wigner Function Dynamics in the Semi-Classical Limit

The crucial feature of the map (1)–(3) is the possibility of introducing a star-product operator [5,66], acting on j-symbols [11]:

$$W_{fg}^j\left(\mathsf{\Theta }\right)=\sum _{j_1,j_2=0,1,\cdots }{L}^{j,j_1{j}_2}\left(W_f^{j_1}\left(\mathsf{\Theta }\right)W_g^{j_2}\left(\mathsf{\Theta }\right)\right),$$

(20)

which is reduced to the standard (local) form [67,68] when the operators $\widehat{f}$ and $\widehat{g}$ are operators from the $se\left(3\right)$ enveloping algebra. The exact form of the star-product operator in general is non-local on the index j and has an involved form, but it is significantly simplified in the limit $j\gg 1$ [12,27],

$$\begin{array}{ccc}\hfill L^{j,j_1{j}_2}& \approx & \mathcal{V}{\int }_0^{2\pi }\frac{d\phi d{\phi }^{\prime }}{{\left(2\pi \right)}^2}{e}^{i\left(j_2-j+{\mathbb{J}}^0\otimes I\right)\phi }{e}^{i\left(j_1-j-I\otimes {\mathbb{J}}^0\right){\phi }^{\prime }},\hfill \end{array}$$

(21)

$$\begin{array}{ccc}\hfill \mathcal{V}& =& e^{-\epsilon {\mathbb{J}}^0\otimes {\mathbb{J}}^0-\frac{\epsilon }{2}\left({\mathbb{J}}^{+}\otimes {\mathbb{J}}^{-}-{\mathbb{J}}^{-}\otimes {\mathbb{J}}^{+}\right)},\hfill \end{array}$$

(22)

where $\epsilon ={(j+1)}^{-1}$,

$${\mathbb{J}}^{\pm }=i{e}^{\mp i\psi }\left[i\frac{\partial }{\partial \theta }\pm cot\theta \frac{\partial }{\partial \psi }\mp \frac{1}{sin\theta }\frac{\partial }{\partial \varphi }\right]{\mathbb{J}}^0=-i\frac{\partial }{\partial \psi }$$

(23)

and the notation $A\otimes B$ means,

$$\left(A\otimes B\right)\left(W_f^{j_1}{W}_g^{j_2}\right)=\left(A{W}_f^{j_1}\right)\left(B{W}_g^{j_2}\right).$$

(24)

Explicitly applying Equation (21) to the symbols of operators $\widehat{f}$ and $\widehat{g}$ and performing an integration and summation, one obtains the following symbolic expression for the symbol of their product

$$W_{fg}^j\left(\mathsf{\Theta }\right)\approx \mathcal{V}\left(W_f^{j+I\otimes {\mathbb{J}}^0}\left(\mathsf{\Theta }\right)W_g^{j-{\mathbb{J}}^0\otimes I}\left(\mathsf{\Theta }\right)\right),$$

(25)

where the operator indices of each symbol are applied to the right or to the left according to (24). The above expression can be further simplified in the limit $j\gg 1$ and considering j as a continuous variable (see Section 4). However, the continuous limits for symbols $W_f^j\left(\mathsf{\Theta }\right)$ are different for even and odd values of the index j due to the parity property (11) and (12). It is convenient to introduce the linear combinations

$$\begin{array}{ccc}\hfill W_f^{j+}\left(\mathsf{\Theta }\right)& =& W_f^j\left(\mathsf{\Theta }\right)+W_f^{j+1}\left(\mathsf{\Theta }\right),\hfill \end{array}$$

(26)

$$\begin{array}{ccc}\hfill W_f^{j-}\left(\mathsf{\Theta }\right)& =& {\left(-1\right)}^j\left(W_f^j\left(\mathsf{\Theta }\right)-W_f^{j+1}\left(\mathsf{\Theta }\right)\right),\hfill \end{array}$$

(27)

which are related through a phase shift,

$$W_f^{j-}(\varphi ,\theta ,\psi )=W_f^{j+}(\varphi ,\theta ,\psi +\pi ).$$

(28)

For instance, $W_{r_k}^{j+}\left(\mathsf{\Theta }\right)=r_k,$ for any (integer) value of the index j. The symbols (26) and (27) become smooth functions of j in the continuous limit, $W_f^{j\pm }\left(\mathsf{\Theta }\right)\to W_f^{\pm }\left(\mathsf{\Theta },j\right)$.

Of particular interests are symbols $W^{j\pm }\left(\mathsf{\Theta }\right)$ with index j distributed in a vicinity, $1\ll \delta \ll j_0,$ of some $j_0$. Actually, smooth and localized functions of $\left(\mathsf{\Theta },j\right)$, with $\delta \sim j_0^{1/2}$ for $j_0\gg 1$ correspond to the so-called semi-classical states. Physically, such states are spread among several $SO\left(3\right)$ invariant subspaces characterized by a large value of spin and localized in angle variables.

The Schrodinger equation

$$i{\partial }_t\widehat{\rho }=[\widehat{H},\widehat{\rho }],$$

$\widehat{H}$ being the Hamiltonian of the system, is mapped into the evolution equations for the Wigner functions

$$i{\partial }_t{W}_{\rho }^j=W_{H\rho }^j-W_{\rho H}^j.$$

(29)

In the continuous limit and for initial semi-classical states, Equation (29) is reduced in the leading order on $j_0$ to the Liouville-type differential equations [12] (see also Appendix A) for $W_{\rho }^{\pm }\left(\mathsf{\Theta },j\right)$,

$${\partial }_t{W}_{\rho }^{\pm }\left(\mathsf{\Theta },j\right)\approx 2\{W_H^{\pm }\left(\mathsf{\Theta },j\right),W_{\rho }^{\pm }\left(\mathsf{\Theta },j\right)\},$$

(30)

where $\{.,.\}$ are the Poisson brackets in the Darboux coordinates $\left(\left(j+1\right)cos\theta ,\varphi \right)$ and $\left(j,\psi \right)$, and $W_H^{\pm }$ are the corresponding symbols of the Hamiltonian. The above equation defines a classical evolution on the symplectic manifold isomorphic to the cotangent bundle $T^{*}{\mathcal{S}}_2$, which corresponds to the co-adjoint orbit of the $SE\left(3\right)$ group fixed by the values of the Casimir operators ${\widehat{\mathbf{r}}}^2=\widehat{I}$ and $\widehat{\mathbf{l}}·\widehat{\mathbf{r}}=0$. Thus, in the semi-classical limit, the Wigner functions $W_{\rho }^{\pm }(\mathsf{\Theta },j)$ can be considered as distributions in a four-dimensional manifold $T^{*}{\mathcal{S}}_2$ and Equation (30) determines “classical trajectories” $({\mathsf{\Theta }}^{cl}\left(t\right),j^{cl}\left(t\right))$ for variable spin systems. A classical observable f can be associated either with $W_f^{+}\left(\mathsf{\Theta },j\right)$ or $W_f^{-}\left(\mathsf{\Theta },j\right)$; however, it is more convenient to choose $W_f^{+}\left(\mathsf{\Theta },j\right)$ due to the relation (28). The explicit form of the Poisson brackets on $T^{*}{\mathcal{S}}_2$ is given in Appendix A, Equation (A6). Strictly speaking, the real expansion parameter, used in transition from (29) to (30), is $j{\left(t\right)}^{-1}$ (for j initially localized close to $j_0\gg 1$). Therefore, the Liouville Equation (30) holds, while, on average over the distribution, $j\left(t\right)\sim {〈{\widehat{\mathbf{l}}}^2〉}^{1/2}\gg 1$.

The evolution of average value of an operator $\widehat{f}$ evaluated according to the overlap relation (8) is convenient to rewrite in terms of symbols $W_f^{j\pm }\left(\mathsf{\Theta }\right)$, as follows:

$$〈\widehat{f}\left(t\right)〉=\sum _{j=0,1,2,\cdots }\frac{j+1}{16{\pi }^2}\int d\mathsf{\Theta }{W}_f^{j+}\left(\mathsf{\Theta }|t\right)\left(W_{\rho }^{j+}\left(\mathsf{\Theta }\right)+{\left(-1\right)}^j{W}_{\rho }^{j-}\left(\mathsf{\Theta }\right)\right),$$

(31)

where $W_f^{j+}\left(\mathsf{\Theta }|t\right)$ is the symbol of the Heisenberg operator $\widehat{f}\left(t\right)$. In the continuous limit, changing summation on j to integration, the contribution of the second term in the above equation becomes negligible. Then, in the spirit of TWA, considering $(\mathsf{\Theta },j)$ as dynamic variables, the average values of physical observables can be estimated according to (8)

$$\langle \widehat{f}\left(t\right)\rangle \approx {\int }_0^{\infty }dj\frac{j+1}{16{\pi }^2}\int d\mathsf{\Theta }{W}_f^{+}({\mathsf{\Theta }}^{cl}\left(t\right),j^{cl}\left(t\right))W_{\rho }^{+}\left(\mathsf{\Theta },j\right),$$

(32)

where the integration is carried out on the initial conditions of the classical trajectories $j^{cl}\left(t\right)=j^{cl}(\mathsf{\Theta },j|t)$, ${\mathsf{\Theta }}^{cl}\left(t\right)={\mathsf{\Theta }}^{cl}(\mathsf{\Theta },j|t)$.

Actually, since it is expected that $W_{\rho }^{+}(\mathsf{\Theta },j)$ is sharply localized at $j\sim j_0\gg 1$ (commonly the width of the semiclassical Wigner distribution is $\sim \sqrt{j_0}$), Equation (32) can be well approximated as

$$\langle \widehat{f}\left(t\right)\rangle \approx {\int }_{-\infty }^{\infty }dj\frac{j+j_0+1}{16{\pi }^2}\int d\mathsf{\Theta }{W}_f^{+}\left(\mathsf{\Theta },j+j_0|t\right)W_{\rho }^{+}(\mathsf{\Theta },j+j_0),$$

(33)

where

$$W_f^{+}\left(\mathsf{\Theta },j+j_0|t\right)\equiv W_f^{+}\left({\mathsf{\Theta }}^{cl}(\mathsf{\Theta },j+j_0|t),j^{cl}(\mathsf{\Theta },j+j_0|t)\right),$$

and $W_{\rho }^{+}(\mathsf{\Theta },j+j_0)$ is now centered at zero in the j axis. In addition, the semiclassical distributions $W_{\rho }^{j\pm }\left(\mathsf{\Theta }\right)$ are approximately normalized

$$\sum _{j=0,1,2,\dots }^{\infty }\frac{j+1}{16{\pi }^2}\int d\mathsf{\Theta }{W}_{\rho }^{j\pm }\left(\mathsf{\Theta }\right)\approx 1.$$

As well as in cases of lower dimensional manifolds, it cannot be expected that propagation along distinguishable classical trajectories (there is no trajectory crossing), originated at every phase-space point of the initial distribution, are able to describe a long-time non-harmonic dynamics [69,70,71]. However, as it will be shown below, there is a natural form to improve the time-validity of TWA for initial semi-classical states of variable spin systems.

It is worth emphasizing that only $W_{\rho }^{\pm }\left(\mathsf{\Theta },j\right)$ combinations satisfy the Liouville evolution Equation (30). Thus, the symbols $W_f^{+}\left(\mathsf{\Theta },j|t\right)$ should be associated with the evolving classical observables according to (31) and (32).

4. Asymptotic Quantization and Discretization Procedure

We start noting that the integral (32) may describe only a destructive dynamic interference corresponding to the initial stage of non-harmonic quantum evolution (collapse time). Such a behavior is due to a continuous superposition of independent classically propagated infinitesimally close fractions of the initial distribution. Several discretization procedures [51,52] have been proposed in order to overcome this problem and to extend the time validity of TWA. Among them, one can mention a semi-classical discretization based on propagating a single trajectory out of each Plank cell in a flat phase-space [44], application of the discrete Wigner function method [72] to a collection of $1/2$-spin systems [45,46,47,48,49,50].

The form of phase-space evaluation of average values in variable spin system (31) suggests an intuitive way for a partial discretization of the initial distribution in the semi-classical limit. In our approximation, the discrete index j, which originally appeared in the expansion (9) for labeling the angular momentum sectors, is now considered as a classical dynamic variable. Then, having applied the star-product (25) to the symbols (26) and (27), we immediately arrive in the continuous limit at the following asymptotic form of the star-product

$$\begin{array}{ccc}\hfill W_{fg}^{\pm }\left(\mathsf{\Theta },j\right)& \approx & e^{-\epsilon {\mathbb{J}}^0\otimes {\mathbb{J}}^0-\frac{\epsilon }{2}\left({\mathbb{J}}^{+}\otimes {\mathbb{J}}^{-}-{\mathbb{J}}^{-}\otimes {\mathbb{J}}^{+}\right)}{e}^{{\partial }_j\otimes {\mathbb{J}}^0-{\mathbb{J}}^0\otimes {\partial }_j}\left[W_f^{\pm }\left(\mathsf{\Theta },j\right)W_g^{\pm }\left(\mathsf{\Theta },j\right)\right]\hfill \end{array}$$

(34)

$$\begin{array}{ccc}& \approx & W_f^{\pm }\left(\mathsf{\Theta },j\right) \ast W_g^{\pm }\left(\mathsf{\Theta },j\right),\hfill \end{array}$$

(35)

where $e^{{\partial }_j\otimes {\mathbb{J}}^0}f\left(j\right)=f(j+I\otimes {\mathbb{J}}^0)$. It is worth noting that the star-product in the form (34) is applicable only to the classical observables, i.e., (26) and (27) symbols.

Strictly speaking, Equation (34) should be applied to functions with the index j distributed in a broad vicinity, $\delta \gg 1,$ of some $j_0\gg 1$, so that $\epsilon \approx {(j_0+1)}^{-1}$ and ${\partial }_j\sim {\delta }^{-1}$. However, a direct application of Equation (34) to the generators of the $se\left(3\right)$ algebra ($\widehat{\mathbf{l}}$, $\widehat{\mathbf{r}}$) shows a good correspondence between the exact commutation relations and their counterparts reconstructed through the asymptotic star-product:

$$\begin{array}{ccc}\hfill W_{[l_k,l_n]}^{+}\left(\mathsf{\Theta },j\right)& =& i{ϵ}_{knm}{W}_{l_m}^{+}\left(\mathsf{\Theta },j\right)+O\left(\epsilon \right)\hfill \\ \hfill W_{[l_k,r_n]}^{+}\left(\mathsf{\Theta },j\right)& =& i{ϵ}_{knm}{W}_{r_m}^{+}\left(\mathsf{\Theta },j\right)+O\left(\epsilon \right),W_{[r_k,r_n]}^{+}\left(\mathsf{\Theta },j\right)=0,\hfill \end{array}$$

where

$$W_{[f,g]}\left(\mathsf{\Theta },j\right)=W_f\left(\mathsf{\Theta },j\right) \ast W_g\left(\mathsf{\Theta },j\right)-W_g\left(\mathsf{\Theta },j\right) \ast W_f\left(\mathsf{\Theta },j\right).$$

(36)

It is worth noting that asymptotically the operator corresponding to the classical observable $(j+j_0)/2$ is conjugated to the operator corresponding to the phase $e^{i\psi }$,

$$W_{[(j+j_0)/2,e^{i\psi }]}^{+}\left(\mathsf{\Theta },j\right)\approx W_{e^{i\psi }}^{+}\left(\mathsf{\Theta },j\right).$$

(37)

In a sense, $(j+j_0)/2$ and $e^{i\psi }$ can be seen as action-angle variables. Here, we consider $\left(j+1\right)/2$ as a physical variable representing the classical angular momentum, as $W_{l_k}^{+}(j,\mathsf{\Theta })\approx \frac{j+1}{2}{n}_k$, according to (16) and (17).

Making use of the asymptotic form of the star-product (34), we can discretize back the variable j following the ideas of deformation quantization. According to the general procedure [66], we solve the eigenvalue equation

$$i{\partial }_{\tau }U=\left(j+j_0\right) \ast U\approx \left(j+j_0-i{\partial }_{\psi }\right)U,$$

$$U\left(\left.j+j_0,\mathsf{\Theta }\right|\tau =0\right)=1,$$

where the expression (34) was employed. A direct expansion of the solution

$$U\left(\left.j+j_0,\mathsf{\Theta }\right|\tau \right)=e^{-i\left(j+j_0\right)\tau },$$

in the Fourier series yields

$$e^{-i\left(j+j_0\right)\tau }=\sum _{L=-\infty }^{\infty }{\mathsf{\Pi }}_L\left(j\right)e^{m\pi i\left(j+j_0\right)}{e}^{-i\left(L+j_0\right)\tau },\left(2m-1\right)\pi <\tau <\left(2m+1\right)\pi ,$$

where

$${\mathsf{\Pi }}_L\left(j\right)=\frac{sin\pi \left(j-L\right)}{\pi \left(j-L\right)}.$$

(38)

Since the Wigner distributions $W_{\rho }^{\pm }\left(\mathsf{\Theta },j+j_0\right)$ have compact supports, with width $\sim \sqrt{j_0}\gg 1$, we can make use of the Whittaker–Shannon–Kotelnikov sampling theorem [73,74,75], approximating

$$\begin{array}{ccc}& & \frac{j+j_0+1}{16{\pi }^2}{W}_f^{+}\left(\mathsf{\Theta },j+j_0|t\right)W_{\rho }^{\pm }\left(\mathsf{\Theta },j+j_0\right)\hfill \\ & \approx & \sum _{L=-\infty }^{\infty }\frac{L+j_0+1}{16{\pi }^2}{W}_f^{+}\left(\mathsf{\Theta },L+j_0|t\right)W_{\rho }^{\pm }\left(\mathsf{\Theta },L+j_0\right){\mathsf{\Pi }}_L\left(j\right).\hfill \end{array}$$

(39)

It is worth noting that, in case of Gaussian function of width $\sim \sqrt{j_0}$, the error of discrete sampling with (38) is of order ∼erfc$\left(\pi \sqrt{j_0}\right)$ [76].

However, a direct discretization (39) of the semi-classical expression (33) is not sufficient for an efficient simulation of the quantum dynamics through phase-space trajectories since the overlap between the classically evolved observable $W_f^{j+}\left(\mathsf{\Theta }|t\right)$ and the branch $W_{\rho }^{j-}\left(\mathsf{\Theta }\right)$ of distribution would be missed. It is worth recalling that $W_{\rho }^{+}\left(\mathsf{\Theta },j\right)$ and $W_{\rho }^{-}\left(\mathsf{\Theta },j\right)$ have maxima at different points of the phase space (shifted in $\pi$ on $\psi$). In order to correct this problem, we rewrite Equation (31)

$$\begin{array}{ccc}\hfill 〈\widehat{f}\left(t\right)〉& =& {\Sigma }_{+}\left(t\right)+\mathcal{P}{\Sigma }_{-}\left(t\right)\hfill \end{array}$$

(40)

$$\begin{array}{ccc}\hfill {\Sigma }_{\pm }\left(t\right)& =& \sum _{j=0,1,2,\cdots }\frac{j+1}{16{\pi }^2}\int d\mathsf{\Theta }{W}_f^{j+}\left(\mathsf{\Theta }|t\right)W_{\rho }^{j\pm }\left(\mathsf{\Theta }\right),\hfill \end{array}$$

(41)

where $\mathcal{P}$ is the parity operator defined according to

$$\mathcal{P}\sum _{j=0,1,2,..}{a}_j=\sum _{j=0,1,2,..}{(-1)}^j{a}_j.$$

Then, considering the semi-classical evolution of $W_f^{+}\left(\mathsf{\Theta },j\right)$ in the continuous limit and applying a subsequent discretization procedure (39)–(41), we arrive at the following discretized version of TWA for variable spin systems,

$$\langle \widehat{f}\left(t\right)\rangle \approx \sum _{L=-\infty }^{\infty }\frac{L+j_0+1}{16{\pi }^2}\int d\mathsf{\Theta }{W}_f^{+}\left(\mathsf{\Theta },L+j_0|t\right)\left(W_{\rho }^{+}\left(\mathsf{\Theta },L+j_0\right)+{(-1)}^L{W}_{\rho }^{-}\left(\mathsf{\Theta },L+j_0\right)\right).$$

(42)

The above equation is just a convolution of the evolving classical observable with a linear combination of the initial distribution $W_{\rho }^{+}$ at $(\varphi ,\theta ,\psi )$ and $(\varphi ,\theta ,\psi +\pi )$, evaluated at the equiseparated points along the “action” variable j. A naive direct discretization of the continuous approximation (32) leads to an incomplete description of the quantum dynamics in the semiclassical limit.

5. Examples

5.1. Rigid Rotor in an External Field

Let us consider the following Hamiltonian governing the evolution of quantum rigid rotor in an external field along the z-axis,

$$\widehat{H}={\widehat{\mathbf{l}}}^2-g{\widehat{z}}^2,$$

(43)

where $\widehat{z}$ is the z-component of the orientation operator (13), and ${\widehat{\mathbf{l}}}^2$ is the square angular momentum operator. The Hamiltonian (43) possesses the $SO\left(3\right)$ symmetry but cannot be reduced to a finite number of spin systems.

The symbol of the Hamiltonian in the continuous limit at the principal order on j has the form

$$W_H^{+}(\mathsf{\Theta },j)\approx \frac{j}{2}\left(\frac{j}{2}+1\right)-g{sin}^2\theta {cos}^2\psi$$

(44)

and leads to the following equations of motion:

$$\begin{array}{ccc}\hfill {\partial }_t\varphi & =& \frac{4g}{j+1}cos\theta {cos}^2\psi \hfill \\ \hfill {\partial }_t\theta & =& -\frac{g}{j+1}sin2\theta sin2\psi \hfill \\ \hfill {\partial }_t\psi & =& j+1-\frac{4g}{j+1}{cos}^2\theta {cos}^2\psi \hfill \\ \hfill {\partial }_{t}j& =& -2g{sin}^2\theta sin2\psi .\hfill \end{array}$$

(45)

As an initial state, we consider a weighted superposition of l-spin coherent states $\left|l;{\vartheta }_0,{\phi }_0\right.〉$

$$\left|\mathsf{\Psi }\left(0\right)\right.〉=\frac{1}{\sqrt{cosh{r}^2}}\sum _{l=0}^{\infty }{e}^{-il{\psi }_0}\frac{r^{2l}}{\sqrt{\left(2l\right)!}}\left|l;{\vartheta }_0,{\phi }_0\right.〉,$$

(46)

where $r^2\gg 1$. The Wigner distribution $W_{\rho }^{+}\left(\mathsf{\Theta },j\right)$ corresponding to the state (46) can be approximated as

$$\begin{array}{ccc}\hfill W_{\rho }^{+}\left(\mathsf{\Theta },j\right)& \approx & \frac{r^{2j}}{\mathsf{\Gamma }\left(j+2\right)cosh{r}^2}\frac{1}{sin\omega /2}{\partial }_{\omega }\left(\frac{{sin}^2\left(\left(j+1\right)\omega /2\right)}{sin\omega /2}\right),\hfill \\ \hfill cos\frac{\omega }{2}& =& \frac{1}{\sqrt{2}}\left(cos\frac{\theta }{2}cos\frac{\varphi +\psi }{2}+sin\frac{\theta }{2}cos\frac{\varphi -\psi }{2}\right),\hfill \end{array}$$

and is localized in $\vartheta \sim {\vartheta }_0=\pi /2$, $\phi \sim {\phi }_0=0$ by construction, with absolute fluctuations $\delta \vartheta \sim \delta \phi \sim r^{-1}$. In addition, $W_{\rho }^{+}\left(\mathsf{\Theta },j\right)$ is localized at $j\sim r^2$ with the fluctuation $\delta j\sim r$. Thus, the phase $\psi$ is also localized with $\delta \psi \sim r^{-1}$, as it is an observable conjugated to j, according to (37), see Figure 1a where the marginal distribution

$$W_{\rho }^{+}\left(\psi ,j\right)=\frac{j+1}{8\pi }{\int }_0^{2\pi }d\varphi {\int }_0^{\pi }sin\theta d\theta W_{\rho }^{+}\left(\mathsf{\Theta },j\right)$$

(47)

is plotted. The same reasoning is applicable to $W_{\rho }^{-}\left(\mathsf{\Theta },j\right)$ due to the relation (28). This means that the state (46) is a semi-classical state on the manifold $T^{*}{\mathcal{S}}_2$.

We have numerically tested the approximation (42) by computing the averages of $\widehat{z}\left(t\right)$, ${\widehat{z}}^2\left(t\right)$, and ${\widehat{\mathbf{l}}}^2\left(t\right)$. For numerical simulations, we have used an adaptive sampling technique that allows us to sample the angular variables in the regions where the initial Wigner function is located, for every fixed value of j. The sampling method is based on an algorithm [77] for numerical integration of multivariate functions and implemented by a modification of the routine cubature [78]. For $r^2=81$, $g=20$, we take on average 7892 points inside $S_3$ sphere for each value of integer $j\in \left[40,160\right]$. The errors obtained in the estimation of the expectation values are lower than $0.0045\%$ at $t=0$. The differential equations were solved by a variable step-size Runge–Kutta method of order $9\left(8\right)$ [79]. The size of the step was adapted to keep the relative errors estimated by the method lower than ${10}^{-6}$.

In Figure 1b–d, we plot the corresponding averages in comparison with the exact calculations and the continuous TWA approximation (33). One can appreciate that Equation (33) coincides with the exact calculations only within the initial collapse, i.e., for times $gt\ll 1$. Conversely, Equation (42) describes very well the evolution of the observables inclusively for much longer intervals that include several revivals, $gt\sim \pi$. For even longer times, $gt\sim r^2$, our approximation starts to deviate from the exact solution, failing to capture the oscillation dephasing and deformations of envelopes of the revivals (although the condition $j\left(t\right)\gg 1$ still holds).

It is worth noting that, in contrast to planar pendulum models, described by periodic Hamiltonians of the form

$$\widehat{H}={\widehat{p}}^2+V\left(\widehat{x}\right),V(x+a)=V\left(x\right),a=const,$$

where $\widehat{p}$ and $\widehat{x}$ are the standard momentum and position operators, a phase-space description of the rigid rotor in an external field is not trivial. For instance, the standard phase-space analysis of the Hamiltonian (43) in $C^2\otimes C^2$, and application of the corresponding TWA [19], faces considerable technical difficulties. In particular, the Schwinger (two-mode) representation of $\widehat{z}$ operator,

$$\widehat{z}=\frac{1}{\sqrt{\widehat{a}{\widehat{a}}^{†}+{\widehat{b}}^{†}\widehat{b}}}\left({\widehat{a}}^{†2}+{\widehat{b}}^{†}\right)\frac{1}{\sqrt{\widehat{a}{\widehat{a}}^{†}+{\widehat{b}}^{†}\widehat{b}}},$$

is quite inconvenient for the $H\left(1\right)\times H\left(1\right)$ phase-space mapping [53,54,55,56], and hides the intrinsic $SO\left(3\right)$ symmetry of the direction operator (13) and (14).

5.2. Spin–Spin Interaction

As another non-trivial example, we consider the following Hamiltonian:

$$\widehat{H}={\widehat{\mathbf{l}}}_1·{\widehat{\mathbf{l}}}_2+\lambda {\widehat{l}}_{z1},$$

(48)

describing an integer spin–spin interaction in the presence of an external non-uniform magnetic field.

Within the framework of our approach, the symbol of the Hamiltonian in the continuous limit has the form

$$W_H^{+}(\mathsf{\Theta },j)=C\left(j\right)+\lambda A\left(j\right)sin\theta cos\psi +\lambda B\left(j\right)cos\theta ,$$

where $A\left(j\right)$, $B\left(j\right)$ and $C\left(j\right)$ are functions of j given in Appendix B.

Taking into account the classical equations of motion on $T^{*}{\mathcal{S}}_2$ (see Appendix B), we compute the evolution of the first spin magnetization $〈{\widehat{l}}_{x1}\left(t\right)〉$ according to the general procedure (42), where the symbol of ${\widehat{l}}_{x1}$ is

$$W_{l_{x1}}^{+}\left(\mathsf{\Theta },j\right)=A\left(j\right)\left(sin\varphi sin\psi -cos\varphi cos\theta cos\psi \right)+B\left(j\right)cos\varphi sin\theta .$$

As the initial state, we consider the product of spin coherent states in x- and y-directions,

$$|\mathsf{\Psi }\left(0\right)\rangle =\left|L_1;{\phi }_1=0,{\vartheta }_1=\pi /2\right.〉\otimes \left|L_2;{\phi }_2=\pi /2,{\vartheta }_2=\pi /2\right.〉.$$

(49)

The general expression for the Wigner function of the state (49) is quite cumbersome, but its localization property at $j\sim 2L_1$ for $L_1\gg L_2$ follows from the marginal distribution

$${\int }_0^{2\pi }\frac{d\psi }{2\pi }{W}_{\rho }^{+}\left(\mathsf{\Theta },j\right)\sim e^{-{\left(j/2-L_1\right)}^2/L_2}{W}_{\rho }\left(\varphi ,\theta ;j\right),$$

where $W_{\rho }\left(\varphi ,\theta ;j\right)=\langle L=j/2;\phi =0,\vartheta =\pi /2|{\widehat{w}}_{L=j/2}(\varphi ,\theta )\left|L=j/2;\phi =0,\vartheta =\pi /2\right.〉$ is the standard $SO\left(3\right)$ Wigner function (19). The localization on the angle $\psi$ follows from its complementarity to the variable j, in the same way as in the rotor case. The marginal distribution (47) corresponding to the state (49) is plotted in Figure 2a.

It is worth noting that this system has $SO\left(3\right)\times SO\left(3\right)$ symmetry (for integer spins), so that the operators acting in the Hilbert space of two-spin system can be mapped into distributions in ${\mathcal{S}}_2\times {\mathcal{S}}_2$ using the following the mapping kernel,

$$\widehat{w}({\mathsf{\Omega }}_1,{\mathsf{\Omega }}_2)={\widehat{w}}_1\left({\mathsf{\Omega }}_1\right)\otimes {\widehat{w}}_2\left({\mathsf{\Omega }}_2\right),$$

(50)

where $\widehat{w}\left(\mathsf{\Omega }\right)$, $\mathsf{\Omega }=\left(\varphi ,\theta \right)$ is defined in (19). In the limit of large spins, the Hamiltonian dynamics can be treated semi-classically [20,21,22,23,27]. The equation of motion for the Wigner function of the whole system ${\tilde{W}}_{\rho }({\mathsf{\Omega }}_1,{\mathsf{\Omega }}_2)$ takes the form:

$$\begin{array}{ccc}\hfill {\partial }_t{\tilde{W}}_{\rho }& \approx & \frac{1}{L_1+1/2}{\{{\tilde{W}}_H,{\tilde{W}}_{\rho }\}}_1+\frac{1}{L_2+1/2}{\{{\tilde{W}}_H,{\tilde{W}}_{\rho }\}}_2,\hfill \\ & & \\ \hfill {\tilde{W}}_H({\mathsf{\Omega }}_1,{\mathsf{\Omega }}_2)& =& \sqrt{L_1(L_1+1)L_2(L_2+1)}{\mathbf{n}}_1·{\mathbf{n}}_2+\lambda \sqrt{L_1(L_1+1)}{n}_{z1},\hfill \end{array}$$

where the canonical variables defining the Poisson brackets ${\{..,..\}}_{1,2}$ are $\left(cos{\theta }_{1,2},{\varphi }_{1,2}\right)$, and the average values are computed according to the standard TWA,

$$\langle \widehat{f}\left(t\right)\rangle \approx \frac{(2L_1+1)(2L_2+1)}{16{\pi }^2}\int d{\mathsf{\Omega }}_{1}d{\mathsf{\Omega }}_2{\tilde{W}}_f(({\mathsf{\Omega }}_1^{cl}\left(t\right),{\mathsf{\Omega }}_2^{cl}\left(t\right)){\tilde{W}}_{\rho }({\mathsf{\Omega }}_1,{\mathsf{\Omega }}_2),$$

(51)

${\mathsf{\Omega }}_1^{cl}\left(t\right)$ and ${\mathsf{\Omega }}_2^{cl}\left(t\right)$ being classical trajectories on ${\mathcal{S}}_2\times {\mathcal{S}}_2$.

For numerical simulations, we used the same adaptive method as in the case of the rotor. For the spin–spin system, $L_1=8$, $L_2=4$, $\lambda =0.2$, $j\in \left[8,24\right],$ the average number of samples is 26,885. The errors in the estimation of the expectation values are lower than $0.0008\%$ at $t=0$.

Comparing the results obtained from (42) and (51), one can observe that the proposed discretization leads to a good coincidence with the exact results significantly beyond the validity of the standard TWA in the framework of Stratonovich–Weyl correspondence. Actually, our approximation describes well the effect of partial revivals produced by the nonlinear term ${\widehat{l}}_1·{\widehat{l}}_2$ at the scale $t\sim 2\pi$, but falling at $t\sim 2L_1$ (independently on the external field coupling constant $\lambda$). The standard TWA breaks down already after the first collapse at $t\lesssim 1$.

6. Conclusions

The semi-classical map (1)–(3) of density matrices of a variable spin system into distributions on four-dimensional symplectic manifold allows for approximating the evolution of such quantum systems in terms of effective classical dynamics on $T^{*}{\mathcal{S}}_2$. The advantage of the map (3) with respect to the standard $SU\left(2\right)/SO\left(3\right)$ case [57] consists of the possibility of a faithful representation of operators whose action is not restricted to a single $SU\left(2\right)$ invariant subspace. In addition, only four Hamilton equations are sufficient to determine the evolution of classical observables for any value of the total angular momentum. However, the simplest Truncated Wigner Approximation suffers from the same intrinsic defects as in the Heisenberg–Weyl and $SU\left(2\right)$ symmetries: it describes well only the short-time evolution of the typical observables of the system (which may include, e.g., generators of $SE\left(3\right)$ group). In order to extend the validity of TWA, while still keeping the idea of classical propagation, we propose to “quantize” back one of the classical dynamic variables, j, which can be considered to some extent as an “action” (actually representing possible values of the classical spin size). We perform such a “quantization” by using the asymptotic form of the star-product within the framework of deformation quantization, leading to a natural discretization of the variable j, and, thus, all of the distributions appearing in the theory, corresponding both to states and observables. A certain subtlety of the proposed method consists of taking into account the parity problem originated from the decomposition of the mapping kernel (4) in the basis of the tensor operators (5). In addition, it results that the obtained discretization of initial distributions with compact support, describing the so-called semi-classical states, is in direct accordance with the famous sampling theorem. This allows the form of calculation of average values to be immediately discretized, basically starting the classical trajectories only at certain points of the initial distribution. The result of such an approach is surprisingly good, as shown in Figure 1 and Figure 2. It is worth noting that the discrete sampling procedure in general is not obvious at all. For instance, applying more sophisticated discretization methods, like, e.g., the adaptive discretization, one obtains much worse results than by following the simple recipe (42). Actually, Equation (42) describes very well all interference effects, such as, e.g., revivals of quantum oscillations that appear due to superpositions of subspaces with different values of the index j, appearing in the exact calculations (8).

The range of applicability of the discretized TWA is considerably longer than the standard semiclassical time ${\tau }_{sem}$, $gt\lesssim j_0^{\beta }{\tau }_{sem}$, $\beta >0$, where $j_0\gg 1$ is the average initial total angular moment. For second degree Hamiltonians, similar to (43) and (48), the leading corrections to the Liouville Equation (30) are of the order $j_0^{-1}$, while the principal term (the Poisson bracket on $T^{*}{\mathcal{S}}_2$) is $\sim j_0$. Actually, the Moyal Equation (29) has in this case the following structure

$${\partial }_t{W}_{\rho }^{}=\left(j_0{\mathcal{L}}_0+j_0^{-1}{\mathcal{L}}_2\right)W_{\rho }^{}+O\left(j_0^{-2}\right),$$

(52)

where ${\mathcal{L}}_0$ is a first order differential operator and ${\mathcal{L}}_2$ contains higher-degree derivatives. Dropping the correction terms $j_0^{-1}{\mathcal{L}}_2$ in the continuous TWA leads to neglecting all of the commutators of order $j_0^0=1$ that appear in the exponent of the formal propagator corresponding to Equation (52). As a result, the description of quantum dynamics occurring in the time-scale $gt\gtrsim 1$ is inaccessible in the semiclassical treatment (32). In practice, the standard TWA [19] breaks down even for shorter time intervals and is unable to describe any genuine quantum effect such as, for example, quantum revivals. It seems that the discretization (42) allows physical processes caused by the interferences between different j-sectors (9) to be emulated until $gt\lesssim j_0$. This is really not surprising, since similar effects take place in almost all nonlinear quantum systems with discrete spectra as a result of a specific composition of some discrete constituents [80]. The difficulty consists, as we have mentioned above, in finding an appropriate discretization of the classical phase-space. In Figure 3a,b, the long time evolution of the rotor (43) and spin observables (48) are shown. One can clearly observe the region of time validity of our approximation.

The present approach is also extendable to half-integer spins, although the calculations become more involved. In such a case, four linear combinations of the Weyl symbols $W_f^j\left(\mathsf{\Theta }\right),$ similar to (26) and (27), should be introduced in order to follow the same procedure as in Section 3 and Section 4. This problem will be considered elsewhere in application to dissipative and pumped down conversion processes.

Recently, a generalized TWA was proposed in [81] for the description of the dynamics of many coupled spins. In this approach, every spin is considered as a discrete variable, i.e., there is no relation to any classical phase-space. Thus, a set of $2\left({\left(2L+1\right)}^2-1\right)$ coupled (first-order) differential equations should be solved even for two interacting spins L. In addition, the rigid rotor evolution in external fields cannot be treated applying the technique [81].

Unfortunately, the map (1)–(4) cannot be used as a faithful (one-to-one map) classical representation of the N-spin, $N\ge 3$ system. However, some global properties of multi-spin systems can be analyzed with our method by using the Schur–Weyl duality [82] and averaging over invariant subspaces of the same dimensions.

Finally, we note that the developed approach can, in principle, be applied to any of the s-parametrized maps [12]. However, the Moyal equation for non self-dual distributions $W_{\rho }^{\left(s\right)}$ contains terms of order one (on the semiclassical parameter), i.e., it has the form

$${\partial }_t{W}_{\rho }^{\left(s\right)}=\left(j_0{\mathcal{L}}_0+s{\mathcal{L}}_1+j_0^{-1}{\mathcal{L}}_2\right)W_{\rho }^{\left(s\right)}+O\left(j_0^{-2}\right),$$

where ${\mathcal{L}}_1$ is a differential operator already containing higher than first degree derivatives. Thus, one may expect that the TWA for $W_{\rho }^{(s\ne 0)}$, which considers the action only of ${\mathcal{L}}_0$, is less precise than for $W_{\rho }^{(s=0)}$.