Positive-P Phase Space Method

• The positive-P phase space method is a doubled coherent-state representation that enables a true positive stochastic sampling of bosonic quantum dynamics.
• It reformulates quantum evolution into a Fokker–Planck equation and corresponding Itô stochastic differential equations, allowing accurate estimation of normally ordered observables.
• Its applications span driven-dissipative quantum optics, many-body simulations, and Bell-state studies, overcoming Hilbert space truncation without approximation.

The positive-P phase space method is a phase-space representation for bosonic quantum systems in which the density operator is expanded over an overcomplete doubled coherent-state basis, so that quantum evolution is recast as a Fokker–Planck equation and, when appropriate, as Itô stochastic differential equations for independent complex variables. Its defining structural feature is the use of two independent complex amplitudes per mode, rather than a single amplitude and its complex conjugate. This doubled phase space permits a true positive distribution for any quantum state and supports exact stochastic sampling of quantum evolution, normally ordered observables, and many multi-time correlations without truncation of the Hilbert space or other approximations to the quantum state (Deuar, 2020, Wüster et al., 2017, Jen, 2012).

1. Formal definition and kernel structure

For an $M$-mode bosonic system, the representation is written as

$$\hat{\rho}=\int d^{4M}\bm{\lambda}\; P_+(\bm{\lambda})\,\hat{\Lambda}(\bm{\lambda}), \qquad \hat{\Lambda}(\bm{\lambda}) = \bigotimes_j \frac{|\alpha_j\rangle \langle \beta_j^*|}{\langle \beta_j^*|\alpha_j\rangle},$$

with $$\bm{\lambda}=\{\alpha_1,\dots,\alpha_M,\beta_1,\dots,\beta_M\}$$. In equivalent notation used elsewhere, the kernel is

$$\hat{\Lambda}(\vec{\alpha},\vec{\alpha}^{+})=\frac{|\vec{\alpha}\rangle\langle (\vec{\alpha}^{+})^{*}|}{\langle (\vec{\alpha}^{+})^{*}|\vec{\alpha}\rangle}.$$

The phase space is doubled because $\alpha_j$ and $\beta_j$ are independent complex variables and are not constrained to be complex conjugates trajectory by trajectory (Deuar, 2020, Rosales-Zárate et al., 2014, Shi et al., 11 Jan 2026).

This doubling is the crucial distinction from the ordinary Glauber–Sudarshan $P$-representation. The positive-P representation exists for all bosonic states, and $P_+(\bm{\lambda})$ can be chosen to be a true positive distribution. The basis is overcomplete and the representation is nonunique, which is operationally useful because different equivalent distributions can have different sampling properties (Deuar, 2020, Rosales-Zárate et al., 2014).

Normally ordered observables are obtained as stochastic moments. In the single-mode form emphasized in driven-dissipative cat-state simulations,

$$\langle (\hat{a}^\dagger)^m \hat{a}^n \rangle = \lim_{s\to\infty}\langle \beta^m \alpha^n \rangle_s,$$

and the full density matrix can in principle be reconstructed from stochastic averages of the kernel, although in practice the method is usually used for observables rather than full state reconstruction (Shi et al., 11 Jan 2026).

2. Operator correspondences and stochastic dynamics

The practical power of the method comes from derivative and derivative-free operator identities for the kernel. In a standard multimode form,

$$\hat{a}_j\hat{\Lambda}=\alpha_j\hat{\Lambda}, \qquad \hat{a}_j^\dagger\hat{\Lambda}=\left[\beta_j+\frac{\partial}{\partial\alpha_j}\right]\hat{\Lambda},$$

$$\hat{\rho}=\int d^{4M}\bm{\lambda}\; P_+(\bm{\lambda})\,\hat{\Lambda}(\bm{\lambda}), \qquad \hat{\Lambda}(\bm{\lambda}) = \bigotimes_j \frac{|\alpha_j\rangle \langle \beta_j^*|}{\langle \beta_j^*|\alpha_j\rangle},$$0

Equivalent single-mode identities are used throughout application papers, for example

