#  Lyapunov spectrum scaling for classical many-body dynamics close to
integrability

Merab Malishava1,2<EMAIL_ADDRESS>Sergej Flach1,2
<EMAIL_ADDRESS>1Center for Theoretical Physics of Complex Systems,
Institute for Basic Science(IBS), Daejeon, Korea, 34126
2Basic Science Program, Korea University of Science and Technology(UST),
Daejeon, Korea, 34113

###### Abstract

We propose a novel framework to characterize the thermalization of many-body
dynamical systems close to integrable limits using the scaling properties of
the full Lyapunov spectrum. We use a classical unitary map model to
investigate macroscopic weakly nonintegrable dynamics beyond the limits set by
the KAM regime. We perform our analysis in two fundamentally distinct long-
range and short-range integrable limits which stem from the type of
nonintegrable perturbations. Long-range limits result in a single parameter
scaling of the Lyapunov spectrum, with the inverse largest Lyapunov exponent
being the only diverging control parameter and the rescaled spectrum
approaching an analytical function. Short-range limits result in a dramatic
slowing down of thermalization which manifests through the rescaled Lyapunov
spectrum approaching a non-analytic function. An additional diverging length
scale controls the exponential suppression of all Lyapunov exponents relative
to the largest one.

††preprint: APS/123-QED

Thermalization is a universal property of the long-time dynamics of generic
nonintegrable many-body systems. Thermal equilibrium is characterized by
stationary distributions and assumes ergodicity and mixing in the phase space
[1]. The thermalization dynamics will in general slow down close to
integrability, and may even cease to be observed [2, 3, 4, 5, 6, 7, 8], which
was also noted in earlier studies of dynamical systems [9, 10, 11, 12]. The
theory of weak nonintegrable perturbations for finite Hamiltonian systems was
pioneered by Kolmogorov in 1954 [13] and later by Arnold [14] and Moser [15].
The corresponding KAM theory demonstrates the violation of the ergodic
hypothesis for sufficiently weak perturbations due to the emergence of a mixed
phase space with a finite fraction of points belonging to regular trajectories
on tori. At a critical strength of the nonintegrable perturbation, all tori
disappear, and the dynamics become fully chaotic allowing for thermalization.
But what exactly is meant by “sufficiently weak” and “critical strength”? As
it turns out the magnitude of the critical perturbation $\tilde{H}$ decays
rapidly with the growth of the number of degrees of freedom. Namely,
$|\tilde{H}|\leq aN^{-b}$ with $b=160$ has been shown as an upper bound for
applicability of KAM theory in lattice systems with short-range interactions
(for example an array of Josephson junctions) [16]. This result suggests that
it is practically impossible to witness quasiperiodic motion suggested by KAM
in macroscopically large systems close to integrability. What is then the
expected behavior of systems with a large number of degrees of freedom in
proximity to an integrable limit? How does one characterize it? Does it have
universality classes? Is there a KAM-like regime for macroscopic models? If
not, what lies beyond the KAM horizon spanned by finite systems?

Figure 1: a) A schematic representation of the unitary circuits map. Large
yellow blocks indicate $\hat{C}$ unitary matrices parametrized by the angle
$\theta$. Small green blocks indicate local nonlinearity generating map
$\hat{G}$ parametrized by the nonlinearity strength $g$. The black arrow on
the right indicates the time flow. b) Control parameter space $\\{\theta,g\\}$
with the highlighted area corresponding to the induced networks. Integrable
limits are reached for $g=0$ (linear evolution of extended normal modes) or
$\theta=0$ (decoupled nonlinear map). Small nonzero $g$ values induce LRNs,
small nonzero $\theta$ values induce SRNs. The network images indicate actions
(filled circles) and couplings induced by the nonintegrable perturbation
(straight lines). Left image - SRN, right image - LRN.  Figure 2: The largest
Lyapunov exponents $\Lambda_{\text{max}}$ in SRN (blue small circles, top) and
LRN (orange small circles, bottom) regime versus the corresponding deviation
from integrable limit $g,\theta$. Solid lines connect the data points and
guide the eye. For the SRN case, the parameter nonlinearity is fixed $g=1.0$,
while for the LRN case the angle $\theta$ is fixed at $0.33\pi$. For both
cases system size $N=200$. The large black circles connected by dashed lines
correspond to data for system size $N=100$.

