text:
# Gauge P representations for quantum-dynamical problems: Removal of boundary terms

## I Introduction

One of the most difficult problems in theoretical physics is also conceptually the simplest. How does one calculate the dynamical time evolution or even the ground state of an interacting many-body quantum system? In essence, this is a natural part of practically any comparison of quantum theory with experiment. The difficulty is that the Hilbert space of all but the most trivial cases can be enormous. This implies that a finite computer is needed to to solve problems that can easily become nearly infinite in dimensionality, if treated using an orthogonal basis expansion.

In this paper, we formally introduce and give examples of techniques for treating general bosonic many-body quantum systems, which we call gauge P representations. These are an extension of the phase-space method called the positive-P representation DG-PosP , and have been recently used in the context of interacting Bose gases ccp2k ; canonical . The advantages of the new technique are the following.

(1) The elimination of certain types of mathematical terms known as boundary-term corrections, which have caused problems in the positive-P representation for over a decade SG-abs ; SS-fail ; GGD-Validity . This is the main focus of the present paper.

(2) Greatly reduced sampling error in computations. Gauge P representations have been used recently to reduce the sampling error in Kerr oscillator simulations ccp2k .

(3) The extension of allowable problems to ’imaginary-time’ canonical ensemble calculations. These problems will be treated elsewhere.

Related extensions to the positive-P representation — although restricted to the scalar interacting Bose gas problem — have also been introduced recently. Different procedures have been introduced by Carusotto, Castin, and Dalibard Paris1 ; Paris2 , and by Plimak, Olsen, and Collett Plimak . These methods implicitly assume the absence of boundary term corrections. This paper unifies and substantially generalizes all these recent advances. It also shows how the gauge method can be used to solve the long-standing problem of boundary-term corrections in the positive P representation. Comparisons to the other methods are given in an Appendix.

Owing to the work of Wilson wilson:74 , and many others Ceperley , we know that large Hilbert space problems can often be treated using stochastic or Monte Carlo techniques for the ground-state, particle masses, and finite-temperature correlations. This is the basis for much work in computational quantum statistical mechanics, and in QCD as well. However, Wilson’s and other related methods are restricted to static or ‘imaginary-time’ calculations, rather than quantum-dynamical problems.

Methods like these that use orthogonal basis sets have not proven useful for quantum dynamics; owing to the notorious phase problem that occurs when trying to sum over families of paths in real-time Feynman path integrals. For this reason, the many-body quantum time-evolution problem is often regarded as inherently insoluble due to its exponential complexity. In fact, it was this very problem that motivated the original proposal of Feynman feynman:82  to develop quantum computers. In these (usually conceptual) devices, the mathematical problem is solved by a physical system consisting of evolving ‘qubits’ or two-state physical devices. Fortunately, this method of doing calculations is not the only one, since no large enough quantum computer exists at presentFactor15 .

Historically, an alternative route is the use of quasi-probability representations of the quantum state, which either implicitly or explicitly make use of a non orthogonal basis. The term quasi-probability is used because there can be no *exact* mapping of all quantum states to a classical phase space with a positive distribution Neumann  that also preserves all the marginal probabilities. These methods include the Wigner Wig-Wigner  (W), Glauber-Sudarshan (P) Gla-P ; Sud-P , and Husimi (Q) Hus-Q ; CG-Q  representations. The classical phase-space representations can be classified according to the operator ordering that stochastic moments correspond to: the W is symmetrically ordered, the Q is anti-normally ordered, while the P representation is normally ordered. Apart from numerous laser physics and quantum optics calculations, these methods have also been used to some extent in quantum statistical mechanics: for example, the theory of BEC phase fluctuations becpf .

None of these methods result in a stochastic time evolution with a positive propagator when there are nonlinearities. To achieve this, a better approach is to use a non-classical phase space of higher dimension. A complex higher-dimensional ‘R representation’ was proposed in Glauber’s seminal paper on coherent state expansions Gla-P . The first probabilistic method of this type was the positive-P representation DG-PosP  (+P), which has proved capable of performing stochastic time-domain quantum calculations in some many-body quantum systems carterdrummond:87 . This uses a basis of coherent states that are not orthogonal, thus allowing freedom of choice in the construction of the representation. The positive-P representation of a quantum state is therefore the most versatile out of a large group of quasi-probability distributions developed to aid quantum mechanical calculations. It has been successfully applied to mesoscopic systems such as quantum solitons carterdrummond:87 ; qsol2 ; qsol3  and the theory of evaporative cooling evapcool , which correctly reproduces the formation of a BEC — as observed in experiment obsex1 ; obsex2 ; obsex3 .

Quasi-probability distributions of this type are computationally superior to direct density matrix methods, which are susceptible to computational complexity blow-up for large Hilbert spaces. Provided certain boundary terms vanish, the usual procedure is to generate a Fokker-Planck equation (which will vary depending on the distribution chosen) from the master equation, and then to convert this to a set of stochastic Langevin equations. For some simple cases, it may even be possible to arrive at appealing results directly from the Fokker-Planck equation (FPE). The resulting stochastic equations can be thought of just as quantum mechanics written in different variables. They have two main advantages over orthogonal basis-state methods, as follows.

First, the whole quantum dynamics can be written exactly in terms of a small number of stochastic equations. In a one-mode case, there is just one complex variable for P and Q and W, and two complex variables for +P. Although a simulation requires us to average over many realizations of the stochastic process, this is often more practical than solving the infinite set of deterministic equations required to solve directly for all the elements of a density matrix. Such an infinite set may be truncated, but this is only a good approximation for a system with few particles, and no more than a few modes.

Second, for a many-mode problem the Hilbert space dimension is $`N=n^M`$ for the case of $`n`$ particles distributed over $`M`$ modes. This gives exponential growth as a function of the number of modes. However, the number of quasi-probability dynamical equations grows only linearly with the number of modes, rather than exponentially in the case of direct methods. Other stochastic methods, known as quantum-trajectory methods, can be used to reduce the $`N^2`$ dimensionality of an $`N\times N`$ density matrix problem to that of the $`N`$-dimensional underlying Hilbert space — but this is clearly insufficient to solve the complexity problem inherent in the exponential growth of the Hilbert space dimension.

There are, however, some caveats when using these distributions. In particular, the vanishing of boundary terms is an important fundamental issue with quasi-probability distributions, and it is this issue that we focus on mostly in this paper. To get an overall picture, consider that once we have a time-evolution problem there are five typical requirements that are encountered in deriving stochastic equations for quasi-probability representations of many-body systems. These requirements occur in closed (unitary evolution) systems, in open systems (in general described by a master equation), or even using a distribution to solve for the canonical ensemble in imaginary time. As such, these requirements are generic to the use of stochastic equations with operator representations:

(1) Positive distribution. A well-behaved positive distributions for all quantum states, including especially the chosen initial condition, is essential for a general algorithm. For example, a number state has a highly singular P distribution, and a W distribution that is negative in some regions of phase space KS-70 , making either distribution impossible to interpret probabilistically for these states. The R distribution is inherently complex. Such problems do not occur for the Q or +P representation — these are positive, and well-behaved for all quantum states DG-PosP .

(2) *Ultraviolet convergence*. While normally-ordered representations are well behaved at large momentum, non-normally-ordered representations of quantum fields — such as the Q or W representations — typically face the problem of ultraviolet divergence in the limit of large momentum cutoff evapcool . This means that almost any observable quantity will involve the simulation of a (nearly) infinitely noisy classical field, leading to diverging standard deviations in two or more space dimensions, even for linear systems. This rules out the Q and W distributions for quantum field simulations in higher than one-dimensional environments.

(3) Second-order derivatives. Only FPEs with second or infinite-order derivatives can be translated into stochastic equations Arnold . Normally-ordered methods such as the P and +P representations can handle most commonly occurring nonlinearities and two-body interactions, with only second-order derivatives. Non-normally ordered representations of quantum fields often lead to third- or higher-order partial derivatives in the Fokker-Planck equation with no stochastic equivalent. For example, the Wigner representation gives such problems for almost any nonlinear term in the master equation.

(4) Positive-definite diffusion. A Fokker-Planck equation must have positive-definite diffusion, to allow simulation with stochastic processes Arnold . When the master equation has nonlinear terms, this does not occur with any of the classical representations. However, the +P representation is guaranteed to always produce positive-definite diffusion DG-PosP , provided no higher derivative terms occur.

(5) Vanishing boundary terms. In the derivation of the Fokker-Planck equations, it is assumed that certain boundary terms arising in partial integration can be neglected. This is not always the case. Boundary terms due to power-law tails can occur when there are moving singularities that can escape to infinity in finite time. In the +P method, such trajectories may cause systematic errors in stochastic averages GGD-Validity , especially for non-integrable dynamical systems. These problems are exponentially suppressed when linear damping rates are increased, but can be large at low damping.

The +P method is often the representation of choice, because it satisfies conditions one to four. Gauge representations (G) combined with stochastic methods to be treated in this paper, share these advantages with the +P representation. However, they can also satisfy the fifth requirement — for an appropriate gauge choice — hence allowing all of the mathematical problems in simulating time evolution to be treated. For this reason, the present paper will focus on solving boundary-term issues encountered with the +P representation for certain nonlinear master equations. The overall picture is summarized in Table 1, as applied to the two-boson anonlinear absorber cases treated here in Sec. IV:

We emphasize that the particular examples treated here have a small particle number and extremely low (or zero) linear damping. As such, they are soluble using other techniques, which allows us to test the accuracy of gauge techniques. Our purpose is to demonstrate the success of the stochastic gauge method in simple cases where boundary terms arise within the +P representation. In this way, we can understand more complex situations where no exact result is known.

We will first derive and describe the stochastic gauge method in Secs. II and III, and subsequently work through two examples: First, solving the boundary-value problem for the driven one- and two-photon absorber in Sec. IV. Second, in Sec. V we will consider the one-mode laser at extremely low power, which exhibits boundary term errors when very non-optimal starting conditions are used. This example will show that gauge methods can also be used to remove errors from this system, but some judgment must be employed to avoid choosing a pathological initial distribution. In the Appendix, we compare the methods derived here with recent related extensions of the positive-P representation by Carusotto and co-workers Paris1 ; Paris2 , Plimak et al. Plimak , and Deuar and Drummond ccp2k .