$$\hat{\rho}=\int d^{4M}\bm{\lambda}\; P_+(\bm{\lambda})\,\hat{\Lambda}(\bm{\lambda}), \qquad \hat{\Lambda}(\bm{\lambda}) = \bigotimes_j \frac{|\alpha_j\rangle \langle \beta_j^*|}{\langle \beta_j^*|\alpha_j\rangle},$$1

together with the derivative identities on the opposite side of the kernel (Deuar, 2020, Shi et al., 11 Jan 2026).

These identities convert the von Neumann or Lindblad master equation into a Fokker–Planck equation. The latter is equivalent to Itô stochastic differential equations of the generic form

$$\hat{\rho}=\int d^{4M}\bm{\lambda}\; P_+(\bm{\lambda})\,\hat{\Lambda}(\bm{\lambda}), \qquad \hat{\Lambda}(\bm{\lambda}) = \bigotimes_j \frac{|\alpha_j\rangle \langle \beta_j^*|}{\langle \beta_j^*|\alpha_j\rangle},$$2

with independent real Gaussian noises satisfying

$$\hat{\rho}=\int d^{4M}\bm{\lambda}\; P_+(\bm{\lambda})\,\hat{\Lambda}(\bm{\lambda}), \qquad \hat{\Lambda}(\bm{\lambda}) = \bigotimes_j \frac{|\alpha_j\rangle \langle \beta_j^*|}{\langle \beta_j^*|\alpha_j\rangle},$$3

In this way, operator evolution is replaced by trajectory evolution, and quantum expectation values are replaced by stochastic averages over trajectories in doubled phase space (Wüster et al., 2017, Deuar, 2020).

A major generalization is the gauge-P representation, in which the density matrix is expanded as

$$\hat{\rho}=\int d^{4M}\bm{\lambda}\; P_+(\bm{\lambda})\,\hat{\Lambda}(\bm{\lambda}), \qquad \hat{\Lambda}(\bm{\lambda}) = \bigotimes_j \frac{|\alpha_j\rangle \langle \beta_j^*|}{\langle \beta_j^*|\alpha_j\rangle},$$4

Here $$\hat{\rho}=\int d^{4M}\bm{\lambda}\; P_+(\bm{\lambda})\,\hat{\Lambda}(\bm{\lambda}), \qquad \hat{\Lambda}(\bm{\lambda}) = \bigotimes_j \frac{|\alpha_j\rangle \langle \beta_j^*|}{\langle \beta_j^*|\alpha_j\rangle},$$5 is a stochastic weight, $$\hat{\rho}=\int d^{4M}\bm{\lambda}\; P_+(\bm{\lambda})\,\hat{\Lambda}(\bm{\lambda}), \qquad \hat{\Lambda}(\bm{\lambda}) = \bigotimes_j \frac{|\alpha_j\rangle \langle \beta_j^*|}{\langle \beta_j^*|\alpha_j\rangle},$$6 is a drift gauge, and $$\hat{\rho}=\int d^{4M}\bm{\lambda}\; P_+(\bm{\lambda})\,\hat{\Lambda}(\bm{\lambda}), \qquad \hat{\Lambda}(\bm{\lambda}) = \bigotimes_j \frac{|\alpha_j\rangle \langle \beta_j^*|}{\langle \beta_j^*|\alpha_j\rangle},$$7 is a diffusion gauge with $$\hat{\rho}=\int d^{4M}\bm{\lambda}\; P_+(\bm{\lambda})\,\hat{\Lambda}(\bm{\lambda}), \qquad \hat{\Lambda}(\bm{\lambda}) = \bigotimes_j \frac{|\alpha_j\rangle \langle \beta_j^*|}{\langle \beta_j^*|\alpha_j\rangle},$$8. The positive-P method is recovered when $$\hat{\rho}=\int d^{4M}\bm{\lambda}\; P_+(\bm{\lambda})\,\hat{\Lambda}(\bm{\lambda}), \qquad \hat{\Lambda}(\bm{\lambda}) = \bigotimes_j \frac{|\alpha_j\rangle \langle \beta_j^*|}{\langle \beta_j^*|\alpha_j\rangle},$$9, $$\bm{\lambda}=\{\alpha_1,\dots,\alpha_M,\beta_1,\dots,\beta_M\}$$0, and $$\bm{\lambda}=\{\alpha_1,\dots,\alpha_M,\beta_1,\dots,\beta_M\}$$1 (Wüster et al., 2017).