To quantify the thermalization of a system one typically chooses a specific
set of observables and studies equipartition and ergodicity for those specific
observables. In their pioneering work Fermi, Pasta, Ulam, and Tsingou
attempted to showcase equipartition using the normal modes of a linear chain
as their choice of observables for a weakly nonlinear chain [17]. In the
absence of nonlinearity, the normal modes are “frozen”, i.e. they become the
actions of the integrable system.

Recent studies attempted to broaden the ergodicity analysis by computing
convergence of finite-time average distributions of observables to their phase
space averages. They revealed that most physical systems belong to two
distinct classes when it comes to thermalization in proximity to integrable
limits [18, 19, 20]. Systems with weak nonlinear perturbations such as FPUT
chains, Josephson junction networks in the limit of small energy density,
discrete nonlinear Schrödinger equations all belong to the class of Long Range
Networks (LRN). On the other hand, a broad range of lattice systems allowing
for proximity to an integrable limit of vanishing lattice coupling belongs to
a class of Short Range Networks (SRN). Examples of SRN include coupled
anharmonic oscillator chains in the limit of weak coupling [18], Josephson
junction chains in the limit of weak Josephson coupling [19], etc.

There are serious limitations of studying thermalization through observable
dynamics. The choice of observables is ambiguous [21, 22], and even for
integrable systems specifically chosen observables show ergodic thermal-like
behavior [23]. Observable dynamics address ergodicity, but not mixing.
However, nonintegrable dynamics are necessarily mixing, show typically
exponential decay of correlations with a macroscopic set of correlation times.

In this Letter, we overcome the above limitations by computing the entire
Lyapunov spectrum [24]. Lyapunov spectra were previously used for diagnosing
phase transitions [25] and energy localization [26]. Here we show that the
scaling properties of the Lyapunov spectrum offer a conceptual novel way for
the description of weakly nonintegrable dynamics in a generic model setup. We
consider a macroscopically large system beyond the limits set by KAM and
characterize thermalization in both SRN and LRN regimes, thus drawing a very
general picture that encapsulates a great number of physically realizable
scenarios and is directly applicable to most weakly nonintegrable classical
systems.

Resolving the entire Lyapunov spectrum for a large system is a numerically
challenging task. It relies on the simultaneous evolution of a large number of
trajectories, [27]. The proximity to integrable limits makes this task even
harder due to an increase of the thermalization times. In view of these
challenges, we need models which possess all physically relevant features to
achieve thermalization and are extremely efficient for the numerical evolution
- unitary maps. The fast, exact, error-free discrete-time evolution is a key
feature of unitary maps which makes them advantageous for heavy numerical
tasks. These properties were on display in recent studies, where discrete
unitary maps were used to achieve record-breaking evolution times for
nonlinear wave-packet spreading tasks [28], Anderson localization [29, 30],
and soliton dynamics [31].

Model – We use classical Unitary Circuit maps. We define a 1D lattice of size
$N$ with one complex component $\psi_{n}$ per site $n$. The classical dynamics
evolves the vector $\vec{\Psi}=\\{\psi_{n}\\}$ in a corresponding phase space
of dimension $2N$ on a deterministic trajectory specified by an initial
condition. The evolution is performed by subsequent applications of the map:

$\displaystyle\hat{U}=\sum_{n\in\mathds{Z}}\hat{G}_{n}\sum_{n\in
2\mathds{Z}+1}\hat{C}_{n,n+1}\sum_{n\in 2\mathds{Z}}\hat{C}_{n,n+1}.$ (1)

The unitary matrices $\hat{C}_{n,n+1}$ are parametrized by the rotation angle
$\theta$ and act as hoppings on pairs of neighboring sites:

$\displaystyle\hat{C}_{n,n+1}\begin{pmatrix}\psi_{n}(t)\\\
\psi_{n+1}(t)\end{pmatrix}=\begin{pmatrix}\cos\theta&\sin\theta\\\
-\sin\theta&\cos\theta\end{pmatrix}\begin{pmatrix}\psi_{n}(t)\\\
\psi_{n+1}(t)\end{pmatrix},$

and the local map $\hat{G}_{n}$ induces nonlinearity:

$\displaystyle\hat{G}_{n}\psi_{n}=e^{ig|\psi_{n}|^{2}}\psi_{n}.$ (3)

