Tr Re Im

# Multi-time correlations in the positive-P, Q, and doubled phase-space
representations

Piotr Deuar<EMAIL_ADDRESS>Institute of Physics, Polish Academy of
Sciences, Aleja Lotników 32/46, 02-668 Warsaw, Poland

###### Abstract

A number of physically intuitive results for the calculation of multi-time
correlations in phase-space representations of quantum mechanics are obtained.
They relate time-dependent stochastic samples to multi-time observables, and
rely on the presence of derivative-free operator identities. In particular,
expressions for time-ordered normal-ordered observables in the positive-P
distribution are derived which replace Heisenberg operators with the bare
time-dependent stochastic variables, confirming extension of earlier such
results for the Glauber-Sudarshan P. Analogous expressions are found for the
anti-normal-ordered case of the doubled phase-space Q representation, along
with conversion rules among doubled phase-space s-ordered representations. The
latter are then shown to be readily exploited to further calculate anti-normal
and mixed-ordered multi-time observables in the positive-P, Wigner, and
doubled-Wigner representations. Which mixed-order observables are amenable and
which are not is indicated, and explicit tallies are given up to 4th order.
Overall, the theory of quantum multi-time observables in phase-space
representations is extended, allowing non-perturbative treatment of many
cases. The accuracy, usability, and scalability of the results to large
systems is demonstrated using stochastic simulations of the unconventional
photon blockade system and a related Bose-Hubbard chain. In addition, a robust
but simple algorithm for integration of stochastic equations for phase-space
samples is provided.

##  1 Introduction

Phase-space representations of quantum mechanics such as the Wigner, P,
positive-P, Q and related approaches are a powerful tool for the study and
understanding of quantum mechanics [Gardiner(1991), Drummond and
Gardiner(1980), Blakie et al.(2008)Blakie, Bradley, Davis, Ballagh, and
Gardiner, Polkovnikov(2010), Lee(1995)]. Their use in recent times has been
directed particularly as a tool in quantum information science (see
[Ferrie(2011)] and [Sperling and Vogel(2020)] for recent reviews) and for the
simulation of large-scale quantum dynamics. Negativity of Wigner, P, and other
phase-space quasi-distributions is a major criterion for quantumness and
closely related to contextuality and nonlocality in quantum mechanics
[Ferrie(2011), Gehrke et al.(2012)Gehrke, Sperling, and Vogel]. The inability
to interpret the Glauber-Sudarshan P [Glauber(1963a), Sudarshan(1963)] and
Wigner distributions in terms of a classical probability density is the
fundamental benchmark for quantum light [Korbicz et al.(2005)Korbicz, Cirac,
Wehr, and Lewenstein, Sperling(2016)]. Quasi-probability representations arise
naturally when looking for hidden variable descriptions and an ontological
model of quantum theory [Ferrie(2011), Spekkens(2008), Lee(1995), Bondar et
al.(2013)Bondar, Cabrera, Zhdanov, and Rabitz]. Wigner and P-distribution
negativity are considered a resource for quantum computation [Albarelli et
al.(2018)Albarelli, Genoni, Paris, and Ferraro, Ferrie(2011)] and closely
related to magic state distillation and quantum advantage [Veitch et
al.(2012)Veitch, Ferrie, Gross, and Emerson]; the Glauber-Sudarshan P
distribution defines a hierarchy of nonclassicalities and nonclassicality
witnesses [Korbicz et al.(2005)Korbicz, Cirac, Wehr, and Lewenstein,
Vogel(2008)] and can be experimentally reconstructed [Agarwal and
Chaturvedi(1994), Kiesel et al.(2008)Kiesel, Vogel, Parigi, Zavatta, and
Bellini]. Moreover, P, Q and Wigner distributions and have been used for
simulations of models important for quantum information topics such as open
spin-qubit systems [Thapliyal et al.(2015)Thapliyal, Banerjee, Pathak, Omkar,
and Ravishankar] and boson sampling [Opanchuk et al.(2018)Opanchuk, Rosales-
Zárate, Reid, and Drummond, Drummond et al.(2016)Drummond, Opanchuk, Rosales-
Zárate, Reid, and Forrester, Opanchuk et al.(2019)Opanchuk, Rosales-Zárate,
Reid, and Drummond].

For simulation of quantum mechanics, the prime advantage of the phase-space
approach is that its computational cost typically grows only linearly with
system size even in interacting systems of many particles. Therefore it
provides a route to the non-perturbative treatment of the quantum dynamics of
large systems. It has been used for solving and simulating a multitude of
problems in various physical fields: e.g. quantum optics [Carter et
al.(1987)Carter, Drummond, Reid, and Shelby, Corney et al.(2008)Corney,
Heersink, Dong, Josse, Drummond, Leuchs, and Andersen, Drummond and
Walls(1980), Corney et al.(1997)Corney, Drummond, and Liebman, Lamprecht et
al.(2002)Lamprecht, Olsen, Drummond, and Ritsch, Olsen et al.(2004)Olsen,
Melo, Dechoum, and Khoury, Dechoum et al.(2004)Dechoum, Drummond, Chaturvedi,
and Reid, Drummond and Opanchuk(2020), Kiesewetter et al.(2014)Kiesewetter,
He, Drummond, and Reid], ultracold atoms [Steel et al.(1998)Steel, Olsen,
Plimak, Drummond, Tan, Collett, Walls, and Graham, Sinatra et
al.(2002)Sinatra, Lobo, and Castin, Carusotto and Castin(2001), Drummond and
Corney(1999), Deuar et al.(2009)Deuar, Sykes, Gangardt, Davis, Drummond, and
Kheruntsyan, Norrie et al.(2005)Norrie, Ballagh, and Gardiner, Deuar and
Drummond(2007), Kheruntsyan et al.(2012)Kheruntsyan, Jaskula, Deuar, Bonneau,
Partridge, Ruaudel, Lopes, Boiron, and Westbrook, Midgley et al.(2009)Midgley,
Wüster, Olsen, Davis, and Kheruntsyan, Świsłocki and Deuar(2016), Lewis-Swan
and Kheruntsyan(2014), Hush et al.(2013)Hush, Szigeti, Carvalho, and Hope,
Lewis-Swan and Kheruntsyan(2015), Wüster et al.(2017)Wüster, Corney, Rost, and
Deuar, Mathey et al.(2014)Mathey, Clark, and Mathey], fermionic systems
[Corney and Drummond(2004), Corney et al.(2006)Corney, Drummond, Heersink,
Josse, Leuchs, and Andersen, Aimi and Imada(2007a), Aimi and Imada(2007b)],
spin systems [Ng et al.(2013)Ng, Sørensen, and Deuar, Barry and
Drummond(2008), Ng and Sørensen(2011)], nuclear physics [Tsukiji et
al.(2016)Tsukiji, Iida, Kunihiro, Ohnishi, and Takahashi], dissipative systems
in condensed matter [Wouters and Savona(2009), Chiocchetta and
Carusotto(2014), Deuar et al.(2021)Deuar, Ferrier, Matuszewski, Orso, and
Szymańska], or cosmology [Opanchuk et al.(2013)Opanchuk, Polkinghorne, Fialko,
Brand, and Drummond]. The vast majority of such calculations so far have
considered only equal time correlations and observables.

Multi-time correlations are also important for answering many physical
questions which cannot necessarily be dealt with by monitoring the time
dependence of equal time observables [Gardiner(1991), Berg et al.(2009)Berg,
Plimak, Polkovnikov, Olsen, Fleischhauer, and Schleich, Polkovnikov(2010),
Lax(1968), Vogel(2008)]. For example, the determination of lifetimes in an
equilibrium or stationary state, response functions, out-of-time-order
correlations (OTOCs), or finding the time resolution required to observe a
transient phenomenon. However their calculation in the phase-space framework
has been restricted because few general results have been available. What is
known is largely restricted to time- and normal-ordered correlations in the
Glauber-Sudarshan P representation[Agarwal and Wolf(1970b), Gardiner(1991)],
time-symmetric ones in the truncated Wigner representation [Berg et
al.(2009)Berg, Plimak, Polkovnikov, Olsen, Fleischhauer, and Schleich, Plimak
and Olsen(2014)], or the corresponding linear response corrections for closed
systems. The Glauber-Sudarshan P representation results provide a particularly
intuitive framework, simply replacing Heisenberg operators
$\widehat{a}^{\dagger}(t)$ and $\widehat{a}(t)$ with time-dependent phase-
space variables in normal-ordered observables such as
$\widehat{a}^{\dagger}(t_{1})\widehat{a}(t_{2})$. They also apply to open
systems [Gardiner(1991), Agarwal(1969b)]. Such ordering corresponds to general
photon counting measurements [Glauber(1963a), Glauber(1963b), Kelley and
Kleiner(1964)]. Still, while scaling with system size is excellent, both above
approaches usually end up being only approximate because one is forced to omit
part of the full quantum mechanics in their numerical implementations
[Gardiner(1991)]. Hence, a broader application of full quantum phase-space
methods to multi-time observables would be advantageous. This is especially
so, since recent years have shown a lot of interest in such quantities, be it
in the study of nonclassicality in the time dimension [Atalaya et
al.(2018)Atalaya, Hacohen-Gourgy, Martin, Siddiqi, and Korotkov, Ringbauer et
al.(2018)Ringbauer, Costa, Goggin, White, and Alessandro, Moreva et
al.(2017)Moreva, Gramegna, Brida, Maccone, and Genovese, Krumm et
al.(2016)Krumm, Sperling, and Vogel], time crystals [Zhang et al.(2017)Zhang,
Hess, Kyprianidis, Becker, Lee, Smith, Pagano, Potirniche, Potter, Vishwanath,
Yao, and Monroe, Else et al.(2016)Else, Bauer, and Nayak, Wilczek(2012),
Syrwid et al.(2017)Syrwid, Zakrzewski, and Sacha], quantum technologies [Liew
and Savona(2010), Ji et al.(2015)Ji, Gladilin, and Wouters, Moczała-Dusanowska
et al.(2019)Moczała-Dusanowska, Dusanowski, Gerhardt, He, Reindl, Rastelli,
Trotta, Gregersen, Höfling, and Schneider] where time correlations and time
resolution are crucial, methods development [Plimak and Olsen(2014), Klobas et
al.(2020)Klobas, Vanicat, Garrahan, and Prosen], and the OTOCs that have
recently found to be important for the study of quantum scrambling [Swingle et
al.(2016)Swingle, Bentsen, Schleier-Smith, and Hayden, Bohrdt et
al.(2017)Bohrdt, Mendl, Endres, and Knap, Gärttner et al.(2018)Gärttner,
Hauke, and Rey], quantum chaos [Maldacena et al.(2016)Maldacena, Shenker, and
Stanford, Shen et al.(2017)Shen, Zhang, Fan, and Zhai] and many-body
localisation [Fan et al.(2017)Fan, Zhang, Shen, and Zhai, He and Lu(2017)].
Therefore, this paper sets out to extend the infrastructure available for
multi-time observables with phase-space methods.

The positive-P representation [Drummond and Gardiner(1980)] does not suffer
from the limitations of the Glauber-Sudarshan P (seen e.g. in Sec. 7.3) due to
its anchoring in doubled phase space. This allows it to incorporate the full
quantum mechanics for most Hamiltonians, including all two-body interactions.
Importantly, as we shall see, the same applies for its differently ordered
analogues like the doubled-Wigner and doubled-Q. While the positive-P is known
to be limited to short times for closed systems due to a noise amplification
instability [Gilchrist et al.(1997)Gilchrist, Gardiner, and Drummond, Deuar
and Drummond(2006)], this abates in dissipative systems [Gilchrist et
al.(1997)Gilchrist, Gardiner, and Drummond, Deuar et al.(2021)Deuar, Ferrier,
Matuszewski, Orso, and Szymańska]. Recent work has shown that a very broad
range of driven dissipative Bose-Hubbard models can be simulated with the
positive-P into and beyond the stationary state [Deuar et al.(2021)Deuar,
Ferrier, Matuszewski, Orso, and Szymańska]. This covers systems of current
interest such as micropillars, transmon qubits, and includes strong quantum
effects such two-photon interference in the unconventional photon blockade
[Liew and Savona(2010), Bamba et al.(2011)Bamba, Imamoğlu, Carusotto, and
Ciuti], making its extension quite timely.

Therefore, after basic background in Sec. 2, the first order of business in
this paper in Sec. 3 is to extend the long known Glauber-Sudarshan P results
on time- and normal-ordered multi-time correlations to the positive-P
representation. This is sort of an obvious extension, but was missing from the
literature, and allows complete quantum calculations. It also lets us explain
the ins and outs and set the stage for the less obvious and broader extensions
that follow:

The first one of those in Sec. 4.1 considers general representations and
reformulates the case of Q representations in an operational form that can be
extended to doubled phase-space. This gives stochastic access to anti-normal
ordered observables such as $\widehat{a}(t_{1})\widehat{a}^{\dagger}(t_{2})$.
It is then shown in Sec. 5 how the Gaussian convolution relationship between
Wigner, P, and Q distributions allows one to easily evaluate anti-normal
ordered correlations in Wigner and P. This is extended also to a wide range of
mixed ordered observables such as
$\widehat{a}^{\dagger}(t_{1})\widehat{a}(t_{2})\widehat{a}^{\dagger}(t_{3})\widehat{a}(t_{4})$
in Sec. 5.4, and further extended to cover doubled phase-space representations
in Sec. 6. In Sec. 7, the use and accuracy of this approach is demonstrated on
the numerical example of the unconventional photon blockade, a system that may
have application to the creation of single-photon sources. Finally, Appendix A
gives details of a robust algorithm for integrating the resulting phase-space
stochastic equations, a matter that has been somewhat neglected in the
literature to date.

##  2 Background

###  2.1 Positive-P representation

Consider an $M$-mode system (modes $1,\dots,j,\dots,M$) with bosonic
annihilation operators $\widehat{a}_{j}$ for the $j$th mode. The density
operator of the system is written in terms of bosonic coherent states of
complex amplitude $\alpha_{j}$ defined relative to the vacuum:

$|\alpha_{j}\rangle_{j}=e^{-|\alpha_{j}|^{2}/2}\,e^{\alpha_{j}\widehat{a}^{\,\dagger}_{j}}|\rm
vac\rangle,$ (1)

using [Drummond and Gardiner(1980)]:

$\displaystyle\widehat{\rho}$ $\displaystyle=$ $\displaystyle\int
d^{4M}\bm{\lambda}\ P_{+}(\bm{\lambda})\,\widehat{\Lambda}(\bm{\lambda})$ (2)
$\displaystyle\widehat{\Lambda}(\bm{\lambda})$ $\displaystyle=$
$\displaystyle\bigotimes_{j}\frac{|\alpha_{j}\rangle_{j}\langle\beta^{*}_{j}|_{j}}{\langle\beta_{j}^{*}|_{j}|\alpha_{j}\rangle_{j}}.$
(3)

The kernel $\widehat{\Lambda}$ involves “ket” $|\alpha_{j}\rangle_{j}$ and
“bra” $\langle\beta^{*}_{j}|_{j}$ states described by separate and independent
variables. Let us define the container variable

$\bm{\lambda}=\left\\{\alpha_{1},\dots,\alpha_{M},\beta_{1},\dots,\beta_{M}\right\\}=\\{\lambda_{\mu}\\}$
(4)

as a shorthand to hold the full configuration information. Individual complex
variables $\lambda_{\mu}$ can be indexed by $\mu=1,\dots,2M$. The _positive-P_
distribution $P_{+}$ is guaranteed real and non-negative [Drummond and
Gardiner(1980)]. This is essential for the utility of the method, which lies
in its ability to represent full quantum mechanics using stochastic
trajectories of the _samples_ $\bm{\lambda}$ of the distribution $P_{+}$.
Notably, the size of these samples scales merely linearly with $M$, allowing
for the simulation of the full quantum dynamics of systems with even millions
of modes [Deuar and Drummond(2007), Kheruntsyan et al.(2012)Kheruntsyan,
Jaskula, Deuar, Bonneau, Partridge, Ruaudel, Lopes, Boiron, and Westbrook].

The distribution (2.1) can be compared to the simpler and more widely known
_Glauber-Sudarshan P representation_ [Glauber(1963a), Sudarshan(1963)], which
uses a single set of coherent state amplitudes and ties the “bra” and ”ket”
amplitudes to be equal. It is written

$\widehat{\rho}=\int d^{2M}{\mathbf{\bm{\alpha}}}\
P(\bm{\alpha})\bigotimes_{j}|\alpha_{j}\rangle_{j}\langle\alpha_{j}|_{j},$ (5)

where the multimode coherent amplitude is
$\bm{\alpha}=\\{\alpha_{1},\dots,\alpha_{M}\\}$. This representation gives a
well behaved representation of only a subset of possible quantum states
[Cahill and Glauber(1969b), Sperling(2016)], though a very useful one that
includes Gaussian density operators [Agarwal and Wolf(1970a)].

For full computational utility of the positive-P method, all quantum
mechanical actions should be able to be written in terms of the samples
$\bm{\lambda}$. Central to this are the following identities:

$\displaystyle\widehat{a}_{j}\widehat{\Lambda}$ $\displaystyle=$
$\displaystyle\alpha_{j}\widehat{\Lambda},$ (6)
$\displaystyle\widehat{a}^{\dagger}_{j}\widehat{\Lambda}$ $\displaystyle=$
$\displaystyle\left[\beta_{j}+\frac{\partial}{\partial\alpha_{j}}\right]\widehat{\Lambda},$
(7) $\displaystyle\widehat{\Lambda}\widehat{a}_{j}$ $\displaystyle=$
$\displaystyle\left[\alpha_{j}+\frac{\partial}{\partial\beta_{j}}\right]\widehat{\Lambda},$
(8) $\displaystyle\widehat{\Lambda}\widehat{a}^{\dagger}_{j}$ $\displaystyle=$
$\displaystyle\beta_{j}\widehat{\Lambda}.$ (9)

A general master equation for the density matrix written in Lindblad form with
reservoir operators $\widehat{R}_{n}$ is

$\frac{\partial\widehat{\rho}}{\partial
t}=-\frac{i}{\hbar}\left[\widehat{H},\widehat{\rho}\right]+\sum_{n}\left[2\widehat{R}_{n}\widehat{\rho}\widehat{R}^{\dagger}_{n}-\widehat{R}^{\dagger}_{n}\widehat{R}_{n}\widehat{\rho}-\widehat{\rho}\widehat{R}^{\dagger}_{n}\widehat{R}_{n}\right].$
(10)

This can be transformed by standard methods [Gardiner(1991)], using identities
(2.1), to a partial differential equation (PDE) for $P_{+}$ of the following
general form:

$\frac{\partial P_{+}}{\partial t}=\left\\{\sum_{n=1}^{n_{\rm
max}}\sum_{\mu_{1},\dots,\mu_{k}}^{2M}\left(\prod_{k=1}^{n}\frac{\partial}{\partial\lambda_{\mu_{k}}}\right)F_{\mu_{1},\dots,\mu_{n}}(\bm{\lambda})\right\\}P_{+},$
(11)

with maximum differential order $n_{\rm max}$ and coefficients $F$ that depend
on the details of (10). Most cases have $n_{\rm max}\leq 2$, in which case
this is a Fokker-Planck equation (FPE):

$\frac{\partial P_{+}}{\partial
t}=\left\\{\sum_{\mu}\frac{\partial}{\partial\lambda_{\mu}}\left[-A_{\mu}(\bm{\lambda})\right]+\sum_{\mu\nu}\frac{\partial^{2}}{\partial\lambda_{\mu}\partial\lambda_{\nu}}\frac{D_{\mu\nu}(\bm{\lambda})}{2}\right\\}P_{+}.$
(12)