Finally, we point out a sixth requirement of containing the growth of sampling error: the averages calculated from the stochastic Langevin equations correspond to quantum mechanical expectation values only in the limit of infinitely many trajectories. Provided boundary terms do not occur, the averages will approach the correct values — within an acceptable sampling error — for sufficiently many trajectories. If this number should increase rapidly with time, the simulation will only be of use for a limited period ccp2k .

The problem of growing sampling error can occur even when there are no boundary terms, and may be regarded as the ultimate frontier in representation theory, just as similar issues dominate the theory of classical chaos. This is less of a fundamental issue, since the sampling error can always be estimated and controlled by increasing the number of trajectories. This is simply a matter of moving to a clustered, parallel computational model, or repeating the calculation many times. Nevertheless, it is of great practical significance. The sampling error problem requires careful gauge optimization, and remains an open area for investigation. An intelligent choice of gauge can often vastly outweigh a brute force computational approach, in terms of sampling error.

## II Gauge Operator Representations

In gauge representations, the density matrix to be computed is expanded in terms of a coherent state basis. For definiteness, we shall focus on the coherent states of the harmonic oscillator, which are useful in expanding Bose fields; but other choices are clearly possible. The expansion kernel is more general than that used in the positive-P representation. In order to define the notation, we start by introducing a set of boson annihilation and creation operators $`\widehat{a}_i`$ , $`\widehat{a}_i^{}`$. The operator $`\widehat{n}_i=\widehat{a}_i^{}\widehat{a}_i`$ is therefore the boson number operator for the $`i`$th mode or site. Boson commutation relations of $`[\widehat{a}_i,\widehat{a}_j^{}]=\delta _{ij}`$ hold for the annihilation and creation operators.

### II.1 Coherent states

If $`𝜶=(\alpha _1,\mathrm{},\alpha _M)`$ is a complex $`M`$-dimensional vector with $`\alpha _i=x_i+iy_i`$, and $`\widehat{𝐚}=(\widehat{a}_1,\mathrm{},\widehat{a}_M)`$ is an $`M`$-dimensional vector of annihilation operators, then the Bargmann coherent state $`𝜶`$ is defined by

$$𝜶=\mathrm{exp}\left[𝜶\widehat{𝐚}^{}\right]|0=\mathrm{exp}\left[|𝜶|^2/2\right]|𝜶,$$
(1)
where $`|𝜶`$ is the usual normalized coherent state which is a simultaneous eigenstate of all the annihilation operators. The inner product of two Bargmann coherent states is

$$𝜷^{}𝜶=\mathrm{exp}\left[𝜶𝜷\right].$$
(2)
It is important to notice here that $`𝜶`$ is an analytic function of the complex vector $`𝜶`$ . The following identities therefore follow immediately:

$`\widehat{a}_i𝜶`$ $`=`$ $`\alpha _i𝜶`$
$`\widehat{a}_i^{}𝜶`$ $`=`$ $`{\displaystyle \frac{}{\alpha _i}}𝜶.`$ (3)
Since $`𝜶`$ is an analytic function, the notation $`/\alpha _i`$ is interpreted here as an analytic derivative, which can be evaluated in either the real or imaginary directions,

$$\frac{}{\alpha _i}𝜶=\frac{}{x_i}𝜶=i\frac{}{y_i}𝜶.$$
(4)
Since the coherent states are an over-complete basis set, any operator can be expanded in more than one way using coherent states. For example, the simplest resolution of the identity operator is

$$\widehat{I}=\frac{1}{\pi ^M}|𝜶𝜶|d^{2M}𝜶.$$
(5)
Thus, introducing a second $`M`$-dimensional vector $`𝜷`$, we can expand any operator $`\widehat{O}`$ directly as

$`\widehat{O}`$ $`=`$ $`{\displaystyle \frac{1}{\pi ^{2M}}}{\displaystyle |𝜶𝜶\left|\widehat{O}\right|𝜷^{\mathbf{}}𝜷^{\mathbf{}}|d^{2M}𝜶d^{2M}𝜷}`$ (6)
$`=`$ $`{\displaystyle O(𝜶,𝜷)|𝜶𝜷^{\mathbf{}}|d^{2M}𝜶d^{2M}𝜷}.`$
Here, we have introduced

$$O(𝜶,𝜷)=\frac{1}{\pi ^{2M}}𝜶\left|\widehat{O}\right|𝜷^{\mathbf{}}.$$
(7)
### II.2 P representations

The possibility of expanding any operator in terms of coherent states leads to the idea that such an expansion can be used to calculate observable properties of a quantum density matrix $`\widehat{\rho }`$ . Historically, this was first proposed by Glauber and Sudarshan Gla-P ; Sud-P , who suggested a diagonal expansion of the form

$$\widehat{\rho }=P(𝜶)|𝜶𝜶|d^{2M}𝜶.$$
(8)
Unlike the direct expansion given above, this has no off-diagonal elements. Surprisingly, expansions of this type always exist, as long as the function $`P(𝜶)`$ is defined to allow highly singular generalized functions and non-positive distributions KS-70 .

As these do not have a stochastic interpretation, the positive-P representation was introduced DG-PosP , which is defined as

$$\widehat{\rho }=P^{(+)}(𝜶,𝜷)\frac{|𝜶𝜷^{\mathbf{}}|}{𝜷^{\mathbf{}}|𝜶}d^{2M}𝜶d^{2M}𝜷$$
(9)
for an $`M`$-mode system.

It is always possible to obtain an explicitly positive-definite distribution of this type DG-PosP , with the definition

$`P^{(+)}(𝜶,𝜷)`$ $`=`$ $`{\displaystyle \frac{1}{(4\pi ^2)^M}}\mathrm{exp}\left[\left|{\displaystyle \frac{𝜶𝜷^{\mathbf{}}}{2}}\right|^2\right]`$ (10)
$`\times {\displaystyle \frac{𝜶+𝜷^{\mathbf{}}}{2}}\left|\widehat{\rho }\right|{\displaystyle \frac{𝜶+𝜷^{\mathbf{}}}{2}}.`$
This form always exists, as do an infinite class of equivalent positive distributions. Even simpler ways to construct the positive-P representation are available in some cases. For example, if the Glauber-Sudarshan representation exists and is positive, then one can simply construct

$$P^{(+)}(𝜶,𝜷)=P(𝜶)\delta ^{2M}(𝜶𝜷^{\mathbf{}}).$$
(11)
The stochastic time evolution of the positive-P distribution does not generally preserve the above compact forms, and may allow less compact positive solutions instead. However, to obtain a time-evolution equation, it is necessary to use partial integration, with the assumption that boundary terms at infinity can be neglected. It is these less compact solutions, occurring during time evolution with a nonlinear Fokker-Planck equation, that lead to power-law tails in the distribution — and hence boundary-term problems caused by the violation of the assumption that these terms vanish.

### II.3 Gauge representations

A technique for constructing an even more general positive distribution is to introduce a quantum complex amplitude $`\mathrm{\Omega }`$, which can be used to absorb the quantum phase factor. This leads to the result that any Hermitian density matrix $`\widehat{\rho }`$ can be expanded in an over-complete basis $`\widehat{\mathrm{\Lambda }}(\stackrel{}{\alpha })`$, where $`\stackrel{}{\alpha }=(\mathrm{\Omega },𝜶,𝜷)`$, and

$`\widehat{\mathrm{\Lambda }}(\stackrel{}{\alpha })`$ $`=`$ $`\mathrm{\Omega }{\displaystyle \frac{𝜶𝜷^{}}{𝜷^{}𝜶}}`$ (12)
$`=`$ $`\mathrm{\Omega }𝜶𝜷^{}\mathrm{exp}\left[𝜶𝜷\right].`$
We define the gauge representation $`G(\stackrel{}{\alpha })`$ as a real, positive function that satisfies the following equation:

$`\widehat{\rho }`$ $`=`$ $`{\displaystyle G(\stackrel{}{\alpha })\left[\widehat{\mathrm{\Lambda }}(\stackrel{}{\alpha })\right]d^{4M+2}\stackrel{}{\alpha }}`$ (13)
$`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle G(\stackrel{}{\alpha })\left[\widehat{\mathrm{\Lambda }}(\stackrel{}{\alpha })+\text{H.c.}\right]d^{4M+2}\stackrel{}{\alpha }}.`$
The last line above follows from the fact that $`\widehat{\rho }`$ is a Hermitian density matrix and $`G(\stackrel{}{\alpha })`$ is real. Here, H.c. is used as an abbreviation for Hermitian conjugate. The use of a complex weight in the above gauge representation is similar to related methods introduced recently for interacting Bose gases Paris1 ; Paris2 , except that we multiply the weight by a normalized (positive-P) projector, in order to simplify the resulting algebra.

As an existence theorem that shows that this representation always exists, consider the complex solution

$$P_0(𝜶,𝜷)=\frac{1}{\pi ^{2M}}𝜶\left|\widehat{\rho }\right|𝜷^{\mathbf{}}𝜷^{\mathbf{}}|𝜶$$
(14)
obtained from Eq. (7), with a phase $`\theta =\mathrm{arg}(P_0)`$ , and simply define

$$G(\stackrel{}{\alpha })=\left|P_0(𝜶,𝜷)\right|\delta ^2(\mathrm{\Omega }\mathrm{exp}[i\theta (𝜶,𝜷)]).$$
(15)
In this type of gauge representation, $`G(\stackrel{}{\alpha })`$ is a positive distribution over a set of Hermitian density-matrix elements $`\widehat{\mathrm{\Lambda }}+\widehat{\mathrm{\Lambda }}^{}`$. It is simple to verify that, by construction

$$\text{Tr}\left(\widehat{\mathrm{\Lambda }}\right)=\mathrm{\Omega }.$$
(16)
For the case of $`\mathrm{\Omega }=1`$, this representation reduces to the positive-P representation, and the kernel $`\widehat{\mathrm{\Lambda }}(\stackrel{}{\alpha })`$ is a projection operator. Since the positive-P representation is a complete representation, it follows that another way to construct the gauge P representation is always available, if one simply defines

$$G(\stackrel{}{\alpha })=P^{(+)}(𝜶,𝜷)\delta ^2(\mathrm{\Omega }1).$$
(17)
As a simple example, a thermal ensemble with $`n_0`$ bosons per mode gives a diagonal P distribution that is Gaussian, so that