3. Observable estimation and multi-time correlations

At equal time, normally ordered moments are obtained directly from the phase-space variables. This directness extends beyond equal-time moments because the positive-P kernel has derivative-free identities on one side for $$\bm{\lambda}=\{\alpha_1,\dots,\alpha_M,\beta_1,\dots,\beta_M\}$$2 and on the other side for $$\bm{\lambda}=\{\alpha_1,\dots,\alpha_M,\beta_1,\dots,\beta_M\}$$3. Those identities allow certain multi-time observables to be estimated using time-dependent samples from the same stochastic trajectory rather than response corrections or auxiliary constructions (Deuar, 2020).

For time-ordered normal-ordered observables, with

$$\bm{\lambda}=\{\alpha_1,\dots,\alpha_M,\beta_1,\dots,\beta_M\}$$4

the central positive-P result is

$$\bm{\lambda}=\{\alpha_1,\dots,\alpha_M,\beta_1,\dots,\beta_M\}$$5

A basic two-time example is

$$\bm{\lambda}=\{\alpha_1,\dots,\alpha_M,\beta_1,\dots,\beta_M\}$$6

The crucial point is that the time history of one trajectory already contains the estimator for the corresponding time-ordered quantum observable (Deuar, 2020).

The same framework supports nonpolynomial observables used in non-Gaussian applications. For a single mode, parity is estimated as

$$\bm{\lambda}=\{\alpha_1,\dots,\alpha_M,\beta_1,\dots,\beta_M\}$$7

and for a multimode system

$$\bm{\lambda}=\{\alpha_1,\dots,\alpha_M,\beta_1,\dots,\beta_M\}$$8

The method is therefore not restricted to low-order moments, although the numerical stability of such estimators can be substantially worse than for ordinary normally ordered correlators (Shi et al., 11 Jan 2026).

The multi-time formalism also clarifies which mixed-order observables are accessible. Early normal-ordered blocks followed by later anti-normal-ordered blocks can be treated by switching from positive-P to doubled-$$\bm{\lambda}=\{\alpha_1,\dots,\alpha_M,\beta_1,\dots,\beta_M\}$$9-type representations, whereas the reverse ordering is not accessible by the same conversion logic because one cannot stochastically convert $$\hat{\Lambda}(\vec{\alpha},\vec{\alpha}^{+})=\frac{|\vec{\alpha}\rangle\langle (\vec{\alpha}^{+})^{*}|}{\langle (\vec{\alpha}^{+})^{*}|\vec{\alpha}\rangle}.$$0-samples back into $$\hat{\Lambda}(\vec{\alpha},\vec{\alpha}^{+})=\frac{|\vec{\alpha}\rangle\langle (\vec{\alpha}^{+})^{*}|}{\langle (\vec{\alpha}^{+})^{*}|\vec{\alpha}\rangle}.$$1-samples. Explicit tallies are given up to fourth order in the cited work (Deuar, 2020).

4. Noise growth, instability, and gauge control

The main practical limitation of positive-P simulations is not formal exactness but noise amplification. In long-range interacting bosonic systems, the useful simulation time is limited by noise arising from interactions, and the paper devoted to that setting emphasizes that the long-range character of the interactions does not further increase the limitations of these methods in the way that it complicates alternatives such as the density matrix renormalisation group (Wüster et al., 2017).

A convenient diagnostic in gauge-P is the characteristic variance

$$\hat{\Lambda}(\vec{\alpha},\vec{\alpha}^{+})=\frac{|\vec{\alpha}\rangle\langle (\vec{\alpha}^{+})^{*}|}{\langle (\vec{\alpha}^{+})^{*}|\vec{\alpha}\rangle}.$$2

When $$\hat{\Lambda}(\vec{\alpha},\vec{\alpha}^{+})=\frac{|\vec{\alpha}\rangle\langle (\vec{\alpha}^{+})^{*}|}{\langle (\vec{\alpha}^{+})^{*}|\vec{\alpha}\rangle}.$$3 grows to around 10, the simulation becomes unreliable. Stochastic gauges are introduced specifically to slow this variance growth without changing physical observables. The practical conclusion reported for nonlocal interactions is that, for small to medium systems, drift gauges are beneficial, whereas for sufficiently large systems it is optimal to use only a diffusion gauge (Wüster et al., 2017).