An exception occurs if irreducible higher order partial derivatives appear,
e.g. due to explicit three-body interactions in the Hamiltonian. Most first
principles models use only two-particle interactions, though.

Standard methods convert an FPE like (12) to the following Ito stochastic
equations of the samples [Gardiner(1991)]:

$\frac{d\lambda_{\mu}}{dt}=A_{\mu}(\bm{\lambda})+\sum_{\sigma}B_{\mu\sigma}(\bm{\lambda})\xi_{\sigma}(t)$
(13)

where

$D_{\mu\nu}=\sum_{\sigma}B_{\mu\sigma}B_{\nu\sigma}$ (14)

i.e. $D=BB^{T}$ in matrix notation, and $\xi_{\sigma}(t)$ are independent real
white noises with variances

$\langle\xi_{\mu}(t)\xi_{\nu}(t^{\prime})\rangle_{\rm
stoch}=\delta(t-t^{\prime})\delta_{\mu\nu}.$ (15)

The notation $\langle\cdot\rangle_{\rm stoch}$ indicates a stochastic average
over the samples. The equations (13) give us quantum mechanical evolution in
terms of the samples, which becomes ever more exact as the number of samples
grows. Quantum expectation values at a given time are evaluated via

$\langle\widehat{a}^{\dagger}_{j_{1}}\cdots\widehat{a}^{\dagger}_{j_{{\mathcal{N}}}}\widehat{a}_{k_{1}}\cdots\widehat{a}_{k_{{\mathcal{M}}}}\rangle=\langle\beta_{j_{1}}\cdots\beta_{j_{{\mathcal{N}}}}\alpha_{k_{1}}\cdots\alpha_{k_{{\mathcal{M}}}}\rangle_{\rm
stoch}.$ (16)

The equivalence can be written more explicitly in terms of ${\mathcal{S}}$
individual samples $\bm{\lambda}^{(u)}$ labelled by $u=1,\dots,{\mathcal{S}}$
as:

$\displaystyle\langle\widehat{a}^{\dagger}_{j_{1}}\cdots\widehat{a}^{\dagger}_{j_{{\mathcal{N}}}}\widehat{a}_{k_{1}}\cdots\widehat{a}_{k_{{\mathcal{M}}}}\rangle$
$\displaystyle\qquad=\lim_{{\mathcal{S}}\to\infty}\frac{1}{{\mathcal{S}}}\sum_{u}\beta^{(u)}_{j_{1}}\cdots\beta^{(u)}_{j_{{\mathcal{N}}}}\alpha^{(u)}_{k_{1}}\cdots\alpha^{(u)}_{k_{{\mathcal{M}}}}.$

The mode labels in the above can be in any combination, provided all
annihilation operators are to the right of all creation operators (which is
called “normal ordering”). Any single-time operator can be expressed as a sum
of normally ordered terms like (16).

###  2.2 Multi-time averages in open systems

The situation is much more complicated when the operators in the expectation
value are not evaluated all at the same time. The root of the difficulty is
that multi-time commutation relations usually depend non-trivially on the full
system dynamics, and a reduction of arbitrary multi-time operators to a
normal-ordered form is not generally possible.

To describe what is meant at the operator level it is first helpful to
introduce the two-sided evolution operator $\breve{V}(t_{1},t_{2})$ such that
evolution by the master equation (10) can be summarised as

$\widehat{\rho}(t_{2})=\breve{V}(t_{2},t_{1})\widehat{\rho}(t_{1}).$ (18)

Notably, the evolution operator has the semigroup property [Spohn(1980)]
$\breve{V}(t_{3},t_{2})\breve{V}(t_{2},t_{1})=\breve{V}(t_{3},t_{1})$. We will
use the convention that two-sided operators, indicated by a breve $\breve{}$ ,
act on everything to their right. Hence

$\breve{V}AB=\breve{V}\left\\{AB\right\\}\neq\left\\{\breve{V}A\right\\}B.$
(19)

Let us also define the _time-ordered form_ via

$\langle\widehat{A}_{1}(t_{1})\widehat{A}_{2}(t_{2})\cdots\widehat{A}_{{{\mathcal{N}}}}(t_{{\mathcal{N}}})\widehat{B}_{1}(s_{1})\widehat{B}_{2}(s_{2})\cdots\widehat{B}_{{\mathcal{M}}}(s_{{\mathcal{M}}})\rangle$
(20)

where the $\widehat{A}_{j}(t)$ and $\widehat{B}_{j}(t)$ are Heisenberg picture
operators, and the times obey

$\displaystyle t_{1}\leq t_{2}\leq$ $\displaystyle\dots$ $\displaystyle\leq
t_{{\mathcal{N}}}$ (21) $\displaystyle s_{1}\geq s_{2}\geq$
$\displaystyle\dots$ $\displaystyle\geq s_{{\mathcal{M}}}.$ (22)

There are ${\mathcal{N}}$ operators $\widehat{A}$ with times labelled
$t_{1},\dots,t_{{\mathcal{N}}}$ increasing to the right, inward, and
${\mathcal{M}}$ operators $\widehat{B}$ with times labelled
$s_{1},\dots,s_{{\mathcal{M}}}$ increasing to the left, also inward. The
location of the inner “meeting point” is arbitrary, and the number of
$\widehat{A}$ or $\widehat{B}$ operators can also be zero. No particular
constraints are imposed on the $\widehat{A}$ and $\widehat{B}$ operators,
except that they should be single time quantities. These operators refer to
measurements made at the respective times $\tau_{r}$, when the density
operator was $\widehat{\rho}(\tau_{r})$. Between measurements the state
evolves according to the master equation (10) (i.e. (18)).

It has been shown [Gardiner(1991)] that multi-time correlations that
correspond to sequences of measurements can always be written in the above
form (20). Therefore this time ordering is not an arbitrary one, and not
particularly restrictive in itself. It has also been shown that a general
time-ordered correlation function (20) obeying (2.2) can be written in a form
that uses the $\breve{V}$. To do so it is necessary to be careful about
operator ordering. Following [Gardiner(1991)], let us order all the times
$t_{p}$ and $s_{q}$ in the correlation in sequence from earliest to latest.
Let us then rename them $\tau_{r}$ so that

$\tau_{1}\leq\tau_{2}\leq\dots\leq\tau_{{\mathcal{R}}-1}\leq\tau_{{\mathcal{R}}}$
(23)

with ${\mathcal{R}}={\mathcal{N}}+{\mathcal{M}}$. We also define the
corresponding two-sided operators $\breve{F}_{r}$

$\breve{F}_{r}\widehat{\rho}=\left\\{\begin{array}[]{l@{\quad}l}\widehat{\rho}\widehat{A}_{p}&\text{if$\tau_r=t_p$}\\\
\widehat{B}_{q}\widehat{\rho}&\text{if$\tau_r=s_q$}\end{array}\right.$ (24)

Then,

$\displaystyle\langle\widehat{A}_{1}(t_{1})\widehat{A}_{2}(t_{2})\cdots\widehat{A}_{{\mathcal{N}}}(t_{{\mathcal{N}}})\widehat{B}_{1}(s_{1})\widehat{B}_{2}(s_{2})\cdots\widehat{B}_{{\mathcal{M}}}(s_{{\mathcal{M}}})\rangle$
(25) $\displaystyle={\rm
Tr}\Big{[}\breve{F}_{{\mathcal{R}}}\breve{V}(\tau_{{\mathcal{R}}},\tau_{{\mathcal{R}}-1})\breve{F}_{{\mathcal{R}}-1}\breve{V}(\tau_{{\mathcal{R}}-1},\tau_{{\mathcal{R}}-2})\cdots$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\cdots\breve{F}_{2}\breve{V}(\tau_{2},\tau_{1})\breve{F}_{1}\widehat{\rho}(\tau_{1}))\Big{]}.\qquad$

###  2.3 Multi-time correlations in the Glauber-Sudarshan P representation

Time-ordered correlations for the Glauber-Sudarshan P representation and a
range of other single-phase space representations like the Husimi Q were first
studied by Agarwal and Wolf in a Greens function framework [Agarwal(1969a)]
extended to the case of open systems [Agarwal(1969b)] and for Hamiltonian
systems in a formal integral form [Agarwal and Wolf(1970b)]. Later, Gardiner
[Gardiner(1991)] used a different more operational approach to derive
equivalent expressions for normally and time-ordered operator averages in the
Glauber-Sudarshan P representation. These have an intuitive form similar to
the single-time stochastic expression (16), as follows:

$\displaystyle\langle\widehat{a}^{\dagger}_{p_{1}}(t_{1})\cdots\widehat{a}^{\dagger}_{p_{{\mathcal{N}}}}(t_{{\mathcal{N}}})\widehat{a}_{q_{1}}(s_{1})\cdots\widehat{a}_{q_{{\mathcal{M}}}}(s_{{\mathcal{M}}})\rangle=\qquad\qquad$
(26)
$\displaystyle\qquad\qquad\langle\alpha^{*}_{p_{1}}(t_{1})\cdots\alpha^{*}_{p_{{\mathcal{N}}}}(t_{{\mathcal{N}}})\alpha_{q_{1}}(s_{1})\cdots\alpha_{q_{{\mathcal{M}}}}(s_{{\mathcal{M}}})\rangle_{\rm
stoch}$

provided the times are ordered according to (2.2). The requirement that the
operators in (26) be normally ordered is an additional constraint on top of
time-ordering (2.2), but one that leads to all operators sorted as they occur
in photo-counting theory [Kelley and Kleiner(1964)]. It covers a very large
subset of the potentially physically interesting correlations. The Gardiner
approach is more amenable to extension to doubled phase space and will be used
in what follows.

This ordering can be contrasted with the _time-symmetric_ ordering for which
straightforward truncated Wigner correspondences for closed system evolution
were found in [Berg et al.(2009)Berg, Plimak, Polkovnikov, Olsen,
Fleischhauer, and Schleich, Plimak and Olsen(2014), Polkovnikov(2010)].
Examples of time-symmetric ordered quantities are
$\tfrac{1}{2}\left[\widehat{a}^{\dagger}(t_{2})\widehat{a}(t_{1})+\widehat{a}(t_{1})\widehat{a}^{\dagger}(t_{2})\right]$
and
$\tfrac{1}{4}\left[\widehat{a}(t_{1})\widehat{a}(t_{2})\widehat{a}^{\dagger}(t_{3})+\widehat{a}(t_{2})\widehat{a}^{\dagger}(t_{3})\widehat{a}(t_{1})+\widehat{a}(t_{1})\widehat{a}^{\dagger}(t_{3})\widehat{a}(t_{2})\right.$
$\left.+\widehat{a}^{\dagger}(t_{3})\widehat{a}(t_{2})\widehat{a}(t_{1})\right]$.

##  3 Time ordered moments in the positive-P representation

###  3.1 Normally ordered observables

To derive an expression like (26) for the _positive_ -P representation, we
will follow Gardiner’s approach [Gardiner(1991)] that was previously used for
the Glauber-Sudarshan P. Smaller steps will, however, be taken here to draw
attention to a few subtleties that will be necessary later.

####  3.1.1 First order correlation function

First consider the correlation

$G^{(1)}(t^{\prime},t)=\langle\widehat{a}^{\dagger}_{j}(t^{\prime})\widehat{a}_{k}(t)\rangle.$
(27)

Comparing to (25), we can identify
$\widehat{A}_{1}=\widehat{a}^{\dagger}_{j}$,
$\widehat{B}_{1}=\widehat{a}_{k}$, and the times $t_{1}=t^{\prime}$,
$s_{1}=t$. For this low order correlation, the time ordering (2.2) sets no
additional conditions. There are two possibilities for $\tau_{1}$, depending
on whether time $t$ or $t^{\prime}$ is later. Consider first $t^{\prime}\geq
t$, so that $\tau_{1}=t$, $\tau_{2}=t^{\prime}$. Using (25) and (24) we have
that

$\displaystyle\langle\widehat{a}^{\dagger}_{j}(t^{\prime})\widehat{a}_{k}(t)\rangle$
$\displaystyle=$
$\displaystyle\TR\left[{\breve{V}(t^{\prime},t)\left\\{\widehat{a}_{k}(t)\widehat{\rho}(t)\right\\}\widehat{a}^{\dagger}_{j}(t^{\prime})}\right]$
(28) $\displaystyle=$
$\displaystyle\TR\left[{\widehat{a}^{\dagger}_{j}(t^{\prime})\breve{V}(t^{\prime},t)\left\\{\widehat{a}_{k}(t)\widehat{\rho}(t)\right\\}}\right]$
(29)

The 2nd line follows from the cyclic property of traces. The $\\{\cdot\\}$ is
kept for now for clarity. We can see that the evolution operator from $t$ to
$t^{\prime}$ acts to the right on all the quantities at time $t$, while the
later-time operator $\widehat{a}^{\dagger}_{j}(t^{\prime})$ acts only on the
evolved quantities to its right. This makes intuitive physical sense. In the
second case of $t^{\prime}<t$, we have $\tau_{1}=t^{\prime}$, $\tau_{2}=t$,
and get that

$\langle\widehat{a}^{\dagger}_{j}(t^{\prime})\widehat{a}_{k}(t)\rangle=\TR\left[{\widehat{a}_{k}(t)\breve{V}(t,t^{\prime})\left\\{\widehat{\rho}(t^{\prime})\widehat{a}^{\dagger}_{j}(t^{\prime})\right\\}}\right].$
(30)

This again has the intuitive form of the operator $\breve{V}$ acting to the
right on all the earlier-time quantities.

Take the first case with $t^{\prime}\geq t$. Expressing the density matrix in
(29) in the positive-P representation (2),

$\displaystyle G^{(1)}$ $\displaystyle=$
$\displaystyle\TR\left[{\widehat{a}^{\dagger}_{j}(t^{\prime})\breve{V}(t^{\prime},t)\left\\{\\!\int
d^{4M}\bm{\lambda}\,P_{+}(\bm{\lambda},t)\,\widehat{a}_{k}(t)\widehat{\Lambda}(\bm{\lambda})\right\\}}\right]$
(31)
$\displaystyle=\TR\left[{\widehat{a}^{\dagger}_{j}(t^{\prime})\breve{V}(t^{\prime},t)\left\\{\\!\int
d^{4M}\bm{\lambda}\,\alpha_{k}P_{+}(\bm{\lambda},t)\,\widehat{\Lambda}(\bm{\lambda})\right\\}}\right].\qquad$

The 2nd line follows from application of (6). We cannot do the same for
$\widehat{a}^{\dagger}_{j}(t^{\prime})$ yet, because the kernel
$\widehat{\Lambda}$ finds itself inside the prior action of the $\breve{V}$
operator.

To deal with this, consider now the action of the evolution operator
$\breve{V}$ on distributions $P_{+}$. i.e. the action of the PDE (11). If we
define the conditional distribution
${\mathcal{P}}(\bm{\lambda},t^{\prime}|{\underline{\bm{\lambda}}},t)$ as the
solution of this PDE at time $t^{\prime}\geq t$ starting from the initial
condition $\delta^{4M}(\bm{\lambda}-{\underline{\bm{\lambda}}})$, i.e. the
“propagator”, then it can be used to formally write

$\displaystyle\breve{V}(t^{\prime},t)\left\\{\\!\int
d^{4M}\bm{\lambda}\,P_{+}(\bm{\lambda},t)\widehat{\Lambda}(\bm{\lambda})\right\\}$
(32) $\displaystyle=$ $\displaystyle\breve{V}(t^{\prime},t)\left\\{\\!\int
d^{4M}\bm{\lambda}\int
d^{4M}{\underline{\bm{\lambda}}}\,\delta^{4M}(\bm{\lambda}-{\underline{\bm{\lambda}}})\,P_{+}({\underline{\bm{\lambda}}},t)\widehat{\Lambda}(\bm{\lambda})\right\\}$
$\displaystyle=$ $\displaystyle\int d^{4M}\bm{\lambda}\int
d^{4M}{\underline{\bm{\lambda}}}\,{\mathcal{P}}(\bm{\lambda},t^{\prime}|{\underline{\bm{\lambda}}},t)P_{+}({\underline{\bm{\lambda}}},t)\widehat{\Lambda}(\bm{\lambda}).$

This now contains no more two-sided operators. Through this convolution,
$P_{+}$ is expressed in ${\underline{\bm{\lambda}}}$ variables which accompany
the earliest time $t$ to aid for later interpretation as part of a joint
probability, while the kernel $\widehat{\Lambda}$ is expressed in the
variables $\bm{\lambda}$ that accompany later times, ready for application of
the next operator identity.

Notice that there are no particularly stringent assumptions about $P_{+}$ for
(32) itself to apply. For example, (32) applies equally well if one replaces
$P_{+}$ with some complex distribution function ${\widetilde{P}}$. This point
will soon be useful. However, there _was_ an assumption that the propagator
${\mathcal{P}}$ is well behaved. This is certainly true if the PDE was of a
Fokker-Planck form (12), and therefore is always justified if we have an exact
mapping of the master equation (10) to stochastic equations (13). However, in
some other cases of third/higher order terms in the PDE, it might not. We will
not be concerned with such cases here.

Now in (31), $\breve{V}$ is acting on a distribution
${\widetilde{P}}(\bm{\lambda},t)=\alpha_{k}P_{+}(\bm{\lambda},t)$. Using (32)
we get

$\displaystyle G^{(1)}=\TR\left[{\widehat{a}^{\dagger}_{j}(t^{\prime})\\!\int
d^{4M}\bm{\lambda}\,\int
d^{4M}{\underline{\bm{\lambda}}}{\mathcal{P}}(\bm{\lambda},t^{\prime}|{\underline{\bm{\lambda}}},t){\widetilde{P}}({\underline{\bm{\lambda}}},t)\,\widehat{\Lambda}(\bm{\lambda})}\right]$
$\displaystyle=\TR\left[{\widehat{a}^{\dagger}_{j}(t^{\prime})\\!\int
d^{4M}\bm{\lambda}\,\int
d^{4M}{\underline{\bm{\lambda}}}\,{\underline{\alpha}}_{k}{\mathcal{P}}(\bm{\lambda},t^{\prime}|{\underline{\bm{\lambda}}},t)P_{+}({\underline{\bm{\lambda}}},t)\,\widehat{\Lambda}(\bm{\lambda})}\right].$

The two-way operator acting on the right that was $\breve{V}$, has now been
gotten rid of, by virtue of being incorporated in the propagator
${\mathcal{P}}$. Therefore, the remaining operator $\widehat{a}^{\dagger}_{j}$
can now be shifted to the right due to the cyclic property of the trace, and
then processed via (9) as so:

$\displaystyle G^{(1)}\\!=\\!\int d^{4M}\bm{\lambda}\int
d^{4M}{\underline{\bm{\lambda}}}{\mathcal{P}}(\bm{\lambda},t^{\prime}|{\underline{\bm{\lambda}}},t){\widetilde{P}}({\underline{\bm{\lambda}}},t)\TR\left[{\widehat{\Lambda}(\bm{\lambda})\widehat{a}^{\dagger}_{j}(t^{\prime})}\right]$
$\displaystyle=\int d^{4M}\bm{\lambda}\,\beta_{j}\int
d^{4M}{\underline{\bm{\lambda}}}\,{\underline{\alpha}}_{k}{\mathcal{P}}(\bm{\lambda},t^{\prime}|{\underline{\bm{\lambda}}},t)P_{+}({\underline{\bm{\lambda}}},t)\TR\left[{\widehat{\Lambda}(\bm{\lambda})}\right]$
$\displaystyle=\int d^{4M}\bm{\lambda}\,\beta_{j}\int
d^{4M}{\underline{\bm{\lambda}}}\,{\underline{\alpha}}_{k}{\mathcal{P}}(\bm{\lambda},t^{\prime}|{\underline{\bm{\lambda}}},t)P_{+}({\underline{\bm{\lambda}}},t).$
(34)

