# Grassmann phase space dynamics of strongly-correlated fermions

Hassan Al-Hamzawi  Alessandro Principi  Leone Di Mauro Villari  Department of
Physics and Astronomy, University of Manchester, Manchester M13 9PL, UK

###### Abstract

We discuss the numerical implementation of two related representations of
fermionic density matrices which have been introduced in Annals of Physics
370, 12 (2016). In both of them, the density matrix is expanded in a basis of
Bargmann coherent states with weights given by the two phase space
distributions. We derive the equations of motion for the distributions when
imaginary time evolution is generated by the Hubbard Hamiltonian. One of them
is a Grassmann Fokker-Planck equation that can be re-cast into a remarkably
simple Itô form involving solely complex variables. In spite of this simple
form, we demonstrate that complications arise in numerically computing the
expectation value of any observable. These are due to exponential growth in
the matrix elements of the stochastic propagator, delicate numerical
sensitivity in performing primitive linear algebra operations, and the re-
appearance of a sign problem.

##  I  Introduction

The study of strongly-correlated systems is a notoriously difficult
mathematical problem since interactions amongst all particles, and hence
nonlinearities of the associated field theory, become fundamental in
determining the state of the system. In this respect, practical applications
favour descriptions that are amenable to numerical simulation. It is therefore
pivotal to find new numerical methods that can overcome the limitations of
existing techniques, such as the sign problem of quantum Monte Carlo
calculations or the substantial overhead in iterant matrix product algorithms.
One way to achieve this is to represent the density matrix by expanding it in
terms of an overcomplete basis, for example given by coherent states. In this
way, states and operators take the form of functions and operations over them,
respectively. Such operations preserve the operator algebra.

The algebra of bosonic operators is based on commutators, therefore states and
operators can be expressed as functions of (commuting) complex variables. This
approach has been extensively studied and has led to numerous formulations,
such as the Glauber and Sudarshan P distribution [1, 2], the Husimi Q
distribution [3], and the Wigner W distribution [4, 5]. These representations
have been widely used to study strongly-interacting and heavily nonlinear
bosonic systems in various fields, from quantum optics [6] to Bose-Einstein
condensation [7] and gravity analogues [8].

The algebra of fermionic operators on the other hand involves anti-commutators
and cannot be represented in terms of complex variables. This complicates the
computational simulation of fermionic systems since it is difficult to
naturally absorb the operator algebra into the description. Calculations
performed using complex variables must explicitly preserve exchange anti-
symmetry. This produces an abundance of terms that have similar moduli and
opposing signs and thus leads to a great deal of cancellation in their
combined contributions (i.e. the infamous sign problem). A finer granularity
and increased computational complexity are hence required to offset the loss
in numerical precision.

Attempts have been made to adapt the techniques of coherent phase space
distributions to fermions, as proposed by Cahill and Glauber [9], and Dalton,
Jeffers, and Barnett [10, 11]. Analogous descriptions have been used to derive
analytical results regarding two-level atoms interacting with electromagnetic
fields [12, 13] as well as in the Jaynes-Cummings model [14]. Corney and
Drummond have also shown that a different but related representation given in
terms of Gaussian states has an equivalent formulation for fermions [15, 16,
17]. The latter has been able to reproduce certain ground state properties of
the Hubbard model [18], albeit with the aid of auxiliary symmetry projectors.

So far, Grassmann phase space formulations have had only very limited
numerical implementations [19] despite being the most natural and intuitive
technique to represent fermionic quantum states as phase space distributions.
These techniques are also very promising since, supposedly, they do not suffer
from some of the shortcomings of conventional methods, such as those mentioned
above. For example, it can be shown that the statistical weights introduced by
the Grassmann phase space distributions are related to the determinants of
what are generally positive-definite matrices. This in turn implies that the
entire formulation is in principle free of the infamous sign problem. Given
these hopeful results, it is fundamental to test how far can these novel
methods go, and to what extent they can enrich the present knowledge of the
phase diagrams of strongly-correlated systems.

This article aims to extend the application of the Grassmann coherent state
phase space developed in Dalton, Jeffers, and Barnett’s work to the simulation
of the Hubbard model, the “gold standard” of strongly-correlated physics. In
fact, the Hubbard model has been the centre of immense theoretical interest
since its inception, due to no shortage of intriguing phenomenology, in
particular its relation to high-temperature superconductivity. Whilst certain
limiting cases have been extensively studied, the full phase diagram remains a
contentious topic with many questions left unanswered [20]. Due to the
aforementioned intrigue, as well as its wide applicability in many-body
quantum field theory, it serves as an ideal test-case to assess the
suitability of any particular choice of representation.

In Sect. II, we provide a brief description of fermionic coherent states, and
a summary of several important results regarding the representation of
physical density matrices as Grassmann phase space distributions. In Sect.
III, a pair of differential equations are derived for the imaginary time
evolution of the grand-canonical ensemble. This is based on a choice of two
available Grassmann representations, one of which results in a differential
equation that has the Fokker-Planck structure. In Sect. IV, the Fokker-Planck
equation is expressed in a stochastic Itô form using only complex numbers, and
the extraction of quantum correlations from the statistical samples is
addressed. Finally, in Sect. V, numerical results are presented in the case of
a $2\times 2$ square lattice. We show that unfortunately, despite the
theoretical validity of the formulation, the stochastic dynamics results in
the exponential growth in the propagator matrix elements. This introduces
numerical instability into our implementation and reflects negatively on the
accuracy and precision of observables.

##  II  Phase Space Representation

In this section we remind the reader of the nature of the basis used to
parameterise the Hilbert space, i.e. the Bargmann coherent states, as well as
of its properties. We introduce fermion creation and annihilation operators
$\hat{a}_{i}^{\dagger}$ and $\hat{a}_{i}$, respectively, where $i=1,\dots,n$
represents a set of quantum numbers. They satisfy the canonical anti-
commutation relations. They also anti-commute with any Grassmann variable
$g_{i}$. In formulae,

$\\{\hat{a}_{i},\hat{a}_{j}^{\dagger}\\}=\delta_{i,j}~{},\\{\hat{a}_{i},\hat{a}_{j}\\}=0,$
(1)

and

$\\{g_{i},\hat{a}_{j}\\}=\\{g_{i},\hat{a}_{j}^{\dagger}\\}=0,$ (2)

where $\\{\cdot,\cdot\\}$ is the anticommutator between pairs of operators or
(Grassmann) variables. The Grassmann variables commute with the vacuum state
$g_{i}\left|0\right\rangle=\left|0\right\rangle g_{i}$. Thus, in the
occupation basis Grassmann variables commute or anti-commute with even or odd
occupancy states respectively. Hermitian conjugation is extended to also act
on any Grassmann coefficients. It reverses the order of the symbols and acts
to conjugate them

$(g_{i}\hat{a}_{j})^{\dagger}=\hat{a}_{j}^{\dagger}g_{i}^{*}.$ (3)

We emphasise that the conjugate of a Grassmann variable is itself an
independent Grassmann variable, see App. A for more details.

###  II.1  Bargmann Coherent States

Bargmann coherent states [10, 11] are defined by the action of the
displacement operator
$\hat{R}(\bm{g})=\exp{(\sum_{i}\hat{a}_{i}^{\dagger}g_{i})}$ on the vacuum
state

$\left|\bm{g}\right\rangle=\hat{R}(\bm{g})\left|0\right\rangle.$ (4)

The bra is then defined as

$\left\langle\bm{g}\right|=\left\langle 0\right|\hat{R}(\bm{g})^{\dagger},$
(5)

where $\hat{R}(\bm{g})^{\dagger}=\exp{(\sum_{i}g_{i}^{*}\hat{a}_{i})}$ and, as
pointed out in the end of the previous section, $g_{i}$ and $g_{i}^{*}$ are
independent variables. By construction they commute with Grassmann variables.
The action of creation and annihilation operators on these states is

$\displaystyle\hat{a}_{i}\left|\bm{g}\right\rangle=g_{i}\left|\bm{g}\right\rangle,~{}\hat{a}_{i}^{\dagger}\left|\bm{g}\right\rangle=-\overrightarrow{\frac{\partial}{\partial
g_{i}}}\left|\bm{g}\right\rangle=\left|\bm{g}\right\rangle\overleftarrow{\frac{\partial}{\partial
g_{i}}},$ (6)
$\displaystyle\left\langle\bm{g}\right|\hat{a}_{i}^{\dagger}=\left\langle\bm{g}\right|g_{i}^{*},~{}\left\langle\bm{g}\right|\hat{a}_{i}=-\left\langle\bm{g}\right|\overleftarrow{\frac{\partial}{\partial
g_{i}^{*}}}=\overrightarrow{\frac{\partial}{\partial
g_{i}^{*}}}\left\langle\bm{g}\right|.$ (7)

Here the arrows make clear the direction in which the derivative operates (see
also App. B.1 for more detail about Grassmann differentiation). Furthermore,
as the conventional coherent states, they satisfy these relations

$\displaystyle\left\langle\bm{u}|\bm{v}\right\rangle=\exp{(\bm{u}^{*}\cdot\bm{v})},$
(8)
$\displaystyle\mathrm{Tr}(\left|\bm{v}\right\rangle\left\langle\bm{u}\right|)=\left\langle-\bm{u}|\bm{v}\right\rangle=\left\langle\bm{u}|-\bm{v}\right\rangle,$
(9) $\displaystyle\hat{I}=\int dg_{n}^{*}dg_{n}\dots dg_{1}^{*}dg_{1}\
\exp{(\bm{g}\cdot\bm{g}^{*})}\left|\bm{g}\right\rangle\left\langle\bm{g}\right|,$
(10)

which in turn imply that the Bargmann states form an overcomplete basis for
the Hilbert space. Note that, since $g_{i}$ and $g_{i}^{*}$ are independent
variables for all $i=1,\dots,n$, the integration in Eq. (10) runs over both
sets of Grassmann variables.

###  II.2  Grassmann Distribution Functions

Dalton, Jeffers, and Barnett [10, 11] demonstrate that any physical fermionic
density matrix can be uniquely expanded as a distribution over Bargmann
coherent states. Using Eq. (10), the density matrix becomes

$\hat{\rho}=\int d^{n}\bm{u}\
d^{n}\bm{v}^{*}B(\bm{u},\bm{v}^{*})\left|\bm{u}\right\rangle\left\langle\bm{v}\right|,$
(11)

where $\bm{u}$ and $\bm{v}$ are two unrelated sets of Grassmann variables (see
also App. B.2 for more detail about Grassmann integration) and

$B(\bm{u},\bm{v}^{*})=\int
d^{n}\bm{u}^{*}d^{n}\bm{v}\left\langle\bm{u}|\hat{\rho}|\bm{v}\right\rangle
e^{\bm{u}\cdot\bm{u}^{*}}e^{\bm{v}\cdot\bm{v}^{*}}.$ (12)

Defining $\bm{g}=\bm{u}$ and $\bm{\tilde{g}}=\bm{v}^{*}$, as well as the pair
of normalised and un-normalised projector operators