$$G_{th}(\stackrel{}{\alpha })\mathrm{exp}\left[\left|𝜶\right|^2/n_0\right]\delta ^{2M}(𝜶𝜷^{\mathbf{}})\delta ^2(\mathrm{\Omega }1).$$
(18)
One advantage of the proposed representation is that it allows more general expansions than the positive-P distribution, and also includes the case of the complex P representation — which has proved useful in solving for non-equilibrium steady-states in quantum systems.

### II.4 Operator identities

The utility of these methods arises when they are used to calculate time (or imaginary time — for which the positive-P distribution cannot be used) evolution of the density matrix. This occurs via a Liouville equation of generic form

$$\frac{}{t}\widehat{\rho }=\widehat{L}(\widehat{\rho }),$$
(19)
where the Liouville superoperator typically involves pre- and post-multiplication of $`\widehat{\rho }`$ by annihilation and creation operators. As an example, the equation for purely unitary time evolution under a Hamiltonian $`\widehat{H}`$ is

$$i\mathrm{}\frac{}{t}\widehat{\rho }=[\widehat{H},\widehat{\rho }].$$
(20)
Effects of the annihilation and creation operators on the projectors are obtained using the results for the actions of operators on the Bargmann states,

$`\widehat{𝐚}\widehat{\mathrm{\Lambda }}(\stackrel{}{\alpha })`$ $`=`$ $`𝜶\widehat{\mathrm{\Lambda }}(\stackrel{}{\alpha })`$
$`\widehat{𝐚}^{}\widehat{\mathrm{\Lambda }}(\stackrel{}{\alpha })`$ $`=`$ $`\left[\mathbf{}_\alpha +𝜷\right]\widehat{\mathrm{\Lambda }}(\stackrel{}{\alpha })`$
$`\widehat{\mathrm{\Lambda }}(\stackrel{}{\alpha })`$ $`=`$ $`\mathrm{\Omega }_\mathrm{\Omega }\widehat{\mathrm{\Lambda }}(\stackrel{}{\alpha }).`$ (21)
For brevity, we use $`\stackrel{}{}=(_\mathrm{\Omega },\mathbf{}_\alpha ,\mathbf{}_\beta )`$ to symbolize either $`\left(_i^x/x_i\right)`$ or $`i\left(_i^y/y_i\right)`$ for each of the $`i=0,\mathrm{},2M`$ complex variables $`\stackrel{}{\alpha }`$. This is possible since $`\widehat{\mathrm{\Lambda }}(\stackrel{}{\alpha })`$ is an analytic function of $`\stackrel{}{\alpha }`$, and an explicit choice of the derivative will be made later.

Using the operator identities given above, the operator equations can be transformed to an integro-differential equation,

$`{\displaystyle \frac{\widehat{\rho }}{t}}`$ $`=`$ $`{\displaystyle G(\stackrel{}{\alpha })\left[_A\widehat{\mathrm{\Lambda }}(\stackrel{}{\alpha })\right]d^{4M+2}\stackrel{}{\alpha }}.`$ (22)
Here the anti-normal ordered notation $`_A`$ indicates an ordering of all the derivative operators to the right. As an example, in the Hamiltonian case, if the original Hamiltonian $`\widehat{H}(\widehat{𝐚},\widehat{𝐚}^{})`$ is normally-ordered (annihilation operators to the right), then

$$_A=\frac{1}{i\mathrm{}}\left[H_A(𝜶,\mathbf{}_\alpha +𝜷)H_A(𝜷,\mathbf{}_\beta +𝜶)\right].$$
(23)
If no terms higher than second order occur, this procedure gives a differential operator with the following general expansion:

$$_A^{(+)}=V+A_j^{(+)}_j+\frac{1}{2}D_{ij}_i_j.$$
(24)
where, to simplify notation, the Latin indices $`i,j,k`$ will from now on be summed over $`i=1,\mathrm{},2M`$, since no derivatives with respect to $`\mathrm{\Omega }`$ are used as yet. V is a term not involving derivative operators with respect to any of the variables in $`\stackrel{}{\alpha }.`$ The drift term $`A_j^{(+)}`$ that is normally found using the positive-P representation is labeled with the superscript $`(+)`$ to identify it.

At this stage, the usual procedure in representation theory is to integrate by parts, provided boundary terms vanish. This gives a normally-ordered differential operator acting on the distribution itself, of form

$$\frac{}{t}G(\stackrel{}{\alpha })=\left[V_jA_j^{(+)}+\frac{1}{2}_i_jD_{ij}\right]G(\stackrel{}{\alpha }).$$
(25)
This type of generalized Fokker-Planck equation can be treated formally using techniques developed by Graham, involving time-symmetric curved-space path integrals Path-int . For computational purposes, we require special choices of the analytic derivatives to obtain a positive-definite diffusion, so that the path integrals have equivalent stochastic equations Arnold . We emphasize here that the equations resulting are quite different to those obtained from the direct insertion of a coherent state identity into a Feynman path integral — which results in severe convergence problems Phase-path . The usual positive-P representation equations are obtained at this stage — provided there is no potential term — and can be transformed to stochastic equations using the techniques described in the following section.

## III Gauge functions

In gauge representations, the time evolution of the representation is modified from the usual positive-P representation equations, by the introduction of a number of arbitrary and freely defined functions on the phase space. This freedom of choice is, of course, not present with an orthogonal basis, and is due to the non-orthogonal nature of a coherent basis set. Although we do not investigate other cases, it is worth noting that a similar gauge freedom is implicitly present whenever a non-orthogonal expansion is used — even if it involves different states from the choice of coherent states made here (e.g., the Fock state wave functions in Refs. Paris1 ; Paris2 ).

### III.1 Diffusion gauges

We first introduce the diffusion gauges, which were implicitly present in the original positive-P representation, but were only recognized recently as allowing improvements in the sampling error. These gauges occur via the non-unique decomposition of the complex diffusion matrix $`D`$, which determines the stochastic correlations in the final equations. Arbitrary functional parameters can therefore be inserted into the final stochastic equations in the noise coefficients, which may lead to further optimization of the simulation. This is because the decomposition of the complex diffusion matrix $`D=BB^T`$ , which is needed to define a stochastic process, does not specify the resulting noise matrix $`B`$ completely.

It has been recently shown by Plimak, Olsen and Collett Plimak  that for the Kerr oscillator using a decomposition different from the obvious diagonal one leads to impressive improvements in the signal-to-noise ratio of the simulation (briefly described in Appendix A.2). This somewhat surprising result leads us to try to quantify the amount of freedom of choice available from this source.

Since $`D=D^T`$, it can always be diagonalized by a complex orthogonal transformation

$$D=O\lambda ^2O^T=B^{(+)}B^{(+)T},$$
(26)
where $`\lambda `$ is the diagonal matrix whose square gives the eigenvalues of $`D`$. Thus $`B^{(+)}=O\lambda `$ can be considered the canonical, or “obvious” choice of decomposition, unique apart from the $`2M`$ signs of the diagonal terms. However, for any orthogonal $`U`$, if $`B^{(+)}`$ is a valid decomposition of $`D`$, then so is the matrix $`B=B^{(+)}U`$. Hence, any matrix in the whole orthogonal family $`B=O\lambda U`$ is a valid decomposition. This can be easily quantified using a basis

$$\sigma _{kl}^{(ij)}=\delta _{ik}\delta _{jl}\delta _{il}\delta _{jk},$$

of the $`M(2M1)`$ independent antisymmetric $`2M\times 2M`$ matrices $`\underset{¯}{\underset{¯}{\sigma }}^{(ij)}`$ . One simply introduces

$$U=\mathrm{exp}\left(\underset{i<j}{}g_{ij}(\stackrel{}{\alpha },t)\underset{¯}{\underset{¯}{\sigma }}^{(ij)}\right).$$
(27)
As an example, for a one-mode case there is one complex gauge function introduced this way, which is $`g^d=g_{12}`$. The resulting transformation is

$`U`$ $`=`$ $`\mathrm{exp}\left(g^d\underset{¯}{\underset{¯}{\sigma }}^{(12)}\right)`$ (28)
$`=`$ $`\mathrm{cos}(g^d)+\underset{¯}{\underset{¯}{\sigma }}^{(12)}\mathrm{sin}(g^d),`$
where the antisymmetric matrix $`\underset{¯}{\underset{¯}{\sigma }}^{(12)}`$ is proportional to a Pauli matrix,

$$\underset{¯}{\underset{¯}{\sigma }}^{(12)}=\left[\begin{array}{ccc}0& & 1\\ 1& & 0\end{array}\right].$$
(29)
Hence, if the noise was diagonal in the canonical form, the transformed (but equivalent) noise matrix becomes

$$B=\left[\begin{array}{ccc}\lambda _{11}\mathrm{cos}(g^d)& & \lambda _{11}\mathrm{sin}(g^d)\\ \lambda _{22}\mathrm{sin}(g^d)& & \lambda _{22}\mathrm{cos}(g^d)\end{array}\right].$$
(30)
Now, the $`2M`$-dimensional (complex) orthogonal matrix family contains $`M(2M1)`$ free complex parameters, so there are $`M(2M1)`$ diffusion gauge functions $`g_{ij}(\stackrel{}{\alpha },t)`$ that one can choose arbitrarily. This represents a large class of specific gauges that can be used directly in simulations, as opposed to the conditions on noise correlations usually given elsewhere Plimak .

As pointed out by Graham Path-int , there is a close similarity between the theory of curved-space metrics, and path integrals with a space-varying diffusion matrix. In the present context, the space is complex, and we have a family of gauges that are generated on taking the matrix square root of the diffusion matrix. We have not yet used this matrix square root, but this decomposition will be applied to obtain positive-definite equations via the choice of analytic derivatives made in the following sections.

The above holds for square noise matrices $`B_s`$, but one is also free to add more noise coefficients in the manner $`B_Q=[B_s,Q]`$. Then

$$B_sB_s^T=\stackrel{~}{D}=DQQ^T,$$
(31)
and all the $`2MW`$ coefficients in the $`2M\times W`$ matrix $`Q`$ are additional arbitrary complex functions. The freedom in $`B_s`$ is the same as before \[i.e. $`M(2M1)`$ independent complex gauge functions\], with the proviso that $`B_s`$ is now given by $`O\stackrel{~}{\lambda }U`$ where the square of $`\stackrel{~}{\lambda }`$ gives the eigenvalues of the modified matrix $`\stackrel{~}{D}`$. The matrix $`B_s`$ would be unchanged if $`QQ^T`$ were set to zero, although this choice of $`Q`$ does not appear to be useful; it just adds extra noise. In general it is not clear whether or not any advantage can be gained by introducing the additional off-square gauge functions contained in $`Q`$.