The last line follows from $\TR\left[{\widehat{\Lambda}}\right]=1$, which is
pre-set by the definition (3).

The quantity ${\mathcal{P}}P_{+}$ is the just the joint probability

$P(\bm{\lambda},t^{\prime};{\underline{\bm{\lambda}}},t)={\mathcal{P}}(\bm{\lambda},t^{\prime}|{\underline{\bm{\lambda}}},t)P_{+}({\underline{\bm{\lambda}}},t)$
(35)

of having configuration ${\underline{\bm{\lambda}}}$ at time $t$ and
configuration $\bm{\lambda}$ at time $t^{\prime}$. (Provided ${\mathcal{P}}$
is well behaved, positive, real, as mentioned before, which is the case for
any model fully described by an FPE (12)). Therefore,

$G^{(1)}=\int d^{4M}\bm{\lambda}\,\int d^{4M}{\underline{\bm{\lambda}}}\
\beta_{j}{\underline{\alpha}}_{k}\
P(\bm{\lambda},t^{\prime};{\underline{\bm{\lambda}}},t).$ (36)

At this stage we can identify probability with stochastic realisations. The
noises $\xi_{\mu}(t)$ introduced during time evolution are independent from
each other, independent at each time step, independent for each sample’s
trajectory. They are also independent of any other random variables used to
produce the initial ensemble
$\\{\bm{\lambda}^{(1)},\dots,\bm{\lambda}^{({\mathcal{S}})}\\}(t)$ that
samples $P_{+}(\bm{\lambda},t)$. Therefore, the evolved configuration
$\bm{\lambda}^{(u)}(t^{\prime})$ at time $t^{\prime}$ depends only on mutually
independent random variables that consist of $\bm{\lambda}^{(u)}(t)$ and the
noise history of the $u$th trajectory. As a result, the combination of initial
configuration $\bm{\lambda}^{(u)}(t)$ and the evolved configuration
$\bm{\lambda}^{(u)}(t^{\prime})$, together form an unbiased sample of the
joint distribution $P(\bm{\lambda},t^{\prime};{\underline{\bm{\lambda}}},t)$.
We arrive then at the final result that

$G^{(1)}(t^{\prime},t)=\langle\widehat{a}^{\dagger}_{j}(t^{\prime})\widehat{a}_{k}(t)\rangle=\langle\beta_{j}(t^{\prime})\alpha_{k}(t)\rangle_{\rm
stoch}.$ (37)

The procedure for the case $t^{\prime}<t$ gives a result analogous to (36):

$G^{(1)}=\int d^{4M}\bm{\lambda}\,\int d^{4M}{\underline{\bm{\lambda}}}\
{\underline{\beta}}_{j}\alpha_{k}\
P(\bm{\lambda},t;{\underline{\bm{\lambda}}},t^{\prime}),$ (38)

which once again leads to (37).

It is especially important to note that the quantity to be averaged comes from
the _time evolution of individual sample trajectories_. Explicitly:

$\langle\widehat{a}^{\dagger}_{j}(t^{\prime})\widehat{a}_{k}(t)\rangle=\lim_{{\mathcal{S}}\to\infty}\frac{1}{{\mathcal{S}}}\sum_{u=1}^{{\mathcal{S}}}\beta^{(u)}_{j}(t^{\prime})\alpha^{(u)}_{k}(t).$
(39)

This allows for very efficient calculations.

####  3.1.2 Higher order correlations

Other time and normal ordered correlations follow a similar pattern. For
example, when $t^{\prime\prime}\geq t^{\prime}\geq t$,

$\displaystyle\langle\widehat{a}_{j}(t^{\prime\prime})\widehat{a}_{k}(t^{\prime})\widehat{a}_{l}(t)\rangle$
$\displaystyle\qquad=\TR\left[{\widehat{a}_{j}(t^{\prime\prime})\breve{V}(t^{\prime\prime},t^{\prime})\left\\{\widehat{a}_{k}(t^{\prime})\breve{V}(t^{\prime},t)\left\\{\widehat{a}_{l}(t)\widehat{\rho}(t)\right\\}\right\\}}\right].$

In this case, following a similar procedure to before, one can act alternately
with (6) on the kernel to extract a factor of $\alpha$, and (32) to convert
the evolution operators to propagators. One finds

$\displaystyle\langle\widehat{a}_{j}(t^{\prime\prime})\widehat{a}_{k}(t^{\prime})\widehat{a}_{l}(t)\rangle$
$\displaystyle=$ $\displaystyle\int
d^{4M}\bm{\lambda}d^{4M}{\underline{\bm{\lambda}}}d^{4M}{\underline{{\underline{\bm{\lambda}}}}}\
\,\alpha_{j}{\underline{\alpha}}_{k}{\underline{{\underline{\alpha}}}}_{l}$
(41)
$\displaystyle\times{\mathcal{P}}(\bm{\lambda},t^{\prime\prime}|{\underline{\bm{\lambda}}},t^{\prime}){\mathcal{P}}({\underline{\bm{\lambda}}},t^{\prime}|{\underline{{\underline{\bm{\lambda}}}}},t)P_{+}({\underline{{\underline{\bm{\lambda}}}}},t).$

Since the times are ordered $t^{\prime\prime}\geq t^{\prime}\geq t$, and the
conditional probabilities follow from the evolution of the FPE,
${\underline{{\underline{\bm{\lambda}}}}}$ are parent variables of the
${\underline{\bm{\lambda}}}$ and so on, and the product of conditional
probabilities is just the joint probability. Hence (41) is

$\int
d^{4M}\bm{\lambda}d^{4M}{\underline{\bm{\lambda}}}d^{4M}{\underline{{\underline{\bm{\lambda}}}}}\
\,\alpha_{j}{\underline{\alpha}}_{k}{\underline{{\underline{\alpha}}}}_{l}\
P(\bm{\lambda},t^{\prime\prime};{\underline{\bm{\lambda}}},t^{\prime};{\underline{{\underline{\bm{\lambda}}}}},t),$
(42)

and

$\langle\widehat{a}_{j}(t^{\prime\prime})\widehat{a}_{k}(t^{\prime})\widehat{a}_{l}(t)\rangle=\langle\alpha_{j}(t^{\prime\prime})\alpha_{k}(t^{\prime})\alpha_{l}(t)\rangle_{\rm
stoch}.$ (43)

Working similarly, using just the 1st and 4th identities in (2.1), one readily
but somewhat cumbersomely finds that the stochastic estimator for the general
time-and-normal-ordered correlation function is

$\displaystyle\langle\widehat{a}^{\dagger}_{p_{1}}(t_{1})\cdots\widehat{a}^{\dagger}_{p_{{\mathcal{N}}}}(t_{{\mathcal{N}}})\widehat{a}_{q_{1}}(s_{1})\cdots\widehat{a}_{q_{{\mathcal{M}}}}(s_{{\mathcal{M}}})\rangle$
(44)
$\displaystyle=\langle\beta_{p_{1}}(t_{1})\cdots\beta_{p_{{\mathcal{N}}}}(t_{{\mathcal{N}}})\alpha_{q_{1}}(s_{1})\cdots\alpha_{q_{{\mathcal{M}}}}(s_{{\mathcal{M}}})\rangle_{\rm
stoch}.\qquad$

Here, of course (2.2) must hold, and the stochastic averaging is over products
constructed using values from the evolution of a single sample, as in (39). It
confirms the suspicion and intuition that the behaviour of the positive-P
representation in this regard should be similar to the earlier expression
(26), for the Glauber-Sudarshan P.

###  3.2 Other ordering in the positive-P

Now to see the limitations of this scheme, consider the anti-normally (but
time-ordered) ordered correlation

${\mathcal{A}}=\langle\widehat{a}_{j}(t^{\prime})\widehat{a}^{\dagger}_{k}(t)\rangle$
(45)

with $t^{\prime}>t$. The first point to make is that we cannot rearrange this
to a normal-ordered form like
$\langle\widehat{a}^{\dagger}_{k}(t)\widehat{a}_{j}(t^{\prime})+\delta_{jk}\rangle$
and then use (44) because generally
$\left[\widehat{a}_{j}(t^{\prime}),\widehat{a}^{\dagger}_{k}(t)\right]\neq\delta_{jk}$
when $t\neq t^{\prime}$. Instead it is some time-dependent operator. Now
applying (25), (45) can be written

${\mathcal{A}}=\langle\widehat{a}_{j}(t^{\prime})\widehat{a}^{\dagger}_{k}(t)\rangle=\TR\left[{\widehat{a}_{j}(t^{\prime})\breve{V}(t^{\prime},t)\left\\{\widehat{a}^{\dagger}_{k}(t)\widehat{\rho}(t)\right\\}}\right].$
(46)

Upon expansion, we will need to act on $\widehat{\Lambda}$ using the 2nd
identity in (2.1), (7), to convert $\widehat{a}^{\dagger}_{k}$ to variable
form. Thus

$\displaystyle{\mathcal{A}}$ $\displaystyle=$ $\displaystyle{\rm
Tr}\Bigg{[}\widehat{a}_{j}(t^{\prime})\breve{V}(t^{\prime},t)\Bigg{\\{}$ (47)
$\displaystyle\int d^{4M}\bm{\lambda}\,P_{+}(\bm{\lambda},t)\
\left[\beta_{k}+\frac{\partial}{\partial\alpha_{k}}\right]\widehat{\Lambda}(\bm{\lambda})\Bigg{\\}}\Bigg{]}.\qquad$

This is not of a form amenable to (32). We can soldier on applying integration
by parts and assuming negligible boundary terms to obtain

$\displaystyle{\mathcal{A}}$ $\displaystyle=$ $\displaystyle{\rm
Tr}\Bigg{[}\widehat{a}_{j}(t^{\prime})\breve{V}(t^{\prime},t)\Bigg{\\{}$
$\displaystyle\\!\int
d^{4M}\bm{\lambda}\,\left[\beta_{k}P_{+}(\bm{\lambda},t)-\frac{\partial
P_{+}(\bm{\lambda},t)}{\partial\alpha_{k}}\right]\widehat{\Lambda}(\bm{\lambda})\Bigg{\\}}\Bigg{]}.$

The matter of whether boundary terms can be discarded has been studied in
depth [Gilchrist et al.(1997)Gilchrist, Gardiner, and Drummond, Deuar and
Drummond(2002), Kinsler and Drummond(1991), Smith and Gardiner(1989),
Carusotto and Castin(2003), Deuar(2005)]. The summary is that one can
determine operationally in a stochastic simulation whether boundary terms are
negligible or not. If deemed negligible, then integration by parts is
justified. In (3.2) we can now identify a distribution
${\widetilde{P}}(\bm{\lambda},t)=[\beta_{k}-\frac{\partial}{\partial\alpha_{k}}]P_{+}(\bm{\lambda},t)$
to act on with (32). Doing so gives:

$\displaystyle{\mathcal{A}}$ $\displaystyle=$ $\displaystyle{\rm
Tr}\Bigg{[}\widehat{a}_{j}(t^{\prime})\\!\int\int
d^{4M}\bm{\lambda}\,d^{4M}{\underline{\bm{\lambda}}}$
$\displaystyle\,{\mathcal{P}}(\bm{\lambda},t^{\prime}|{\underline{\bm{\lambda}}},t)\left[{\underline{\beta}}_{k}P_{+}({\underline{\bm{\lambda}}},t)-\frac{\partial
P_{+}({\underline{\bm{\lambda}}},t)}{\partial{\underline{\alpha}}_{k}}\,\right]\widehat{\Lambda}(\bm{\lambda})\Bigg{]}.$

The second operator $\widehat{a}_{j}(t^{\prime})$ can now act on
$\widehat{\Lambda}$ from the left. One obtains

$\displaystyle{\mathcal{A}}$ $\displaystyle=$ $\displaystyle\int\int
d^{4M}\bm{\lambda}\,d^{4M}{\underline{\bm{\lambda}}}\,$
$\displaystyle{\mathcal{P}}(\bm{\lambda},t^{\prime}|{\underline{\bm{\lambda}}},t)\alpha_{j}\left[{\underline{\beta}}_{k}P_{+}({\underline{\bm{\lambda}}},t)-\frac{\partial
P_{+}({\underline{\bm{\lambda}}},t)}{\partial{\underline{\alpha}}_{k}}\,\right].$

While this is formally acceptable (given those negligible boundary terms), and
could be used for some analytic work in small systems such as demonstrated in
[Gardiner(1991)], unfortunately the derivative of $P_{+}$ is not amenable to
interpretation in terms of stochastic samples. At least not the direct samples
we are investigating here. It may be partially treatable using the quantum
jump and response theory approach previously applied to truncated Wigner [Berg
et al.(2009)Berg, Plimak, Polkovnikov, Olsen, Fleischhauer, and Schleich,
Polkovnikov(2010), Plimak and Olsen(2014)], which is a topic for another time.

However – in Sections 5 and 6, a different direct way to evaluate anti-normal
ordered observables such as
${\mathcal{A}}=\langle\widehat{a}_{j}(t^{\prime})\widehat{a}^{\dagger}_{k}(t)\rangle$
will be demonstrated.

##  4 Other phase-space representations

In the derivations of Sec. 3.1, one can see that the crucial aspect for
obtaining a stochastically useful expression is to use only those identities
which do not contain derivatives111Strictly speaking, a small exception to
this appears if the derivative appears only at the final time when the only
remaining operator is $\widehat{\Lambda}$ since
Tr$[\frac{\partial}{\partial\lambda_{\mu}}\widehat{\Lambda}]=\frac{\partial}{\partial\lambda_{\mu}}{\rm
Tr}[\widehat{\Lambda}]=0$ removes any awkward terms. Such a case can, however,
also be treated by an identity without any derivatives after using the cyclic
property of the trace at the right step.. This suggests that convenient
expressions for multi-time correlations similar to (44) will be obtainable
whenever the operators in the correlation can be converted to phase-space
variables without resorting to identities with derivatives.

However, it happens that such identities without derivatives are not
particularly abundant in other phase-space representations. For example, the
Wigner representation [Wigner(1932), Moyal(1949), Polkovnikov(2010)] has
derivatives in all identities, as does its dimension-doubled analogue [Plimak
et al.(2001)Plimak, Olsen, Fleischhauer, and Collett, Hoffmann et
al.(2008)Hoffmann, Corney, and Drummond]. Notably, while the trace of the
symmetric form in the Wigner representation corresponds to $|\alpha|^{2}$ with
no corrections: ${\rm
Tr}\left[\tfrac{1}{2}\left(\widehat{a}^{\dagger}\widehat{a}+\widehat{a}\widehat{a}^{\dagger}\right)\widehat{\Lambda}_{\rm
Wig}\right]=|\alpha|^{2}$, this does not remove derivatives in the
corresponding identity222Here, $\widehat{\Lambda}_{\rm
Wig}=\widehat{\Lambda}_{s}$ from (4.1) with $s=0$.. The best that appears to
be achievable in this way is

$\frac{\widehat{a}^{\dagger}\widehat{a}+\widehat{a}\widehat{a}^{\dagger}}{2}\widehat{\Lambda}_{\rm
Wig}=\left[|\alpha|^{2}+\frac{1}{4}\frac{\partial^{2}}{\partial\alpha\partial\alpha^{*}}\right]\widehat{\Lambda}_{\rm
Wig}.$ (51)

Hence, the path-integral and time-symmetric approach [Berg et al.(2009)Berg,
Plimak, Polkovnikov, Olsen, Fleischhauer, and Schleich] is more suited to the
Wigner representation. Also the phase-space representations developed for spin
systems [Barry and Drummond(2008), Mandt et al.(2015)Mandt, Sadri, Houck, and
Türeci, Ng et al.(2013)Ng, Sørensen, and Deuar], contain derivatives for all
identities.

One notable exception is the Q representation, which admits derivative-free
identities similar to (2.1) for _anti_ -normally ordered operators, and so is
a good candidate for convenient phase-space expressions. Time-ordered anti-
normal operators occur for example in the theory of photon detectors that
operate via emission rather than absorption of photons [Mandel(1966)]. Formal
integral expressions for this kind of correlation were also provided in
[Agarwal and Wolf(1970b)]. Below, the operational stochastic expressions are
found using the Gardiner approach.

###  4.1 The case of the Q representation

The Husimi Q representation [Husimi(1940)] is defined as

$Q(\bm{\alpha})=\frac{1}{\pi^{M}}\langle\bm{\alpha}|\,\widehat{\rho}\,|\bm{\alpha}\rangle$
(52)

and is positive for any $\widehat{\rho}$. Due to the Q distribution being
defined in this explicit way, rather than the implicit form (2.1), observable
expressions in the Q distribution have traditionally been found by simply
expanding the trace:

$\TR\left[{\widehat{O}\,\widehat{\rho}\,}\right]=\frac{1}{\pi^{M}}\int
d^{2M}\bm{\alpha}\
\langle\bm{\alpha}|\widehat{O}\,\widehat{\rho}\,|\bm{\alpha}\rangle,$ (53)

and applying the eigenvalue equation for coherent states

$\widehat{a}_{j}|\bm{\alpha}\rangle=\alpha_{j}|\bm{\alpha}\rangle.$ (54)

For anti-normal ordered correlations at equal times, the cyclic property of
traces gives
$\langle\widehat{a}_{j_{1}}\cdots\widehat{a}_{j_{{\mathcal{N}}}}\widehat{a}^{\dagger}_{k_{1}}\cdots\widehat{a}^{\dagger}_{k_{{\mathcal{M}}}}\rangle=\TR\left[{\widehat{a}^{\dagger}_{k_{1}}\cdots\widehat{a}^{\dagger}_{k_{{\mathcal{M}}}}\widehat{\rho}\
\widehat{a}_{j_{1}}\cdots\widehat{a}_{j_{{\mathcal{N}}}}}\right]$ which
immediately leads to

$\langle\widehat{a}_{j_{1}}\cdots\widehat{a}_{j_{{\mathcal{N}}}}\widehat{a}^{\dagger}_{k_{1}}\cdots\widehat{a}^{\dagger}_{k_{{\mathcal{M}}}}\rangle=\langle\alpha_{j_{1}}\cdots\alpha_{j_{{\mathcal{N}}}}\alpha^{*}_{k_{1}}\cdots\alpha^{*}_{k_{{\mathcal{M}}}}\rangle_{\rm
stoch}.$ (55)

However, this traditional approach fails with time-ordered correlations.
Whatever way one orders the operators, working this way on expression (25)
will lead to a form

$\displaystyle\frac{1}{\pi^{M}}\int
d^{2M}\bm{\alpha}\langle\bm{\alpha}|\Bigg{[}\breve{F}_{R}$
$\displaystyle\breve{V}(\tau_{R},\tau_{R-1})\left\\{\breve{F}_{R-1}\breve{V}(\tau_{R-1},\tau_{R-2})\big{\\{}\cdots\breve{F}_{1}\widehat{\rho}(\tau_{1})\big{\\}}\right\\}\Bigg{]}|\bm{\alpha}\rangle$

in which it is the _latest time_ operator $\breve{F}_{R}(\tau_{R})$ that acts
on the outer coherent states. One can convert $\breve{F}_{R}$ to phase-space
variables using (54), and move the state vectors $|\bm{\alpha}\rangle$ or
$\langle\bm{\alpha}|$ closer to the density matrix and the form (52). However,
in the next putative step, there is no clear way to convert the evolution
operator $\breve{V}(\tau_{R},\tau_{R-1})$ to phase space form. Moreover, when
there is phase-space diffusion, the propagator ${\mathcal{P}}$ is well-behaved
only in the forward time direction, so there is no way to act with
${\mathcal{P}}$ on the outer state vectors and variables which correspond to
later times than the inner ones.