$\hat{\Lambda}^{(n)}(\bm{g},\bm{\tilde{g}})=\frac{\left|\bm{g}\right\rangle\left\langle\bm{\tilde{g}}^{*}\right|}{\mathrm{Tr}(\left|\bm{g}\right\rangle\left\langle\bm{\tilde{g}}^{*}\right|)},\
\hat{\Lambda}(\bm{g},\bm{\tilde{g}})=\left|\bm{g}\right\rangle\left\langle\bm{\tilde{g}}^{*}\right|,$
(13)

it is possible to re-express the density matrix in two equivalent forms

$\displaystyle\hat{\rho}$ $\displaystyle=\int d^{n}\bm{g}\
d^{n}\bm{\tilde{g}}\
B(\bm{g},\bm{\tilde{g}})\hat{\Lambda}(\bm{g},\bm{\tilde{g}})$ (14)
$\displaystyle=\int d^{n}\bm{g}\ d^{n}\bm{\tilde{g}}\
P(\bm{g},\bm{\tilde{g}})\hat{\Lambda}^{(n)}(\bm{g},\bm{\tilde{g}}).$ (15)

Both distribution functions are Grassmann even, with the $B$-distribution
expressed as the Grassmann integral in Eq. (12) and related to
$P$-distribution through

$P(\bm{g},\bm{\tilde{g}})=B(\bm{g},\bm{\tilde{g}})e^{\bm{g}\cdot\bm{\tilde{g}}}.$
(16)

The expectation value of any normally ordered product of creation and
annihilation operators can be evaluated through a Grassmann integral weighted
by the $P$-distribution function as

$\displaystyle\begin{aligned}
\mathrm{Tr}[\hat{a}_{\mu_{1}}^{\dagger}\dots\hat{a}_{\mu_{p}}^{\dagger}\hat{a}_{\nu_{q}}\dots\hat{a}_{\nu_{1}}\hat{\rho}]&=\int
d^{n}\bm{g}\ d^{n}\bm{\tilde{g}}\ P(\bm{g},\bm{\tilde{g}})\\\ &\times
g_{\nu_{q}}\dots
g_{\nu_{1}}\tilde{g}_{\mu_{1}}\dots\tilde{g}_{\mu_{p}}.\end{aligned}$ (17)

Finally, via the integral expansion of the density matrix given in Eqs. (14)
and (15), the relations stated in Eqs. (6) and (7), and the integration by
parts formulae found in Eqs. (78) and (79), it is possible to show the
following operator to phase space mappings

$\displaystyle\begin{gathered}\hat{a}_{i}\hat{\rho}\to g_{i}P,\\\
\hat{\rho}\hat{a}_{i}\to
P\left(\overleftarrow{\frac{\partial}{\partial\tilde{g}_{i}}}-g_{i}\right),\\\
\hat{a}_{i}^{\dagger}\hat{\rho}\to\left(\overrightarrow{\frac{\partial}{\partial
g_{i}}}-\tilde{g}_{i}\right)P,\\\ \hat{\rho}\hat{a}_{i}^{\dagger}\to
P\tilde{g}_{i},\end{gathered}$ (22)

for the $P$-distribution and

$\displaystyle\begin{gathered}\hat{a}_{i}\hat{\rho}\to g_{i}B,\\\
\hat{\rho}\hat{a}_{i}\to
B\overleftarrow{\frac{\partial}{\partial\tilde{g}_{i}}},\\\
\hat{a}_{i}^{\dagger}\hat{\rho}\to\overrightarrow{\frac{\partial}{\partial
g_{i}}}B,\\\ \hat{\rho}\hat{a}_{i}^{\dagger}\to B\tilde{g}_{i},\end{gathered}$
(27)

for the $B$-distribution. In the case of sequences of several operators, the
previous correspondences are applied consecutively in order of proximity to
the density matrix.

##  III  Imaginary Time Evolution

The Hubbard model [21, 22, 23] is one of the most established and well studied
formulations of interacting fermions on a lattice. It is described by the
Hamiltonian

$\hat{H}=-\sum\limits_{i,j,\sigma}\bm{\tau}_{i,j}\hat{a}_{i,\sigma}^{\dagger}\hat{a}_{j,\sigma}+U\sum\limits_{k}\hat{a}_{k,\uparrow}^{\dagger}\hat{a}_{k,\uparrow}\hat{a}_{k,\downarrow}^{\dagger}\hat{a}_{k,\downarrow},$
(28)

where $\bm{\tau}$ is the hopping matrix. Its elements represent hopping
amongst nearest neighbours, next-to-nearest neighbours, and so on. In this
paper we will consider only the former. We parametrise the hopping between any
pair of neighbouring sites $(i,j)$ with the amplitude $t$. In Eq. (28),
$\sigma$ denotes the particle spin and takes on the symbolic values
$\uparrow\downarrow$, while $U$ is the on-site interaction potential. For
later factorisation convenience, it is possible to re-cast the second term of
the Hamiltonian in the form of a normally ordered product as

$\hat{H}_{I}=\sum\limits_{i,j,\sigma,\sigma^{\prime}}\bm{U}_{i\sigma,j\sigma^{\prime}}\hat{a}_{i,\sigma}^{\dagger}\hat{a}_{j,\sigma^{\prime}}^{\dagger}\hat{a}_{j,\sigma^{\prime}}\hat{a}_{i,\sigma},$
(29)

where

$\bm{U}_{i\sigma,j\sigma^{\prime}}=-\frac{|U|}{2}\sum\limits_{k}\delta_{i,k}\kappa(\sigma)\delta_{j,k}\kappa(\sigma^{\prime}),$
(30)

and

$\kappa(\sigma)=\begin{cases}1&\sigma=\uparrow\\\
-\mathrm{sign}(U)&\sigma=\downarrow\end{cases}.$ (31)

To derive the equations of motion satisfied by the $P$\- and
$B$-distributions, we note that the density matrix for the grand-canonical
ensemble satisfies the following Matsubara differential equation

$\frac{d}{d\beta}\hat{\rho}=-\frac{1}{2}\\{\hat{H}-\mu\hat{N},\hat{\rho}\\},$
(32)

where $\beta^{-1}=k_{B}T$. Henceforth, we absorb the chemical potential $\mu$
into the diagonal elements of the hopping matrix $\bm{\tau}$ to simplify our
notation.

Appendix C shows how the correspondences of Eqs. (22) can be utilised to
derive the differential equation describing the dynamics of the
$P$-distribution. The key result is given in Eq. (82), which represents a
family of coupled differential equations relating the complex valued
coefficients of the multinomial associated with $P$. The presence of terms
with a mismatched number of Grassmann multiplications and differentiations
hinders the ability to re-express the differential equation in any form
amenable to numerical integration.

The $B$-distribution on the other hand has simpler correspondences, as shown
in Eqs. (27). Using them, it is possible to show that it satisfies a relation
that has the general structure of a Fokker-Planck equation, though the
coefficients are Grassmann variables or products thereof [see Eq. (83) derived
in App. C]. With summation implied over all the present indices, the
differential equation is re-written in a compact form as

$\displaystyle\begin{gathered}\frac{d}{d\beta}B=-\left(BA_{i\sigma}^{\mu}\right)\overleftarrow{\frac{\partial}{\partial
g_{i\sigma}^{\mu}}}+\frac{1}{2}\left(BD_{i\sigma,j\sigma^{\prime}}^{\mu,\nu}\right)\overleftarrow{\frac{\partial}{\partial
g_{j\sigma^{\prime}}^{\nu}}}\overleftarrow{\frac{\partial}{\partial
g_{i\sigma}^{\mu}}},\\\
A_{i\sigma}^{\mu}=-\frac{1}{2}\sum\limits_{j,\sigma^{\prime},\nu}\delta_{\sigma,\sigma^{\prime}}^{\mu,\nu}\left(\delta^{\nu,1}\bm{\tau}_{i,j}+\delta^{\nu,2}\bm{\tau}_{i,j}^{*}\right)g_{j\sigma^{\prime}}^{\nu},\\\
D_{i\sigma,j\sigma^{\prime}}^{\mu,\nu}=\frac{|U|}{2}\delta^{\mu,\nu}\sum\limits_{k}\delta_{i,k}g_{i\sigma}^{\mu}\kappa(\sigma)\delta_{j,k}g_{j\sigma^{\prime}}^{\nu}\kappa(\sigma^{\prime}),\end{gathered}$
(36)

where the notation $g^{\mu}$ has been introduced to distinguish the two sets
of Grassmann variables, i.e. $g^{1}=\bm{g}$ and $g^{2}=\bm{\tilde{g}}$.
Despite removing the mismatch in the number of Grassmann multiplications and
differentiations, Eq. (36) viewed directly still represents a family of
coupled differential equations. However, this current form presents the
opportunity for a stochastic transformation.

Considering any arbitrary Grassmann functions $F(\bm{g},\bm{\tilde{g}})$ and
$H(\bm{g},\bm{\tilde{g}})=F(\bm{g},\bm{\tilde{g}})e^{\bm{g}\cdot\bm{\tilde{g}}}$
then their phase space distribution expectations are related through

$\displaystyle\begin{aligned} \overline{F(\bm{g},\bm{\tilde{g}})}|_{P}&=\int
d^{n}\bm{g}\ d^{n}\bm{\tilde{g}}\
F(\bm{g},\bm{\tilde{g}})P(\bm{g},\bm{\tilde{g}},\beta)\\\ &=\int d^{n}\bm{g}\
d^{n}\bm{\tilde{g}}\
F(\bm{g},\bm{\tilde{g}})B(\bm{g},\bm{\tilde{g}},\beta)e^{\bm{g}\cdot\bm{\tilde{g}}}\\\
&=\overline{H(\bm{g},\bm{\tilde{g}})}|_{B}.\end{aligned}$ (37)

Taking the derivative of Eq. (37) with respect to $\beta$ (which only acts on
the distribution function), and using the Fokker-Planck relation in Eq. (36)
we obtain

$\displaystyle\begin{aligned}
\frac{d}{d\beta}\overline{H(\bm{g},\bm{\tilde{g}})}|_{B}&=\int d^{n}\bm{g}\
d^{n}\bm{\tilde{g}}\
H(\bm{g},\bm{\tilde{g}})\frac{d}{d\beta}B(\bm{g},\bm{\tilde{g}},\beta)\\\
&=\int d^{n}\bm{g}\ d^{n}\bm{\tilde{g}}\
[-\sum\limits_{i,\sigma,\mu}H\overleftarrow{\frac{\partial}{\partial
g_{i\sigma}^{\mu}}}A_{i\sigma}^{\mu}\\\
&+\frac{1}{2}\sum\limits_{i,j,\sigma,\sigma^{\prime},\mu,\nu}H\overleftarrow{\frac{\partial}{\partial
g_{j\sigma^{\prime}}^{\nu}}}\overleftarrow{\frac{\partial}{\partial
g_{i\sigma}^{\mu}}}D_{i\sigma,j\sigma^{\prime}}^{\mu,\nu}]B.\end{aligned}$
(38)

Here, we integrated by parts using the standard rules of Grassmann calculus,
which are provided in Eqs. (68) and (79) for convenience. Eq. (38) is a useful
relation to derive stochastic evolution in the following section.