If $`B`$ is given a functional form dependent on the phase-space variables, it may lead to additional terms in the Stratonovich form of the equations, which are considered later in this section. In this situation one must be careful not to introduce additional boundary-term errors arising from an excessively rapid growth of the noise gauges.

There is a subtlety here which one must take some care with. The complex noise matrix $`B`$ is not the matrix that usually appears in the theory of stochastic equations. Instead, this matrix is subsequently transformed into an ‘equivalent’ stochastic form, by taking advantage of the analyticity of the Bargmann states. This means that the effect of the diffusion gauges on the final equations also makes use of the non-uniqueness of the coherent basis set itself.

### III.2 Drift gauges

While the diffusion gauges can control sampling error due to the correlations of noise terms, they cannot eliminate boundary terms due to singular trajectories in the drift equations. The extra variable $`\mathrm{\Omega }`$ allows the $`_\mathrm{\Omega }`$ identity to be used to convert any potential term $`V`$ to a derivative term, and also to introduce a stochastic gauge to stabilize the resulting drift equations. This defines an infinite class of formally equivalent Fokker-Planck equations, in a similar way to related procedures in QED and QCD. To demonstrate this, we introduce $`2M`$ arbitrary complex drift gauge functions $`𝐠=[g_i(\stackrel{}{\alpha },t)]`$, to give a new differential operator $`_{GA}`$ whose form differs from the original $`_A^{(+)}`$ by terms that vanish identically when applied to the kernel $`\widehat{\mathrm{\Lambda }}(\stackrel{}{\alpha })`$,

$$_{GA}=_A^{(+)}+\left[V+\frac{1}{2}𝐠𝐠\mathrm{\Omega }_\mathrm{\Omega }+g_kB_{jk}_j\right]\left[\mathrm{\Omega }_\mathrm{\Omega }1\right].$$
(32)
The total differential operator $`_{GA}`$ has an anti-normal Fokker-Planck form. Extending the drift and diffusion matrices to include the extra variable $`\mathrm{\Omega }`$, we can write this — summing repeated $`a,b,c`$ indices over $`a=0,\mathrm{},2M`$ — as

$$_{GA}=\left[A_a_a+\frac{1}{2}D_{ab}_a_b\right].$$
(33)
The *total* complex drift vector is $`\stackrel{}{A}=(A_0,A_1,\mathrm{},A_{2M})`$; where

$`A_0`$ $`=`$ $`\mathrm{\Omega }V`$
$`A_j`$ $`=`$ $`A_j^{(+)}g_kB_{jk}.`$ (34)
The new diffusion matrix $`\underset{¯}{\underset{¯}{D}}`$ with elements $`D_{ab}`$ is not diagonal, but it can be factorized. Explicitly, it is now a square $`(2M+1)\times (2M+1)`$ complex matrix, given by

$`\underset{¯}{\underset{¯}{D}}`$ $`=`$ $`\left[\begin{array}{cc}\mathrm{\Omega }^2\mathrm{𝐠𝐠}^T& \mathrm{\Omega }𝐠B^T\\ B𝐠^T\mathrm{\Omega }& BB^T\end{array}\right]`$ (37)
$`=`$ $`\left[\begin{array}{cc}0& \mathrm{\Omega }𝐠\\ 0& B\end{array}\right]\left[\begin{array}{cc}0& 0\\ \mathrm{\Omega }𝐠^T& B^T\end{array}\right]=\underset{¯}{\underset{¯}{B}}\underset{¯}{\underset{¯}{B}}^T.`$ (42)
Thus, we now have a new stochastic noise matrix with one added dimension,

$$\underset{¯}{\underset{¯}{B}}=\left[\begin{array}{cc}0& \mathrm{\Omega }𝐠\\ 0& B\end{array}\right].$$
(43)
The operator (32) was chosen to give this form for $`\underset{¯}{\underset{¯}{B}}`$, so that the only change in noise is for the $`\mathrm{\Omega }`$ variable.

### III.3 Positive-definite diffusion

It is always possible to transform these second-derivative terms into a positive semi-definite diffusion operator on a real space, which is a necessary requirement for a stochastic equation. When $`\underset{¯}{\underset{¯}{D}}=\underset{¯}{\underset{¯}{B}}\underset{¯}{\underset{¯}{B}}^T`$, divide $`\underset{¯}{\underset{¯}{B}}=\underset{¯}{\underset{¯}{B}}^x+i\underset{¯}{\underset{¯}{B}}^y`$ into its real and imaginary parts. A similar procedure is followed for $`\stackrel{}{A}`$.

Recalling that the original kernel was analytic, thus allowing for more than one choice of derivatives, the choice for $`_a`$ can now be made definite by choosing it so that the resulting drift and diffusion terms are always real,

$`A_a_a`$ $``$ $`A_a^x_a^x+A_a^y_a^y,`$ (44)
$`D_{ab}_a_b`$ $``$ $`B_{ac}^xB_{bc}^x_a^x_b^x+B_{ac}^yB_{bc}^x_a^y_b^x+(xy).`$
Hence, the gauge differential operator can now be written explicitly as

$$_{GA}=\left[\stackrel{~}{A}_\mu _\mu +\frac{1}{2}\stackrel{~}{D}_{\mu \nu }_\mu _\nu \right],$$
(45)
where the indices $`\mu ,\nu `$ cover the $`(4M+2)`$-dimensional phase-space of the real and imaginary parts of $`\stackrel{}{\alpha }`$ , so that $`\stackrel{~}{\alpha }=(\stackrel{}{x},\stackrel{}{y})`$, and $`_\mu =/\stackrel{~}{\alpha }_\mu `$. The diffusion matrix $`\underset{¯}{\overset{~}{\underset{¯}{D}}}=\underset{¯}{\overset{~}{\underset{¯}{B}}}\underset{¯}{\overset{~}{\underset{¯}{B}}}^T`$ is now positive semi-definite, since, by construction

$$\underset{¯}{\overset{~}{\underset{¯}{B}}}=\left[\begin{array}{cc}0& \underset{¯}{\underset{¯}{B}}^x\\ 0& \underset{¯}{\underset{¯}{B}}^y\end{array}\right].$$
(46)
so that the diffusion matrix is the square of a real matrix — explicitly,

$$\underset{¯}{\overset{~}{\underset{¯}{D}}}=\left[\begin{array}{cc}0& \underset{¯}{\underset{¯}{B}}^x\\ 0& \underset{¯}{\underset{¯}{B}}^y\end{array}\right]\left[\begin{array}{cc}0& 0\\ \left(\underset{¯}{\underset{¯}{B}}^x\right)^T& \left(\underset{¯}{\underset{¯}{B}}^y\right)^T\end{array}\right].$$
(47)
As $`_{GA}`$ is now explicitly real as well as positive-definite by construction, it can be applied to the Hermitian conjugate kernel as well, resulting in the final time-evolution equation,

$`{\displaystyle \frac{\widehat{\rho }}{t}}`$ $`=`$ $`{\displaystyle G(\stackrel{~}{\alpha })\left[_{GA}\widehat{\mathrm{\Lambda }}(\stackrel{~}{\alpha })\right]d^{\mathrm{\hspace{0.17em}4}M+2}\stackrel{~}{\alpha }}.`$ (48)
On integrating by parts, provided boundary terms vanish, at least one solution will satisfy the following (normally-ordered) positive-definite Fokker-Planck equation — with the differential operators on the left, each acting on all terms to the right,

$$\frac{G}{t}=_{GN}G\left[_\mu \stackrel{~}{A}_\mu +\frac{1}{2}_\mu _\nu \stackrel{~}{D}_{\mu \nu }\right]G.$$
(49)
This implies that we have an equivalent set of Ito stochastic differential equations available, with $`2M`$ *real* Gaussian noises $`dW_i`$ , which are

$`d\mathrm{\Omega }`$ $`=`$ $`\mathrm{\Omega }(Vdt+g_kdW_k)`$
$`d\alpha _j`$ $`=`$ $`(A_j^{(+)}g_kB_{jk})dt+B_{jk}dW_k.`$ (50)
The noises obey $`dW_idW_j=\delta _{ij}dt`$, and are uncorrelated between time steps.

Numerical simulations are usually done in the Stratonovich calculus, due to superior convergence properties MD-semimp , so the equivalent *complex* Stratonovich equation allows us to write efficient algorithms,

$`d\alpha _a`$ $`=`$ $`dx_a+idy_a`$ (51)
$`=`$ $`\left[A_a{\displaystyle \frac{1}{2}}\left(B_{bk}_b\right)B_{ak}\right]dt+B_{ak}dW_k,`$
where $`\left(B_{bk}_b\right)\left(B_{bk}^x_b^x+B_{bk}^y_b^y\right)`$. The derivative terms above are the Stratonovich correction in the drift, corresponding to related terms obtained in curved-space path integrals.

These gauge terms are now utilized to stabilize coherent-state paths entering into highly non-classical regions of phase space. This allows one to benefit from the over-completeness of coherent states, in reducing the sampling error and eliminating boundary terms.

### III.4 Moments

The procedure for calculating observable moments is slightly different for the gauge representation than for the positive P. Any moment can be written in terms of the normally ordered operator products $`\widehat{a}^n\widehat{a}^m`$, and their expectation values are given by

$$\widehat{a}^n\widehat{a}^m_{\text{quant}}=\frac{\beta ^n\alpha ^m\mathrm{\Omega }+(\alpha ^n\beta ^m\mathrm{\Omega })^{}_{\text{stoch}}}{\mathrm{\Omega }+\mathrm{\Omega }^{}_{\text{stoch}}}.$$
(52)
which differs from the positive-P situation whenever $`\mathrm{\Omega }`$ differs from unity.

The average norm $`\mathrm{\Omega }`$ is always preserved if there is no potential term ($`V=0`$ ), since the resulting equation for the weight variable is

$$d\mathrm{\Omega }=\mathrm{\Omega }g_kdW_k.$$
(53)
The decorrelation property of Ito equations Arnold  then implies that

$$d\mathrm{\Omega }=\mathrm{\Omega }g_kdW_k=0.$$
(54)
### III.5 Gauge properties