Therefore, we proceed in a non-traditional way, using an implicit form similar
to what was done for the positive-P distribution in Sec. 3.1. The Q
representation is the $s\to-1$ limiting case of the family of s-ordered
representations $W_{s}(\bm{\alpha})$ studied by Cahill and Glauber [Cahill and
Glauber(1969a), Cahill and Glauber(1969b)] (the Glauber-Sudarshan P and Wigner
correspond to $s=1$ and $s=0$, respectively). All these distributions can be
written using coherent displacement operators

$\widehat{D}_{j}(\alpha_{j})=e^{\alpha_{j}\widehat{a}^{\,\dagger}_{j}-\alpha_{j}^{*}\widehat{a}_{j}},\
;\ \widehat{D}(\bm{\alpha})=\prod_{j}\widehat{D}_{j}(\alpha_{j}),$ (57)

in an implicit form similar to (2.1):

$\displaystyle\widehat{\rho}$ $\displaystyle=$ $\displaystyle\int
d^{2M}\bm{\alpha}\,W_{s}(\bm{\alpha})\widehat{\Lambda}_{s}(\bm{\alpha})$ (58)
$\displaystyle\widehat{\Lambda}_{s}(\bm{\alpha})$ $\displaystyle=$
$\displaystyle\prod_{j}\widehat{D}_{j}(\alpha_{j})\widehat{T}_{j}(0,-s)\widehat{D}_{j}(-\alpha_{j}),$
(59)

with the base operator [Cahill and Glauber(1969b)]

$\widehat{T}_{j}(0,-s)=\frac{2}{1+s}\left(\frac{s-1}{1+s}\right)^{\widehat{a}^{\,\dagger}_{j}\widehat{a}_{j}}.$
(60)

One has $\TR\left[{\widehat{\Lambda}_{s}}\right]=1$ and the operator
identities

$\displaystyle\widehat{a}_{j}\widehat{\Lambda}_{s}$ $\displaystyle=$
$\displaystyle\left[\alpha_{j}-\frac{1-s}{2}\frac{\partial}{\partial\alpha^{*}_{j}}\right]\widehat{\Lambda}_{s},$
(61) $\displaystyle\widehat{a}^{\dagger}_{j}\widehat{\Lambda}_{s}$
$\displaystyle=$
$\displaystyle\left[\alpha^{*}_{j}+\frac{1+s}{2}\frac{\partial}{\partial\alpha_{j}}\right]\widehat{\Lambda}_{s},$
(62) $\displaystyle\widehat{\Lambda}_{s}\widehat{a}_{j}$ $\displaystyle=$
$\displaystyle\left[\alpha_{j}+\frac{1+s}{2}\frac{\partial}{\partial\alpha^{*}_{j}}\right]\widehat{\Lambda}_{s},$
(63) $\displaystyle\widehat{\Lambda}_{s}\widehat{a}^{\dagger}_{j}$
$\displaystyle=$
$\displaystyle\left[\alpha^{*}_{j}-\frac{1-s}{2}\frac{\partial}{\partial\alpha_{j}}\right]\widehat{\Lambda}_{s},$
(64)

which can be verified by equating the LHS and RHS when $\widehat{T}(0,-s)$ is
expanded in number states. Therefore, in the limit $s\to-1$ corresponding to
the Q representation with $\widehat{\Lambda}_{-1}=\widehat{\Lambda}_{Q}$,

$\displaystyle\widehat{a}_{j}\widehat{\Lambda}_{Q}$ $\displaystyle=$
$\displaystyle\left[\alpha_{j}-\frac{\partial}{\partial\alpha^{*}_{j}}\right]\widehat{\Lambda}_{Q},$
(65) $\displaystyle\widehat{a}^{\dagger}_{j}\widehat{\Lambda}_{Q}$
$\displaystyle=$ $\displaystyle\alpha^{*}_{j}\widehat{\Lambda}_{Q},$ (66)
$\displaystyle\widehat{\Lambda}_{Q}\widehat{a}_{j}$ $\displaystyle=$
$\displaystyle\alpha_{j}\widehat{\Lambda}_{Q},$ (67)
$\displaystyle\widehat{\Lambda}_{Q}\widehat{a}^{\dagger}_{j}$ $\displaystyle=$
$\displaystyle\left[\alpha^{*}_{j}-\frac{\partial}{\partial\alpha_{j}}\right]\widehat{\Lambda}_{Q}.$
(68)

We see then that multi-time correlations in which only the orderings
$\widehat{a}^{\dagger}_{j}\widehat{\Lambda}_{Q}$ and
$\widehat{\Lambda}_{Q}\widehat{a}_{j}$ appear could have the capacity to
correspond to simple stochastic expressions. This implies anti-normal ordered
moments since by the cyclic property of traces

$\TR\left[{\widehat{a}^{\dagger}_{q_{1}}\cdots\widehat{a}^{\dagger}_{q_{{\mathcal{M}}}}\widehat{\rho}\,\widehat{a}_{p_{1}}\cdots\widehat{a}_{p_{{\mathcal{N}}}}}\right]=\langle\widehat{a}_{p_{1}}\cdots\widehat{a}_{p_{{\mathcal{N}}}}\widehat{a}^{\dagger}_{q_{1}}\cdots\widehat{a}^{\dagger}_{q_{{\mathcal{M}}}}\rangle.$
(69)

The implicit form (4.1) for $\widehat{\rho}$ allows us to proceed in a similar
fashion to Sec. 3.1. Consider then the anti-normal and time-ordered
correlation

${\mathcal{A}}^{\prime}=\langle\widehat{a}_{p_{1}}(t_{1})\cdots\widehat{a}_{p_{{\mathcal{N}}}}(t_{{\mathcal{N}}})\widehat{a}^{\dagger}_{q_{1}}(s_{1})\cdots\widehat{a}^{\dagger}_{q_{{\mathcal{M}}}}(s_{{\mathcal{M}}})\rangle$
(70)

with times obeying (2.2). A potentially awkward issue is that the kernel
$\widehat{\Lambda}_{s}$ is not very well bounded in the $s\to-1$ limit, with
$\langle\bm{\beta}|\widehat{\Lambda}_{s}(\bm{\alpha})|\bm{\beta}\rangle=\frac{2}{\varepsilon}\exp\left[-2|\bm{\alpha}-\bm{\beta}|^{2}/\varepsilon\right]$
where $s=\varepsilon-1$. To gauge whether this is a problem, we will work
using infinitesimal $\varepsilon>0$ in which case

$\widehat{\Lambda}_{s}\to\widehat{\Lambda}_{Q}+{\mathcal{O}}(\varepsilon).$
(71)

Using the form (25) on (69) with now $B=\widehat{a}^{\dagger}$ and
$A=\widehat{a}$, one obtains

$\displaystyle{\mathcal{A}}^{\prime}={\mathcal{O}}(\varepsilon)+\int
d^{2M}\bm{\alpha}\,{\rm
Tr}\Bigg{[}\breve{F}_{R}\breve{V}(\tau_{R},\tau_{R-1})\Bigg{\\{}$
$\displaystyle\cdots\breve{F}_{3}\breve{V}(\tau_{3},\tau_{2})\Big{\\{}\breve{F}_{2}\breve{V}(\tau_{2},\tau_{1})\big{\\{}W_{s}(\bm{\alpha},\tau_{1})\breve{F}_{1}\widehat{\Lambda}_{Q}(\bm{\alpha})\big{\\}}\Big{\\}}\Bigg{\\}}\Bigg{]}.$

The part in the inner $\\{\ \\}$ brackets will be either

$\begin{array}[]{l@{\qquad\qquad}l}W_{s}(\bm{\alpha},t_{1})\widehat{\Lambda}_{Q}\widehat{a}_{p_{1}}&\text{if$t_1\leq
s_{{\mathcal{M}}}$}\\\
W_{s}(\bm{\alpha},s_{{\mathcal{M}}})\widehat{a}^{\dagger}_{q_{{\mathcal{M}}}}\widehat{\Lambda}_{Q}&\text{otherwise.}\end{array}$
(73)

Both cases allow us to use derivative-free identities (67) or (66), to replace
(73) with ${\widetilde{P}}_{s}(\bm{\alpha},\tau_{1})\widehat{\Lambda}_{Q}$
where

$\begin{array}[]{l@{\quad\qquad}l}{\widetilde{P}}_{s}=\alpha_{p_{1}}W_{s}(\bm{\alpha},s_{{\mathcal{M}}})&\text{if$t_1\leq
s_{{\mathcal{M}}}$}\\\
{\widetilde{P}}_{s}=\alpha^{*}_{q_{{\mathcal{M}}}}W_{s}(\bm{\alpha},t_{1})&\text{otherwise.}\end{array}$
(74)

Now assuming an acceptable propagator
${\mathcal{P}}_{s}(\bm{\alpha},\tau_{2}|{\underline{\bm{\alpha}}},\tau_{1})$
exists for the evolution of the s-ordered distribution $W_{s}$, use of (32)
leads to

$\displaystyle{\mathcal{A}}^{\prime}={\mathcal{O}}(\varepsilon)+\int
d^{2M}\bm{\alpha}\,\int d^{2M}{\underline{\bm{\alpha}}}\,{\rm
Tr}\Bigg{[}\breve{F}_{R}\breve{V}(\tau_{R},\tau_{R-1})\Bigg{\\{}$ (75)
$\displaystyle\cdots\breve{F}_{3}\breve{V}(\tau_{3},\tau_{2})\Big{\\{}{\mathcal{P}}_{s}(\bm{\alpha},\tau_{2}|{\underline{\bm{\alpha}}},\tau_{1}){\widetilde{P}}_{s}({\underline{\bm{\alpha}}},\tau_{1})\breve{F}_{2}\widehat{\Lambda}_{Q}\Big{\\}}\Bigg{\\}}\Bigg{]},\qquad$

in which the objects related to the first time interval have been fully
converted to phase-space quantities. Proceeding in this fashion with
increasing time for the remaining $\widehat{a}$ and $\widehat{a}^{\dagger}$
operators, one arrives at

$\displaystyle{\mathcal{A}}^{\prime}=\bigg{[}\int
d^{2M}\bm{\alpha}(\tau_{1})\,\int d^{2M}\bm{\alpha}(\tau_{2})\,\cdots\,\int
d^{2M}\bm{\alpha}(\tau_{R})$ (76)
$\displaystyle\qquad\qquad\qquad\alpha_{p_{1}}(t_{1})\cdots\alpha_{p_{{\mathcal{N}}}}(t_{{\mathcal{N}}})\alpha^{*}_{q_{1}}(s_{1})\cdots\alpha^{*}_{q_{{\mathcal{M}}}}(s_{{\mathcal{M}}})\times$
$\displaystyle{\mathcal{P}}_{s}(\bm{\alpha}(\tau_{R}),\tau_{R}|\bm{\alpha}(\tau_{R-1}),\tau_{R-1})\cdots{\mathcal{P}}_{s}(\bm{\alpha}(\tau_{2}),\tau_{2}|\bm{\alpha}(\tau_{1}),\tau_{1})$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\times
W_{s}(\bm{\alpha}(\tau_{1}),\tau_{1})\bigg{]}+{\mathcal{O}}(\varepsilon).$

Variables $\bm{\alpha},{\underline{\bm{\alpha}}}$ etc were relabelled to
$\bm{\alpha}(\tau_{R}),\bm{\alpha}(\tau_{R-1})$.

In the limit $\varepsilon\to 0$ that we are considering,
$W_{s}(\bm{\alpha})\to Q(\bm{\alpha})$, which is real non-negative. To be
consistent, the propagator ${\mathcal{P}}_{-1}$ must also be real nonnegative.
When the master equation (10) is faithfully reproduced by a Fokker-Planck
equation for $Q(\bm{\alpha})$, this will be the case. Then, the
${\mathcal{P}}_{s}\dots{\mathcal{P}}_{s}W_{s}$ factors can be interpreted
similarly to (35) as the joint probability of samples
$\bm{\alpha}(\tau_{1}),\bm{\alpha}(\tau_{2}),\dots$ at successive times. With
that, we arrive at the hoped for result that a time ordered (as (2.2)) and
_anti_ -normally ordered correlation is evaluated as

$\displaystyle\langle\widehat{a}_{p_{1}}(t_{1})\cdots\widehat{a}_{p_{{\mathcal{N}}}}(t_{{\mathcal{N}}})\widehat{a}^{\dagger}_{q_{1}}(s_{1})\cdots\widehat{a}^{\dagger}_{q_{{\mathcal{M}}}}(s_{{\mathcal{M}}})\rangle$
(77)
$\displaystyle=\langle\alpha_{p_{1}}(t_{1})\cdots\alpha_{p_{{\mathcal{N}}}}(t_{{\mathcal{N}}})\alpha^{*}_{q_{1}}(s_{1})\cdots\alpha^{*}_{q_{{\mathcal{M}}}}(s_{{\mathcal{M}}})\rangle_{\rm
stoch}\qquad$

using _Q representation_ samples $\alpha_{j}$. Stochastic averaging is over
products constructed using values from the evolution of a single sample, like
in (39).

##  5 Evaluation by conversion to the Q representation

###  5.1 Conversion between samples of s-ordered distributions

The s-ordered distributions $W_{s}$ introduced in (4.1) are mutually related
by [Cahill and Glauber(1969b)]

$\displaystyle W_{s}(\bm{\alpha}^{\prime})$
$\displaystyle=\left(\frac{2}{s_{0}-s}\right)^{M}\int\frac{d^{2M}\bm{\alpha}}{\pi^{M}}\exp\left[-\frac{2|\bm{\alpha}^{\prime}-\bm{\alpha}|^{2}}{s_{0}-s}\right]W_{s_{0}}(\bm{\alpha}).$

in the sense that $W_{s}$ and $W_{s_{0}}$ represent the same quantum density
matrix $\widehat{\rho}$. When $s_{0}>s$, this is a Gaussian convolution of the
more normally-ordered distribution $W_{s_{0}}$, and reflects the well known
property that Q distributions ($s=-1$) are broader and more smoothed than
Wigner ($s=0$), which are in turn broader than Glauber-Sudarshan P
distributions ($s=1$) for the same state. Importantly for us here, this means
that if we have samples $\bm{\alpha}$ of a more normally ordered distribution,
we can easily also obtain samples $\bm{\alpha}^{\prime}$ of the less normally
ordered distributions simply by adding Gaussian noise. The prescription is

$\alpha^{\prime}_{j}=\alpha_{j}+\sqrt{\frac{s_{0}-s}{2}}\ \zeta_{j}$ (79)

for each mode $j$, with each $\zeta_{j}$ a complex random variable of variance
1:

$\langle\zeta_{j}\rangle_{\rm stoch}=0\ ;\
\langle\zeta_{j}\zeta_{k}\rangle_{\rm stoch}=0\ ;\
\langle\zeta^{*}_{j}\zeta_{k}\rangle_{\rm stoch}=1.$ (80)

In particular, converting samples of a P distribution to samples of Q requires

$\alpha^{\prime}_{j}=\alpha_{j}+\zeta_{j}.$ (81)

Converting samples of a Wigner distribution (if it is nonnegative to begin
with) to samples of Q can be done with

$\alpha^{\prime}_{j}=\alpha^{\rm Wig}_{j}+\frac{\zeta_{j}}{\sqrt{2}}.$ (82)

###  5.2 Evaluation of anti-normal ordered moments starting from P and Wigner
representations

Therefore, if one has samples $\bm{\alpha}(t_{0})$ of a P or Wigner
representation up to a time $t_{0}$ (e.g. from a prior stochastic evolution),
the prescription (79) can be used to convert them to samples
$\bm{\alpha}^{\prime}(t_{0})$ of the Q representation at that time. The noises
$\zeta_{j}$ are generated just once at this time. Subsequent evolution
according to the Q stochastic equations then leads to Q samples
$\bm{\alpha}^{\prime}(t)$ at later times $t>t_{0}$. These can be directly used
in expression (77) to evaluate anti-normal ordered multi-time observables.

This is potentially a little less convenient that remaining in one
distribution throughout because the conversion time $t_{0}$ has to be chosen
before starting a simulation, and the anti-normal ordered operators cannot
extend to times before $t_{0}$. It does preserve the principal advantages of
phase-space simulation, though: intuitive and computationally tractable
expressions for observables, stochasticity, gentle scaling with system size.
The evolution equations (13) are usually of similar form in all s-ordered
representations, apart from simplification at special $s$ values.

One cannot convert the other way with this procedure (e.g. from Wigner to P),
so normally ordered multi-time correlations cannot be extracted this way from
samples of Wigner or Q distributions.

###  5.3 Mixed-order moments

The above discussion suggests a way that some mixed-order multi-time
correlations that contain both normal and anti-normal factors could be
evaluated. Suppose early time operators ($t\leq t_{0}$) are normally ordered
(can be evaluated using the P distribution), while later time operators
($t>t_{0}$) are anti-normally ordered (can be evaluated using the Q
representation). A switch according to (81) could be made at $t_{0}$ after
evaluating any normally-ordered factors, and the Q distribution evolved and
used at later times to obtain the remaining inner factors with $t>t_{0}$. Let
us check in detail whether this is feasible.

Consider the time-ordered but neither anti- or normally ordered correlation

$\displaystyle{\mathcal{A}}^{\prime\prime}=\langle\widehat{a}^{\dagger}_{p_{1}}(t_{1})\widehat{a}_{p_{2}}(t_{2})\widehat{a}^{\dagger}_{q_{1}}(s_{1})\widehat{a}_{q_{2}}(s_{2})\rangle$
$\displaystyle=\int d^{2M}\bm{\alpha}\,{\rm
Tr}\bigg{[}\,\widehat{a}^{\dagger}_{q_{1}}\breve{V}(s_{1},t_{2})\bigg{\\{}$
(83)
$\displaystyle\breve{V}(t_{2},s_{2})\Big{\\{}\widehat{a}_{q_{2}}\breve{V}(s_{2},t_{1})\big{\\{}P(\bm{\alpha},t_{1})\widehat{\Lambda}_{1}(\bm{\alpha})\widehat{a}^{\dagger}_{p_{1}}\big{\\}}\Big{\\}}\widehat{a}_{p_{2}}\bigg{\\}}\bigg{]}$

whose times satisfy (2.2), and $t_{1}\leq s_{2}\leq t_{2}<s_{1}$ (for
instance). The initial expansion of $\widehat{\rho}$ is made in the $P$
representation, where $P(\bm{\alpha})=W_{1}(\bm{\alpha})$ and
$\widehat{\Lambda}_{1}(\bm{\alpha})=|\bm{\alpha}\rangle\langle\bm{\alpha}|$.
The first two operators $\widehat{a}^{\dagger}_{p_{1}}$,
$\widehat{a}_{q_{2}}$, and first two $\breve{V}$ convert easily to phase-space
variables via (6) and (9), giving

$\displaystyle{\mathcal{A}}^{\prime\prime}=\iiint
d^{2M}\bm{\alpha}\,d^{2M}{\underline{\bm{\alpha}}}\,d^{2M}{\underline{{\underline{\bm{\alpha}}}}}\,{\rm
Tr}\Bigg{[}\ \widehat{a}^{\dagger}_{q_{1}}\breve{V}(s_{1},t_{2})\Bigg{\\{}$
(84)
$\displaystyle{\mathcal{P}}_{1}(\bm{\alpha},t_{2}|{\underline{{\underline{\bm{\alpha}}}}},s_{2}){\mathcal{P}}_{1}({\underline{{\underline{\bm{\alpha}}}}},s_{2}|{\underline{\bm{\alpha}}},t_{1})P({\underline{\bm{\alpha}}},t_{1})\,{\underline{{\underline{\alpha}}}}_{q_{2}}\,{\underline{\alpha}}^{*}_{p_{1}}\widehat{\Lambda}_{1}(\bm{\alpha})\widehat{a}_{p_{2}}\Bigg{\\}}\Bigg{]}.$