##  IV  Stochastic Sampling

Instead of considering the expectation value of a fixed Grassmann function
$H(\bm{g},\bm{\tilde{g}})$ and an evolving distribution
$B(\bm{g},\bm{\tilde{g}},\beta)$, stochastic calculus is employed to
facilitate sampling. The Grassmann variables themselves are allowed to evolve
under a static distribution $B(\bm{g},\bm{\tilde{g}})$ which does not depend
explicitly on $\beta$ (though it can depend on it through $\bm{g}$ and
$\bm{\tilde{g}}$). The evolution of such variables, which we denote as
$g_{i\sigma}^{\mu}(\beta)$, is such that it ensures that Eq. (38) continues to
hold in the mean.

Expanding the change of Eq. (37) to second order using the rules of Grassmann
calculus recalled in Eq. (73) for convenience, we obtain

$\displaystyle\begin{aligned} \overline{\delta H}&=\int d^{n}\bm{g}\
d^{n}\bm{\tilde{g}}\
(H(\beta+\delta\beta)-H(\beta))B(\bm{g},\bm{\tilde{g}})\\\ &=\int d^{n}\bm{g}\
d^{n}\bm{\tilde{g}}\
[\sum\limits_{i,\sigma,\mu}H\overleftarrow{\frac{\partial}{\partial
g_{i\sigma}^{\mu}}}\delta g_{i\sigma}^{\mu}\\\
&+\frac{1}{2}\sum\limits_{i,j,\sigma,\sigma^{\prime},\mu,\nu}H\overleftarrow{\frac{\partial}{\partial
g_{j\sigma^{\prime}}^{\nu}}}\overleftarrow{\frac{\partial}{\partial
g_{i\sigma}^{\mu}}}\delta g_{i\sigma}^{\mu}\delta
g_{j\sigma^{\prime}}^{\nu}]B(\bm{g},\bm{\tilde{g}}).\end{aligned}$ (39)

We assume a Langevin form for the stochastic differential equation satisfied
by the Grassmann variables, i.e.

$\delta
g_{i\sigma}^{\mu}=\mathcal{A}_{i\sigma}^{\mu}\delta\beta+\sum\limits_{k}\mathcal{B}_{i\sigma}^{\mu,k}\delta
W_{k}+\mathcal{O}(\delta\beta^{2}),$ (40)

where $\delta W_{k}$ are independent Wiener increments [24], i.e. they are
normally distributed random variables with zero mean and variance
$\delta\beta$. Thus the averages of one and two Grassmann variable variations
satisfy

$\displaystyle\begin{gathered}\left\langle\delta
g_{i\sigma}^{\mu}\right\rangle_{\text{Stochastic}}=\mathcal{A}_{i\sigma}^{\mu}\delta\beta+\mathcal{O}(\delta\beta^{2}),\\\
\left\langle\delta g_{i\sigma}^{\mu}\delta
g_{j\sigma^{\prime}}^{\nu}\right\rangle_{\text{Stochastic}}=\sum\limits_{k}\mathcal{B}_{i\sigma}^{\mu,k}\mathcal{B}_{j\sigma^{\prime}}^{\nu,k}\delta\beta+\mathcal{O}(\delta\beta^{2}),\end{gathered}$
(43)

where $\left\langle\dots\right\rangle_{\text{Stochastic}}$ is the average over
noise realisations. Then, assuming that

$A_{i\sigma}^{\mu}=-\mathcal{A}_{i\sigma}^{\mu},\
D_{i\sigma,j\sigma^{\prime}}^{\mu,\nu}=\sum\limits_{k}\mathcal{B}_{i\sigma}^{\mu,k}\mathcal{B}_{j\sigma^{\prime}}^{\nu,k},$
(44)

it ensures that Eq. (39) reduces to Eq. (38) in the mean and to leading order
in $\delta\beta$.

Combining Eq. (36) with Eqs. (40) and (44), the infinitesimal evolution is
finally expressed in terms of the following linear operators

$\displaystyle\begin{gathered}\bm{g}_{\uparrow}(\beta+\delta\beta)=\left(I+\frac{\delta\beta}{2}\bm{\tau}+\sqrt{\frac{|U|}{2}}R\right)\bm{g}_{\uparrow}(\beta),\\\
\bm{g}_{\downarrow}(\beta+\delta\beta)=\left(I+\frac{\delta\beta}{2}\bm{\tau}-\mathrm{sign}(U)\sqrt{\frac{|U|}{2}}R\right)\bm{g}_{\downarrow}(\beta),\\\
\bm{\tilde{g}}_{\uparrow}(\beta+\delta\beta)=\left(I+\frac{\delta\beta}{2}\bm{\tau}^{*}+\sqrt{\frac{|U|}{2}}\tilde{R}\right)\bm{\tilde{g}}_{\uparrow}(\beta),\\\
\bm{\tilde{g}}_{\downarrow}(\beta+\delta\beta)=\left(I+\frac{\delta\beta}{2}\bm{\tau}^{*}-\mathrm{sign}(U)\sqrt{\frac{|U|}{2}}\tilde{R}\right)\bm{\tilde{g}}_{\downarrow}(\beta),\end{gathered}$
(49)

where $I$ is the identity matrix, $\bm{\tau}$ is the hopping matrix (which
also includes the chemical potential), and $R$ and $\tilde{R}$ are both
diagonal matrices whose entries are independent Wiener increments. Integrating
Eq. (49) over the finite interval between $\beta_{0}$ and $\beta$, we obtain

$\bm{g}^{\mu}(\beta)=U^{\mu}(\beta;\beta_{0})\bm{g}^{\mu}(\beta_{0}),$ (50)

where $U^{\mu}(\beta;\beta_{0})$ is the imaginary time- ($\beta$-)ordered
product of the operators above.

Having established a relation for the dynamics of the Grassmann variables, the
static phase space distribution remains to be determined. Choosing initial
conditions at $\beta_{0}=0$ such that $\hat{\rho}\propto\hat{I}$ and
evaluating the Grassmann integral in Eq. (11), the $B$-distribution acquires
the form

$B(\bm{g}_{0},\bm{\tilde{g}}_{0})\propto\exp{(\bm{g}_{0}\cdot\bm{\tilde{g}}_{0})}.$
(51)

where $\bm{g}_{0}=(\bm{g}_{0\uparrow}^{T},\bm{g}_{0\downarrow}^{T})^{T}$ and
$\bm{\tilde{g}}_{0}=(\bm{\tilde{g}}_{0\uparrow}^{T},\bm{\tilde{g}}_{0\downarrow}^{T})^{T}$.

Figure 1: Ensemble statistics of $10^{5}$ samples on a $2\times 2$ lattice
with $U=4t$, $\mu=0$. The stochastic process realisations were generated via
Euler-Maruyama integration at a step size of $\delta\beta t=10^{-4}$. (a)
Exponential growth of the modulus of the statistical weights. (b) Fluctuations
in the sign of the statistical weights. (c) Expectation value of the
Hamiltonian as extrapolated from the ensemble correlations.

We can now introduce the characteristic function to obtain the normally
ordered correlations. Utilising two further sets of Grassmann variables
$\bm{u}$ and $\bm{v}$, the characteristic function is defined to be

$\chi(\bm{u},\bm{v})=\mathrm{Tr}\left[\exp{\left(i\sum\limits_{n}v_{n}\hat{a}_{n}\right)}\hat{\rho}\exp{\left(i\sum\limits_{m}\hat{a}_{m}^{\dagger}u_{m}\right)}\right].$
(52)

Recalling Eqs. (16) and (17), Eq. (52) reduces to

$\displaystyle\begin{aligned} \chi(\bm{u},\bm{v})&=\int d^{n}\bm{g}\
d^{n}\bm{\tilde{g}}\
B(\bm{g},\bm{\tilde{g}})e^{(\bm{g}\cdot\bm{\tilde{g}}+i\bm{v}\cdot\bm{g}+i\bm{\tilde{g}}\cdot\bm{u})}\\\
&=\overline{\exp{(\bm{g}\cdot\bm{\tilde{g}}+i\bm{v}\cdot\bm{g}+i\bm{\tilde{g}}\cdot\bm{u})}}|_{B}.\end{aligned}$
(53)

Hence, correlations can be obtained as derivatives of the characteristic
function. The Grassmann variables in Eq. (53) are evaluated at imaginary time
$\beta$. Using Eq. (50), we can rewrite Eq. (53) in terms of Grassmann
variables at imaginary time $\beta_{0}$ as

$\displaystyle\begin{gathered}\chi(\bm{u},\bm{v})=\overline{\exp{(\bm{g}\cdot\bm{\tilde{g}}+i\bm{v}\cdot\bm{g}+i\bm{\tilde{g}}\cdot\bm{u})}}|_{\beta}\\\
=\int d^{n}\bm{g}_{0}\ d^{n}\bm{g}_{0}\
e^{({\bm{g}}_{0}^{T}U^{T}\tilde{U}\bm{\tilde{g}}_{0}+i\bm{v}^{T}U\bm{g}_{0}+i{\bm{\tilde{g}}}_{0}^{T}\tilde{U}^{T}\bm{u})}B(\bm{g}_{0},\bm{\tilde{g}}_{0}),\end{gathered}$
(56)

with $B(\bm{g}_{0},\bm{\tilde{g}}_{0})$ being the exponential in Eq. (51).
This has the form of a Gaussian-type integral over Grassmann variables and can
be shown to give [25]

$\chi(\bm{u},\bm{v})\propto\mathrm{det}(I+U^{T}\tilde{U})e^{(\bm{u}^{T}\tilde{U}(I+U^{T}\tilde{U})^{-1}U^{T}\bm{v})}.$
(57)

Integration of Eq. (50) produces a statistical ensemble of trajectories, each
associated with its own set of moments. The dynamics, chosen to preserve Eq.
(38), ensures that the ensemble average converges to the true physical
correlation.

Let the realisations of the stochastic process be indexed by $\alpha$.
Identifying $Z^{\alpha}=\mathrm{det}(I+U^{T}_{\alpha}\tilde{U}_{\alpha})$ and
$M^{\alpha}=\tilde{U}_{\alpha}(I+U^{T}_{\alpha}\tilde{U}_{\alpha})^{-1}U^{T}_{\alpha}$,
where $U_{\alpha}$ and $\tilde{U}_{\alpha}$ are the evolution operators for a
given realisation, the first order correlations can be expressed as

$\langle\hat{a}_{i,\sigma}^{\dagger}\hat{a}_{j,\sigma^{\prime}}\rangle=\frac{\mathrm{Tr}[\hat{a}_{i,\sigma}^{\dagger}\hat{a}_{j,\sigma^{\prime}}\hat{\rho}]}{\mathrm{Tr}[\hat{\rho}]}=\frac{\sum_{\alpha}Z^{\alpha}M^{\alpha}_{i\sigma,j\sigma^{\prime}}}{\sum_{\alpha}Z^{\alpha}},$
(58)