A distinct numerical pathology is the appearance of diverging or runaway trajectories. In an $$\hat{\Lambda}(\vec{\alpha},\vec{\alpha}^{+})=\frac{|\vec{\alpha}\rangle\langle (\vec{\alpha}^{+})^{*}|}{\langle (\vec{\alpha}^{+})^{*}|\vec{\alpha}\rangle}.$$4 coherent-state formulation, a discrete analogue of the positive-P representation is used to address this by evolving trajectories in continuous phase space, projecting trajectories that become too large onto equivalent discrete phase-space points, and then continuing the evolution from the discrete kernel. Because the current kernel can be re-expressed as a convex combination of discrete kernels, this resampling suppresses runaway behavior without changing the exact ensemble average in the large-sample limit (Žunkovič, 2015).

Boundary-term errors provide another major limitation. In multimode cat-state simulations, parity is identified as the most unstable observable: it can show unexpected decay even when photon numbers, coherent amplitudes, and spatial correlations remain accurate, and spikes in parity error bars are indicators of boundary-term errors. In unstable runs, reconstruction of the Wigner function at the origin loses the central interference fringes, suggesting drift toward the classical mixture

$$\hat{\Lambda}(\vec{\alpha},\vec{\alpha}^{+})=\frac{|\vec{\alpha}\rangle\langle (\vec{\alpha}^{+})^{*}|}{\langle (\vec{\alpha}^{+})^{*}|\vec{\alpha}\rangle}.$$5

rather than a coherent cat superposition. In that particular multimode linear model, stochastic gauge techniques do not significantly improve stability and may fail to reduce noise (Shi et al., 11 Jan 2026).

5. Applications and generalizations

The positive-P method has been applied across driven-dissipative quantum optics, many-body bosonic dynamics, Bell-correlation simulations, and symmetry-adapted phase-space constructions. In the cited literature it is used both as a foundational exact representation and as a scalable computational framework for systems beyond direct density-matrix reach.

System or setting	Role of positive-P	Representative source
Cascade atomic ensemble	Simulates fluctuation-initiated superradiant signal and idler emission; computes intensities and signal-idler correlations	(Jen, 2012)
Bell-violating photonic states	Samples Bell states and parametric down-conversion dynamics probabilistically but exactly	(Rosales-Zárate et al., 2014)
Long-range interacting bosons	Provides positive-P and gauge-P trajectory simulations without Hilbert-space truncation	(Wüster et al., 2017)
Multimode cat-state arrays	Derives exact SDEs and simulates transient dynamics for networks up to $$\hat{\Lambda}(\vec{\alpha},\vec{\alpha}^{+})=\frac{|\vec{\alpha}\rangle\langle (\vec{\alpha}^{+})^{*}|}{\langle (\vec{\alpha}^{+})^{*}|\vec{\alpha}\rangle}.$$6 sites	(Shi et al., 11 Jan 2026)
Multi-time driven-dissipative lattices	Computes two-time and higher multi-time correlators; demonstrated on unconventional photon blockade and a 32-site Bose–Hubbard chain	(Deuar, 2020)
Discrete $$\hat{\Lambda}(\vec{\alpha},\vec{\alpha}^{+})=\frac{|\vec{\alpha}\rangle\langle (\vec{\alpha}^{+})^{*}|}{\langle (\vec{\alpha}^{+})^{*}|\vec{\alpha}\rangle}.$$7 phase space	Embeds positive-P as an $$\hat{\Lambda}(\vec{\alpha},\vec{\alpha}^{+})=\frac{|\vec{\alpha}\rangle\langle (\vec{\alpha}^{+})^{*}|}{\langle (\vec{\alpha}^{+})^{*}|\vec{\alpha}\rangle}.$$8 kernel and introduces a discrete analogue for runaway-trajectory control	(Žunkovič, 2015)
Symmetry-aware matrix phase space	Recovers ordinary positive-P when symmetry projection is absent; reduces variance when symmetry sectors are built into the basis	(Drummond et al., 11 Jun 2026)
Majorana fermions	Supplies a fermionic Positive-P/Q-like phase-space analogue based on Gaussian operators and real antisymmetric covariance variables	(Joseph et al., 2017)