Further work by this route is now closed because
$\widehat{\Lambda}_{1}\widehat{a}_{p_{2}}$ involves (8) and derivatives.
Instead, another result that follows from from [Cahill and Glauber(1969a)]
allows us to convert kernels:

$\displaystyle\widehat{\Lambda}_{s}(\bm{\alpha})$
$\displaystyle=\left(\frac{2}{s-s_{0}}\right)^{M}\int\frac{d^{2M}\bm{\alpha}^{\prime}}{\pi^{M}}\exp\left[-\frac{2|\bm{\alpha}^{\prime}-\bm{\alpha}|^{2}}{s-s_{0}}\right]\widehat{\Lambda}_{s_{0}}(\bm{\alpha}^{\prime}).$

Taking $s=1$ and $s_{0}=-1$, to move to a Q representation kernel, one finds

$\displaystyle\widehat{\Lambda}_{1}(\bm{\alpha})=\int\frac{d^{2M}\bm{\zeta}}{\pi^{M}}\exp\left[-|\bm{\zeta}|^{2}\right]\widehat{\Lambda}_{Q}(\bm{\alpha}^{\prime}).$
(86)

where $\bm{\zeta}=\bm{\alpha}^{\prime}-\bm{\alpha}=\\{\zeta_{j}\\}$ has
exactly the properties of the noise in (80), and
$\bm{\alpha}^{\prime}=\bm{\alpha}+\bm{\zeta}$ is given by (81). After
substituting (86) into (84), applying (67) to
$\widehat{\Lambda}_{Q}(\bm{\alpha}^{\prime})\widehat{a}_{p_{2}}$, and defining
the distribution

$\displaystyle{\widetilde{P}}(\bm{\alpha}^{\prime},t_{2})$ $\displaystyle=$
$\displaystyle\alpha^{\prime}_{p_{2}}\\!\\!\int\frac{d^{2M}\bm{\alpha}}{\pi^{M}}e^{-|\bm{\alpha}^{\prime}-\bm{\alpha}|^{2}}\\!\\!\\!\iint\\!\\!d^{2M}{\underline{\bm{\alpha}}}\,d^{2M}{\underline{{\underline{\bm{\alpha}}}}}\qquad$
$\displaystyle{\mathcal{P}}_{1}(\bm{\alpha},t_{2}|{\underline{{\underline{\bm{\alpha}}}}},s_{2}){\mathcal{P}}_{1}({\underline{{\underline{\bm{\alpha}}}}},s_{2}|{\underline{\bm{\alpha}}},t_{1})P({\underline{\bm{\alpha}}},t_{1})\,{\underline{{\underline{\alpha}}}}_{q_{2}}\,{\underline{\alpha}}^{*}_{p_{1}}.$

we have

${\mathcal{A}}^{\prime\prime}=\int d^{2M}\bm{\alpha}^{\prime}\,{\rm
Tr}\Bigg{[}\
\widehat{a}^{\dagger}_{q_{1}}\breve{V}(s_{1},t_{2})\Bigg{\\{}{\widetilde{P}}(\bm{\alpha}^{\prime},t_{2})\,\widehat{\Lambda}_{Q}(\bm{\alpha}^{\prime})\Bigg{\\}}\Bigg{]}.$
(88)

Now using (32), and taking care with variable labels:

$\displaystyle{\mathcal{A}}^{\prime\prime}$ $\displaystyle=$
$\displaystyle\iint
d^{2M}\bm{\alpha}^{\prime}\,d^{2M}{\underline{\bm{\alpha}}}^{\prime}\
{\mathcal{P}}_{-1}(\bm{\alpha}^{\prime},s_{1}|{\underline{\bm{\alpha}}}^{\prime},t_{2}){\widetilde{P}}({\underline{\bm{\alpha}}}^{\prime},t_{2})\,$
(89) $\displaystyle\qquad\qquad\qquad\qquad\times{\rm
Tr}\left[\widehat{a}^{\dagger}_{q_{1}}\widehat{\Lambda}_{Q}(\bm{\alpha}^{\prime})\right].$

After the variable change caused by (32), now

${\underline{\bm{\alpha}}}^{\prime}=\bm{\alpha}+\bm{\zeta}.$ (90)

Notice also the $-1$ label on the latest propagator in (89), indicating that
it is according to the Q representation equations in the time interval
$(t_{2},s_{2}]$. Finally, applying (66), one arrives at:

$\displaystyle{\mathcal{A}}^{\prime\prime}=\int
d^{2M}{\underline{\bm{\alpha}}}\,d^{2M}{\underline{{\underline{\bm{\alpha}}}}}\,d^{2M}\bm{\alpha}\,d^{2M}\bm{\zeta}\,\left(\frac{e^{-|\bm{\zeta}|^{2}}}{\pi^{M}}\right)\,d^{2M}\bm{\alpha}^{\prime}$
$\displaystyle\qquad\times\
{\underline{{\underline{\alpha}}}}_{q_{2}}\,{\underline{\alpha}}^{*}_{p_{1}}(\alpha_{p_{2}}+\zeta_{p_{2}})\alpha^{\prime*}_{q_{1}}{\mathcal{P}}_{-1}(\bm{\alpha}^{\prime},s_{1}|\bm{\alpha}+\zeta,t_{2})$
$\displaystyle\qquad\times{\mathcal{P}}_{1}(\bm{\alpha},t_{2}|{\underline{{\underline{\bm{\alpha}}}}},s_{2}){\mathcal{P}}_{1}({\underline{{\underline{\bm{\alpha}}}}},s_{2}|{\underline{\bm{\alpha}}},t_{1})P({\underline{\bm{\alpha}}},t_{1}).$
(91)

This encodes the following sequence of operations:

  1. 1.

Start with initial samples ${\underline{\bm{\alpha}}}$ at $t_{1}$.

  2. 2.

Propagate P representation equations to $s_{2}$ obtaining samples
${\underline{{\underline{\bm{\alpha}}}}}$.

  3. 3.

Propagate P representation equations to $t_{2}$ obtaining samples
$\bm{\alpha}$.

  4. 4.

Add Gaussian noise as per (90), to get samples
${\underline{\bm{\alpha}}}^{\prime}$.

  5. 5.

Propagate Q representation equations to $s_{1}$ obtaining samples
$\bm{\alpha}^{\prime}$.

Along the way, samples are collected to use in the final stochastic average.
The ${\mathcal{P}}_{s}$ factors in (5.3) together with the Gaussian factor in
the top line form the joint probability
$P(\bm{\alpha}^{\prime},s_{1};{\underline{\bm{\alpha}}}^{\prime},t_{2};{\underline{{\underline{\bm{\alpha}}}}},s_{2};{\underline{\bm{\alpha}}},t_{1})$.
Therefore, the final estimator for the observable
${\mathcal{A}}^{\prime\prime}$ is

$\displaystyle{\mathcal{A}}^{\prime\prime}=\langle\widehat{a}^{\dagger}_{p_{1}}(t_{1})\widehat{a}_{p_{2}}(t_{2})\widehat{a}^{\dagger}_{q_{1}}(s_{1})\widehat{a}_{q_{2}}(s_{2})\rangle=$
$\displaystyle\langle\alpha_{p_{1}}^{*}(t_{1})\alpha^{\prime}_{p_{2}}(t_{2})\alpha^{\prime*}_{q_{1}}(s_{1})\alpha_{q_{2}}(s_{2})\rangle_{\rm
stoch}.$

Primed variables $\alpha^{\prime}$ are samples of the Q distribution, while
un-primed ones $\alpha$ are samples of the P distribution. This all matches
intuitively with the evolution and the expectation that anti-normally ordered
elements will use Q distribution samples and normally-ordered elements samples
of the P distribution.

###  5.4 Most general ordering case

The widest generalisation of this procedure to other (time-ordered) cases is
as follows: If there is an early time set of normally ordered operators, on
either side of the correlation, it can be dealt with by sampling the P
distribution according to the replacements

$\widehat{a}_{j}\to\alpha_{j}\qquad;\qquad\widehat{a}^{\dagger}_{j}\to\alpha_{j}^{*}.$
(93)

Once this avenue becomes exhausted, one adds noise via (81) to convert
$\alpha$’s to Q distribution samples $\alpha^{\prime}$. If the remaining later
time inner factors are anti-normally ordered, they can then be dealt with
using the replacement

$\widehat{a}_{j}\to\alpha^{\prime}_{j}\qquad;\qquad\widehat{a}^{\dagger}_{j}\to\alpha_{j}^{\prime*}.$
(94)

The above covers both fully normal and fully anti-normal ordered products as
special cases. On the other hand, the case of an early time anti-normally
ordered block and a later time normally ordered block containing several times
is not amenable to this approach because one cannot stochastically convert Q
samples to P samples.

A large number of cases can also be reduced to a form amenable to this
procedure by the use of the commutator $[\widehat{a},\widehat{a}^{\dagger}]=1$
for equal time factors. Also, simulations starting in a Wigner representation
can be used for evaluation of outer symmetric ordered parts, and then a switch
can be made to the Q representation via (82) to evaluate any remaining inner
anti-normal ordered parts.

###  5.5 Correlation function coverage

Table 1 counts the number of distinct 1, 2, and 3 annihilation/creation
operator products that can be evaluated using the various methodology
described above. All one and two-factor combinations can be evaluated (though
4/12 of the latter require use of the Q representation). For third order
correlations, almost all time-ordered cases can be evaluated (74 out of 80).
The majority require use of the Q representation or a shift from a P to Q
representation as described in Sec. 5.4 to work. There are 6 exceptions that
cannot be evaluated:
$\langle\widehat{a}^{\dagger}(t_{2})\widehat{a}^{\dagger}(t_{3})\widehat{a}^{\dagger}(t_{1})\rangle$,
$\langle\widehat{a}(t_{2})\widehat{a}^{\dagger}(t_{3})\widehat{a}(t_{1})\rangle$,
$\langle\widehat{a}(t_{1})\widehat{a}(t_{3})\widehat{a}(t_{2})\rangle$,
$\langle\widehat{a}(t_{1})\widehat{a}^{\dagger}(t_{3})\widehat{a}(t_{2})\rangle$,
$\langle\widehat{a}(t_{3})\widehat{a}(t_{2})\widehat{a}^{\dagger}(t_{1})\rangle$,
$\langle\widehat{a}^{\dagger}(t_{3})\widehat{a}(t_{2})\widehat{a}^{\dagger}(t_{1})\rangle$,
where $t_{1}<t_{2}<t_{3}$ is assumed. These all share the feature that the
earliest factor already requires the Q representation, while the next operator
in time requires the P representation to which one cannot return. There are
also a number of correlations that are not time ordered, for which the
earliest time $t_{1}$ is on the middle operator, and these are not possible to
evaluate according to the schemes presented here.

Order | 1st | 2nd | 3rd
---|---|---|---
(number of operators) | order | order | order
Total permutations | 2 | 12 | 104
single time correlations | 2 | 4 | 8
multi-time accessible |  |  |
with P representation | – | 4 | 22
additional accessible |  |  |
with Q representation | – | 4 | 22
additional accessible |  |  |
with mixed order (Sec. 5.4) | – | – | 22
Total doable | 2 | 12 | 74
time ordered not doable | – | – | 6
Not time ordered, not doable | – | – | 24
Table 1: A tally of $\widehat{a}$, $\widehat{a}^{\dagger}$ product
permutations that can/cannot be evaluated with the various approaches
discussed. The general form considered is
$\langle\widehat{A}(t_{a})\widehat{B}(t_{b})\widehat{C}(t_{c})\rangle$, where
$\widehat{A},\widehat{B},\widehat{C}$ can be either of $\widehat{a}$ or
$\widehat{a}^{\dagger}$ (same mode), and the time arguments can take up to
three distinct times $t_{1}<t_{2}<t_{3}$. Permutations with the same time
topology (e.g. $\widehat{A}(t_{1})\widehat{B}(t_{1})\widehat{C}(t_{2})$ and
$\widehat{A}(t_{2})\widehat{B}(t_{2})\widehat{C}(t_{3})$) are counted only
once.

The greatest interest in multi-time correlations usually concerns those
involving two times (say $t=0$ and $t=\tau>0$). Examples are counting
correlations like
$\langle\widehat{a}^{\dagger}(0)\widehat{a}^{\dagger}(\tau)\widehat{a}(\tau)\widehat{a}(0)\rangle$,
and pair correlations such as
$\langle\widehat{a}^{\dagger}(0)\widehat{a}^{\dagger}(0)\widehat{a}(\tau)\widehat{a}(\tau)\rangle$.
The case count for these is summarised by Table 2. A we can see, all 160 kinds
of time ordered four-operator products of this form can be evaluated,
including atypical combinations such as
$\langle\widehat{a}(0)\widehat{a}^{\dagger}(\tau)\widehat{a}^{\dagger}(\tau)\widehat{a}(0)\rangle$,
but very many require the Q representation. (72 out of the 160 accessible ones
require the use of the mid-simulation switching to the Q representation
described in Sec. 5.4). The only correlations that are inaccessible are the
non time ordered ones such as e.g.
$\langle\widehat{a}^{\dagger}(\tau)\widehat{a}^{\dagger}(0)\widehat{a}(0)\widehat{a}(\tau)\rangle$.

Order | 2nd | 3rd | 4th
---|---|---|---
(number of operators) | order | order | order
Total permutations | 12 | 56 | 240
single time correlations | 4 | 8 | 16
multi-time accessible |  |  |
with P representation | 4 | 14 | 36
additional accessible |  |  |
with Q representation | 4 | 14 | 36
additional accessible |  |  |
with mixed order (Sec. 5.4) | – | 12 | 72
Total doable | 12 | 48 | 160
time ordered not doable | – | – | –
Not time ordered, not doable | – | 8 | 80
Table 2: A tally of $\widehat{a}$, $\widehat{a}^{\dagger}$ products involving
up to four operators, evaluated at one of two times. The general form
considered is
$\langle\widehat{A}(t_{a})\widehat{B}(t_{b})\widehat{C}(t_{c})\widehat{D}(t_{d})\rangle$,
where $\widehat{A},\widehat{B},\widehat{C},\widehat{D}$ can be either of
$\widehat{a}$ or $\widehat{a}^{\dagger}$ (same mode), and the time arguments
can take up to two distinct times $t=0$ and $t=\tau>0$.

The change of distribution can introduce additional restrictions. In
particular, Wigner and Q distributions have a tendency to involve higher order
derivatives in the PDE (11) than P distributions, so that a standard diffusion
process no longer captures the full quantum dynamics. An example are two-
photon losses with operators $\widehat{R}=\widehat{a}^{2}$ for which P
distributions produce only 2nd order derivatives but Q distributions 4th order
terms, and Wigner distributions 3rd order ones. This precludes fully accurate
calculation of correlation functions such as
$\langle\widehat{a}^{\dagger}(0)\widehat{a}(t)\widehat{a}(t)\widehat{a}^{\dagger}(0)\rangle=\langle\alpha^{*}(0)\alpha^{\prime}(t)\alpha^{\prime}(t)\alpha^{\prime*}(0)\rangle_{\rm
stoch}$ requiring Q evolution after the changeover, but not cases where the
change to a Q distribution is only needed at the final time such as
$\langle\widehat{a}^{\dagger}(0)\widehat{a}(t)\widehat{a}^{\dagger}(t)\widehat{a}(0)\rangle=\langle\alpha^{*}(0)\alpha^{\prime}(t)\alpha^{\prime*}(t)\alpha(0)\rangle_{\rm
stoch}$. Stochastic techniques for dealing with higher order PDE terms have
been investigated [Jumarie(1999), Olsen et al.(2002)Olsen, Plimak, and
Fleischhauer, Jumarie(2005), Plimak et al.(2003)Plimak, Fleischhauer, Olsen,
and Collett] particularly in [Drummond(2014)] for doubled phase space, though
attempts to date have shown strong time and stability limitations.

##  6 Q representations and s-ordering in doubled phase space

To use the mechanisms described in Sec. 5 in the doubled phase-space
representations that give more complete coverage of quantum mechanics, we need
also to consider the doubled phase space Q representation and how to switch to
it from the positive-P.

###  6.1 Doubled phase space s-ordered representations

The s-ordered representations were first generalised to doubled phase space by
de Oliveira [de Oliveira(1992)]. Later studies [Plimak et al.(2001)Plimak,
Olsen, Fleischhauer, and Collett, Hoffmann et al.(2008)Hoffmann, Corney, and
Drummond, Deuar(2009)] used a different normalisation that is closer to the
original single-phase space formulation (57) of Cahill and Glauber [Cahill and
Glauber(1969a)]. It will also be used here. The explicit expansion of the
density matrix was given in [Deuar(2009)] as333In the supplemental material
therein.:

$\displaystyle\widehat{\rho}$ $\displaystyle=$ $\displaystyle\int
d^{4M}\bm{\lambda}\,W^{+}_{s}(\bm{\lambda})\widehat{\Lambda}^{+}_{s}(\bm{\lambda})$
(95) $\displaystyle\widehat{\Lambda}^{+}_{s}(\bm{\lambda})$ $\displaystyle=$
$\displaystyle\prod_{j}\widehat{d}_{j}(\alpha_{j},\beta_{j})\widehat{T}_{j}(0,-s)\widehat{d}_{j}(-\alpha_{j},-\beta_{j}),\qquad$
(96)

where the displacement-like operator $\widehat{d}$ is

$\widehat{d}_{j}(\alpha_{j},\beta_{j})=e^{\alpha_{j}\widehat{a}^{\dagger}_{j}-\beta_{j}\widehat{a}_{j}}\quad;\quad\widehat{d}_{j}(\alpha,\alpha^{*})=\widehat{D}_{j}(\alpha),$
(97)

obtained by the replacement $\alpha^{*}\to\beta$ in (57). To distinguish from
single phase space, the superscript + is used where necessary. The usual
properties $\TR\left[{\widehat{\Lambda}^{+}_{s}}\right]=1$ and
$\widehat{d}^{-1}(\alpha,\beta)=\widehat{d}(-\alpha,-\beta)$ continue to
apply. The operator identities are now

$\displaystyle\widehat{a}_{j}\widehat{\Lambda}^{+}_{s}$ $\displaystyle=$
$\displaystyle\left[\alpha_{j}-\frac{1-s}{2}\frac{\partial}{\partial\beta_{j}}\right]\widehat{\Lambda}^{+}_{s},$
(98) $\displaystyle\widehat{a}^{\dagger}_{j}\widehat{\Lambda}^{+}_{s}$
$\displaystyle=$
$\displaystyle\left[\beta_{j}+\frac{1+s}{2}\frac{\partial}{\partial\alpha_{j}}\right]\widehat{\Lambda}^{+}_{s},$
(99) $\displaystyle\widehat{\Lambda}^{+}_{s}\widehat{a}_{j}$ $\displaystyle=$
$\displaystyle\left[\alpha_{j}+\frac{1+s}{2}\frac{\partial}{\partial\beta_{j}}\right]\widehat{\Lambda}^{+}_{s},$
(100) $\displaystyle\widehat{\Lambda}^{+}_{s}\widehat{a}^{\dagger}_{j}$
$\displaystyle=$
$\displaystyle\left[\beta_{j}-\frac{1-s}{2}\frac{\partial}{\partial\alpha_{j}}\right]\widehat{\Lambda}^{+}_{s}.$
(101)