We turn briefly here to the question of gauge classification and properties. Just as in QED, the over-complete nature of the coherent-state expansion means that many equivalent, stable gauges exist. However, they may not be equivalent in terms of boundary terms. These are determined by the tails of the distribution function, which depends intimately on the gauge chosen for the time evolution. It is essential that the distribution tails are sufficiently bounded to eliminate boundary terms arising in partial integration. It is sufficient to bound tails better than any inverse power law, for which it is conjectured to require (as a necessary condition) that all deterministic trajectories are bounded over any finite time interval GGD-Validity . This issue is discussed in greater detail below, and in Ref. gauge\_paper .

The main criteria for a useful gauge are the elimination of boundary terms and the reduction of sampling error. However, there is an enlarged space of variables for the Fokker-Planck equation here. For this reason, it is possible to stabilize trajectories in the usual positive-P phase space, while introducing new gauge-induced boundary terms in the $`\mathrm{\Omega }`$ space. When it comes to the formation of boundary terms, the phase of $`\mathrm{\Omega }`$ is generally innocuous provided the gauge is periodic in this variable, but the gauge distribution must be strongly bounded as $`|\mathrm{\Omega }|\mathrm{}`$ to prevent new boundary terms from arising.

We can classify gauges according to their real or imaginary nature, and their functional dependence; which can be on just the phase-space variables, just the quantum phase, or on both. This gives rise to nine gauge types, depending on the following criteria.

#### III.5.1 Gauge complexity

Gauges are in general complex functions, which leads to the following classification of gauge complexity:

1. Real gauge
2. Imaginary gauge
3. Complex gauge

In general, we find that trajectories can be stabilized by real, imaginary or complex gauges, provided they have some $`(𝜶,𝜷)`$ phase-space dependence.

It is worthwhile to note that the imaginary and real parts of the gauges affect the behavior of sampling error differently. In the Ito calculus, the evolution of the weight $`\mathrm{\Omega }`$ due to the gauges is simply $`d\mathrm{\Omega }=\mathrm{\Omega }g_kdW_k`$. Typically, i.e., when there are no significant correlations between the phase of $`𝜶`$ ( or $`𝜷`$) and $`\mathrm{\Omega }`$, the weight factor appearing in moment calculations is just approximately $`\text{Re}[\mathrm{\Omega }]`$. As a general rule, sampling errors are partially due to stochastic fluctuations in the phase-space trajectories, and partially due to stochastic fluctuations in the weight function. Thus there is a trade-off; a gauge that is strongly stabilizing may reduce phase-space fluctuations at the expense of increased weight variance, and vice versa.

To understand the different types of gauges in somewhat greater detail, we consider the evolution of the weight variance for real and imaginary gauges, in a simple case where gauge and weight are decorrelated, with $`\mathrm{\Omega }=1`$ initially. Let $`\mathrm{\Omega }=\mathrm{\Omega }^{}+i\mathrm{\Omega }^{\prime \prime }`$ and $`g_k=g_k^{}+ig_k^{\prime \prime }`$ , then

$`d\mathrm{\Omega }^{}`$ $`=`$ $`\left(\mathrm{\Omega }^{}g_k^{}\mathrm{\Omega }^{\prime \prime }g_k^{\prime \prime }\right)dW_k,`$
$`d\mathrm{\Omega }^{\prime \prime }`$ $`=`$ $`\left(\mathrm{\Omega }^{}g_k^{\prime \prime }+\mathrm{\Omega }^{\prime \prime }g_k^{}\right)dW_k.`$ (55)
If we consider the evolution of the squares of these terms, the Ito rules of stochastic calculus give

$`d[\mathrm{\Omega }^{}]^2`$ $`=`$ $`\left(\mathrm{\Omega }^{}g_k^{}\mathrm{\Omega }^{\prime \prime }g_k^{\prime \prime }\right)^2dt,`$
$`d[\mathrm{\Omega }^{\prime \prime }]^2`$ $`=`$ $`\left(\mathrm{\Omega }^{}g_k^{\prime \prime }+\mathrm{\Omega }^{\prime \prime }g_k^{}\right)^2dt.`$ (56)
Suppose for simplicity that the $`g_k`$ and $`\mathrm{\Omega }`$ are approximately uncorrelated, then we have two cases to consider.

1. Real gauge:

$`d[\mathrm{\Omega }^{}]^2`$ $`=`$ $`[\mathrm{\Omega }^{}]^2d\tau ,`$ (57)
where $`d\tau =g_kg_kdt`$. This initially leads to linear growth in the variance, and hence in the sampling error. The real part of the gauge will cause noise directly in $`\mathrm{\Omega }^{}`$, producing asymetric spreading in $`\mathrm{\Omega }^{}`$, which can lead to a few rare very highly weighted trajectories for times $`\tau 1`$. The effect of the real gauge may become misleading once the distribution becomes highly skewed, as the rare trajectories that are important for moment calculations may be missed if the sample is too small. At long times, if $`g_kg_k`$ is constant and uncorrelated with $`\mathrm{\Omega }`$, then the growth becomes exponential, with $`[\mathrm{\Omega }^{}]^2=e^\tau `$.
2. Imaginary gauge:

$`d[\mathrm{\Omega }^{}]^2`$ $`=`$ $`[\mathrm{\Omega }^{\prime \prime }]^2d\tau ,`$
$`d[\mathrm{\Omega }^{\prime \prime }]^2`$ $`=`$ $`[\mathrm{\Omega }^{}]^2d\tau .`$ (58)
where $`d\tau =g_k^{\prime \prime }g_k^{\prime \prime }dt`$. This leads initially to quadratic growth in the variance of $`\mathrm{\Omega }^{}`$, and hence a slower growth in the sampling error. If $`g_kg_k`$ is constant and remains uncorrelated with $`\mathrm{\Omega }`$, then the growth is given by $`[\mathrm{\Omega }^{}]^2=\mathrm{cosh}(\tau )`$, $`[\mathrm{\Omega }^{\prime \prime }]^2=\mathrm{sinh}(\tau )`$. An imaginary gauge will cause mutual canceling of trajectories that have weights of randomly positive and negative sign once $`\tau \pi `$. This can also have deleterious effects for small samples, if the average sample weight becomes negative — of course, this cannot be true over the entire stochastic population.

The generic behavior is more complex than in the examples given above, due to correlations between the gauge and the normalization.

Clearly any type of gauge tends to cause growth in the norm variance. However, there is an exception to this rule: the norm-preserving gauges. This class of gauges is of special interest as they generate trajectories having an invariant normalization, so that $`\text{Re}[d\mathrm{\Omega }]0`$. From the equation for the norm variance, Eq. (56), it follows that a necessary and sufficient condition for a norm-preserving gauge is that $`\mathrm{\Omega }^{}g_k^{}=\mathrm{\Omega }^{\prime \prime }g_k^{\prime \prime }`$. If $`\mathrm{\Omega }^{}=1`$ initially, this implies that $`g_k=i\mathrm{\Omega }^{}f_k=i(1i\mathrm{\Omega }^{\prime \prime })f_k`$ , where $`f_k`$ is a real function. Unless $`g_k=0`$, norm-preserving gauges are generally functions of both the phase-space variables and the weight $`\mathrm{\Omega }`$ . A preliminary study of these gauges has shown that these gauges can greatly reduce sampling error, although gauge-induced boundary terms are also possible ccp2k , depending on the choice of $`f_k`$.

#### III.5.2 Functional dependence

From the above analysis, we see that gauges can functionally depend on any phase-space variable, as well as the generalized quantum phase variable or weight $`\mathrm{\Omega }`$ . This leads to three functional types:

1. Autonomous (depends on $`\mathrm{\Omega }`$ only)
2. Space dependent (depends on phase-space only)
3. Mixed (depends on all components of $`\stackrel{}{\alpha }`$ including $`\mathrm{\Omega }`$ )

Autonomous gauges appear to be the least useful since they do not affect $`𝜶`$ or $`𝜷`$ behavior, but gauges of either purely space-dependent or mixed type can be used.

A possible caveat with mixed gauges is that they may be much harder to analyze, as two-way couplings will occur between the normal phase-space variables $`𝜶`$, $`𝜷`$ and the weight.

## IV Nonlinear Absorber case

The nonlinear absorber is an example of a nonlinear master equation that can give either correct or incorrect results when treated with the usual positive-P representation methods, if the boundary terms are ignored. Generally, problems only arise when the linear damping has exceptionally small values or the number of bosons per mode is small (see Fig. 2), so this is not a practical problem in optics. However, for other physical systems such as a BEC this may be significant. It is a well-studied case, and a detailed treatment can be found in Ref. GGD-Validity . It also has the merit that exact solutions can be readily found using other means. By analyzing this example we can ensure that the modifications to the drift equations obtained from gauge terms, do eliminate boundary terms and give correct results.

Consider a cavity mode driven by coherent radiation, and damped by a zero temperature bath that causes both one and two photon losses. We have scaled time so that the rate of two-photon loss is unity. Without this nonlinear process, nothing unusual happens. The scaled one-photon loss rate is $`\gamma `$, and $`\epsilon `$ is the scaled (complex) driving field amplitude. The master equation is

$`{\displaystyle \frac{\widehat{\rho }}{t}}`$ $`=`$ $`[\epsilon \widehat{a}^{}\epsilon ^{}\widehat{a},\widehat{\rho }]+{\displaystyle \frac{\gamma }{2}}(2\widehat{a}\widehat{\rho }\widehat{a}^{}\widehat{a}^{}\widehat{a}\widehat{\rho }\widehat{\rho }\widehat{a}^{}\widehat{a})`$ (59)
$`+{\displaystyle \frac{1}{2}}(2\widehat{a}^2\widehat{\rho }\widehat{a}^2\widehat{a}^2\widehat{a}^2\widehat{\rho }\widehat{\rho }\widehat{a}^2\widehat{a}^2).`$
Following the treatment of Sec. II, we arrive at the gauge representation Stratonovich stochastic equations

$`d\alpha `$ $`=`$ $`[\epsilon ^{}\alpha (\alpha \beta +ig+(\gamma 1)/2)]dt+i\alpha dW,`$
$`d\beta `$ $`=`$ $`[\epsilon ^{}\beta (\alpha \beta +i\overline{g}+(\gamma 1)/2)]dt+i\beta d\overline{W},`$
$`d\mathrm{\Omega }`$ $`=`$ $`S_\mathrm{\Omega }dt+\mathrm{\Omega }\left[gdW+\overline{g}d\overline{W}\right].`$ (60)
Here $`S_\mathrm{\Omega }dt`$ is the appropriate Stratonovich correction term \[given by the derivative terms in Eq. (51) \], which depends on the particular gauges chosen.