with similar relations for the higher order moments obtained by expanding Eq.
(52) and Eq. (57). Since Wiener increments at different imaginary times are
uncorrelated, the average $\langle
I+U^{T}_{\alpha}\tilde{U}_{\alpha}\rangle_{\text{Stochastic}}$ is a positive-
definite matrix. This suggests that, excluding outliers, the statistical
weights appearing above ought to generally remain positive throughout the
ensemble.

##  V  Analysis and Results

Figure 2: Ensemble statistics of $10^{5}$ samples on a $2\times 2$ lattice
with $U=4t$, $\mu=2t$. The stochastic process realisations were generated via
Euler-Maruyama integration at a step size of $\delta\beta t=10^{-4}$. To
enhance readability, the (blue) scaled bars for $\beta t\geq 2.5$ represents
half the base-10 logarithm of the calculated error.

With all the necessary frameworks established, we now apply the formalism to
the case of a $2\times 2$ square lattice. This system is of a sufficiently
small size that the stochastic sampling may be contrasted with the results of
an exact diagonalisation.

We first discuss the results obtained in the symmetric configuration with
$\mu=0$. Before addressing the behaviour of physical quantities, we highlight
that Eq. (49) implies that the propagator matrix elements grow asymptotically
in an exponential manner. This is illustrated in Fig. 1(a), where we plot the
median statistical weight of the trajectories $Z^{\alpha}$ as a function of
$\beta t$ ($t$ being the hopping amplitude). This highlights that the
reliability of the results could be particularly dependent on the precision of
the underlying computer arithmetic system. Thus, the evaluation of
determinants as well as the inversion of matrices can become numerically
unstable operations. As a consequence, $Z^{\alpha}$ cannot be guaranteed to
have a fixed sign for all values of $\beta t$. We plot the average sign of
$Z^{\alpha}$ in Fig. 1(b), which is found to remain positive for all values of
$\beta t$, though fluctuations do begin to arise briefly before the simulation
cutoff. Fig. 1(c) showcases the expectation value of the Hamiltonian. Some
deviations from the exact treatment are observed in addition to a steadily
increasing uncertainty as $\beta$ increases.

The quality of the results drastically change when applied at half-filling
($\mu=2t$). The rate of exponential growth demonstrated in Fig. 2(a) is
substantially steeper than the previous configuration. In addition, we find
that $Z^{\alpha}$ has equal probability to be either positive or negative
beyond a certain value of $\beta t$. This is shown in Fig. 2(b), where the
mean sign of $Z^{\alpha}$ is zero above $\beta t\gtrsim 12$. While for small
values of $\beta t$ the statistical weight is positive for all trajectories,
its sign gets rapidly randomised as $\beta t$ increases. Thus a sign problem
re-emerges. Both the aforementioned issues seem to be inherent attributes of
this description of the system. Employing more sophisticated integration
techniques and increasing the lattice size not only failed to alleviate the
instability but rather made it more pronounced. We note that while
$Z^{\alpha}$ cannot be guaranteed to remain of fixed sign indefinitely, the
fact that $\langle
I+U^{T}_{\alpha}\tilde{U}_{\alpha}\rangle_{\text{Stochastic}}$ is positive
definite suggests that the sign problem might itself be a symptom of numerical
instability in the evaluation of the determinant rather than an independent
issue entirely.

In order to systematically account for these shortcomings and their influence
over the accuracy and precision of the calculations, the correlations are re-
defined to be

$\langle\hat{a}_{i,\sigma}^{\dagger}\hat{a}_{j,\sigma^{\prime}}\rangle=\left.\frac{\sum_{\alpha}Z^{\alpha}M^{\alpha}_{i\sigma,j\sigma^{\prime}}}{\sum_{\alpha}|Z^{\alpha}|}\middle/\frac{\sum_{\alpha}Z^{\alpha}}{\sum_{\alpha}|Z^{\alpha}|}\right.,$
(59)

where both the numerator and the denominator are to be evaluated separately.

This modified expectation value of the Hamiltonian is shown in Fig. 2(c).
There, we find that the average energy starts fluctuating at a relatively
small value of $\beta t$ due to the steady loss of precision in matrix
inversion combined with the strongly correlated dynamics. Furthermore, the
reliability of the results rapidly deteriorates as shown by the larger and
larger error bars and as expected from the prior observations. We thus
conclude that, for the most interesting cases, the Grassmann phase space
representation of the Hubbard model, though mathematically valid, falls short
of providing the level of numerical stability desirable for practical
applications.

##  VI  Conclusions

The theory of Grassmann phase space distributions was introduced and applied
to the grand-canonical ensemble of the Hubbard model. Two differential
equations were obtained to describe the dynamics of two related phase space
distributions in imaginary time. The $B$-distribution produced a Fokker-Planck
differential equation which can be re-cast into a stochastic evolution of
Grassmann variables.

Several critical limitations were found in performing numerical calculations
in this formulation. Firstly, the matrix elements of the stochastic propagator
were shown to grow exponentially. Secondly, the statistical weight associated
to each trajectory exhibited severe fluctuations in sign that produced
catastrophic cancellation. The uncertainty these findings introduced to
ensemble correlations was most critical in the non-perturbative parameter
regimes, where other numerical methods also face similar convergence issues.
Hence, our results show that any extension of the phase diagram beyond what is
presently known is not possible, at least within our current implementation.

Despite the numerical fragility encountered in this investigation, the
Grassmann phase space theory still demonstrates the possibility of mapping
Hamiltonians containing two-body interactions into stochastic differential
equations that are fitting for computational simulation. Its application is
deemed most suitable in circumstances involving a small number of interacting
modes. These are cases typically encountered in the quantum optics of two-
level systems, when the initial distribution is considerably restrictive in
its non-zero correlations, thereby sidestepping most of the highlighted
issues, or when the dynamics of the $P$-distribution renders it directly
amenable to linear stochastic methods.

Akin to similar probabilistic techniques, this approach achieves its best
accuracy at high temperatures (small imaginary times) before the development
of numerical instabilities. It is however necessary to ensure that the spectra
of the operators in Eq. (49) remains bounded when taking the thermodynamic
limit. Otherwise, the exponential run-off would present a problem even in such
a semi-classical regime.

##  Appendix A  Grassmann Algebra

The set of $n$ Grassmann variables, or Grassmann numbers, $g_{i}$, where
$i=1,\dots,n$, are defined to be the generators of the Grassmann algebra, or
exterior algebra, of an $n$ dimensional vector space. They compose
associatively and distributively, but satisfy the anti-commutation relation

$\\{g_{i},g_{j}\\}=0.$ (60)

The most general function of Grassmann variables is a multinomial of the form

$F(\bm{g})=f+\sum\limits_{1\leq i\leq n}f_{i}g_{i}+\sum\limits_{1\leq i,j\leq
n}f_{i,j}g_{i}g_{j}+~{}\dots,$ (61)

wherein all coefficients are elements of the underlying field, which is taken
to be the complex numbers. This definition fails to be unique unless further
conditions are specified to accommodate the anti-commutativity. This can be
achieved by restricting the summations to follow some prescribed ordering,
such as $i<j<\dots$, or by demanding total anti-symmetry in the coefficients
themselves so that $f_{i,j,\dots}=\varepsilon_{i,j,\dots}f_{1,2,\dots}$.

Multinomials in general neither commute nor anti-commute with Grassmann
variables or with each other. However, they may always be expressed as
$F=F_{E}+F_{O}$, where $F_{E}g_{i}=g_{i}F_{E}$ and contains all the even
summands and $F_{O}g_{i}=-g_{i}F_{O}$ and contains all the odd summands.

Further defined is another set of Grassmann variables $g_{i}^{*}$ taken as
conjugate companions of the former. The two sets are independent [26, 27] and
anti-commute amidst themselves as well as between each other. They map via
conjugation which is an involution that acts accordingly

$(g_{i}g_{j}\dots g_{u}g_{v})^{*}=g_{v}^{*}g_{u}^{*}\dots g_{j}^{*}g_{i}^{*},\
(g_{i}^{*})^{*}=g_{i}.$ (62)

##  Appendix B  Grassmann Calculus

It is possible to extend the notions of calculus to deal with Grassmann
quantities. However, their anti-commuting nature makes it is difficult to
assign the usual interpretations associated with calculus over the real or
complex numbers to such relations. What is obtained is a formal system wherein
differentiation and integration act equivalently rather than as inverses of
one another.

###  B.1  Grassmann Differentiation

Grassmann differentiation may act either from the left or from the right,
unlike regular differentiation which acts identically in both cases, under
these rules

$\displaystyle\begin{gathered}\overrightarrow{\frac{\partial}{\partial
g_{i}}}g_{j}=\delta_{i,j}~{},(\overrightarrow{\frac{\partial}{\partial
g_{i}}}g_{j}F)=(\overrightarrow{\frac{\partial}{\partial
g_{i}}}g_{j})F-g_{j}(\overrightarrow{\frac{\partial}{\partial g_{i}}}F),\\\
g_{j}\overleftarrow{\frac{\partial}{\partial
g_{i}}}=\delta_{i,j}~{},(Fg_{j}\overleftarrow{\frac{\partial}{\partial
g_{i}}})=F(g_{j}\overleftarrow{\frac{\partial}{\partial
g_{i}}})-(F\overleftarrow{\frac{\partial}{\partial g_{i}}})g_{j},\\\
\overrightarrow{\frac{\partial}{\partial
g_{i}}}1=1\overleftarrow{\frac{\partial}{\partial g_{i}}}=0.\end{gathered}$
(66)

Informally stated, the anti-commutative property of the Grassmann variables is
used to move $g_{i}$ to either the leftmost or rightmost position in the
expression and then it cancels with the derivative operator. Alternatively,
the derivative operator itself may be treated as anti-commuting with the
Grassmann variables until it is allowed to reach its associated $g_{i}$ and
cancel it.

This definition of Grassmann differentiation implies the following useful
results, in which $F_{X}$ is either even or odd with $\mathrm{sign}(X)$ being
$+1$ or $-1$ respectively

$\overrightarrow{\frac{\partial}{\partial
g_{i}}}(\overrightarrow{\frac{\partial}{\partial
g_{j}}}F)=-\overrightarrow{\frac{\partial}{\partial
g_{j}}}(\overrightarrow{\frac{\partial}{\partial g_{i}}}F),$ (67)
$(F\overleftarrow{\frac{\partial}{\partial
g_{j}}})\overleftarrow{\frac{\partial}{\partial
g_{i}}}=-(F\overleftarrow{\frac{\partial}{\partial
g_{i}}})\overleftarrow{\frac{\partial}{\partial g_{j}}},$ (68)
$(\overrightarrow{\frac{\partial}{\partial
g_{i}}}F_{X}H)=(\overrightarrow{\frac{\partial}{\partial
g_{i}}}F_{X})H+\mathrm{sign}(X)F_{X}(\overrightarrow{\frac{\partial}{\partial
g_{i}}}H),$ (69) $(HF_{X})\overleftarrow{\frac{\partial}{\partial
g_{i}}}=H(F_{X}\overleftarrow{\frac{\partial}{\partial
g_{i}}})+\mathrm{sign}(X)(H\overleftarrow{\frac{\partial}{\partial
g_{i}}})F_{X},$ (70) $\overrightarrow{\frac{\partial}{\partial
g_{i}}}F_{X}=-\mathrm{sign}(X)F_{X}\overleftarrow{\frac{\partial}{\partial
g_{i}}}.$ (71)

