Stochastic phase-space simulation of multimode cat states via the positive-P representation

We present a comprehensive study of the transient dynamics of multimode Schrödinger cat states in dissipatively coupled resonator arrays using the positive-P phase-space method. By employing the positive-P representation, we derive the exact stochastic differential equations governing the system’s dynamics, enabling the simulation of system sizes significantly larger than those accessible via direct master equation simulation. We demonstrate the utility of this method by simulating transient dynamics for networks up to N=21 sites. Furthermore, we critically examine the method’s usefulness and limitations, specifically highlighting the computational instability encountered when estimating the state parity in the systems. Our results provide a pathway for scalable simulations of non-Gaussian states in large open quantum systems.

💡 Research Summary

The paper addresses a pressing challenge in the simulation of large open quantum systems that host non‑Gaussian resources such as multimode Schrödinger cat states. Conventional master‑equation approaches quickly become intractable because the Hilbert space dimension grows exponentially with the number of resonators (N). To overcome this bottleneck, the authors employ the positive‑P phase‑space representation, which maps the density matrix onto a probability distribution over pairs of complex variables (αj, βj) for each mode. By converting the Lindblad master equation into a Fokker‑Planck equation and then into Itô stochastic differential equations (SDEs), they obtain an exact set of equations (Eqs. 3a‑3b) that can be integrated numerically using ensembles of stochastic trajectories.

The key advantage of this approach is that the computational cost scales linearly with the number of modes and the number of trajectories, allowing the authors to simulate arrays up to N = 21 resonators—far beyond the N ≈ 3 limit of direct master‑equation integration. The SDEs contain independent white‑noise terms ξj(t) and \tilde{ξ}j(t), which are realized as Gaussian random numbers with variance 1/Δt. Normal‑ordered observables are obtained directly from ensemble averages of products of α and β (Eq. 4), and the full density matrix can in principle be reconstructed from the kernel operator (Eq. 5), although this is computationally expensive.

The authors first benchmark the method on a single‑mode system (N = 1, γ = 0). They examine four observables: average photon number n, coherent amplitude ζ, second‑order correlation g₂, and parity Π. By scanning the parameter space of single‑photon loss κ₁ and two‑photon loss κ₂, they identify three regimes. In the “orange” region (large κ₂) trajectories become unstable before any transient cat state forms. In the “blue” region (moderate κ₂) n, ζ, and g₂ agree with exact master‑equation results, but Π decays unexpectedly, indicating loss of parity information due to boundary‑term errors inherent to the positive‑P method. In the “green” region (dominant κ₁) all observables are stable, but the system relaxes to a trivial vacuum or heavily damped state, limiting its relevance for cat‑state physics.

Next, the method is extended to multimode arrays (N = 3) with strong non‑local dissipation γ ≫ κ₁, κ₂ (the Zeno limit). Comparisons with quantum‑trajectory simulations performed using QuTiP’s mcsolve show excellent agreement for local photon numbers, coherent amplitudes, and both first‑order (g₁) and second‑order (g₂) spatial correlation functions. However, the global parity of the multimode cat state again exhibits a rapid, unphysical decay, while local parity (Π_loc) remains accurate—a stability that the authors attribute to the stabilizing effect of the large γ term on the stochastic trajectories.

A major contribution of the work is the demonstration that positive‑P enables the study of spatial correlations in large arrays (N = 21) with realistic parameters (κ₁/ε = 10⁻³, κ₂/ε = 0.2). In the absence of non‑local dissipation (γ = 0) each resonator independently forms a cat state, leading to vanishing off‑diagonal elements of g₁ and Poissonian statistics (g₂ ≈ 1). Introducing weak non‑local dissipation (γ/ε = 0.2) yields a finite coherence length: g₁ displays an exponential decay away from the diagonal, while g₂ remains close to unity but shows a subtle ridge indicating weak inter‑site correlations. These results illustrate how engineered dissipation can mediate controlled propagation of quantum coherence across a network of resonators.

The authors also explore mitigation strategies for the observed instability. They implement drift‑gauge modifications—additional deterministic terms designed to suppress the growth of trajectories that lead to boundary‑term errors. While these gauges improve stability in the single‑mode case, they provide only marginal benefit for multimode systems, especially when γ is small. The persistent difficulty in estimating parity reflects a fundamental limitation: parity involves the exponential of the product αβ, which is highly sensitive to the tails of the stochastic distribution. Consequently, even with large ensembles (up to 10⁶ trajectories) the parity estimator suffers from large statistical fluctuations and systematic decay.

In summary, the paper establishes the positive‑P representation as a powerful tool for simulating the transient dynamics of large, dissipatively coupled resonator arrays that host multimode cat states. It delivers accurate predictions for photon numbers, coherent amplitudes, and spatial correlation functions across system sizes previously inaccessible to exact methods. At the same time, it clearly delineates the method’s shortcomings—most notably the unreliable estimation of parity and the associated boundary‑term errors. The work points toward future directions, such as stochastic‑gauge extensions, hybrid complex‑P/positive‑P schemes, or alternative phase‑space techniques, to overcome these hurdles and enable fully reliable simulations of non‑Gaussian quantum resources in scalable quantum hardware.