With no gauge ($`g=\overline{g}=0`$), the positive-P Stratonovich equations are recovered,

$`d\alpha `$ $`=`$ $`[\epsilon \alpha (\alpha \beta +\{\gamma 1\}/2)]dt+i\alpha dW,`$
$`d\beta `$ $`=`$ $`[\epsilon ^{}\beta (\alpha \beta +\{\gamma 1\}/2)]dt+i\beta d\overline{W}.`$ (61)
We will concentrate on the various simplifications of this model, which correspond to existing literature, and simpler analysis.

### IV.1 Relevance to many-body problems

The nonlinearity seen here can occur directly in the form of a nonlinear collisional damping term in a many-body system, so that it can be referred to generically as ‘two-boson absorption’. This type of damping is common both to nonlinear photonic and atomic interactions.

It is of nearly the same form as for an ‘imaginary-time’ thermal equilibrium state calculation for the usual model of an alkali-metal Bose gas or BEC BECH . There, for example, the interaction energy between identical bosons of mass $`m`$ and $`s`$-wave scattering length $`a_s`$ in $`D`$-dimensional space is given by

$$\widehat{H}=\frac{2\pi \mathrm{}^2a_s}{m}d^D𝐱\widehat{\psi }^2(𝐱)\widehat{\psi }^2(𝐱),$$
(62)
provided that $`a_s`$ is much smaller than other characteristic lengths of the system (which is usually the case). The master equation for an imaginary-time calculation is

$$\frac{\widehat{\rho }_e}{\tau }=\frac{1}{2}\{\widehat{H}\mu \widehat{N},\widehat{\rho }_e\}_+,$$
(63)
where $`\widehat{\rho }_e`$ is the thermal canonical ensemble density matrix, $`\mu `$ is the chemical potential, $`N`$ is the number operator for the entire system, and $`\tau =1/k_BT`$ is an inverse temperature. Apart from the fact that it is not trace-preserving, this is a nonlinearity very similar to that occurring in the nonlinear absorber master equation.

While boundary-term discrepancies only occur with this nonlinearity for low occupations per mode (see also Fig. 2), for a many-mode system at finite temperature one expects a large number of modes to have just such a low occupation. Thus, it is important to check that boundary terms are indeed eliminated. Note that the gauge representation simulation is efficient over a wide range of occupation numbers. See, for example, Fig. 3. More details of applications to both real and imaginary time many-body systems with many modes will be given elsewhere.

### IV.2 Two-boson absorber

In its simplest form, corresponding to $`\gamma =\epsilon =0`$, only two-boson absorption takes place. We expect that for a state $`|\psi =_nc_n|n`$ all even boson number components will decay to vacuum, and all odd-numbered components will decay to $`|1`$, leaving a mixture of vacuum and one-boson states at long times.

The positive-P representation has been found to give erroneous results SG-abs ; CM-abs ; MC-abs ; GGD-using  due to the existence of moving singularities GGD-Validity , which cause power-law tails in the distribution leading to boundary terms. The moment usually concentrated on in this system is the number of bosons $`\widehat{n}=\widehat{a}^{}\widehat{a}`$, which corresponds to the statistical average of $`n=\alpha \beta `$ in the positive-P representation. This has a convenient closed equation (Stratonovich),

$$dn=n(n+i\stackrel{~}{g}1/2)d\tau +indW^+$$
(64)
with $`dW^+=(dW+d\overline{W})`$, $`\tau =2t`$, and $`\stackrel{~}{g}=(g+\overline{g})/2`$.

Let us examine the behavior of the above equation, when $`\stackrel{~}{g}=0`$, i.e., in the standard, un-gauged formulation. The deterministic part of the evolution has a repellor at $`n=0`$, and an attractor at $`n=\frac{1}{2}`$. The noise is finite, and of standard deviation $`\sqrt{dt/2}`$ at the attractor. We can see that the deterministic part of the evolution has a single trajectory of measure zero which can escape to infinity along the negative real axis,

$$\alpha =\beta =\frac{1}{\sqrt{\tau _0\tau }},$$
(65)
where $`\tau _0=1/\alpha (0)^2=1/n(0)`$. This moving singularity is known to cause the power law behavior of the Fokker-Planck solution at large $`|n|`$, which means that integration by parts is not in fact valid — which leads to incorrect results.

Indeed, it can be easily seen that in the steady-state limit, all trajectories in a simulation will head toward $`n=\frac{1}{2}`$, making $`lim_t\mathrm{}\widehat{n}=\frac{1}{2}`$. Quantum mechanics, however, predicts that if we start from a state $`\widehat{\rho }_0`$, the steady state will be

$$\underset{t\mathrm{}}{lim}\widehat{n}=\underset{j=0}{\overset{\mathrm{}}{}}1+2j\left|\widehat{\rho }_0\right|1+2j.$$
(66)
For a coherent state $`|\alpha _0`$ input, say, this will be

$$\underset{t\mathrm{}}{lim}\widehat{n}=\frac{1}{2}\left(1e^{2|\alpha _0|^2}\right).$$
(67)
Thus we can expect that the positive P simulation will give correct results only when $`e^{|\alpha _0|^2}1`$.

To correct the problem we have to change the phase-space topology in some way to prevent the occurrence of moving singularities. We have found that a good gauge for a two-boson absorber nonlinearity in general is

$$g=\overline{g}=\stackrel{~}{g}=i(n|n|).$$
(68)
This replaces the $`n^2`$ term in Eq. (64), which may become repulsive from zero, with $`n|n|`$ which is always a restoring force, and so never leads to super-exponential escape.

With the gauge (68), the Stratonovich equations become

$`dn`$ $`=`$ $`n(|n|1/2)d\tau +indW^+,`$ (69)
$`d\mathrm{\Omega }`$ $`=`$ $`\mathrm{\Omega }\left\{[n+(n|n|)^2]d\tau /2+i(n|n|)dW^+\right\}.`$
Phase-space trajectories have changed now, but since it has all come from the same master equation, it still describes the same system. Consider the equations for the polar decomposition of $`n=re^{i\varphi }`$,

$`dr`$ $`=`$ $`r(r1/2)d\tau ,`$
$`d\varphi `$ $`=`$ $`dW^+.`$ (70)
This is exact, and shows that now we have an attractor on the circle $`|n|=\frac{1}{2}`$, and a repellor at $`n=0`$, with free phase diffusion in the tangential direction. Once trajectories reach the attractor, only phase diffusion occurs. Some more complicated evolution is occurring in the $`\mathrm{\Omega }`$ variable. In any case, there are now no moving singularities anywhere in the phase space, and simulations correspond exactly to quantum mechanics.

Figure 1 compares results for a truncated number-state basis calculation, a positive-P calculation, and a “circular” gauge (68) calculation for an initial coherent state of $`\alpha _0=1/\sqrt{2}`$. Figure 2 compares steady-state values for exact, positive-P, and gauge calculations for various initial coherent states in a wide range. It is seen that the gauge calculation is correct to within the small errors due to finite sample size.

### IV.3 One- and two-boson absorber

If we now turn on the one-boson decay as well, but still do not have any driving, we expect that all states will decay to the vacuum on two time scales $`1`$ and $`1/\gamma `$. If $`\gamma 1`$, nothing interesting happens, however if $`\gamma 1`$, we should first see a rapid decay to a mixture of vacuum and one-boson states due to the two-boson process, and then a slow decay of the one-boson state to the vacuum on a time scale of $`\tau 2/\gamma `$.

In this case the positive-P equations display different behavior depending on whether $`\gamma `$ is above or below the threshold $`\gamma =1`$. Below threshold, we have an attractor at $`n=(1\gamma )/2`$, and a repellor at $`n=0`$, while above threshold, the attractor is at $`n=0`$, and the repellor at $`n=(\gamma 1)/2`$. In either case, there is a singular trajectory along the negative real axis, which can cause boundary-term errors. It turns out that the steady state calculated this way is erroneous while $`\gamma <1`$, and there are transient boundary-term errors while $`\gamma <2`$ SG-abs . The false steady state below threshold lies at the location of the attractor: $`(1\gamma )/2`$.

Let us try to fix this problem using the same circular gauge (68) as before. The equation for $`r`$ is now

$$dr=r(r[1\gamma ]/2)d\tau ,$$
(71)
while the $`\varphi `$ and $`\mathrm{\Omega }`$ evolution is unchanged. So, above threshold we are left with only an attractor at $`n=0`$, while below threshold we have a repellor at $`n=0`$ surrounded by an attracting circle at $`r=(1\gamma )/2`$. This phase space again has no moving singularities.

The results of simulations for the parameter $`\gamma =0.1`$ are shown in Fig. 3. The gauge simulation tracks the exact results. We have chosen $`\gamma 1`$ so that a system with two widely differing time scales is tested. The circular gauge avoids the false results of the positive-P simulation. Note also that the gauge simulation remains efficient for a wide range of occupation numbers — from $`\widehat{n}1001`$, where the positive P is also accurate, to $`\widehat{n}0.11`$ where it is totally incorrect.

### IV.4 Driven two-boson absorber

The other type of situation to consider is when we have a driving field as well as two-boson damping. In these considerations we have set the one-boson damping rate to zero ($`\gamma =0`$), since this process never causes any of the simulation problems anyway, but leaving it out simplifies analysis. Failure of the positive-P representation method has been found in this limit as well SS-fail , and is evident in Fig. 4. The equation for $`n`$ is no longer stand-alone in this case, and we must simulate all three complex variables as in Eq. (60), the $`\mathrm{\Omega }`$ equation being the same as in the undriven case (69).

A treatment of the singular trajectory problem with the same circular gauge (68) leads again to correct results, as seen in Fig. 4.

## V The single-mode laser

Let us now consider the second quantum system for which systematic errors have been seen with the positive-P representation. We will see that the problem here is somewhat different than in the previous case. The difference is that for two-boson damping, boundary-term errors occur even when we choose an optimal (i.e., compact) initial distribution to represent our starting state, whereas here systematic errors occur only for unreasonably broad initial distributions. Nevertheless, since normally it is assumed that the initial condition can be of arbitrary breadth it is instructive to investigate how this problem can be tackled with stochastic gauge methods.