The classical unitary circuit dynamics is schematically represented in Fig. 1.

The map dynamics of classical unitary circuits is mixing and therefore
ergodic. The map possesses two distinct LRN and SRN integrable limits. We
checked that the thermalization properties of observables in unitary circuits
are in line with previous observations for Hamiltonian systems [18, 19].

Figure 3: Renormalized Lyapunov spectrum $\\{\bar{\Lambda}_{i}\\}$ for SRN in
proximity to the integrable limit. The nonlinearity strength is fixed at
$g=1.0$. Angle $\theta$ varies from $10^{-1}$ (blue, top) to $5\cdot 10^{-4}$
(red, bottom). The inset shows the coefficient $1/\beta$ of the exponential
decay (see also Fig. 4) of the curves as a function of $\theta$. System size
$N=200$.

SRN integrable limit – We consider the limiting case $\theta=0$ with a fixed
nonzero value of nonlinearity strength $g$. In this setup the matrices
$\hat{C}_{n,n+1}$ become unity, thus decoupling the sites. The unitary
evolution applies a nonlinear norm dependent phase shift at each site:

$\displaystyle\hat{U}_{\text{int}}^{\text{SRN}}\vec{\Psi}=\sum_{n\in\mathds{Z}}e^{ig|\psi_{n}|^{2}}\psi_{n}.$
(4)

The system turns integrable with the norm at each site $|\psi_{n}|$ being a
constant of motion. By introducing a weak deviation from the limit
$0<\theta\ll 1$ one induces a network with next-to-nearest-neighbor hopping -
an SRN (see Supplemental Material for more details
[[SeetheSupplementaryMaterialatthe[URL]]SM]). We schematically represent the
parameter space and corresponding network in Fig. 1(b).

LRN integrable limit – Vanishing nonlinearity strength $g=0$ results in a
linear evolution with corresponding eigenvalue problem
$e^{i\omega}\vec{\Psi}(t)=\hat{U}_{\text{int}}^{\text{LRN}}\vec{\Psi}(t)$ (see
details in Supplemental Material [32]). The evolution corresponding to this
integrable limit is given by:

$\displaystyle\hat{U}_{\text{int}}^{\text{LRN}}\vec{\Psi}=\sum_{k}\hat{u}_{k}\vec{\psi}_{k},$
(5)

where $\vec{\psi}_{k}$ are the normal modes of the system and $u_{k}$ is the
evolution map in reciprocal space corresponding to the wavenumber $k$. In this
limit the absolute values of the normal mode amplitudes
$|c_{k}|=|\vec{\psi}_{k}^{\dagger}\cdot\vec{\Psi}(t)|$ are the constants of
motion. The deviation from this limit $g\neq 0$ induces an all-to-all coupling
among the normal modes of the system which respects translational invariance
through selection rules [32]. This by definition constitutes an LRN. The green
region in the control parameter space in Fig. 1(b) corresponds to that LRN
with the schematic representation of the network sketched right to it.

Figure 4: Renormalized Lyapunov spectrum $\\{\bar{\Lambda}_{i}\\}$ for SRN in
log scale (corresponding to the parameters and data of Fig. 3). The dashed
line is to guide the eye for the fit of exponential decay.

We compute the Lyapunov spectrum of Unitary Circuits in proximity to
integrable limits in order to resolve the entire set of characteristic time
scales. We follow the evolution of a set of orthogonal tangent vectors
$\\{\vec{w}_{i}\\}$ in the $2N$ dimensional phase space and compute the
increment $\gamma_{i}(t)=|\vec{w}_{i}(t)|$. Details on the approach can be
found in Section IV of the supplemental material [32]. For each vector we
compute the transient value $X_{i}(t)=1/t\sum_{\tau}^{t}\log\gamma(\tau)$
which in the infinite time limit turns into the Lyapunov characteristic
exponent (LCE) $\Lambda_{i}=\lim_{t\rightarrow\infty}X_{i}(t)$. The
numerically computed LCEs are the values of $X(t)$ extracted at the last step
of the dynamics. The LCEs are ordered from largest to smallest value upon
incrementing the index $i$. Due to the symplectic nature of the map the
spectrum is symmetric with LCEs coming in pairs
$\Lambda_{i}=-\Lambda_{2N-i+1}$. Norm conservation ensures two vanishing LCEs
$\Lambda_{N}=\Lambda_{N+1}=0$ [33]. According to the numerical setup, however,
it is impossible to achieve exact $\Lambda_{i}=0$ values, with bounds on the
smallest computed Lyapunov exponents $\Lambda_{\text{min}}\sim 1/t$.