Several application-specific conclusions recur. In superradiant cascade emission, increasing optical depth shortens the signal-idler correlation time scale, implying a broader spectral window needed to store or retrieve the idler pulse (Jen, 2012). In Bell simulations, the method reproduces Bell violations because the stochastic variables are not restricted to the eigenvalue range of the measured observables, so the representation is not a hidden-variable model in Bell’s sense (Rosales-Zárate et al., 2014). In dissipatively coupled resonator arrays, the positive-P method extends simulations to system sizes significantly larger than those accessible via direct master-equation simulation and, in principle, to $$\hat{\Lambda}(\vec{\alpha},\vec{\alpha}^{+})=\frac{|\vec{\alpha}\rangle\langle (\vec{\alpha}^{+})^{*}|}{\langle (\vec{\alpha}^{+})^{*}|\vec{\alpha}\rangle}.$$9 sites, while remaining especially effective for normally ordered observables and correlation functions (Shi et al., 11 Jan 2026).

The method also serves as a platform for structurally richer representations. Matrix phase-space expands the kernel to include an $\alpha_j$0 stochastic matrix of symmetry weights, and for bosons the matrix-P representation reduces to the gauge-P for normal ordering and $\alpha_j$1, while $\alpha_j$2 gives the ordinary Positive-P expansion. A key existence theorem in that framework proves that for any bosonic density matrix, a positive semidefinite matrix-P representation exists (Drummond et al., 11 Jun 2026).

6. Relation to $\alpha_j$3, $\alpha_j$4, Wigner, and other positive phase-space constructions

The positive-P method is closely related to other phase-space representations but should not be conflated with them. The ordinary $\alpha_j$5-representation writes

$\alpha_j$6

which is classically suggestive but is not always a regular positive function and can be singular. The Husimi $\alpha_j$7-function is positive and smooth, but it is broadened. By contrast, the positive-P representation keeps a positive stochastic interpretation by enlarging the phase space through off-diagonal coherent-state kernels and independent complex variables (Abe et al., 2014, Deuar, 2020).

The distinction from Wigner-based constructions is sharper. Exact Wigner dynamics on ordinary $\alpha_j$8 phase space is not generally the density of a positive stochastic process; for non-quadratic potentials the Wigner–Moyal generator contains higher odd derivatives and signed momentum-transfer terms that are not the Fokker–Planck generator of a positive Brownian diffusion. In that setting, the exact Wigner function must be reconstructed as a weighted empirical measure,

$\alpha_j$9

where the carrier process is positive but the weights may be signed. Positive-P addresses the positivity problem differently: it preserves a positive stochastic density by moving to a doubled or enlarged phase space rather than by insisting on a positive density for the exact Wigner function on ordinary phase space (Limkumnerd et al., 7 May 2026).

It is also distinct from variational programs that seek a single positive probability density in ordinary phase space that is as close as possible to the Wigner function. One such construction defines a positive normalized density $\beta_j$0 closest to $\beta_j$1 in a quadratic metric and interprets the minimum deviation as a measure of quantumness. That approach is explicitly not about stochastic simulation or phase-space doubling; it is about optimal classicalization of a given Wigner distribution (Roy et al., 2013).

A common misconception is that probabilistic sampling with positive-P implies a hidden-variable model. The Bell-violation simulations make the contrary point explicit: the sampled variables are quasi-observables and may assume values outside the normal eigenvalue range, in parallel with weak values. The representation remains exactly equivalent to quantum mechanics because it reproduces the density operator and the required operator moments, not because its stochastic variables correspond to measurement outcomes (Rosales-Zárate et al., 2014).

Taken together, these comparisons place the positive-P phase space method in a precise niche. It is an exact, doubled-phase-space stochastic representation for bosonic quantum dynamics; it generalizes naturally to gauge-P, discrete $\beta_j$2, and symmetry-projected matrix forms; it supports direct estimators for large classes of equal-time and multi-time observables; and its principal obstacles are not representational but numerical, namely noise growth, boundary terms, and instability in observables that are especially sensitive to rare trajectories.