The $s\to 1$ limit of all the above gives the positive-P representation, $s=0$
the doubled-Wigner representation of [Hoffmann et al.(2008)Hoffmann, Corney,
and Drummond], and the limit $s\to-1$ a doubled phase space analogue to the Q
representation (“doubled-Q”). Eqs. (6.1) are equivalent to the correspondences
found in [de Oliveira(1992)]. An advantage of the doubled-Q and doubled Wigner
representations relative to their single-phase space analogues is that all 2nd
order derivative terms in the PDE can be made positive-definite and converted
fully to stochastic equations via the same analytic kernel trick as for the
positive-P representation. For example, spontaneous emission in a Schwinger
boson system need not be amputated like in [Huber et al.(2021)Huber, Kirton,
and Rabl] for the Wigner representation.

The kernel transform between different orderings in doubled phase-space is
found to be

$\displaystyle\widehat{\Lambda}^{+}_{s}(\bm{\alpha},\bm{\beta})$
$\displaystyle=\left(\frac{2}{s-s_{0}}\right)^{M}\\!\\!\\!\\!\int\frac{d^{2M}\bm{\zeta}}{\pi^{M}}\exp\left[-\frac{2|\bm{\zeta}|^{2}}{s-s_{0}}\right]\widehat{\Lambda}^{+}_{s_{0}}(\bm{\alpha}+\bm{\zeta},\bm{\beta}+\bm{\zeta}^{*}),$

which is easily verified with the help of

$\widehat{T}_{j}(0,-s)=\frac{1}{\pi}\int d^{2}\gamma
e^{-\tfrac{1}{2}s|\gamma|^{2}}\widehat{D}_{j}(\gamma)\\\ $ (103)

from [Cahill and Glauber(1969a)], the easy to show identities

$\displaystyle\widehat{D}_{j}(\alpha)\widehat{D}_{j}(\gamma)$ $\displaystyle=$
$\displaystyle\widehat{D}_{j}(\alpha+\gamma)e^{\tfrac{1}{2}(\alpha\gamma^{*}-\alpha^{*}\gamma)}$
(104)
$\displaystyle\widehat{d}_{j}(\alpha,\beta)\widehat{d}_{j}(\alpha^{\prime},\beta^{\prime})$
$\displaystyle=$
$\displaystyle\widehat{d}_{j}(\alpha+\alpha^{\prime},\beta+\beta^{\prime})e^{\tfrac{1}{2}(\alpha\beta^{\prime}-\beta\alpha^{\prime})}\\!,\qquad$
(105)

and Gaussian integrals. The distribution transform

$\displaystyle W^{+}_{s}(\bm{\alpha},\bm{\beta})$
$\displaystyle=\left(\frac{2}{s_{0}-s}\right)^{M}\\!\\!\\!\\!\int\frac{d^{2M}\bm{\zeta}}{\pi^{M}}\exp\left[-\frac{2|\bm{\zeta}|^{2}}{s_{0}-s}\right]W^{+}_{s_{0}}(\bm{\alpha}+\bm{\zeta},\bm{\beta}+\bm{\zeta}^{*}).$

is readily found by equating two (95) expansions of $\widehat{\rho}$ which
have different $s$, and applying (6.1) to the one with higher $s$ ($=s_{0}$).

###  6.2 Non normally-ordered correlations in the positive-P representation

The ideas from Sec. 5 can be used to treat non normally-ordered correlations
in the positive-P representation, which – unlike the Glauber-Sudarshan P – is
applicable for all quantum states and systems. With the tools for the doubled
s-ordered phase space described in Sec. 6.1, derivation of the expressions for
mixed-time expectation values follows the same path as in Secs. 5.1 and 5.3,
but some care needs to be taken to incorporate the doubled phase-space. We
obtain that:

(1) To shift samples from a more to a less normally ordered doubled-phase
space representation, one adds the following noise:

$\alpha^{\prime}_{j}=\alpha_{j}+\sqrt{\frac{s_{0}-s}{2}}\ \zeta_{j}\ ;\
\beta^{\prime}_{j}=\beta_{j}+\sqrt{\frac{s_{0}-s}{2}}\ \zeta^{*}_{j}.$ (107)

Notably – the same noise for $\alpha$ and $\beta^{*}$, which was not obvious
_a priori_. The prefactor is $1$ for positive-P to “doubled”-Q, and
$1/\sqrt{2}$ for positive-P to doubled-Wigner.

(2) Anti-normal ordered mixed-time correlations are evaluated as

$\displaystyle\langle\widehat{a}_{p_{1}}(t_{1})\cdots\widehat{a}_{p_{{\mathcal{N}}}}(t_{{\mathcal{N}}})\widehat{a}^{\dagger}_{q_{1}}(s_{1})\cdots\widehat{a}^{\dagger}_{q_{{\mathcal{M}}}}(s_{{\mathcal{M}}})\rangle$
(108)
$\displaystyle=\langle\alpha^{\prime}_{p_{1}}(t_{1})\cdots\alpha^{\prime}_{p_{{\mathcal{N}}}}(t_{{\mathcal{N}}})\beta^{\prime}_{q_{1}}(s_{1})\cdots\beta^{\prime}_{q_{{\mathcal{M}}}}(s_{{\mathcal{M}}})\rangle_{\rm
stoch}.\qquad$

where $\alpha^{\prime}$ and $\beta^{\prime}$ are doubled-Q representation
samples, possibly created via (107) (with $s_{0}-s=2$) from initial positive-P
samples, and later evolved via the appropriate doubled-Q representation
evolution equations.

(3) For mixed ordering, the procedure in Sec. 5.4 follows with the same
structure, except that the correspondences are

$\widehat{a}_{j}\to\alpha_{j}\qquad;\qquad\widehat{a}^{\dagger}_{j}\to\beta_{j}.$
(109)

in the positive-P and

$\widehat{a}_{j}\to\alpha^{\prime}_{j}\qquad;\qquad\widehat{a}^{\dagger}_{j}\to\beta^{\prime}_{j}.$
(110)

in the doubled-Q. For example, the expression for the correlation
${\mathcal{A}}^{\prime\prime}$ (5.3) using samples starting in the positive-P
is

${\mathcal{A}}^{\prime\prime}=\langle\beta_{p_{1}}(t_{1})\alpha^{\prime}_{p_{2}}(t_{2})\beta^{\prime}_{q_{1}}(s_{1})\alpha_{q_{2}}(s_{2})\rangle_{\rm
stoch}$ (111)

(4) The correlation tallies of Sec. 5.5 apply without change to the doubled
phase space representations.

##  7 Test example: unconventional photon blockade

###  7.1 The system

Consider the two-site Bose-Hubbard Hamiltonian

$\displaystyle\widehat{H}$ $\displaystyle=$
$\displaystyle\sum_{j=1,2}\widehat{a}^{\dagger}_{j}\left[-\Delta+\frac{U}{2}\widehat{a}^{\dagger}_{j}\widehat{a}_{j}\right]\widehat{a}_{j}$
(112)
$\displaystyle+F\left[\widehat{a}^{\dagger}_{1}+\widehat{a}_{1}\right]-J\left[\widehat{a}^{\dagger}_{1}\widehat{a}_{2}+\widehat{a}^{\dagger}_{2}\widehat{a}_{1}\right]$

with standard annihilation operators $\widehat{a}_{j}$ for modes $j$ using
units $\hbar=m=1$. This describes two sites with (real) tunnelling $J$, local
on-site interaction constant $U$, detuning $\Delta$, and a coherent drive $F$
(real) only on the first site. There is also a decay process with rate
$\gamma$ that is treated by describing the evolution of the system with the
master equation

$\displaystyle\frac{\partial\widehat{\rho}}{\partial t}$ $\displaystyle=$
$\displaystyle-i\left[\widehat{H},\widehat{\rho}\right]+\frac{\gamma{\overline{N}}}{2}\sum_{j}\left[2\widehat{a}^{\dagger}_{j}\widehat{\rho}\widehat{a}_{j}-\widehat{a}_{j}\widehat{a}^{\dagger}_{j}\widehat{\rho}-\widehat{\rho}^{\dagger}\widehat{a}_{j}\widehat{a}_{j}\right]$
(113)
$\displaystyle+\frac{\gamma({\overline{N}}+1)}{2}\sum_{j}\left[2\widehat{a}_{j}\widehat{\rho}\widehat{a}^{\dagger}_{j}-\widehat{a}^{\dagger}_{j}\widehat{a}_{j}\widehat{\rho}-\widehat{\rho}\widehat{a}^{\dagger}_{j}\widehat{a}_{j}\right].$

To test the performance a bit more beyond the standard model, we have also
addded a thermal bath with mean occupation ${\overline{N}}$. Such a system can
be realised using e.g. micropillars [Goblot et al.(2019)Goblot, Rauer,
Vicentini, Le Boité, Galopin, Lemaître, Le Gratiet, Harouri, Sagnes, Ravets,
Ciuti, Amo, and Bloch, Schneider et al.(2016)Schneider, Winkler, Fraser, Kamp,
Yamamoto, Ostrovskaya, and Höfling] or transmon qubits [Schmidt and
Koch(2013), Houck et al.(2012)Houck, Türeci, and Koch]. The nontrivial feature
here is a two-boson destructive interference effect between photons injected
by the drive, and other photons that have been previously injected, tunnelled
to site 2 and then back, returning with a relative phase of $\pi$ [Bamba et
al.(2011)Bamba, Imamoğlu, Carusotto, and Ciuti]. The result of this is that in
the steady state only single photons can be present at the driven site 1,
providing possibly an avenue to create single-photon sources [Liew and
Savona(2010)]. The lack of double occupation is evidenced in a single-time two
body correlation function $g_{11}=\langle\widehat{a}_{1}^{\dagger
2}\widehat{a}_{1}^{2}\rangle/\langle\widehat{a}^{\dagger}_{1}\widehat{a}_{1}\rangle^{2}$,
which is very close to zero. However, for practical application, it is of
particular interest to find out how large a time mismatch between measured
photons can be accommodated without significantly increasing $g^{(2)}_{11}$
from this low level. If it is too short, then the single-photon source will
have a too short active time for practical applications. Therefore the
correlation function of particular interest is

$g_{11}(\tau)=\frac{\langle\widehat{a}^{\dagger}_{1}(t)\widehat{a}^{\dagger}_{1}(t+\tau)\widehat{a}_{1}(t+\tau)\widehat{a}_{1}(t)\rangle}{n_{1}^{2}}$
(114)

in the stationary state, with delay time $\tau$. The mean occupation is
$n_{1}=\langle\widehat{a}^{\dagger}_{1}\widehat{a}_{1}\rangle$.

The Ito stochastic evolution equations in the positive-P ($s=1$) and doubled-Q
($s=-1$) representations are

$\displaystyle\frac{d\alpha_{j}}{dt}$ $\displaystyle=$
$\displaystyle\left[-iU(\alpha_{j}\beta_{j}+s-1)-\frac{\gamma}{2}+\sqrt{-isU}\,\xi_{j}(t)\right]\alpha_{j}-iF_{j}$
(115) $\displaystyle-i\sum_{k}H^{\rm
sp}_{jk}\alpha_{k}+\sqrt{\gamma\left({\overline{N}}+\tfrac{1-s}{2}\right)}\,\eta_{j}(t)$
$\displaystyle\frac{d\beta_{j}}{dt}$ $\displaystyle=$
$\displaystyle\left[iU(\alpha_{j}\beta_{j}+s-1)-\frac{\gamma}{2}+\sqrt{isU}\,{\widetilde{\xi}}_{j}(t)\right]\beta_{j}+iF_{j}$
(116) $\displaystyle+i\sum_{k}H^{\rm
sp}_{jk}\beta_{k}+\sqrt{\gamma\left({\overline{N}}+\tfrac{1-s}{2}\right)}\,\eta^{*}_{j}(t).$

These were derived in [Deuar et al.(2021)Deuar, Ferrier, Matuszewski, Orso,
and Szymańska] for the positive-P. Here the matrix elements of the one-
particle Hamiltonian are $H^{\rm sp}_{jj}=-\Delta$, $H^{\rm
sp}_{12}=H_{21}=-J$, the drives are $F_{1}=F$, $F_{2}=0$, while $\xi_{j}$ and
${\widetilde{\xi}}_{j}$ are delta-correlated independent white real noises
with variance

$\langle\xi_{j}(t)\xi_{j^{\prime}}(t^{\prime})\rangle=\delta_{jj^{\prime}}\delta(t-t^{\prime})$
(117)

and $\eta_{j}$ are similarly correlated independent complex noises:

$\langle\eta_{j}(t)\eta^{*}_{j^{\prime}}(t^{\prime})\rangle=\delta_{jj^{\prime}}\delta(t-t^{\prime});\qquad\langle\eta_{j}(t)\eta_{j^{\prime}}(t^{\prime})\rangle=0.$
(118)

Appendix A.3 gives some detail on generation of the noise. In some cases,
calculations with these stochastic equations can already be faster than brute
force calculations directly with the density matrix $\widehat{\rho}$ in a
suitably truncated number state basis. For large systems, they are the only
tractable way to access full quantum mechanics. Appendix A gives details of an
integration algorithm that is robust to the multiplicative noise appearing in
(7.1), and was used for the simulations reported here and elsewhere [Deuar and
Pietraszewicz(2019), Deuar et al.(2021)Deuar, Ferrier, Matuszewski, Orso, and
Szymańska, Ross et al.(2021)Ross, Deuar, Shin, Thomas, Henson, Hodgman, and
Truscott, Świsłocki and Deuar(2016), Deuar(2016), Kheruntsyan et
al.(2012)Kheruntsyan, Jaskula, Deuar, Bonneau, Partridge, Ruaudel, Lopes,
Boiron, and Westbrook, Deuar et al.(2014)Deuar, Jaskula, Bonneau,
Krachmalnicoff, Boiron, Westbrook, and Kheruntsyan, Deuar et al.(2011)Deuar,
Chwedeńczuk, Trippenbach, and Ziń, Deuar et al.(2013)Deuar, Wasak, Ziń,
Chwedeńczuk, and Trippenbach, Ng et al.(2013)Ng, Sørensen, and Deuar,
Deuar(2009)].

###  7.2 Two-time photon-photon correlations

Consider first the strong antibunching parameters studied in [Bamba et
al.(2011)Bamba, Imamoğlu, Carusotto, and Ciuti] and [Deuar et al.(2021)Deuar,
Ferrier, Matuszewski, Orso, and Szymańska]: $U=0.0856$, $J=3$,
$\Delta=-0.275$, $\gamma=1$, $F=0.01$. We will study the correlations
primarily in the stationary state. To this end, a positive-P simulation is
initialised in the vacuum $\bm{\alpha}=\bm{\beta}=0$, and evolved up to
$t=30$. The stationary state is attained after $t\gtrsim 15$, and we calculate
multi time correlation functions from times $t_{0}=20$ to $t_{0}+\tau$ using
the available samples. (${\mathcal{S}}=2^{16}$ in all cases). Uncertainty in
predictions is calculated using sub-ensemble averaging, as explained in
Appendix A.4. These errors are shown as triple lines in all the plots.

The two-time photon-photon correlations between one photon measured at site 1
at time $t_{0}$ and the other at site $j$ after a delay time of $\tau$ are

$g_{1,j}(\tau)=\frac{\langle\widehat{a}^{\dagger}_{1}(t_{0})\widehat{a}^{\dagger}_{j}(t_{0}+\tau)\widehat{a}_{j}(t_{0}+\tau)\widehat{a}_{1}(t_{0})\rangle}{n_{1}(t_{0})n_{j}(t_{0}+\tau)},$
(119)

where $n_{j}(t)=\langle\widehat{a}^{\dagger}(t)\widehat{a}(t)\rangle={\rm
Re}\langle\alpha_{j}(t)\beta_{j}(t)\rangle_{\rm stoch}$ is the mean occupation
of site $j$. In the positive-P representation the correlation (119) is
calculated via

$g_{1,j}(\tau)=\frac{{\rm
Re}\left\langle\alpha_{1}(t_{0})\beta_{1}(t_{0})\alpha_{j}(t_{0}+\tau)\beta_{j}(t_{0}+\tau)\right\rangle_{\rm
stoch}}{n_{1}(t_{0})n_{j}(t_{0}+\tau)}.$ (120)

We can take the real parts above, because the imaginary parts must converge to
zero in the ${\mathcal{S}}\to\infty$ limit. This behaviour is verified in the
simulations.

Figure 1: Time duration of antibunching in the unconventional photon blockade
system, $F=0.01$. Top: clean system (${\overline{N}}=0$), strong antibunching.
Bottom: degradation with background thermal intensity ${\overline{N}}>0$.
Black dashed lines: direct solution of the master equation. Triple lines show
$1\sigma$ statistical uncertainty.

Fig. 1 (top) shows the multi-time local photon correlation (119) at site 1
(i.e. $j=1$) for the clean case (${\overline{N}}=0$), and verifies that the
result perfectly matches the exact brute force solution. The desired anti-
correlation dip is seen around $\tau=0$, along with characteristic
oscillations out to delay times of about $\tau=5$. The bottom panel shows how
the anti-correlation dip degrades when the system is linked to a particle
reservoir (growing ${\overline{N}}$). Notably the dip timescale does not
change appreciably as the minimum correlation rises. A 10% remnant correlation
which might still be acceptable for applications is found when
${\overline{N}}=10^{-8}=0.026n_{1}$. This sets a limit on how much background
photon reservoir is acceptable.

Figure 2: Breakdown of the Glauber-Sudarshan P representation: shown in pink;
Positive-P simulation: yellow; exact: black dashed. $F=0.01$.

###  7.3 Breakdown of the Glauber-Sudarshan P

This system also shows the breakdown of the single-phase-space P
representation very clearly. The evolution equation is (115) with $s=1$, but
it is a correct representation of the quantum FPE (12) only when the diffusion
matrix $D_{\mu\nu}$ in the FPE has all nonnegative eigenvalues
[Gardiner(1991)]. Here, the diffusion matrix for real, imaginary parts of
$\alpha_{j}=\alpha_{jr}+i\alpha_{ji}$ has elements $D_{jr,jr}=\tfrac{U}{2}{\rm
Im}(\alpha_{j}^{2})+\tfrac{\gamma{\overline{N}}}{2}$,
$D_{ji,ji}=-\tfrac{U}{2}{\rm
Im}(\alpha_{j}^{2})+\tfrac{\gamma{\overline{N}}}{2}$,
$D_{jr,ji}=D_{ji,jr}=-\tfrac{U}{2}{\rm Re}(\alpha_{j}^{2})$. with eigenvalues
$\lambda_{j,\pm}=\tfrac{\gamma{\overline{N}}}{2}\pm\tfrac{U}{2}|\alpha_{j}|^{2}$.
These only become non-negative once $\gamma{\overline{N}}>U|\alpha_{j}|^{2}$,
i.e. $\gamma{\overline{N}}\gtrsim Un_{j}$. The question then is: for what
parameters does the evolution remain well described? When ${\overline{N}}=0$
The antibunched mode has mean occupation $n_{1}=3.87\times 10^{-7}$, naively
suggesting ${\overline{N}}\sim 5\times 10^{-8}=1.5Un_{1}/\gamma$ to already be
a value for which the description is good. However, we can see in Fig. 2 that
it does not give correct results at all. This is because the problem lies in
the nr. 2 mode with $n_{2}=1.07\times 10^{-5}$. Taking a far larger
${\overline{N}}=2\times 10^{-5}$ (in which case $\gamma{\overline{N}}\sim
10Un_{1},7Un_{2}$) gives almost correct results with the Glauber-Sudarshan P
(though still not fully), but of course antibunching in mode 1 is long gone
for such a relatively high thermal noise level.

This is an indication that skimping on full quantum effects by trying
approximate semiclassical methods is not a good strategy for this kind of
system.