The evolution of the phase space vector $\vec{\Psi}$ is obtained from
subsequent applications of the map $\hat{U}$. We use periodic boundary
conditions $\psi_{N+1}=\psi_{1}$. The initial conditions for the amplitudes of
the local complex components are drawn from an exponential distribution
$p(x)=e^{-x}$, while their phases were generated as uncorrelated and random
numbers chosen uniformly from the interval $[0,2\pi]$. The state vector is
then uniformly rescaled such that the norm density
$\frac{1}{N}\sum|\psi_{n}|^{2}=1$. The largest integration time varied between
$t_{\text{max}}=10^{8}$ and $t_{\text{max}}=10^{9}$. We have performed
computations for a set of initial conditions to ensure the independence of
results on the choice of initial state. Unless stated otherwise the system
size is set to $N=200$.

Figure 5: Renormalized Lyapunov spectrum $\\{\bar{\Lambda}_{i}\\}$ for LRN in
proximity to the integrable limit. The angle $\theta=0.33\pi$ is fixed. The
deviation from integrable limit $g$ varies from $10^{-1}$ (blue) to $10^{-3}$
(purple). The dashed line is to show the asymptotic curve as $g\rightarrow 0$.
In the inset we showcase the asymptotic curve as parameter $\theta$ is varied.
For all cases system size $N=200$.

First, we show the dependence of the largest Lyapunov exponent $\Lambda_{max}$
on the distance $g$ or $\theta$ to the integrable limit in Fig. 2 for both
networks. Both curves show a dependence which might resemble a power law
$\Lambda_{max}\sim g^{\nu}$ and $\Lambda_{max}\sim\theta^{\mu}$ with
$\nu\approx 1/2$ and $\mu\approx 3/2$. Remarkably the SRN case shows a much
slower diminishing of $\Lambda_{max}$ upon approaching the integrable limit as
compared to the LRN. This is similar to the study of a Hamiltonian system
dynamics [18, 34, 35]. Our data in Fig.2 are obtained for two different system
sizes $N=100,200$ and show very good agreement, therefore we can exclude
finite-size corrections. We now proceed to the analysis of the entire Lyapunov
spectrum. In Fig. 3 and Fig. 5 we show the renormalized Lyapunov spectrum
$\bar{\Lambda}_{i}=\Lambda_{i}/\Lambda_{\text{max}}$ for SRN and LRN
respectively. The index of Lyapunov exponents is rescaled $\rho=i/N$ so that
all positive LCEs $\bar{\Lambda}(\rho)$ correspond to $\rho\in[0,1]$. We
notice a dramatic qualitative difference between the two regimes. For the LRN
case the renormalized Lyapunov spectrum $\bar{\Lambda}(\rho)$ converges to a
limiting smooth curve for $g\rightarrow 0$. For the SRN instead that curve
vanishes in an exponential way. We will explain these observations in detail
below.

For the SRN an increasing number of Lyapunov exponents seems to be vanishing
upon approaching the integrable limit as seen in Fig. 3. We replot the same
spectrum in log scale in Fig. 4 and notice an exponential decay of the
renormalized spectrum:

$\displaystyle\Lambda_{\rho}^{\text{SRN}}=\Lambda_{\text{max}}e^{-\rho/\beta}.$
(6)

We fit the exponential decay and plot the exponent $1/\beta$ versus $\theta$
in the inset in Fig. 3. We observe that the exponent is rapidly diverging upon
approaching the integrable limit such that $\beta(\theta\rightarrow
0)\rightarrow 0$. The entire Lyapunov spectrum of the SRN is therefore
characterized by two scaling parameters - the largest Lyapunov exponent
$\Lambda_{max}$ which is an inverse time scale, and the parameter $\beta$
which is an inverse length scale. This result explains and agrees with
previous studies on dynamical glass in Hamiltonian systems [18, 19] where the
largest Lyapunov exponent stems from local resonances with rapidly increasing
distance between them upon approaching the integrable limit. Our results show
that the Lyapunov spectrum contains the quantitative scaling parameters of
that dynamical glass theory.