Though it might seem arbitrary, as Grassmann variables lack a natural notion
of scale, it proves useful to consider the deformation of a Grassmann function
$\delta F$ under the transformation $\bm{g}\to\bm{g}+\bm{\delta g}$. The
$\delta g_{i}$ themselves are Grassmann variables, or even more generally any
Grassmann odd function. An expression reminiscent of a Taylor series is
obtained, where careful attention ought to be paid to ensure the correct
ordering of the terms in these relations

$\displaystyle\delta F=\sum\limits_{1\leq i\leq n}\delta
g_{i}\overrightarrow{\frac{\partial}{\partial
g_{i}}}F+\frac{1}{2!}\sum\limits_{1\leq i,j\leq n}\delta g_{i}\delta
g_{j}\overrightarrow{\frac{\partial}{\partial
g_{j}}}\overrightarrow{\frac{\partial}{\partial g_{i}}}F+~{}\dots,$ (72)
$\displaystyle\delta F=\sum\limits_{1\leq i\leq
n}F\overleftarrow{\frac{\partial}{\partial g_{i}}}\delta
g_{i}+\frac{1}{2!}\sum\limits_{1\leq i,j\leq
n}F\overleftarrow{\frac{\partial}{\partial
g_{i}}}\overleftarrow{\frac{\partial}{\partial g_{j}}}\delta g_{j}\delta
g_{i}+~{}\dots\ .$ (73)

###  B.2  Grassmann Integration

Grassmann integration too may act from either direction, though the present
focus is on left integration as it ties more closely with the phase space
theory. It obeys the following rules

$\displaystyle\begin{gathered}\int dg_{i}\ g_{j}=\delta_{i,j}~{},\int dg_{i}\
1=0,\\\ \left(\int dg_{i}\ g_{j}F\right)=\left(\int dg_{i}\
g_{j}\right)F-g_{j}\left(\int dg_{i}\ F\right),\end{gathered}$ (76)

where the parentheses indicate which terms are being included in the
integrand. These are entirely analogous to the rules governing left
differentiation in Eqs. (66). Differentials may be treated as anti-commuting
between themselves and the Grassmann variables, and nested integrals are
evaluated from the innermost expression

$\int dg_{i}\left(\int dg_{j}\ F\right)=-\int dg_{j}\left(\int dg_{i}\
F\right),$ (77)

maintaining consistency with the anti-commutativity stated in Eq. (67).

Furthermore, as the derivative of a Grassmann function with respect to the
Grassmann variable $g_{i}$ does not contain any instance of said variable,
integrating Eq. (69) and Eq. (70) hence gives

$\int dg_{i}\ (\overrightarrow{\frac{\partial}{\partial
g_{i}}}F_{X})H=-\mathrm{sign}(X)\int dg_{i}\
F_{X}(\overrightarrow{\frac{\partial}{\partial g_{i}}}H),$ (78) $\int dg_{i}\
H(F_{X}\overleftarrow{\frac{\partial}{\partial g_{i}}})=-\mathrm{sign}(X)\int
dg_{i}\ (H\overleftarrow{\frac{\partial}{\partial g_{i}}})F_{X}.$ (79)

##  Appendix C  Phase Space Dynamics

To show how the differential equations arise, the two components of the
Hamiltonian in Eq. (28) are considered in turn and the correspondences in Eqs.
(22) are applied. The kinetic part produces

$\displaystyle\sum\limits_{i,j,\sigma}\bm{\tau}_{i,j}\hat{a}_{i,\sigma}^{\dagger}\hat{a}_{j,\sigma}\hat{\rho}\to\sum\limits_{i,j,\sigma}\bm{\tau}_{i,j}(\overrightarrow{\frac{\partial}{\partial
g_{i\sigma}}}-\tilde{g}_{i\sigma})g_{j\sigma}P=\sum\limits_{i,j,\sigma}\bm{\tau}_{i,j}(Pg_{j\sigma}\overleftarrow{\frac{\partial}{\partial
g_{i\sigma}}}-P\tilde{g}_{i\sigma}g_{j\sigma}),$
$\displaystyle\sum\limits_{i,j,\sigma}\bm{\tau}_{i,j}\hat{\rho}\hat{a}_{i,\sigma}^{\dagger}\hat{a}_{j,\sigma}\to\sum\limits_{i,j,\sigma}\bm{\tau}_{i,j}P\tilde{g}_{i\sigma}(\overleftarrow{\frac{\partial}{\partial\tilde{g}_{j\sigma}}}-g_{j\sigma})=\sum\limits_{i,j,\sigma}\bm{\tau}_{i,j}(P\tilde{g}_{i\sigma}\overleftarrow{\frac{\partial}{\partial\tilde{g}_{j\sigma}}}-P\tilde{g}_{i\sigma}g_{j\sigma}),$

where the even nature of the $P$-distribution has been used to apply the
derivative swap in Eq. (71). Hence, taking the sum of these two terms gives

$\sum\limits_{i,j,\sigma}\bm{\tau}_{i,j}(Pg_{j\sigma}\overleftarrow{\frac{\partial}{\partial
g_{i\sigma}}}+P\tilde{g}_{i\sigma}\overleftarrow{\frac{\partial}{\partial\tilde{g}_{j\sigma}}}-2P\tilde{g}_{i\sigma}g_{j\sigma}).$
(80)

The interaction potential follows, it leads to some prolonged expressions that
are possible to collapse down

$\displaystyle\sum\limits_{i,j,\sigma,\sigma^{\prime}}\bm{U}_{i\sigma,j\sigma^{\prime}}\hat{a}_{i,\sigma}^{\dagger}\hat{a}_{j,\sigma^{\prime}}^{\dagger}\hat{a}_{j,\sigma^{\prime}}\hat{a}_{i,\sigma}\hat{\rho}\to\sum\limits_{i,j,\sigma,\sigma^{\prime}}\bm{U}_{i\sigma,j\sigma^{\prime}}(\overrightarrow{\frac{\partial}{\partial
g_{i\sigma}}}-\tilde{g}_{i\sigma})(\overrightarrow{\frac{\partial}{\partial
g_{j\sigma^{\prime}}}}-\tilde{g}_{j\sigma^{\prime}})g_{j\sigma^{\prime}}g_{i\sigma}P,$
$\displaystyle\sum\limits_{i,j,\sigma,\sigma^{\prime}}\bm{U}_{i\sigma,j\sigma^{\prime}}\hat{\rho}\hat{a}_{i,\sigma}^{\dagger}\hat{a}_{j,\sigma^{\prime}}^{\dagger}\hat{a}_{j,\sigma^{\prime}}\hat{a}_{i,\sigma}\to\sum\limits_{i,j,\sigma,\sigma^{\prime}}\bm{U}_{i\sigma,j\sigma^{\prime}}P\tilde{g}_{i\sigma}\tilde{g}_{j\sigma^{\prime}}(\overleftarrow{\frac{\partial}{\partial\tilde{g}_{j\sigma^{\prime}}}}-g_{j\sigma^{\prime}})(\overleftarrow{\frac{\partial}{\partial\tilde{g}_{i\sigma}}}-g_{i\sigma}),$

considering the cross terms, it is found after re-arranging that

$\displaystyle\bm{U}_{i\sigma,j\sigma^{\prime}}\tilde{g}_{i\sigma}\overrightarrow{\frac{\partial}{\partial
g_{j\sigma^{\prime}}}}g_{j\sigma^{\prime}}g_{i\sigma}P=\bm{U}_{i\sigma,j\sigma^{\prime}}\overrightarrow{\frac{\partial}{\partial
g_{j\sigma^{\prime}}}}\tilde{g}_{i\sigma}g_{i\sigma}g_{j\sigma^{\prime}}P,$
$\displaystyle\bm{U}_{i\sigma,j\sigma^{\prime}}P\tilde{g}_{i\sigma}\tilde{g}_{j\sigma^{\prime}}\overleftarrow{\frac{\partial}{\partial\tilde{g}_{j\sigma^{\prime}}}}g_{i\sigma}=\bm{U}_{i\sigma,j\sigma^{\prime}}P\tilde{g}_{j\sigma^{\prime}}\tilde{g}_{i\sigma}g_{i\sigma}\overleftarrow{\frac{\partial}{\partial\tilde{g}_{j\sigma^{\prime}}}},$

then noting
$\bm{U}_{i\sigma,j\sigma^{\prime}}=\bm{U}_{j\sigma^{\prime},i\sigma}$,
switching the direction of the derivatives as per Eq. (71), and combining all
the terms

$\sum\limits_{i,j,\sigma,\sigma^{\prime}}\bm{U}_{i\sigma,j\sigma^{\prime}}(Pg_{i\sigma}g_{j\sigma^{\prime}}\overleftarrow{\frac{\partial}{\partial
g_{j\sigma^{\prime}}}}\overleftarrow{\frac{\partial}{\partial
g_{i\sigma}}}+P\tilde{g}_{i\sigma}\tilde{g}_{j\sigma^{\prime}}\overleftarrow{\frac{\partial}{\partial\tilde{g}_{j\sigma^{\prime}}}}\overleftarrow{\frac{\partial}{\partial\tilde{g}_{i\sigma}}}-2P\tilde{g}_{j\sigma^{\prime}}g_{j\sigma^{\prime}}g_{i\sigma}\overleftarrow{\frac{\partial}{\partial
g_{i\sigma}}}-2P\tilde{g}_{i\sigma}\tilde{g}_{j\sigma^{\prime}}g_{j\sigma^{\prime}}\overleftarrow{\frac{\partial}{\partial\tilde{g}_{i\sigma}}}+2P\tilde{g}_{i\sigma}\tilde{g}_{j\sigma^{\prime}}g_{j\sigma^{\prime}}g_{i\sigma}).$
(81)

Finally, bringing Eqs. (80) and (81) together as prescribed in Eq. (32) gives

$\displaystyle\frac{d}{d\beta}P=-\frac{1}{2}$
$\displaystyle[-\sum\limits_{i,j,\sigma}\bm{\tau}_{i,j}(Pg_{j\sigma}\overleftarrow{\frac{\partial}{\partial
g_{i\sigma}}}+P\tilde{g}_{i\sigma}\overleftarrow{\frac{\partial}{\partial\tilde{g}_{j\sigma}}}-2P\tilde{g}_{i\sigma}g_{j\sigma})+\sum\limits_{i,j,\sigma,\sigma^{\prime}}\bm{U}_{i\sigma,j\sigma^{\prime}}(Pg_{i\sigma}g_{j\sigma^{\prime}}\overleftarrow{\frac{\partial}{\partial
g_{j\sigma^{\prime}}}}\overleftarrow{\frac{\partial}{\partial g_{i\sigma}}}$
(82)
$\displaystyle+P\tilde{g}_{i\sigma}\tilde{g}_{j\sigma^{\prime}}\overleftarrow{\frac{\partial}{\partial\tilde{g}_{j\sigma^{\prime}}}}\overleftarrow{\frac{\partial}{\partial\tilde{g}_{i\sigma}}}-2P\tilde{g}_{j\sigma^{\prime}}g_{j\sigma^{\prime}}g_{i\sigma}\overleftarrow{\frac{\partial}{\partial
g_{i\sigma}}}-2P\tilde{g}_{i\sigma}\tilde{g}_{j\sigma^{\prime}}g_{j\sigma^{\prime}}\overleftarrow{\frac{\partial}{\partial\tilde{g}_{i\sigma}}}+2P\tilde{g}_{i\sigma}\tilde{g}_{j\sigma^{\prime}}g_{j\sigma^{\prime}}g_{i\sigma})].$