###  7.4 Differently ordered correlations

Figure 3: Estimation in positive-P and doubled-Q simulations: occupation of
mode 1, when $F=3$, after switching from the positive-P to the doubled-Q
representation at $t=t_{0}=20$.

To test how the prescriptions developed in Secs. 4.1 and 5 work, we use a
different driving, $F=3$, which generates larger occupations (in the
stationary state $n_{1}\approx 0.043$ and $n_{2}\approx 0.98$) and as a result
more interesting anti-normal and mixed-order correlations. The doubled-Q
simulations are also too noisy to get useful predictions for the parameters in
Fig. 1 because Q distributions have a width of ${\mathcal{O}}(1)$ even in
vacuum. This is a certain limitation to the Q distribution approach. The
difference in sampling accuracy can be nicely seen in Fig. 3 which shows the
mean and uncertainty of $n_{1}$ during the simulation used to generate Figs.
4-5. The samples are switched from positive-P to doubled-Q at $t=t_{0}=20$.

Figure 4: Anti-normal ordered moments. Strongly pumped case ($F_{1}=3$), mode
2. Shown are ${\mathcal{G}}_{a}$ and the real and imaginary parts of
${\mathcal{G}}_{b}$ as per (121). Exact results: dashed black.  Figure 5:
Moments that are neither normal nor anti-normal. $F=3$ like in Fig. 4. Shown
are ${\mathcal{G}}_{c}$ and the real and imaginary parts of
${\mathcal{G}}_{d}$ as per (7.4). Exact results: dashed black lines.

Arguably the primary practical interest lies in cases with four factors and
two times – corresponding roughly to either detection of single particles at
two times, or creation and destruction of pairs. Fig. 4 shows the simulation
of some anti-normally ordered correlations of this kind (not normalised),
evaluated using the doubled-Q representation:

$\displaystyle{\mathcal{G}}_{a}$ $\displaystyle=$
$\displaystyle\langle\widehat{a}_{2}(t_{0})\widehat{a}_{2}(t_{0}+\tau)\widehat{a}^{\dagger}_{2}(t_{0}+\tau)\widehat{a}^{\dagger}_{2}(t_{0})\rangle$
(121) $\displaystyle=$ $\displaystyle{\rm
Re}\langle\alpha^{\prime}_{2}(t_{0})\alpha^{\prime}_{2}(t_{0}+\tau)\beta^{\prime}_{2}(t_{0}+\tau)\beta^{\prime}_{2}(t_{0})\rangle_{\rm
stoch}$ $\displaystyle{\mathcal{G}}_{b}$ $\displaystyle=$
$\displaystyle\langle\widehat{a}_{2}(t_{0}+\tau)^{2}[\widehat{a}^{\dagger}_{2}(t_{0})]^{2}\rangle=\langle\alpha^{\prime}_{2}(t_{0}+\tau)^{2}\beta^{\prime}_{2}(t_{0})^{2}\rangle_{\rm
stoch}$

The latter ${\mathcal{G}}_{b}$ is an anomalous correlation that has both real
and imaginary components. Similarly, Fig. 5 shows simulations of a number of
mixed-order correlations which absolutely require a swapping from positive-P
to doubled-Q representation using the procedure of Sec. 5.4:

$\displaystyle{\mathcal{G}}_{c}$ $\displaystyle=$
$\displaystyle\langle\widehat{a}_{2}(t_{0}+\tau)\widehat{a}^{\dagger}_{2}(t_{0}+\tau)\widehat{a}^{\dagger}_{2}(t_{0})\widehat{a}_{2}(t_{0})\rangle$
$\displaystyle=$ $\displaystyle{\rm
Re}\langle\alpha^{\prime}_{2}(t_{0}+\tau)\beta^{\prime}_{2}(t_{0}+\tau)\beta^{\prime}_{2}(t_{0})\alpha_{2}(t_{0})\rangle_{\rm
stoch}$ $\displaystyle{\mathcal{G}}_{d}$ $\displaystyle=$
$\displaystyle\langle\widehat{a}_{2}(t_{0})[\widehat{a}^{\dagger}_{2}(t_{0}+\tau)]^{2}\widehat{a}_{2}(t_{0})\rangle$
(122) $\displaystyle=$
$\displaystyle\langle\alpha^{\prime}_{2}(t_{0})\beta^{\prime}_{2}(t_{0}+\tau)^{2}\alpha_{2}(t_{0})\rangle_{\rm
stoch}$

In all cases, primed variables are evaluated in the doubled-Q representation,
un-primed in the positive-P.

Notably, in both figures the stochastic simulations perfectly and very
accurately agree with brute force calculations using the density matrix. This
is despite the regime being one which is very poorly treated by approximate
semiclassical methods (occupations are ${\mathcal{O}}(1)$ or smaller). This is
strong evidence that the intuitive expressions and approach laid out in the
previous sections is appropriate.

###  7.5 Correlation dynamics in the stationary state

As a larger-sized example, consider the same Hamiltonian and master equation
as (112) and (113) (take ${\overline{N}}=0$), but with a longer chain of 32
sites, which is already beyond or at the limit of the capabilities of
alternative methods such as brute force or corner space renormalisation
[Finazzi et al.(2015)Finazzi, Le Boité, Storme, Baksic, and Ciuti, Casteels et
al.(2016)Casteels, Rota, Storme, and Ciuti]. The tunnelling term in (112) now
becomes
$-J\sum_{j=1}^{31}\left[\widehat{a}^{\dagger}_{j}\widehat{a}_{j+1}+\widehat{a}^{\dagger}_{j+1}\widehat{a}_{j}\right],$
and a driving of $F=3$ remains only at the $j=1$ site. While the stationary
state does not show any time-dependent change in expectation values, we should
expect to still see signatures of transport in its multi-time correlations as
a function of distance and delay time if the method is good.

Figure 6: Spreading of correlations in a 32 site chain, in the stationary
state. $F=3$. The correlation is given by (119).

Fig. 6 shows such a calculation in the 32 site system, plotting the normalised
photon-photon correlation $g_{1,j}(\tau)$ with spatial separation $j-1$ and
delay time $\tau$, as defined in (119). The correlations are evaluated from
$t_{0}=20$. Despite this being the stationary state, a very clear anti-
correlation signal can be seen moving steadily with delay time. Its speed is
approximately $2J$, twice the tunnelling rate. Correlation waves often travel
at twice the characteristic speed for single particles [Cheneau et
al.(2012)Cheneau, Barmettler, Poletti, Endres, Schauß, Fukuhara, Gross, Bloch,
Kollath, and Kuhr, Läuchli and Kollath(2008), Pietraszewicz et
al.(2019)Pietraszewicz, Stobińska, and Deuar], so this is not unexpected.
However, the speed is clearly not related to the superfluid speed of sound,
here $\sim\sqrt{Un_{j}}$, which is far lower.

##  8 Summary

The known framework for evaluation of multi-time observables from the work of
Gardiner in the P representation [Gardiner(1991)] has been extended here to
include firstly: the much more widely applicable positive-P representation
(44); Secondly the Q (77), doubled-Q and other single and double phase space
s-ordered representations; Thirdly: other orderings such as anti-normal and
mixed-ordered observable products (Sec. 5, and especially (5.3) and the
algorithm of Sec. 5.4). These results allow the evaluation of a very wide
range of quantum multi-time observables in bosonic systems (classified in
Tables. 1-2) in a way that contains the full quantum mechanics, and is
scalable to large systems. Systems with localised interactions require
computational effort proportional to $M$ or $M\log M$ with the number of modes
or sites, $M$.

While out-of-time correlations (OTOCs) are not directly accessible through the
mechanism outlined here, time reversal schemes like those used in experiment
[Gärttner et al.(2017)Gärttner, Bohnet, Safavi-Naini, Wall, Bollinger, and
Rey] or theory [Dowling et al.(2005)Dowling, Drummond, Davis, and Deuar,
Swingle et al.(2016)Swingle, Bentsen, Schleier-Smith, and Hayden, Bohrdt et
al.(2017)Bohrdt, Mendl, Endres, and Knap] could be attempted by changing the
signs of constants in the equation of motion. When combined with the
techniques developed here, this would give access to a wide range of
information about quantum chaos [Bohrdt et al.(2017)Bohrdt, Mendl, Endres, and
Knap, Maldacena et al.(2016)Maldacena, Shenker, and Stanford, Foini and
Kurchan(2019)], many body localisation [Fan et al.(2017)Fan, Zhang, Shen, and
Zhai, He and Lu(2017)], and quantum phase transitions [Shen et al.(2017)Shen,
Zhang, Fan, and Zhai] in larger systems than were accessible to date [Bohrdt
et al.(2017)Bohrdt, Mendl, Endres, and Knap, Gärttner et al.(2017)Gärttner,
Bohnet, Safavi-Naini, Wall, Bollinger, and Rey].

Along the way, a number of additional results were obtained regarding
conversion formulae between different orderings in the doubled-phase space
representations ((6.1) and (6.1)) and their stochastic samples (107). A clear
case of the breakdown of the Glauber-Sudarshan P representation was seen in
Fig. 2. Also, in Sec. 7 some results regarding the correlation functions
(Figs. 4 and 5), susceptibility to background thermal density (Fig. 1) and
signal speed (Fig. 6) in the unconventional photon blockade system were
obtained.

Finally, Appendix A details a convenient algorithm for simulation of phase-
space stochastic equations, an item that has been hard to find in the
literature in the past.

I am grateful to Michał Matuszewski, Marzena Szymańska, and Alex Ferrier for
discussions on topics leading to this paper, and to Wouter Vestraelen for
discussions of numerics. This research was supported by the National Science
Centre (Poland) grant No. 2018/31/B/ST2/01871.

##  Appendix A Numerical techniques for stochastic simulation of phase-space
trajectories

This appendix details how the stochastic equations for the examples in Sec. 7
were integrated, and describes a robust and general algorithm that is
particularly applicable to stochastic trajectories generated by phase-space
descriptions of quantum systems. The bulk of the method is based on the semi-
implicit midstep algorithm described in [Drummond and Mortimer(1991), Werner
and Drummond(1997)], adding some modified propagators as per [Deuar and
Pietraszewicz(2019)] which allow for more efficient treatment of exponential
and dominant deterministic processes. Some practical elements that are usually
skipped over in the literature are also pointed out. Example textbooks for the
broad topic of stochastic integration include [Kloeden and Platen(1992), Ikeda
and Watanabe(1988)].

###  A.1 Integration algorithm

Let us introduce the following notation. The Ito stochastic equations for a
set of variables $\vec{v}$ with elements $v_{j}$ indexed by $j$ will be
written

$\frac{dv_{j}}{dt}={\mathcal{D}}_{j}(\vec{v},t)=A_{j}(\vec{v},t)+X_{j}(\vec{v},t)$
(123)

where $X_{j}$ are noise terms with zero mean $\langle X_{j}\rangle_{\rm
stoch}=0$, whereas $A_{j}$ are deterministic “drift terms” that contain no
explicit noise contribution, apart from the statistical spread of the
variables $v_{j}$ themselves. At the level of the Ito equations (123) we also
require that the noise is not correlated in time $\langle
X_{j}(t)X_{j^{\prime}}(t^{\prime})\rangle_{\rm
stoch}\propto\delta(t-t^{\prime})$, i.e. proportional to Wiener increments.
For example, in (7.1):

$\displaystyle X_{\alpha_{j}}$ $\displaystyle=$
$\displaystyle\alpha_{j}\sqrt{-isU}\,\xi_{j}+\sqrt{\gamma\left({\overline{N}}+\frac{1-s}{2}\right)}\,\eta_{j},\qquad$
(124) $\displaystyle X_{\beta_{j}}$ $\displaystyle=$
$\displaystyle\beta_{j}\sqrt{isU}\,{\widetilde{\xi}}_{j}+\sqrt{\gamma\left({\overline{N}}+\frac{1-s}{2}\right)}\,\eta^{*}_{j}.$
(125)

For various reasons, the standard mathematical environments and numerical
packages tend not to perform well when integrating stochastic equations444A
notable exception is the XMDS package [Dennis et al.(2013)Dennis, Hope, and
Johnsson]., especially ones that contain variable-dependent diffusion
coefficients such as in (A.1) or the prototypical Kubo oscillator [Kubo et
al.(1985)Kubo, Toda, and Hashitsume] with multiplicative noise. Phase space
methods that treat the full quantum dynamics also typically contain such
terms. Once one leaves the simplest Ito-Euler algorithm which requires very
small step sizes, two of the recurring issues are that advanced integration
algorithms usually assume sufficiently continuous derivatives (completely
violated by white noise) or try to adapt the timestep in real time (which can
introduce systematic errors when noise is involved). A third difficulty is the
need for autocorrelation corrections (seen Sec. A.2) which become much more
difficult to derive as the complexity of the algorithm grows.

The semi-implicit midpoint approach [Drummond(1983), Smith and Gardiner(1989)]
has been shown to be robust against spurious numerical instability as
carefully analysed in [Drummond and Mortimer(1991), Werner and
Drummond(1997)], but remains simple enough for any autocorrelation corrections
to be readily calculated.

Suppose we are to advance time by $\Delta t$, starting from variables
$v_{j}^{(0)}$ at time $t$ to $v^{\prime}_{j}$ at $t+\Delta t$. Define a
constituent substep or propagator that advances times by $\tau$ and variables
to $v_{j}^{\rm
step}\left[\vec{v}_{j}^{(0)},{\mathcal{D}}(\vec{v}^{\,(0)},t),\tau\right]$
which depends in order on the starting variables, the form of the derivative,
and the timestep $\tau$. The simplest choice for this propagator is the Ito-
Euler form

$v_{j}^{\rm
step}\left[v_{j}^{(0)},{\mathcal{D}}(\vec{v}^{\,(0)},t),\tau\right]=v_{j}^{(0)}+{\mathcal{D}}(\vec{v}^{\,(0)},t)\,\tau$
(126)

which just uses the values at the beginning of the step for evaluating the
derivative, as per the definition of Ito stochastic calculus.

The semi-implicit midpoint algorithm can be viewed as a container procedure
that calculates the final variables $v^{\prime}_{j}$ via a series of
iterations [Werner and Drummond(1997)]

$\displaystyle v^{(1)}_{j}$ $\displaystyle=$ $\displaystyle v^{\rm
step}_{j}\left[v_{j}^{(0)},{\mathcal{D}}(\vec{v}^{\,(0)},t),\Delta
t/2\right],$ $\displaystyle v^{(2)}_{j}$ $\displaystyle=$ $\displaystyle
v^{\rm step}_{j}\left[v_{j}^{(0)},{\mathcal{D}}(\vec{v}^{\,(1)},t),\Delta
t/2\right],$ (127) $\displaystyle\cdots$ $\displaystyle v^{(n)}_{j}$
$\displaystyle=$ $\displaystyle v^{\rm
step}_{j}\left[v_{j}^{(0)},{\mathcal{D}}(\vec{v}^{\,(n-1)},t),\Delta
t/2\right],$ $\displaystyle v^{\prime}_{j}(t+\Delta t)$ $\displaystyle=$
$\displaystyle v^{\rm
step}_{j}\left[v_{j}^{(0)},{\mathcal{D}}(\vec{v}^{\,(n)},t+\Delta t/2),\Delta
t\right].\qquad$ (128)

The numbered iterations always start from the initial point, but iteratively
find a self-consistent estimate of the midpoint value of the derivatives. The
last iteration (128) uses this midpoint value and midstep time to make the
final advancement of the variables. Usually, a single midpoint iteration (i.e.
$n=1$) is sufficient. When using the basic Ito-Euler propagator (126) at time
argument $t+\Delta t/2$ in (128), the final step is accurate to
${\mathcal{O}}(\Delta t^{2})$ in the deterministic parts $A_{j}$, and the
step’s mean to order ${\mathcal{O}}(\Delta t)$ in the noise terms $X_{j}$
[Drummond and Mortimer(1991)], provided the autocorrelation correction is
incorporated into $v^{\rm step}_{j}$, as explained in Sec. A.2. This is so-
called “weak” stochastic convergence to ${\mathcal{O}}(\Delta t)$ in the
terminology of [Drummond and Mortimer(1991)] because the accuracy for single
trajectories is ${\mathcal{O}}(\sqrt{\Delta t})$. “Strong” stochastic
convergence to ${\mathcal{O}}(\Delta t)$ for single trajectories requires a
more involved and more inconvenient autocorrelation correction than that
presented in Sec. A.2. Nevertheless, the “weak” level of accuracy in the noise
term is actually quite high already. For example, even 4th or 9th order Runge-
Kutta implementations made available are generally only accurate to
${\mathcal{O}}(\sqrt{\Delta t})$ in this regard [Dennis et al.(2013)Dennis,
Hope, and Johnsson].

Moreover, in practice, one can often easily but substantially improve on the
Euler step (126). For example, if $\Delta$, $\gamma$, or $U$ are large in
(7.1), the evolution is primarily exponential. Another typical case is in
continuum models when the kinetic energy or external potential energy
contributions are dominant and trivial to integrate [Deuar and
Pietraszewicz(2019)]. A pragmatic and flexible propagator choice is to
separate the evolution into parts that will be treated exponentially (E),
strictly linearly (L) and the rest (R):

${\mathcal{D}}_{j}(\vec{v})={\mathcal{D}}_{j}^{E}(\vec{v})\,v_{j}+{\mathcal{D}}_{j}^{R}(\vec{v})+{\mathcal{D}}_{j}^{L}(\vec{v}),$
(129)

and solve the equation for the E and R parts

$\frac{dv_{j}}{d\tau}={\mathcal{D}}_{j}^{E}\,v_{j}+{\mathcal{D}}_{j}^{R}$
(130)

as if the coefficients were constant [Deuar and Pietraszewicz(2019)]. This
gives the full propagator

$\displaystyle v_{j}^{\rm
step}\left[v_{j}^{(0)},{\mathcal{D}}(\vec{v}^{\,(0)},t),\tau\right]=$ (131)
$\displaystyle
v_{j}^{(0)}\,e^{\tau\,{\mathcal{D}}_{j}^{E}(\vec{v}^{\,(0)},t)}+\left(e^{\tau\,{\mathcal{D}}_{j}^{E}(\vec{v}^{\,(0)},t)}-1\right)\
\frac{{\mathcal{D}}^{R}_{j}(\vec{v}^{\,(0)},t)}{{\mathcal{D}}^{E}_{j}(\vec{v}^{\,(0)},t)}$
$\displaystyle+{\mathcal{D}}^{L}_{j}(\vec{v}^{\,(0)},t)\ \tau.$

Then, the leading processes are often integrated exactly, with only higher
order corrections needed to be worked on by the midpoint algorithm. This is
particularly efficient when the stochastic terms are a perturbation on the
dominant exponential deterministic evolution. It is usually optimal to place
all non-exponential parts into the “R” part, but some special cases can
require avoidance of the nonlinear solution [Deuar and Pietraszewicz(2019)].
Both $A$ drift and $X$ noise parts can enter the coefficients as one chooses.

In the case of the calculations in this paper using (7.1),