In contrast, the LRN spectrum is characterized by single parameter scaling.
The renormalized Lyapunov spectrum approaches a smooth limiting curve
$\bar{\Lambda}(\rho)$ as seen in Fig. 5. We compute the limiting curves by a
linear fit of $\bar{\Lambda}_{\rho}(g)$ at each value of $\rho_{j}$. Thus in
the LRN regime, the final form of the spectrum is given by:

$\displaystyle\Lambda_{\rho}^{\text{LRN}}=\Lambda_{\text{max}}f(\rho,\theta)$
(7)

The limiting curves for different values of $\theta$ are plotted in the inset
of Fig.5 and show little if any variation. It appears that the limiting curve
$f(\rho)$ is universal for all LRN parameter choices.

To further characterize the chaotic dynamics and showcase the difference
between SR and LR networks we compute the Kolmogorov-Sinai entropy
$K_{\text{KS}}=\int_{0}^{1}\bar{\Lambda}_{\rho}d\rho$. In the SRN case, from
eq. (6) follows

$\displaystyle
K_{\text{KS}}^{\text{SRN}}=\Lambda_{\text{max}}\beta(1-e^{-1/\beta})\;.$ (8)

Therefore the renormalized Kolmogorov-Sinai entropy
$k_{\text{KS}}=K_{\text{KS}}/\Lambda_{max}$ will tend to zero in the
integrable limit $k_{\text{KS}}\approx\beta$.

In the LRN regime the integral over the asymptotic function $f(\rho,\theta)$
(see eq. (7)) will lead to finite values of the renormalized KS entropy
$k_{\text{KS}}=\int_{0}^{1}f(\rho)d\rho>0$ at the very integrable limit.

To conclude, we identified the Lyapunov spectrum as a universal characteristic
descriptor of the complex phase space dynamics of a macroscopic system in
proximity to an integrable limit. The limit is characterized by a macroscopic
number of conserved actions. We identify two classes of nonintegrable
perturbation networks - short and long-range ones. Long-range networks are
characterized by a single parameter scaling of the Lyapunov spectrum - knowing
the largest Lyapunov exponent allows to reconstruct the entire spectrum.
Consequently all Lyapunov exponents scale as the largest one upon approaching
the integrable limit. Typical long-range networks are realized with
translationally invariant lattice systems in the limit of weak nonlinearity.
In that case, the actions correspond to normal modes extended over the entire
real space. Nonintegrable perturbations will typically couple them all. On the
other side, short-range networks are characterized by a two-parameter scaling.
In addition to the largest Lyapunov exponent, a diverging length scale results
in a suppression of the renormalized Lyapunov spectrum upon approaching the
integrable limit. Typical short-range networks are realized with lattice
systems and local (short-range) nonlinearities in the limit of weak coupling.
Interestingly the short-range network case appears to include disordered
systems as well. We, therefore, expect that disordered systems with weak
short-range nonlinearities will correspond to the SRN universality class.
Quantizing the classical dynamics could lead to many-body localization in the
case of short-range networks, as opposed to long-range networks.

Acknowledgments: This work was supported by the Institute for Basic Science
(Project number: IBS-R024-D1).