The other differential equation is found by noting that that the Hubbard
Hamiltonian in Eq. (28) is number conserving and that the $B$-distribution
correspondences in Eqs. (27) only change the total Grassmann order by either
$+1$ or $-1$. Consequentially, omitting the non-zero order terms in Eq. (82)
then directly gives

$\frac{d}{d\beta}B=-\frac{1}{2}[-\sum\limits_{i,j,\sigma}\bm{\tau}_{i,j}(Bg_{j\sigma}\overleftarrow{\frac{\partial}{\partial
g_{i\sigma}}}+B\tilde{g}_{i\sigma}\overleftarrow{\frac{\partial}{\partial\tilde{g}_{j\sigma}}})+\sum\limits_{i,j,\sigma,\sigma^{\prime}}\bm{U}_{i\sigma,j\sigma^{\prime}}(Bg_{i\sigma}g_{j\sigma^{\prime}}\overleftarrow{\frac{\partial}{\partial
g_{j\sigma^{\prime}}}}\overleftarrow{\frac{\partial}{\partial
g_{i\sigma}}}+B\tilde{g}_{i\sigma}\tilde{g}_{j\sigma^{\prime}}\overleftarrow{\frac{\partial}{\partial\tilde{g}_{j\sigma^{\prime}}}}\overleftarrow{\frac{\partial}{\partial\tilde{g}_{i\sigma}}})].$
(83)

## References

  * Glauber [1963] R. J. Glauber, Coherent and incoherent states of the radiation field, Phys. Rev. 131, 2766 (1963).
  * Sudarshan [1963] E. C. G. Sudarshan, Equivalence of semiclassical and quantum mechanical descriptions of statistical light beams, Phys. Rev. Lett. 10, 277 (1963).
  * HUSIMI [1940] K. HUSIMI, Some formal properties of the density matrix, Proceedings of the Physico-Mathematical Society of Japan. 3rd Series 22, 264 (1940).
  * Wigner [1932] E. Wigner, On the quantum correction for thermodynamic equilibrium, Phys. Rev. 40, 749 (1932).
  * Moyal [1949] J. E. Moyal, Quantum mechanics as a statistical theory, Mathematical Proceedings of the Cambridge Philosophical Society 45, 99–124 (1949).
  * Drummond and Carter [1987] P. D. Drummond and S. J. Carter, Quantum-field theory of squeezing in solitons, J. Opt. Soc. Am. B 4, 1565 (1987).
  * Opanchuk  _et al._ [2012] B. Opanchuk, M. Egorov, S. Hoffmann, A. I. Sidorov, and P. D. Drummond, Quantum noise in three-dimensional bec interferometry, Europhysics Letters 97, 50003 (2012).
  * Carusotto  _et al._ [2008] I. Carusotto, S. Fagnocchi, A. Recati, R. Balbinot, and A. Fabbri, Numerical observation of hawking radiation from acoustic black holes in atomic bose–einstein condensates, New Journal of Physics 10, 103001 (2008).
  * Cahill and Glauber [1999] K. E. Cahill and R. J. Glauber, Density operators for fermions, Phys. Rev. A 59, 1538 (1999).
  * Dalton  _et al._ [2014] B. J. Dalton, J. Jeffers, and S. M. Barnett,  _Phase Space Methods for Degenerate Quantum Gases_ (Oxford University Press, 2014).
  * Dalton  _et al._ [2016] B. Dalton, J. Jeffers, and S. Barnett, Grassmann phase space methods for fermions. i. mode theory, Annals of Physics 370, 12 (2016).
  * Anastopoulos and Hu [2000] C. Anastopoulos and B. L. Hu, Two-level atom-field interaction: Exact master equations for non-markovian dynamics, decoherence, and relaxation, Phys. Rev. A 62, 033821 (2000).
  * Shresta  _et al._ [2005] S. Shresta, C. Anastopoulos, A. Dragulescu, and B. L. Hu, Non-markovian qubit dynamics in a thermal field bath: Relaxation, decoherence, and entanglement, Phys. Rev. A 71, 022109 (2005).
  * Dalton  _et al._ [2013] B. Dalton, B. Garraway, J. Jeffers, and S. Barnett, Grassmann phase space theory and the jaynes–cummings model, Annals of Physics 334, 100 (2013).
  * Corney and Drummond [2004] J. F. Corney and P. D. Drummond, Gaussian quantum monte carlo methods for fermions and bosons, Phys. Rev. Lett. 93, 260401 (2004).
  * Corney and Drummond [2005] J. F. Corney and P. D. Drummond, Gaussian operator bases for correlated fermions, Journal of Physics A: Mathematical and General 39, 269 (2005).
  * Corney and Drummond [2006] J. F. Corney and P. D. Drummond, Gaussian phase-space representations for fermions, Phys. Rev. B 73, 125112 (2006).
  * Assaad  _et al._ [2005] F. F. Assaad, P. Werner, P. Corboz, E. Gull, and M. Troyer, Symmetry projection schemes for gaussian monte carlo methods, Phys. Rev. B 72, 224518 (2005).
  * Plimak  _et al._ [2001] L. I. Plimak, M. J. Collett, and M. K. Olsen, Langevin equations for interacting fermions and cooper-like pairing in trapped one-dimensional fermions, Phys. Rev. A 64, 063409 (2001).
  * Arovas  _et al._ [2022] D. P. Arovas, E. Berg, S. A. Kivelson, and S. Raghu, The hubbard model, Annual Review of Condensed Matter Physics 13, 239 (2022).
  * Hubbard [1963] J. Hubbard, Electron correlations in narrow energy bands, Proceedings of the Royal Society of London. Series A, Mathematical and Physical Sciences 276, 238 (1963).
  * Gutzwiller [1963] M. C. Gutzwiller, Effect of correlation on the ferromagnetism of transition metals, Phys. Rev. Lett. 10, 159 (1963).
  * Kanamori [1963] J. Kanamori, Electron Correlation and Ferromagnetism of Transition Metals, Progress of Theoretical Physics 30, 275 (1963).
  * Gardiner [2009] C. W. Gardiner,  _Stochastic methods : a handbook for the natural and social sciences_ , fourth edition. ed., Springer series in synergetics (Unnumbered) (Springer, Berlin, 2009).
  * Zinn-Justin [2002] J. Zinn-Justin,  _Quantum Field Theory and Critical Phenomena_ (Oxford University Press, 2002).
  * Berezin [1966] F. A. Berezin,  _The method of second quantization_ , Pure and applied physics, v. 24 (Academic Press, New York, 1966).
  * Altland and Simons [2010] A. Altland and B. D. Simons,  _Condensed Matter Field Theory_, 2nd ed. (Cambridge University Press, 2010).


# An analytical approach to the Mean–Return-Time Phase of isotropic stochastic
oscillators

Konstantin Holzhausen  Bernstein Center for Computational Neuroscience Berlin,
Philippstr. 13, Haus 2, 10115 Berlin, Germany  Physics Department of Humboldt
University Berlin, Newtonstr. 15, 12489 Berlin, Germany  Peter J. Thomas
Department of Mathematics, Applied Mathematics and Statistics, 212 Yost Hall,
Case Western Reserve University, 10900 Euclid Avenue, Cleveland, Ohio, USA
Benjamin Lindner  Bernstein Center for Computational Neuroscience Berlin,
Philippstr. 13, Haus 2, 10115 Berlin, Germany  Physics Department of Humboldt
University Berlin, Newtonstr. 15, 12489 Berlin, Germany

###### Abstract

One notion of phase for stochastic oscillators is based on the mean return-
time (MRT): a set of points represents a certain phase if the mean time to
return from any point in this set to this set after one rotation is equal to
the mean rotation period of the oscillator (irrespective of the starting
point). For this so far only algorithmically defined phase, we derive here
analytical expressions for the important class of isotropic stochastic
oscillators. This allows us to evaluate cases from the literature explicitly
and to study the behavior of the MRT phase in the limits of strong noise. We
also use the same formalism to show that lines of constant return time
variance (instead of constant mean return time) can be defined, and that they
in general differ from the MRT-isochrons.

Oscillatory processes encountered in physics and biology often come along with
a substantial component of randomness and are then termed _stochastic
oscillations_. Intracellular calcium waves SkuKet08 , neural action potentials
WalAlv06 , rhythmic fluctuations in lasing intensity McNWie88 , and population
oscillations of prey-predator systems McKNew05  all show pronounced
fluctuations in their amplitude as well as in the timing of the cycles.

In order to find a reduced description of such random oscillations, in recent
years efforts have been made to generalize the phase concept for deterministic
oscillators Win74 ; Guc75  (which has also been applied to systems with weak
noise, see e.g. FreNei00 ; FreSch03 ; ZhoBur13 ; MaKli14 ) to the stochastic
case SchPik13  (see also ThoLin14 ; GiaPoqSha2018 ; AmiHolSri2019 ; EngKue2021
). The first proposal, the mean–return-time (MRT) phase by Schwabedal and
Pikovsky SchPik13 , is an intuitive generalization of the stroboscopic return-
time-phase definition for deterministic oscillators: in the two-dimensional
case, with any starting point on the line of equal phase (isochron), the mean
return-time to the very same line after one rotation is equal to the mean
rotation period of the oscillator. Schwabedal and Pikovsky suggested an
algorithmic procedure to deform a line until (experimental or simulated) data
satisfy this criterion. Cao et al. CaoLin20  showed more recently, that the
MRT phase for planar white-Gaussian noise driven oscillators obeys a partial
differential equation with an unusual jump condition.

Based on CaoLin20  we present an analytical solution for the MRT phase for the
important class of isotropic stochastic oscillators. The solution is given in
terms of quadratures and is exploited to explore the limit of strong noise and
the role of boundaries for so-far purely numerically treated examples from the
literature. Finally, our approach allows us to go beyond the familiar notion
of isochrons (lines of constant mean rotation time) to introduce a new
construct, the _iso-variance_ lines (lines of constant variance of the
rotation time). We demonstrate that the iso-variance lines strongly differ
from the isochrons.

Derivation of the isochron expression.– We consider a class of rotationally
invariant planar stochastic oscillators (also called isotropic stochastic
oscillators) that obey

$\displaystyle\dot{\rho}$
$\displaystyle=g(\rho)+q_{\rho}(\rho)\,\xi_{\rho}(t)$ (1)
$\displaystyle\dot{\phi}$
$\displaystyle=f(\rho)+q_{\phi}(\rho)\,\xi_{\phi}(t)$