$\displaystyle{\mathcal{D}}_{\alpha_{j}}^{E}$ $\displaystyle=$ $\displaystyle-
iU(\alpha_{j}\beta_{j}+s-1)-\frac{\gamma}{2}+\sqrt{-isU}\,\xi_{j}-iH^{\rm
sp}_{jj}+C^{E}_{\alpha_{j}},$ $\displaystyle{\mathcal{D}}_{\alpha_{j}}^{R}$
$\displaystyle=$ $\displaystyle-iF_{j}-i\sum_{k\neq j}H^{\rm
sp}_{jk}\alpha_{k}+\sqrt{\gamma\left({\overline{N}}+\frac{1-s}{2}\right)}\,\eta_{j},$
$\displaystyle{\mathcal{D}}_{\beta_{j}}^{E}$ $\displaystyle=$ $\displaystyle
iU(\alpha_{j}\beta_{j}+s-1)-\frac{\gamma}{2}+\sqrt{isU}\,{\widetilde{\xi}}_{j}+iH^{\rm
sp}_{jj}+C^{E}_{\beta_{j}},$ $\displaystyle{\mathcal{D}}_{\beta_{j}}^{R}$
$\displaystyle=$ $\displaystyle iF_{j}+i\sum_{k\neq j}H^{\rm
sp}_{jk}\beta_{k}+\sqrt{\gamma\left({\overline{N}}+\frac{1-s}{2}\right)}\,\eta^{*}_{j},$

and ${\mathcal{D}}^{L}_{j}=0$ were used for the propagator (131) and combined
with a single iteration ($n=1$) of the midpoint algorithm (A.1). Note the
presence of the autocorrelation corrections

$\displaystyle C^{E}_{\alpha_{j}}$ $\displaystyle=$
$\displaystyle\frac{isU}{2}$ (133) $\displaystyle C^{E}_{\beta_{j}}$
$\displaystyle=$ $\displaystyle-\frac{isU}{2}$ (134)

as explained in Sec. A.2.

Importantly – the form of the corrections (A.1) above arises when the noises
$\xi$, $\eta$ etc. appearing in the $X_{j}$ are calculated just once at the
beginning of the entire procedure for making the full $\Delta t$ step. Other
noise input can change the expressions (A.1).

The actual timestep $\Delta t$ that must be used is constrained by the need
for the coefficients ${\mathcal{D}}_{j}$ to change little over a single
timestep $\Delta t$. A good practical rule of thumb is to require

$\frac{\left|{\mathcal{D}}^{(\nu)}(\vec{v},t)\,\Delta
t\right|}{\left|v_{j}^{(0)}\right|}\lesssim 0.1,$ (135)

for each term ${\mathcal{D}}^{(\nu)}$ in
${\mathcal{D}}=\sum_{\nu}{\mathcal{D}}^{(\nu)}$. One caveat is that one should
take care not to estimate ${\mathcal{D}}(\vec{v}^{\,(0)},t)$ on the run using
noise or even variable values from individual trajectories. Doing so can, and
often does, introduce autocorrelations between noise values and time step
lengths, leading to systematic errors. For treating noise terms, it is useful
to use the RMS estimate
$X_{j}\approx\sqrt{\langle|X_{j}|^{2}\rangle}\propto\sqrt{\Delta t}$ based on
expected mean variable values $\langle v_{j}\rangle$ for timestep estimates.

Finally, for completeness and general applicability beyond the examples of
Sec. 7, it can be advantageous to treat widely differing processes using a
split-step approach [Werner and Drummond(1997)]. In particular, the treatment
of kinetic energy and tunnelling is much bettered for many-mode systems. A
split step allows to take advantage of the fact that kinetic energy and
tunnelling are diagonal in k-space and can be done exactly there. k-space is
reached using a discrete Fast Fourier Transform that has computational cost of
only $M\log M$ for an M-site system [Frigo and Johnson(2005)]. The split-step
is applied in three parts:
(i) First evolve terms diagonal in k-space for $\Delta t/2$;
(ii) then change variables to x-space and evolve the remaining terms there for
$\Delta t$;
(iii) finally return to k-space for another k-space evolution of $\Delta t/2$.

For cases like (7.1) in which the main mean-field evolution has a Gross-
Pitaevskii form, the split step method is symplectic (conserves energy) and
has been shown to have an accuracy of ${\mathcal{O}}(\Delta t)^{2}$ for the
long time solution [Javanainen and Ruostekoski(2006)] provided that the latest
available copy of the field is input as $\vec{v}^{\,(0)}$ after each Fourier
transform [Feit et al.(1982)Feit, Fleck, and Steiger]. Within each k- or
x-space split step, the midpoint method (A.1) can employed if the coefficients
${\mathcal{D}}$ of the substep are variable-dependent, or otherwise if the
coefficients are constant just a plain propagator substep (131) is done. Thus,
one has the split-step algorithm as a container for midpoint iterations, which
themselves are containers for the base propagators. Long-range interactions
(eg. dipole-dipole) can also often be treated efficiently through Fourier
transforms [Wüster et al.(2017)Wüster, Corney, Rost, and Deuar].

Sumarizing, the algorithm used for the simulations of Sec. 7 consists of
(A.1), (128), (131)–(A.1).

###  A.2 Timestep autocorrelation correction

Realisations of Ito stochastic equations (123) must at the least produce the
required average and noise variance to lowest order. That is, writing

$\Delta v_{j}=v^{\prime}_{j}(t+\Delta t)-v_{j}^{(0)},$ (136)

one needs

$\displaystyle\langle\Delta v_{j}\rangle_{\rm stoch}$ $\displaystyle=$
$\displaystyle\left\langle A_{j}(\vec{v}^{\,(0)})\Delta t\right\rangle_{\rm
stoch}+{\mathcal{O}}(\Delta t)^{3/2}$ $\displaystyle\\!\\!\\!\langle\Delta
v_{j}\Delta v_{k}\rangle_{\rm stoch}$ $\displaystyle=$
$\displaystyle\left\langle X_{j}(\vec{v}^{\,(0)})X_{k}(\vec{v}^{\,(0)})\Delta
t^{2}\right\rangle_{\rm stoch}\\!\\!\\!\\!\\!\\!+{\mathcal{O}}(\Delta
t)^{3/2}.$

In orders of $\Delta t$, the drift is $A_{j}\Delta t\sim{\mathcal{O}}(\Delta
t)$ and the noise $X_{j}\Delta t\sim{\mathcal{O}}(\sqrt{\Delta t})$. The time-
dependence of coefficients $A$ and $X$ does not enter at this order of $\Delta
t$. A simple Ito-Euler timestep $v_{j}^{(0)}\to
v_{j}^{(0)}+A_{j}(\vec{v}^{\,(0)})\,\Delta t+X_{j}(\vec{v}^{\,(0)})\,\Delta t$
meets both conditions (A.2).

However, the midstep algorithm with an Ito-Euler propagator (126) can
introduce an extra correlation of order $\Delta t$ by virtue of using the same
noise for both the halfstep of $\Delta t/2$ and then for the final full step.
The first halfstep advance using (126) produces

$\displaystyle\Delta v_{j}^{(1)}$ $\displaystyle=$ $\displaystyle
v_{j}^{(1)}-v_{j}^{(0)}=A_{j}(\vec{v}^{\,(0)})\frac{\Delta
t}{2}+X_{j}(\vec{v}^{\,(0)})\,\frac{\Delta t}{2},$ (138) $\displaystyle\sim$
$\displaystyle{\mathcal{O}}(\sqrt{\Delta t})$

Note that the $X_{j}$ themselves are ${\mathcal{O}}(\frac{1}{\sqrt{\Delta
t}})$, which is important below. The next step in (A.1) involves coefficients
$A_{j}(\vec{v}^{\,(1)})$ and $X_{j}(\vec{v}^{\,(1)})$. Taylor expanding around
$\vec{v}^{\,(0)}$,

$X_{j}(\vec{v}^{\,(1)})=X_{j}(\vec{v}^{\,(0)})+\sum_{k}\frac{dX_{j}}{dv_{k}}(\vec{v}^{\,(0)})\Delta
v_{k}^{(1)}+{\mathcal{O}}(\sqrt{\Delta t}),$ (139)

and analogously for $A_{j}$. Using (138) and discarding terms of subleading
order ${\mathcal{O}}(\sqrt{\Delta t})$ in the new coefficients one finds

$\displaystyle X_{j}(\vec{v}^{\,(1)})$ $\displaystyle=$ $\displaystyle
X_{j}(\vec{v}^{\,(0)})$
$\displaystyle+\sum_{k}\frac{dX_{j}}{dv_{k}}(\vec{v}^{\,(0)})\,X_{j}(\vec{v}^{\,(0)})\frac{\Delta
t}{2}+{\mathcal{O}}(\sqrt{\Delta t}),$ $\displaystyle A_{j}(\vec{v}^{\,(1)})$
$\displaystyle=$ $\displaystyle
A_{j}(\vec{v}^{\,(0)})+{\mathcal{O}}(\sqrt{\Delta t}).$ (140)

Notice the appearance of a term in $X_{j}(\vec{v}^{\,(1)})$ of the same order
as the original coefficients $X_{j}(\vec{v}^{\,(0)})$. Assuming first just one
$n=1$ iteration, the final step (128) becomes

$\displaystyle v^{\prime}$ $\displaystyle=$ $\displaystyle
v^{(0)}_{j}+A_{j}(\vec{v}^{\,(0)})\Delta t+X_{j}(\vec{v}^{\,(0)})\,\Delta t$
$\displaystyle+\sum_{k}\frac{dX_{j}}{dv_{k}}(\vec{v}^{\,(0)})\,X_{j}(\vec{v}^{\,(0)})\frac{(\Delta
t)^{2}}{2}+{\mathcal{O}}(\Delta t)^{3/2}.$

It turns out that iterations of $n>1$ do not affect the expression (A.2) at
${\mathcal{O}}(\Delta t)$. The new term in (A.2) breaks the equivalence of the
algorithm (A.2) because it is of the same order as the drift $A_{j}\Delta t$.

To restore equivalence, one can add a correction $C_{j}$ to the drift used in
the algorithm as per

$A_{j}\to A_{j}+C_{j}.$ (142)

in either the “strongly convergent” form

$C_{j}^{\rm strong}=-\frac{\Delta
t}{2}\sum_{k}\frac{dX_{j}}{dv_{k}}(\vec{v})\,X_{j}(\vec{v})$ (143)

or the “weakly convergent” averaged form $C_{j}=\left\langle C_{j}^{\rm
strong}\right\rangle_{\rm stoch}$. For the case of Ito-Euler propagator in the
midstep algorithm,

$C_{j}=C_{j}^{\rm Strat}=-\frac{\Delta
t}{2}\left\langle\sum_{k}\frac{dX_{j}}{dv_{k}}(\vec{v})\,X_{j}(\vec{v})\right\rangle_{\rm
stoch}.$ (144)

The strong forms (143) have different values for each stochastic realization
[Drummond and Mortimer(1991)] and obtain ${\mathcal{O}}(\Delta t)$ timestep
accuracy for the noise terms, while the weak forms use a pre-calculated mean
(144). The weak variant is more commonly used in practice and sufficient to
obtain overall stable and accurate integration. It is what is applied in (A.1)
and the calculations of Sec. 7.

The weak autocorrelation correction $C_{j}$ has often been dubbed a
“Stratonovich correction” because the form (144) is identical to the
correction used to move between Ito and Stratonovich stochastic calculus.
However, this is a misnomer because the match is merely a coincidence that
occurs when an Ito-Euler propagator is used with the midpoint algorithm. For
the case of the partly exponential propagator (131), the same procedure as
above leads to the following, different, form of the autocorrelation
corrections. When $n\geq 1$ in the midstep algorithm:

$C_{j}=C_{j}^{\rm Strat}+\frac{\Delta t}{2}\left\langle
X^{E}_{j}(\vec{v})X^{L}_{j}(\vec{v})\right\rangle_{\rm stoch},$ (145)

or when one uses no midstepping ($n=0$):

$C_{j}=-\frac{\Delta t}{2}\left\langle
X^{E}_{j}(\vec{v})X^{E}_{j}(\vec{v})v_{j}+X^{E}_{j}(\vec{v})X^{R}_{j}(\vec{v})\right\rangle_{\rm
stoch}.$ (146)

As it happens, in the case of the equations (7.1), the corrections are the
same whether using an Ito-Euler propagator ((A.1)) or the partially
exponential one (145)-(146), but that is not always the case.

###  A.3 Noise implementation

Since the underlying stochastic equations are defined in the infinitesimal
limit $\Delta t\to 0$, in principle any implementation of the noise that has
zero mean and satisfies the variance conditions (117) and (118) in the $\Delta
t\to 0$ limit is correct. By the central limit theorem, after several
timesteps the effective distribution of the sum of the noises will always
converge to a Gaussian one. In practice, however, it is usually desirable to
use explicitly Gaussian distributed noises in the implementation. Doing so
already accurately depicts the limiting distribution for each $\Delta t$ step
and avoids a potential need to reduce timestep further to capture the right
noise distribution. For real noises $\xi$, one generates fresh Gaussian random
variables of variance $1/\Delta t$ for each $j$ at the beginning of each
timestep. Generating new ones for each time step ensures the
$\delta(t-t^{\prime})$ limit as $\Delta t\to 0$. The Box-Muller algorithm [Box
and Muller(1958)] is a simple way to obtain two such independent Gaussian
variables from two uniformly distributed noises, whereas a variety of more
efficient though less transparent algorithms are also known [Bell(1968),
Brent(1974)]. In our case Box-Muller was used, while the underlying uniformly
distributed noise was generated using the SFMT fast Mersenne twister method
[Saito and Matsumoto(2008)]
555http://www.math.sci.hiroshima-u.ac.jp/m-mat/MT/SFMT/http://www.math.sci.hiroshima-u.ac.jp/m-mat/MT/SFMT/.
Complex noise $\eta$ is simply constructed as
$\eta=(\xi_{r}+i\xi_{i})/\sqrt{2}$ using two real Gaussian noises $\xi_{i,r}$.

When using efficient random number generation such as [Saito and
Matsumoto(2008)] it turns out that the creation of Gaussian noises from
uniform ones is the most computationally costly step. The Brent algorithm
[Brent(1974)] alleviates this appreciably, while a simple alternative route
applicable at least for stochastic simulations – due to their central limit
properties – is to produce binomially distributed noise instead. One adds
$n_{B}$ uniform noises $r_{n}$ on $[0,1]$ as per
$\xi\approx\sqrt{\frac{12}{n_{B}\Delta
t}}\sum_{n=1}^{n_{B}}(r_{n}-\tfrac{1}{2})$. In practice $n_{B}=4$ or $n_{B}=3$
is already sufficient.

###  A.4 Error estimation

Statistical uncertainty in quantities ${\mathcal{Q}}$ calculated from the
simulations can be robustly estimated via sub-ensemble averaging [Werner and
Drummond(1997)] even when the distribution of samples is unknown. The full
ensemble of, say ${\mathcal{S}}$, realisations is divided into a smaller
number $u$ of sub-ensembles. Auxiliary subensemble means
${\overline{{\mathcal{Q}}}}_{(i=1,\dots,u)}$ are calculated for each sub-
ensemble individually. Then by the central limit theorem, the $1\sigma$
uncertainty in the full ensemble prediction is

$\Delta{\mathcal{Q}}=\sqrt{\frac{{\rm
var}\left[{\overline{{\mathcal{Q}}}}_{(i)}\right]}{u-1}}.$ (147)

In practice $u\lesssim\sqrt{{\mathcal{S}}}$ is useful. In the simulations
reported in this paper, $u=32$ was used.

A separate issue is testing the time discretization accuracy in stochastic
equations. In a simple implementation where noise is generated sequentially,
changing timestep changes also the noise history, making comparison of single
trajectories with different $\Delta t$ useless for determination of accuracy.
Yet, comparing the ensemble means for different timesteps can be onerous in
large systems. A solution is to compare two runs of a single realisation with
identical underlying noise sources [Drummond and Mortimer(1991)]: One with
time step $\Delta t$ and noises $\xi_{j}^{(i)}$ generated sequentially for
time steps numbered $i=1,\dots$. The other with integration timestep $2\Delta
t$ but noises generated by adding pairs of noises from the first simulation:
$\xi_{j}^{(j)\prime}=\xi_{j}^{(2j-1)}+\xi_{j}^{(2j)}$. This allows separation
of the time discretization error from the stochastic randomness.

## References

  * [Agarwal(1969a)] G. S. Agarwal.  Phase-space analysis of time-correlation functions.  _Phys. Rev._ , 177:400–407, 1969a.  doi: https://doi.org/10.1103/PhysRev.177.400.
  * [Agarwal(1969b)] G. S. Agarwal.  Master equations in phase-space formulation of quantum optics.  _Phys. Rev._ , 178:2025–2035, 1969b.  doi: https://doi.org/10.1103/PhysRev.178.2025.
  * [Agarwal and Chaturvedi(1994)] G. S. Agarwal and S. Chaturvedi.  Scheme to measure the positive $P$ distribution.  _Phys. Rev. A_ , 49:R665–R667, 1994.  doi: https://doi.org/10.1103/PhysRevA.49.R665.
  * [Agarwal and Wolf(1970a)] G. S. Agarwal and E. Wolf.  Calculus for functions of noncommuting operators and general phase-space methods in quantum mechanics. II. quantum mechanics in phase space.  _Phys. Rev. D_ , 2:2187–2205, 1970a.  doi: https://doi.org/10.1103/PhysRevD.2.2187.
  * [Agarwal and Wolf(1970b)] G. S. Agarwal and E. Wolf.  Calculus for functions of noncommuting operators and general phase-space methods in quantum mechanics. III. a generalized Wick theorem and multitime mapping.  _Phys. Rev. D_ , 2:2206–2225, 1970b.  doi: https://doi.org/10.1103/PhysRevD.2.2206.
  * [Aimi and Imada(2007a)] Takeshi Aimi and Masatoshi Imada.  Gaussian-Basis Monte Carlo Method for Numerical Study on Ground States of Itinerant and Strongly Correlated Electron Systems.  _J. Phys. Soc. Jpn._ , 76:084709, 2007a.  doi: https://doi.org/10.1143/JPSJ.76.084709.
  * [Aimi and Imada(2007b)] Takeshi Aimi and Masatoshi Imada.  Does Simple Two-Dimensional Hubbard Model Account for High-$T_{c}$ Superconductivity in Copper Oxides?  _J. Phys. Soc. Jpn._ , 76:113708, 2007b.  doi: https://doi.org/10.1143/JPSJ.76.113708.
  * [Albarelli et al.(2018)Albarelli, Genoni, Paris, and Ferraro] Francesco Albarelli, Marco G. Genoni, Matteo G. A. Paris, and Alessandro Ferraro.  Resource theory of quantum non-Gaussianity and Wigner negativity.  _Phys. Rev. A_ , 98:052350, 2018.  doi: https://doi.org/10.1103/PhysRevA.98.052350.
  * [Atalaya et al.(2018)Atalaya, Hacohen-Gourgy, Martin, Siddiqi, and Korotkov] Juan Atalaya, Shay Hacohen-Gourgy, Leigh S. Martin, Irfan Siddiqi, and Alexander N. Korotkov.  Multitime correlators in continuous measurement of qubit observables.  _Phys. Rev. A_ , 97:020104, 2018.  doi: https://doi.org/10.1103/PhysRevA.97.020104.
  * [Bamba et al.(2011)Bamba, Imamoğlu, Carusotto, and Ciuti] Motoaki Bamba, Atac Imamoğlu, Iacopo Carusotto, and Cristiano Ciuti.  Origin of strong photon antibunching in weakly nonlinear photonic molecules.  _Phys. Rev. A_ , 83:021802, 2011.  doi: https://doi.org/10.1103/PhysRevA.83.021802.
  * [Barry and Drummond(2008)] D. W. Barry and P. D. Drummond.  Qubit phase space: SU$(n)$ coherent-state $P$ representations.  _Phys. Rev. A_ , 78:052108, 2008.  doi: https://doi.org/10.1103/PhysRevA.78.052108.
  * [Bell(1968)] J. R. Bell.  Algorithm 334, normal random deviates.  _Communications of the ACM_ , 11(7):498, 1968\.  doi: https://doi.org/10.1145/363397.363547.