We have found that stochastic gauges can be used to increase the allowable breadth to include all reasonable starting conditions, but once one tries to increase the initial spread too much, it becomes unlikely that any gauge will remove systematic errors, without introducing too much sampling (i.e. random) error instead.

### V.1 The laser model

Ito stochastic differential equations for a simple photonic or atomic laser model that can be derived from the positive-P distribution are SS-fail ; GGD-Validity

$`d\stackrel{~}{\alpha }`$ $`=`$ $`(G\stackrel{~}{\alpha }\stackrel{~}{\beta })\stackrel{~}{\alpha }d\tau +\sqrt{Q}d\eta ,`$
$`d\stackrel{~}{\beta }`$ $`=`$ $`(G\stackrel{~}{\alpha }\stackrel{~}{\beta })\stackrel{~}{\beta }d\tau +\sqrt{Q}d\eta ^{}`$ (72)
in appropriate scaled variables, with the complex Gaussian noise $`d\eta `$ obeying $`d\eta d\eta ^{}=2d\tau `$. In terms of physical parameters, we have

$`\stackrel{~}{\alpha }`$ $`=`$ $`\alpha /\sqrt{𝒩}`$
$`\stackrel{~}{\beta }`$ $`=`$ $`\beta /\sqrt{𝒩},`$ (73)
where $`\tau `$ is the scaled time, and $`𝒩1`$ is a scaling parameter that equals the number of gain atoms in a simple photonic laser model. Both $`G`$ the gain parameter and $`QG/𝒩`$, the noise parameter, are real and positive.

Since this time we are again interested in the (scaled) boson number $`\stackrel{~}{n}=\stackrel{~}{\alpha }\stackrel{~}{\beta }=\widehat{n}/𝒩`$, its evolution can be written as a closed equation

$$d\stackrel{~}{n}=2(\stackrel{~}{n}a)(\stackrel{~}{n}b)d\tau +2\sqrt{Q\stackrel{~}{n}}dW,$$
(74)
where now the real Gaussian noise obeys $`dWdW=d\tau `$, and the deterministic stationary points in the Stratonovich calculus are

$`a`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left(G+\sqrt{G^2+2Q}\right),`$
$`b`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left(G\sqrt{G^2+2Q}\right).`$ (75)
We find that the stationary point at $`a`$ is an attractor, and at $`b`$ we have a repellor. Defining $`\mathrm{\Delta }=b\stackrel{~}{n}`$, we get

$$d\mathrm{\Delta }=2\mathrm{\Delta }(\mathrm{\Delta }+\sqrt{G^2+2Q})+\text{noise},$$
(76)
which shows that we again have a singular trajectory escaping to infinity in finite time along the negative real axis for $`\stackrel{~}{n}<b`$.

### V.2 Initial conditions

Let us consider the usual case of vacuum initial conditions. A vacuum can be represented by

$$P^{(+)}(\stackrel{~}{\alpha },\stackrel{~}{\beta })=\delta (\stackrel{~}{\alpha })\delta (\stackrel{~}{\beta }),$$
(77)
but also by Gaussian distributions of any variance $`\sigma _0^2`$, around the above,

$$P^{(+)}(\stackrel{~}{\alpha },\stackrel{~}{\beta })=\frac{1}{4\pi ^2\sigma _0^4}\mathrm{exp}\left\{\frac{|\stackrel{~}{\alpha }|^2+|\stackrel{~}{\beta }|^2}{2\sigma _0^2}\right\}.$$
(78)
Note: the distribution of $`\stackrel{~}{n}`$ is non-Gaussian, but has a standard deviation of $`\sigma _{\stackrel{~}{n}}\sqrt{2}\sigma _0^2`$ in both the real and imaginary directions.

It has been found by Schack and Schenzle SS-fail  that for the single-mode laser model, a positive-P simulation of pumping from a vacuum will give correct answers if the usual $`\delta `$-function initial condition (77) is used, but will have systematic errors if the initial condition used has a sufficiently large variance (see Fig. 5). We emphasize here that this is not a real problem in practical cases, as the variance required to cause systematic errors is typically extremely large, once the scaling needed to obtain the usual (approximate) laser model is taken into account.

This can be understood because if we have a sufficiently broad initial distribution, the region of phase space that includes the singular trajectory will be explored by the distribution. Even if initially $`\sigma _{\stackrel{~}{n}}|b|`$, the region $`\stackrel{~}{n}<b`$ may be subsequently explored due to the presence of the noise terms.

Apart from the obvious $`\delta `$-function initial condition, one might want to try the canonical distribution of Eq. (10), which is a standard positive-P representation construction DG-PosP . It will not cause problems as its variance is $`\sigma _0^2=1/𝒩`$, which for any realistic case will be very small (i.e., $`\sigma _{\stackrel{~}{n}}|b|`$). Schack and Schenzle discovered anomalous results when they chose $`\sigma _0^2=1`$, due to an erroneous procedure of scaling the equations — while not scaling the canonical initial condition in $`\alpha `$. Nevertheless, since any $`\sigma _0`$ is supposed to represent the same state, insight into what can be achieved using gauge methods is gained if we analyze the systematic errors for such a relatively large $`\sigma _0`$.

### V.3 Gauge corrections

The Fokker-Planck equation corresponding to Eq. (72) is

$$\frac{P}{\tau }=\left\{\frac{}{\stackrel{~}{\alpha }}[\stackrel{~}{n}G]\stackrel{~}{\alpha }+\frac{}{\stackrel{~}{\beta }}[\stackrel{~}{n}G]\stackrel{~}{\beta }+2Q\frac{^2}{\stackrel{~}{\alpha }\stackrel{~}{\beta }}\right\}P.$$
(79)
We now introduce gauges using the same method as in Sec. II. This leads to the Ito stochastic equations

$`d\stackrel{~}{\alpha }`$ $`=`$ $`\stackrel{~}{\alpha }(G\stackrel{~}{n})d\tau \sqrt{Q}(g+i\overline{g})d\tau +\sqrt{Q}d\eta ,`$
$`d\stackrel{~}{\beta }`$ $`=`$ $`\stackrel{~}{\beta }(G\stackrel{~}{n})d\tau \sqrt{Q}(gi\overline{g})d\tau +\sqrt{Q}d\eta ^{},`$
$`d\mathrm{\Omega }`$ $`=`$ $`\mathrm{\Omega }[(gi\overline{g})d\eta +(g+i\overline{g})d\eta ^{}]/2.`$ (80)
It is convenient to define a transformed gauge function $`\stackrel{~}{g}`$, which is also arbitrary, such that

$`g={\displaystyle \frac{(\stackrel{~}{\alpha }+\stackrel{~}{\beta })\stackrel{~}{g}}{2\sqrt{Q}}}`$ $`,`$
$`\overline{g}={\displaystyle \frac{(\stackrel{~}{\alpha }\stackrel{~}{\beta })\stackrel{~}{g}}{2i\sqrt{Q}}}.`$ (81)
Changing to $`\stackrel{~}{n}`$ and $`\mathrm{\Theta }=\mathrm{ln}(\mathrm{\Omega })`$ variables we obtain the Stratonovich equation

$`d\stackrel{~}{n}`$ $`=`$ $`2\stackrel{~}{n}(G\stackrel{~}{n}\stackrel{~}{g})d\tau +Qd\tau +2\sqrt{Q\stackrel{~}{n}}dW,`$
$`d\mathrm{\Theta }`$ $`=`$ $`{\displaystyle \frac{\stackrel{~}{n}\stackrel{~}{g}^2}{2Q}}d\tau +S_\mathrm{\Theta }d\tau +\stackrel{~}{g}\sqrt{{\displaystyle \frac{\stackrel{~}{n}}{Q}}}dW,`$ (82)
with $`S_\mathrm{\Theta }dt`$ being the appropriate Stratonovich correction \[given by the derivative terms in Eq. (51) \] for a particular gauge function $`\stackrel{~}{g}`$.

### V.4 Correcting for the moving singularities

Consider the deterministic evolution of the real part, $`\stackrel{~}{n}_x`$, of $`\stackrel{~}{n}=\stackrel{~}{n}_x+i\stackrel{~}{n}_y`$,

$$d\stackrel{~}{n}_x=2\stackrel{~}{n}_x^2+2G\stackrel{~}{n}_x+Q+2\stackrel{~}{n}_y^22\stackrel{~}{n}_x\text{Re}[\stackrel{~}{g}]+2\stackrel{~}{n}_y\text{Im}[\stackrel{~}{g}].$$
(83)
The moving singularity is due to the $`2\stackrel{~}{n}_x^2`$ leading term for negative values of $`\stackrel{~}{n}_x`$. We now consider criteria for choosing the drift gauges as follows.

(1) It is desirable to keep the gauge terms to a minimum because whenever they act the weights of trajectories become more randomized — see Sec. III.5.1. Thus, let us restrict ourselves to functions $`\stackrel{~}{g}`$ that are only nonzero for $`\stackrel{~}{n}_x<0`$.

(2) This immediately leads to another restriction on $`\stackrel{~}{g}`$: To be able to use the efficient numerical algorithms in the Stratonovich calculus, we must be able to calculate the correction term $`S_\mathrm{\Theta }`$, which depends on derivatives of $`\stackrel{~}{g}\sqrt{\stackrel{~}{n}/Q}`$. This immediately suggests that $`\stackrel{~}{g}`$ must always be continuous, hence, in particular, $`lim_{n_x0}(\stackrel{~}{g})=0`$. For ease of analysis, let us start with a simple form for the gauge, $`\stackrel{~}{g}=c\lambda \stackrel{~}{n}_x+\lambda _y\stackrel{~}{n}_y`$. This restriction immediately implies $`c=\lambda _y=0`$, hence