## References

  * Huang [1987] K. Huang,  _Statistical Mechanics_ , 2nd ed. (John Wiley & Sons, 1987).
  * Gogolin  _et al._ [2011] C. Gogolin, M. P. Müller, and J. Eisert, Physical review letters 106, 040401 (2011).
  * Rigol [2009] M. Rigol, Physical review letters 103, 100403 (2009).
  * Campbell  _et al._ [2005] D. K. Campbell, P. Rosenau, and G. M. Zaslavsky, Chaos: An Interdisciplinary Journal of Nonlinear Science 15, 015101 (2005).
  * Gaveau and Schulman [2015] B. Gaveau and L. S. Schulman, The European Physical Journal Special Topics 224, 891 (2015).
  * Bouchaud [1992] J.-P. Bouchaud, Journal de Physique I 2, 1705 (1992).
  * Bel and Barkai [2006] G. Bel and E. Barkai, EPL (Europhysics Letters) 74, 15 (2006).
  * Bel and Barkai [2005] G. Bel and E. Barkai, Physical review letters 94, 240602 (2005).
  * Ford [1992] J. Ford, Physics Reports 213, 271 (1992).
  * Zabusky and Deem [1967] N. J. Zabusky and G. S. Deem, Journal of Computational Physics 2, 126 (1967).
  * Toda [1967] M. Toda, Journal of the Physical Society of Japan 22, 431 (1967), https://doi.org/10.1143/JPSJ.22.431 .
  * Hénon and Heiles [1964] M. Hénon and C. Heiles, The Astronomical Journal 69, 73 (1964).
  * Kolmogorov [1954] A. N. Kolmogorov, in  _Dokl. Akad. Nauk SSSR_ , Vol. 98 (1954) pp. 527–530.
  * Arnold [2009] V. I. Arnold, Collected Works: Representations of Functions, Celestial Mechanics and KAM Theory, 1957–1965 , 267 (2009).
  * Moser [1962] J. Moser, Nachr. Akad. Wiss. Göttingen, II , 1 (1962).
  * Wayne [1984] C. E. Wayne, Communications in mathematical physics 96, 311 (1984).
  * Fermi  _et al._ [1955] E. Fermi, P. Pasta, S. Ulam, and M. Tsingou,  _STUDIES OF THE NONLINEAR PROBLEMS_ , Tech. Rep. (Los Alamos Scientific Lab., N. Mex., 1955).
  * Danieli  _et al._ [2019] C. Danieli, T. Mithun, Y. Kati, D. K. Campbell, and S. Flach, Phys. Rev. E 100, 032217 (2019).
  * Mithun  _et al._ [2019] T. Mithun, C. Danieli, Y. Kati, and S. Flach, Physical review letters 122, 054102 (2019).
  * Mithun  _et al._ [2021] T. Mithun, C. Danieli, M. V. Fistul, B. L. Altshuler, and S. Flach, Phys. Rev. E 104, 014218 (2021).
  * Goldfriend and Kurchan [2019] T. Goldfriend and J. Kurchan, Physical Review E 99, 022146 (2019).
  * Ganapa  _et al._ [2020] S. Ganapa, A. Apte, and A. Dhar, Journal of Statistical Physics 180, 1010 (2020).
  * Baldovin  _et al._ [2021] M. Baldovin, A. Vulpiani, and G. Gradenigo, Journal of Statistical Physics 183, 1 (2021).
  * Oseledets [1968] V. I. Oseledets, Trudy Moskovskogo Matematicheskogo Obshchestva 19, 179 (1968).
  * de Wijn  _et al._ [2015] A. S. de Wijn, B. Hess, and B. V. Fine, Physical Review E 92, 062929 (2015).
  * Iubini and Politi [2021] S. Iubini and A. Politi, Chaos, Solitons & Fractals 147, 110954 (2021).
  * Benettin  _et al._ [1980] G. Benettin, L. Galgani, A. Giorgilli, and J.-M. Strelcyn, Meccanica 15, 9 (1980).
  * Vakulchyk  _et al._ [2019] I. Vakulchyk, M. V. Fistul, and S. Flach, Physical review letters 122, 040501 (2019).
  * Vakulchyk  _et al._ [2017] I. Vakulchyk, M. V. Fistul, P. Qin, and S. Flach, Phys. Rev. B 96, 144204 (2017).
  * Malishava  _et al._ [2020] M. Malishava, I. Vakulchyk, M. Fistul, and S. Flach, Physical Review B 101, 144201 (2020).
  * Vakulchyk  _et al._ [2018] I. Vakulchyk, M. Fistul, Y. Zolotaryuk, and S. Flach, Chaos: An Interdisciplinary Journal of Nonlinear Science 28, 123104 (2018).
  * [32] .
  * Skokos [2010] C. Skokos, The lyapunov characteristic exponents and their computation, in  _Dynamics of Small Solar System Bodies and Exoplanets_, edited by J. J. Souchay and R. Dvorak (Springer Berlin Heidelberg, Berlin, Heidelberg, 2010) pp. 63–135.
  * Mulansky  _et al._ [2011] M. Mulansky, K. Ahnert, A. Pikovsky, and D. Shepelyansky, Journal of Statistical Physics 145, 1256 (2011).
  * De Wijn  _et al._ [2013] A. De Wijn, B. Hess, and B. Fine, Journal of Physics A: Mathematical and Theoretical 46, 254012 (2013).