in polar coordinates, where $\xi_{\rho,\phi}$ denote white Gaussian noise
sources with
$\left\langle\xi_{i}(t)\,\xi_{j}(t^{\prime})\right\rangle=\delta_{i,j}\,\delta(t-t^{\prime}),\;i,j\in\\{\rho,\,\phi\\}$.
If the noise is multiplicative ($q_{\rho}(\rho)\neq\text{const}$) we interpret
it in the sense of Ito Gar85 . We impose reflecting boundaries at an inner
circle with radius $\rho_{-}$ and an outer circle $\rho_{+}$ such that the
dynamics is restricted to $\rho_{-}<\rho(t)<\rho_{+}$. We assume that there is
a mean rotation around the origin, the rotational sense of which may change
depending on the parameters.

Because none of the functions on the right-hand side depend explicitely on the
phase $\phi$, the stochastic dynamics is rotationally invariant, which is why
we refer to eq. (1) as an isotropic oscillator. In particular, we expect that
the isochrons for different phases $\psi\in[0,2\pi]$ will be rotationally
invariant as well, and can be expressed via
$\varphi(\rho)=\phi_{\text{I}}(\rho)+\psi$. In the following we calculate
$\phi_{\text{I}}(\rho)$.

According to CaoLin20 , the isochrons are the level sets of the mean-first-
passage-time function $T(\rho,\,\phi)$ that is the solution to the partial
differential equation

$\mathcal{L}^{\dagger}\,T(\rho,\,\phi)=\left[g\partial_{\rho}+f\partial_{\phi}+\frac{q_{\rho}^{2}}{2}\partial_{\rho}^{2}+\frac{q_{\phi}^{2}}{2}\partial_{\phi}^{2}\right]T(\rho,\,\phi)=-1\,$
(2)

satisfying the periodic-plus-jump condition at an arbitrary connection between
the inner and outer boundaries (here we choose for simplicity the line
$\phi=0$ and consider a counterclockwise rotation):

$T(\rho,\,2\pi)+\overline{T}=T(\rho,0)$ (3)

and adjoint reflecting boundary conditions at $\rho_{\pm}$. $\overline{T}$
denotes the mean period, the time it takes one realization of the process on
average to perform one revolution around the origin. Following CaoLin20 , we
assume $f,g,q_{\rho},q_{\phi}\in\mathcal{C}^{2}$ as well as
$q_{\rho}^{2}(\rho)>0$ and $q_{\phi}^{2}(\rho)>0$ for all
$\rho\in[\rho_{-},\,\rho_{+}]$. The choice of polar coordinates translates the
rotational invariance of eq. (1) into a translational invariance in the
angular coordinate $\phi$. Correspondingly, a symmetry-adapted solution is of
the general form

$T(\rho,\,\phi)=-\frac{\overline{T}}{2\,\pi}\,\phi+T_{\rho}(\rho)$ (4)

that also fulfills the periodic-plus-jump condition. This symmetry-adapted
ansatz reduces the number of dimensions in eq. (2), giving an ordinary
differential equation for $T_{\rho}$ in $\rho$, which can be solved explicitly
in terms of quadratures. The isochron $\phi_{\textrm{I}}$ as a special level
set of $T(\rho,\phi)$ can then be parametrized according to
$\phi_{\textrm{I}}(\rho)=2\pi T_{\rho}(\rho)/\overline{T}$, leading to our
main result:

$\phi_{\textrm{I}}(\rho)\\!=\\!2\\!\int\limits_{\rho_{-}}^{\rho}\\!dq\\!\int\limits_{\rho_{-}}^{q}\\!du\frac{f(u)-\overline{\omega}}{q_{\rho}^{2}(u)}\exp\left[{-2\int_{u}^{q}dv\frac{g(v)}{q_{\rho}^{2}(v)}}\right].$
(5)

The mean rotation frequency $\overline{\omega}$ (essentially, the inverse mean
rotation period $2\pi/\bar{T}$) is given by

$\overline{\omega}=\frac{2\pi}{\overline{T}}=\frac{\int_{\rho_{-}}^{\rho_{+}}d\rho\,f(\rho)/q_{\rho}^{2}(\rho)\,e^{-2\int_{\rho}^{\rho_{+}}d\rho^{\prime}q_{\rho}^{-2}(\rho^{\prime})\,g(\rho^{\prime})}}{\int_{\rho_{-}}^{\rho_{+}}d\rho\,q_{\rho}^{-2}(\rho)\,e^{-2\int_{\rho}^{\rho_{+}}d\rho^{\prime}q_{\rho}^{-2}(\rho^{\prime})\,g(\rho^{\prime})}}.$
(6)

We can make a few conclusions from the analytical result eq. (5) without any
numerical evaluation: i) $q_{\phi}$, i.e. the strength of phase diffusion,
does not enter the isochron parametrization at all (all oscillators that only
differ with respect to $q_{\phi}$ have the same isochrons); ii) the slope of
the phase $d\phi/d\rho$ is given by the difference between the local phase
progression speed $f(\rho)$ and the average speed $\overline{\omega}$ weighted
with the (normalized) steady state density of the radius and averaged over the
interval $(\rho_{-},\rho)$; iii) if the overall strength of the radial noise
increases without bounds and we keep the two boundary values $\rho_{\pm}$
fixed at non-vanishing finite values, the integrand in eq. (5) approaches
zero, i.e. the isochron $\phi_{\textrm{I}}$ does not depend on the radius
anymore and isochrons approach the shape of spokes of a wheel.

_Newby-Schwemmer oscillator.–_ This is given by NewSch14

$\displaystyle\dot{\rho}$
$\displaystyle=-\gamma\,\rho\,(\rho^{2}-1)+\frac{D}{\rho}+\sqrt{2\,D}\,\xi_{\rho}(t),$
(7) $\displaystyle\dot{\phi}$
$\displaystyle=\omega+\omega\,\gamma\,c\,(1-\rho)^{2}+\frac{\sqrt{2\,D}}{\rho}\,\xi_{\phi}(t).$

The sense of rotation is counterclockwise on the limit cycle at $\rho=1$ but,
for $c<0$, turns to clockwise when the trajectory deviates sufficiently both
towards the inside and the outside. With noise, the occupation probability
w.r.t. $\rho$ determines whether the net rotation goes clockwise or
counterclockwise NewSch14 , and at the transition point $D=D_{t}$ between both
regimes, the mean rotation time diverges (see Fig. 1c)

Figure 1: MRT-Isochrons of the Newby-Schwemmer model. (a): Isochron
parametrization $\phi_{\textrm{I}}(\rho)$, eq. (5). (b): Isochrons (colored),
reflecting boundaries (dotted) and limit cycle (solid black). (c): Mean period
$\overline{T}$ vs $D$ according to theory (lines) and simulations (circles).
(d): Mean return-time $\left\langle T(\rho_{0})\right\rangle$ to a curve vs
initial position $\rho_{0}$ on this curve, for the curve being an isochron
(circles, different colors represent different $D$) or a spoke (triangles,
here $D=0.1$).

The deterministic isochron (dark green line in Fig. 1a and b) has a hook-like
shape, because points away from the limit cycle need some head start (net
rotation is counterclockwise). For non-vanishing noise this principal shape is
maintained: for $D<D_{t}$ the net rotation is still counterclockwise, and
points off the limit cycle need a headstart; for $D>D_{t}$ these points need a
“later” start (net rotation is clockwise), thus the isochron has the same
curvature. In addition, with increasing noise intensity $D$ the isochrons
become flatter, approaching the spoke of a wheel (dashed black line) as
predicted by our general conclusion (iii) above.

We checked the constant-MRT property by performing extensive simulations with
ensembles of stochastic trajectories for a certain noise intensity $D$
starting on the corresponding isochron. Indeed, the MRT for a given noise
intensity from and to the isochron do not depend on the specific starting
point along the isochron (circles in Fig. 1d). In marked contrast, choosing
(for a comparatively small noise intensity) as initial set and target line the
spoke of a wheel, the mean return times _do_ depend on the initial position
(triangles in Fig. 1d).

Guckenheimer-Schwabedal-Pikovsky oscillator.– A classic example of a
deterministic system, in which the definition of phase is problematic, is one
with two limit cycles with distinct rotation frequencies (this example goes
back to Guckenheimer Guc75 ). As shown numerically by Schwabedal and Pikovsky
SchPik13 , with noise, transitions occur between the attraction domains of the
two limit cycles and thus the MRT phase uniquely defines one phase for the
entire system. This system is given by

$\displaystyle\dot{\rho}$
$\displaystyle=\rho\,(1-\rho)(3-\rho)(c-\rho)+\frac{\sigma^{2}\,\rho}{2}+\sigma\,\rho\xi(t)$
(8) $\displaystyle\dot{\phi}$
$\displaystyle=\omega+\delta\,(\rho-2)-(1-\rho)\,(3-\rho)$

Figure 2: MRT-Isochrons of the Guckenheimer-Schwabedal-Pikovsky oscillator.
Panels as in Fig. 1 but for eq. (8). Parameters: $\omega=3,c=1.8$,
$\rho_{-}=0.5,\rho_{+}=4.0$.

The unstable limit cycle at $\rho=c$ separates the basins of attraction at
$\rho=1$ and $\rho=3$. With noise, the unstable limit cycle can be surmounted,
smoothing the discontinuity, and thus resulting in one stochastic isochron
that is continuous over the whole domain connecting both basins of attraction
(see Fig. 2, a,b).

While MRTs for a given noise intensity remain constant for initial positions
along the corresponding isochron, they vary significantly if the initial set
and target line is a ray emanating from the origin (cf. Fig. 2d). For $D=0$,
the return time suffers a jump by $2\,\delta$ between the basins of attraction
(i.e. for an initial point $\rho=c$), reflecting the above mentioned fact that
there is no unique phase for the deterministic system. The stronger this jump
is, the more twisted the stochastic MRT isochron becomes (cf. lines for
$\delta=0.01$ and $\delta=0.2$ in Fig. 2a,b); $\delta$ may also change how the
mean rotation period depends on the noise level (cf. Fig. 2c) We note that for
strongly increasing noise, the shape of the isochron flattens as again
expected from our general consideration.

Figure 3: MRT-isochrons of the monomial model. Presentation in the panels is
as in Fig. 1 but for eq. (9) with $n=1,\rho_{-}=0,\rho_{+}=\infty$.
Simulations in (d) are, however, performed with $\rho_{-}=10^{-4}$
$\rho_{+}=4.0$ consequently, the MRT-property is only obeyed for initial
points far from the outer boundary.

A monomial model.– According to what we have found so far, we may expect, that
in the strong-noise limit (until now, the less explored regime of stochastic
oscillators), the isochrons always approach a radial spoke (conclusion (iii)
above). This, however, is only true for systems with finite boundaries. To
show this, we consider a particularly simple model system with monomial
frequency dependence, for which the integrals in eq. (5) can be expressed by
elementary functions. The model obeys

$\displaystyle\dot{\rho}$
$\displaystyle=-\alpha+\frac{D}{\rho}+\sqrt{2\,D}\,\xi_{\rho}(t)$ (9)
$\displaystyle\dot{\phi}$
$\displaystyle=-\beta\,\rho^{n}+\frac{\sqrt{2\,D}}{\rho}\,\xi_{\phi}(t)$