$$\stackrel{~}{g}=\{\begin{array}{cc}\lambda \text{Re}[\stackrel{~}{n}]& \text{if }\text{Re}[\stackrel{~}{n}]<0\hfill \\ 0& \text{if }\text{Re}[\stackrel{~}{n}]0\hfill \end{array},$$
(84)
and $`S_\mathrm{\Theta }=\lambda (\text{Re}[\stackrel{~}{n}]+\stackrel{~}{n}+|\stackrel{~}{n}|)/2.`$ when $`\text{Re}[\stackrel{~}{n}]<0`$, zero otherwise.

(3) The next necessary condition, to remove moving singularities, is that the $`2\stackrel{~}{n}_x^2`$ term is canceled, hence:

$$\lambda 1.$$
(85)
(4) Now, if $`\lambda =1`$ there are no systematic errors, but the sampling error very quickly obscures everything because $`n_x`$ still heads to $`\mathrm{}`$ exponentially due to the $`2G\stackrel{~}{n}_x`$ term. This takes it into regions of everincreasing $`|\stackrel{~}{g}|`$, and weights quickly become randomized. For slightly larger parameters $`\lambda `$, the $`\stackrel{~}{n}_x`$ evolution takes trajectories to a point lying far into the negative $`n_x`$ region where the two leading terms balance. Here the trajectories sit, and quickly accumulate weight noise. It is clear that for an optimum simulation all stationary points of $`\stackrel{~}{n}_x`$ in the nonzero gauge region must be removed. In this system this condition is

$$\lambda >1+\frac{G^2}{2Q}.$$
(86)
An example has been plotted in Fig. 5 where we have parameters $`G=1`$, $`Q=0.25`$ (leading to $`a1.1124`$ and $`b0.1124`$ ). We are considering an initial condition of $`\sigma _0^2=0.1`$, which is already much larger than the canonical variance for physically likely parameters. Typical values of $`\stackrel{~}{n}`$ initially will be of order $`\sigma _{\stackrel{~}{n}}0.14|b|`$ here. A good choice of gauge has $`\lambda =4`$. The use of this gauge clearly restores the correct results.

### V.5 Non-optimal initial conditions

As we increase the spread of the initial distribution beyond $`\sigma _{\stackrel{~}{n}}|b|`$, it becomes increasingly difficult to find a gauge that will give reasonable simulations. (For example we have tried a wide variety of what seemed like promising gauges for $`\sigma _0^2=0.3`$, with the above values of parameters $`Q`$ and $`G`$, and none have come close to success). The problem is that while we can remove systematic errors, large random noise appears and obscures whatever we are trying to calculate.

Trajectories that start off at a value of $`\stackrel{~}{n}`$ lying significantly beyond $`b`$ require a lot of modification to their subsequent evolution to (1) stop them from escaping to $`\mathrm{}`$ and (2) move them out of the gauged region of phase space so that they do not accumulate excessive weight noise. If there are many of these, the trade-off between the gauge size and length of time spent in the gauged region does not give much benefit anymore. Nevertheless, one may be sure that if this is the case, results will at worst be noisy and unusable, rather than being systematically incorrect.

We stress again that this whole matter of non-optimal initial conditions is not a major hurdle to dynamical simulations because a compact starting distribution is generally found very easily.

## VI Conclusions

The positive-P representation is well suited to complex quantum mechanical problems, such as many-body systems, but has been known for about a decade to have systematic errors in some cases of its use — due to non-vanishing boundary terms. The gauge P representation, a variant on the usual positive-P representation, can be used to eliminate boundary terms and consequently all the systematic errors that were encountered previously. It can also reduce sampling error in a simulation, and allows ‘imaginary time’ calculations of thermal equilibrium states. The fact that correct results are immediately obtained in every case where systematic errors were found with the positive-P method, is strong evidence that these previous problems were indeed due to boundary terms caused by moving singularities in the analytically continued deterministic equations. Of course, boundary terms can occur for other reasons (for example, if the noise term grows too rapidly with radius), so caution is still needed in the gauge choice.

The technique appears to be broadly applicable, and only requires the recognition of what instabilities in the stochastic equations could lead to problems. It does not require detailed knowledge of what the boundary terms are, provided instabilities are removed. However, we remark here that the general specification of necessary and sufficient conditions to eliminate boundary terms remains an open problem, and clearly requires growth restrictions on the gauge terms, both in phase space and quantum-amplitude space. Care is also required with the choice of the gauge and initial distribution. However, using unsuitable gauges or initial conditions may only lead to large sampling errors, not systematic errors, provided the gauge is chosen to eliminate boundary corrections in the first place. Sampling error then allows for a confident assessment of the magnitude of inaccuracies in a simulation, which can be supplemented by numerical analysis of the distribution tails.

The main conclusion we come to is that this method does, in the cases studied, provide a complete solution to the problem of simulation of a many-body quantum system in phase space, under conditions where previous direct simulation techniques were not practicable. All known technical requirements on the path to obtaining a stochastically equivalent description to quantum mechanics, which is applicable to both large and small particle numbers, have been satisfied by this method. For this reason, we believe that gauge simulations can be used to simulate many quantum systems without systematic errors when carrying out more difficult calculations, where no exact result is known.

These conclusions must be supplemented by the detailed study of relevant gauges for particular quantum systems. We note, however, that the mathematical techniques employed here for generating stochastic gauges, may well be useful for other representations as well as the gauge P representation described here.

###### Acknowledgements.

Numerical calculations were carried out using open software from the XMDS project xmds . Thanks to the Australian Research Council and the Alexander von Humboldt-Stiftung for providing research support. We thank Damian Pope for suggesting the name stochastic gauges.

## Appendix A Other extensions of the positive-P representation

### A.1 The work of Carusotto, Castin, and Dalibard

Recently, Carusotto, Castin, and Dalibard Paris1 ; Paris2  (CCD) have made related extensions to the positive-P representation. These were derived for the particular case of an interacting scalar Bose gas, and led to a number of conditions for an Ito stochastic evolution to be equivalent to a master equation evolution.

It can be shown quite simply that the equations (50) generated by the gauge P representation for this Hamiltonian satisfy the CCD conditions. We conjecture that these provide the most general possible solution to the stochastic problem posed by these authors. In particular, $`db=\mathrm{\Pi }[g_kdW_k\overline{N}(\varphi _1dB_2^{}+\varphi _2^{}dB_1)]`$, using the above paper’s formalism. Our methods can also treat a much larger class of Hamiltonians and master equations than considered in the CCD treatment.

In Ref. Paris1  systematic errors due to boundary terms were not considered. However, evolutions satisfying "exactness" conditions derived using the same procedure can contain such errors.

As an example, following the CCD procedure Paris1  for a one-mode two-boson absorber master equation, as in Eq. (59) with $`\gamma =\epsilon =0`$, one arrives at the conditions

$`dB_1dB_2^{}`$ $`=`$ $`0,`$
$`dB_a^2`$ $`=`$ $`\varphi _a^2,`$
$`F_1`$ $`=`$ $`dbdB_1/\mathrm{\Pi },`$ (87)
$`F_2`$ $`=`$ $`db^{}dB_2/\mathrm{\Pi }^{},`$
$`f`$ $`=`$ $`\mathrm{\Pi }(\overline{N}\varphi _1\varphi _2^{})^2,`$ (88)
where (referring back to the notation in this present paper),

$`d\varphi _1=`$ $`d\alpha /\sqrt{\overline{N}}=`$ $`F_1dt+dB_1,`$
$`d\varphi _2=`$ $`d\beta ^{}/\sqrt{\overline{N}}=`$ $`F_2dt+dB_2,`$
$`d\mathrm{\Pi }=`$ $`d[\mathrm{\Omega }e^{\varphi _1\varphi _2^{}\overline{N}}]=`$ $`fdt+db.`$ (89)
It can be seen that the positive-P equations (IV) satisfy these conditions, while producing the erroneous evolution seen in Fig. 1. In summary, the methods of the CCD paper do not obviate the need to choose gauges that eliminate boundary terms.

### A.2 Noise optimization by Plimak, Olsen, and Collett

In Ref. Plimak , Plimak, Olsen and Collett have found that for some systems (the Kerr oscillator $`\widehat{H}=\omega _0\widehat{a}^{}\widehat{a}+\kappa \widehat{a}^2\widehat{a}^2/2`$, in particular), the most obvious (diagonal) choice of noise matrix $`B`$ may not be the optimal one.

For example, for the above Hamiltonian, one finds that the diffusion matrix (in $`\alpha ,\beta `$) variables is

$$D=i\kappa \left[\begin{array}{ccc}\alpha ^2& & 0\\ 0& & \beta ^2\end{array}\right]=BB^T.$$
(90)
Following the procedure in Eq. (30), an equivalent but broader choice of noise matrix $`B`$ can be any of

$$B=\sqrt{i\kappa }\left[\begin{array}{ccc}i\alpha \mathrm{cos}(g)& & i\alpha \mathrm{sin}(g)\\ \beta \mathrm{sin}(g)& & \beta \mathrm{cos}(g)\end{array}\right],$$
(91)
with the usual diagonal decomposition given by $`g=0`$.

However, in Ref. Plimak  it was found that for a positive-P simulation, different decompositions with nonzero constant $`g`$ gave the lowest sampling error for coherent state initial conditions. In their notation, they introduce $`\sqrt{A+1}=\sqrt{2}\mathrm{cos}(g)`$, and consider the case of real $`A1`$ (i.e., imaginary $`g`$ ) only.

### A.3 Stochastic gauges for the Kerr oscillator

In Ref. ccp2k , the sampling error in a Kerr oscillator simulation — equivalent to a one-mode BEC model, apart from linear terms — was reduced substantially by using a representation similar to the gauge P representation formally introduced here. The basic differences were the following.

1. Instead of a complex gauge $`\mathrm{\Omega }`$, a phase factor $`e^{i\theta }`$ with a real $`\theta `$ variable, was used.
2. The normalization with respect to the behavior of $`\theta `$ was carried out explicitly inside the kernel, rather than post-simulation in the moments as in Eq. (52).

This type of representation is a norm-preserving gauge P representation, as discussed earlier. A parametrized family of gauges led to stable trajectories (as opposed to the large sampling error present with a positive-P simulation). However, some systematic errors were seen due to boundary terms. These boundary terms occurred because of the stochastic growth of the gauge term in $`\mathrm{\Omega }`$ space, when $`\theta `$ approached $`\pm \pi /2`$. With the gauge P representation introduced in this paper, a wide range of gauges do not lead to any systematic errors gauge\_paper , provided gauge growth is controlled.

We note here that the norm-preserving gauges have the property that, in the present notation, $`g_k=i[1i\mathrm{\Omega }^{\prime \prime }]f_k`$ . However, while the growth of $`\mathrm{\Omega }^{}`$ is stabilized, there is growth in the variance of $`\mathrm{\Omega }^{\prime \prime }`$. This means that the function $`f_k`$ must behave as a decreasing function of $`\mathrm{\Omega }^{\prime \prime }`$ in order to ensure that the distribution is bounded sufficiently in the weight-function space to avoid finite boundary terms. The detailed requirements and conditions for this type of gauge will be treated elsewhere.