with positive parameters $\alpha,\beta$ and integer $n$, resulting in
clockwise rotations with a speed that increases as a power law with increasing
radius. We are specifically interested in the limit in which the inner
boundary approaches the origin and the outer boundary goes to infinity. In
this limit, the MRT and the isochrons can be expressed by elementary functions
Hol21 ; particularly simple and striking is the expression for $n=1$:

$\lim_{\begin{subarray}{c}\rho_{-}\rightarrow 0\\\
\rho_{+}\rightarrow\infty\\\
\alpha>0\end{subarray}}\phi(\rho)=\frac{\beta}{\alpha}\,\frac{\rho^{2}}{2},\quad\lim_{\begin{subarray}{c}\rho_{-}\rightarrow
0\\\ \rho_{+}\rightarrow\infty\\\
\alpha>0\end{subarray}}\overline{T}=\frac{\alpha}{2\beta D}$ (10)

Hence, quite surprisingly, the MRT phase is completely independent of the
noise level but clearly different from a simple ray (see Fig. 3a,b).
Specifically in the limit $D\to\infty$ the isochrons do not converge to the
spokes of wheel, which is due to the absence of finite boundaries: for
stronger and stronger noise, the main share of probability moves to larger and
larger radii and the mean rotation time drops strongly (see Fig. 3c).

We can test how important the infinite boundary condition is by running
stochastic simulations with initial set and target line with the isochron eq.
(10) but imposing a reflecting boundary at a finite value $\rho_{+}<\infty$
(for this setting, eq. (10) is _not_ the exact isochron). For all noise
intensities, the MRT is flat as a function of the initial radius except for a
finite region close to the outer boundary (see Fig. 3d).

Figure 4: Iso-variance curves for the Newby-Schwemmer oscillator.(a):
Isochronal curve and curve of constant return-time variance of the Newby-
Schwemmer model in terms of angle variable vs radius. (b) Same curves as in
(a), in cartesian coordinates. (c): Mean return time to isochrons, iso-
variance line, or spoke as function of the initial point on the respective
line. (d) Variances of the Return times $\left\langle\Delta^{2}\right\rangle>$
as function of initial position on the respective lines.

Iso-variance curves.– The same approach used for the isochrons permits to
calculate the lines of constant variance, that is, the lines for which the
return time to the line after one rotation possesses the same variance
irrespective of the starting point. It turns out that these lines are the
contour lines of the variance function, obeying

$L^{\dagger}\Delta=q_{\rho}^{2}\,\left(\partial_{\rho}\,T(\rho,\,\phi)\right)^{2}+q_{\phi}^{2}\,\left(\partial_{\phi}\,T(\rho,\,\phi)\right)^{2}$
(11)

with a jump condition
$\Delta(\rho,\,\phi)=\Delta(\rho,\,\phi+2\,\pi)+\overline{\Delta}$ which can
be solved for the isotropic oscillator by a similar ansatz
$\Delta(\rho)=\Delta_{\rho}(\rho)-\phi\bar{\Delta}/(2\pi)$ resulting in the
parametrization of the iso-variance curve

$\displaystyle\phi_{\Delta}(\rho)$
$\displaystyle=2\int\limits_{\rho_{-}}^{\rho}\\!\\!dq\int\limits_{\rho_{-}}^{q}\\!du\left[\frac{2\pi}{\overline{\Delta}}\left(q_{\rho}^{2}(u)\left(\partial_{u}T_{\rho}(u)\right)^{2}+\frac{q_{\phi}^{2}(u)\,\overline{T}^{2}}{4\pi^{2}}\right)\right.$
(12)
$\displaystyle\left.\frac{}{}-f(u)\right]\,\exp{\left(-2\int\limits_{u}^{q}dv\,\frac{g(v)}{q_{\rho}^{2}(v)}\right)}\Big{/}q_{\rho}^{2}(u)$

Here the variance of the rotation time $\bar{\Delta}$ is given by

$\bar{\Delta}=-2\pi\frac{\int\limits_{\rho_{-}}^{\rho_{+}}d\rho\,\left[\left(\partial_{\rho}T_{\rho}(\rho)\right)^{2}+\frac{q_{\rho}^{2}}{q_{\phi}^{2}}\frac{\overline{T}^{2}}{4\pi^{2}}\right]\,\exp\left(-2\int\limits_{\rho}^{\rho_{+}}du\,\frac{g(u)}{q_{\rho}^{2}(u)}\right)}{\int\limits_{\rho_{-}}^{\rho_{+}}d\rho\,f(\rho)\exp\left(-2\int\limits_{\rho}^{\rho_{+}}du\,\frac{g(u)}{q_{\rho}^{2}(u)}\right)/q^{2}_{\rho}(\rho)}$
(13)

We have evaluated and tested these formulas for the Newby-Schwemmer oscillator
eq. (7) and also compare the resulting iso-variance lines to the isochrons and
the spoke of wheel in Fig. 4. We observe that isochron and iso-variance curve
do not agree in this case: the iso-variance curve (pink line in Fig. 4a,b) is
significantly less twisted than the isochron (green line). In panels c,d we
test the defining properties of the two lines: the MRT is approximately
independent of the starting point only for the isochron but neither for the
iso-variance curve nor the spoke (c); the variance of the return time is
approximately flat as a function of the starting point for the iso-variance
line but neither for the isochron nor the spoke.

_Conclusions.–_ We have found the analytical mapping for Schwabedal and
Pikovsky’s mean return-time (MRT) phase for the important class of planar
isotropic stochastic oscillators and have tested it for three examples. We
have seen that for systems with finite boundaries, the phase description in
the strong noise limit always yields the geometric phase (spokes of a wheel);
for small to moderate noise, the curvature of the MRT phase isochrons reflect
the stochastic interplay between radial and angular dynamics. We also
demonstrated that constancy of the first cumulant (the MRT) over the isochron
does not imply constancy of higher cumulants, for instance, the variance of
the return time for different initial positions on the isochron. We could find
an equation and solve it for the iso-variance curve: in our example the
resulting line was quite different from the MRT isochron.

We hope that our results can be used to address important open issues in the
analysis of stochastic oscillators: i) the comparison to the asymptotic phase
ThoLin14  in some analytically tractable cases (see e.g. ThoLin19 ); ii) the
oscillator response to external perturbations via a generalized phase-response
curve GalErm05 ; ErmGal07 ; iii the analysis of variability in neuronal firing
patterns PuTho2021 ; iv) the treatment of networks of strongly stochastic
oscillators ErSau06 .

## References

  * (1) A. Skupin, H. Kettenmann, U. Winkler, M. Wartenberg, H. Sauer, S. C. Tovey, C. W. Taylor, and M. Falcke.  How does intracellular ca2+ oscillate: By chance or by the clock?  Biophys. J., 94:2404, 2008.
  * (2) J. T. Walter, K. Alvina, M. D. Womack, C. Chevez, and K. Khodakhah.  Decreases in the precision of Purkinje cell pacemaking cause cerebellar dysfunction and ataxia.  Nat. Neurosci, 9:389, 2006.
  * (3) B. McNamara, K. Wiesenfeld, and R. Roy.  Observation of stochastic resonance in a ring laser.  Phys. Rev. Lett., 60:2626, 1988.
  * (4) A. J. McKane and T. J. Newman.  Predator-prey cycles from resonant amplification of demographic stochasticity.  Phys. Rev. Lett., 94(21):218102, 2005.
  * (5) AT Winfree.  Patterns of phase compromise in biological cycles.  Journal of Mathematical Biology, 1(1):73–93, 1974.
  * (6) J. Guckenheimer.  Isochrons and Phaseless Sets.  J. Math. Biol., 1(3):259, 1975.
  * (7) J. A. Freund, A. B. Neiman, and L. Schimansky-Geier.  Analytic description of noise-induced phase synchronization.  Europhys. Lett., 50:8, 2000.
  * (8) J. A. Freund, L. Schimansky-Geier, and P. Hänggi.  Frequency and phase synchronization in stochastic systems.  Chaos, 13:225, 2003.
  * (9) P. Zhou, S. D. Burton, N. Urban, and G. B. Ermentrout.  Impact of neuronal heterogeneity on correlated colored-noise-induced synchronization.  Front. Comput. Neurosci., 7:113, 2013.
  * (10) R. Ma, G. S. Klindt, I. H. Riedel-Kruse, F. Jülicher, and B. M. Friedrich.  Active Phase and Amplitude Fluctuations of Flagellar Beating.  Phys. Rev. Lett., 113:048101, 2014.
  * (11) J. Schwabedal and A. Pikovsky.  Phase Description of Stochastic Oscillations.  Phys. Rev. Lett., 110:204102, 2013.
  * (12) P. J. Thomas and B. Lindner.  Asymptotic phase of stochastic oscillators.  Phys. Rev. Lett., 113:254101, 2014.
  * (13) G. Giacomin, C. Poquet, and A. Shapira.  Small noise and long time phase diffusion in stochastic limit cycle oscillators.  J. Diff. Equ., 264:1019, 2018.
  * (14) Z. Aminzare, Philip P. Holmes, and V. Srivastava.  On phase reduction and time period of noisy oscillators.  In 2019 IEEE 58th conference on decision and control (CDC), page 4717. IEEE, 2019.
  * (15) M. Engel and Ch. Kühn.  A random dynamical systems perspective on isochronicity for stochastic oscillations.  Commun. Math. Phys., 386:1603, 2021.
  * (16) A.Cao, B. Lindner, and P. J. Thomas.  A partial differential equation for the mean—first–return-time phase of planar stochastic oscillators.  SIAM J. App. Math., 80:422, 2020.
  * (17) C. W. Gardiner.  Handbook of Stochastic Methods.  Springer-Verlag, Berlin, 1985.
  * (18) J.M. Newby and M.A. Schwemmer.  Effects of Moderate Noise on a Limit Cycle Oscillator: Counterrotation and Bistability.  Phys. Rev. Lett., 112:114101, 2014.
  * (19) K. Holzhausen.  An analytic approach to the mean-first-return-time phase of planar isotropic stochastic oscillators.  Master’s thesis, Humboldt Universität zu Berlin, 2021.
  * (20) P. J. Thomas and B. Lindner.  Phase descriptions of a multidimensional Ornstein-Uhlenbeck process.  Phys. Rev. E, 99:062221, 2019.
  * (21) R. Galán, G. Ermentrout, and N. Urban.  Efficient estimation of phase-resetting curves in real neurons and its significance for neural-network modeling.  Phys. Rev. Lett., 94:158101, 2005.
  * (22) G. B. Ermentrout, R. F. Galán, and N. N. Urban.  Relating neural dynamics to neural coding.  Phys. Rev. Lett., 99:248103, 2007.
  * (23) S. Pu and P. J. Thomas.  Resolving molecular contributions of ion channel noise to interspike interval variability through stochastic shielding.  Biol. Cyb., 115:267, 2021.
  * (24) B. Ermentrout and D. Saunders.  Phase resetting and coupling of noisy neural oscillators.  J. Comput. Neurosci., 20:179, 2006.