text:
# Quantum-field-theoretical techniques for stochastic representation of quantum problems

## I Introduction

There is a well-known duality between Fokker-Planck equations (FPE) and stochastic differential equations (SDE), which goes back as far as Einstein’s and Langevin’s theories of Brownian motion (see Risken’s book  for a detailed discussion of FP equation and related issues). Langevin equations or, more generally, stochastic differential equations, have long been a successful computational tool in quantum stochastics , allowing for the numerical stochastic integration of systems for which analytical solution would be, at best, extremely difficult. There are long-standing rules, especially in quantum optics, for beginning with a particular system Hamiltonian and mapping this onto stochastic equations of motion for the field variables. However, this process depends on there being an FP equation equivalent to the master equation for the Hamiltonian in question; that is the partial differential equation for the appropriate probability or pseudoprobability distribution must contain derivatives of no higher than second order. This is the content of Pawula’s theorem . It restricts the method to a certain class of problems which, although containing many interesting cases, is not exhaustive. A generalised FP equation with higher order derivatives has no mapping onto a stochastic differential equation. As descriptions in terms of FP equations proper are the exception rather than the rule, it is of clear interest to extend present methods to allow for at least the numerical modelling of processes involving higher order noises.

For purposes of numerical simulations, time may always be regarded as discretised, so the development of a stochastic difference equation approximating a generalised FP will generally be sufficient. Since Pawula’s theorem only applies in the continuous time limit, this creates the loophole we need. What we will show here is that discretised stochastic equations may be devised for generalised FPEs, which, in a manner reminiscent of the positive-P representation of quantum optics, require a doubled phase-space. In agreement with Pawula’s theorem, the method we use to develop these noises has no natural extension to the continuous time limit. This notwithstanding, any average corresponding to a physical quantity is expected to be correct in this limit (whereas the rest of the averages diverge). The question whether a stochastic process in a certain generalised mathematical sense can be defined corresponding to our methods is a subject for futher investigation. This is not, however, a problem for practical computer simulations, as these are necessarily performed on a discrete time grid.

As in the well-known positive-P representation , nonphysical averages in our method determine the sampling noise. The fact that they are expected to diverge in the continuous time limit means that sampling noise grows as the time-grid spacing descreases. This is in contrast to the properties of the Wiener process, where sampling noise is independent of the grid spacing for a given sample size. In practice, however, this distinction between the Gaussian and higher-order noises is quantitative rather than qualitative, because for a given computing time the sample size necessarily decreases with the grid spacing. In either case the growth of sampling noise prevents the grid being made too fine.

As a specific example, in this paper we consider a quantum oscillator with Kerr nonlinearity. It is well known that the Wigner function of this system obeys a generalised FP equation with 3rd order derivatives. Our goal is to show that the latter is (approximately) equivalent to a system of stochastic difference equations. This equivalence is not straighforward but rather bears a lot of similarity to the positive-P representation known in quantum optics . We shall therefore call these equations a positive-W representation of the nonlinear quantum oscillator.

The positive-P representation originates with the single-time P-distribution, yet it allows one to calculate a much wider class of quantum averages , namely, averages of multi-time, time-normally ordered  operator products. The positive-P representation may thus be regarded as a constructive mapping of a quantum problem to a classical stochastic problem. Its characterictic property is that time-normal averages are mapped directly on classical averages. In a much similar manner, the positive-W equations, which originate with the single-time W-distribution, allow one to calculate multi-time, time-Wigner ordered operator averages. This new type of operator ordering is introduced in this paper; to the best of our knowledge it has never been discussed before. It generalises the single-time symmetric operator ordering, in the same way as the time-normal ordering generalises the single-time normal ordering. The positive-W representation is hence a constructive mapping of a quantum problem onto a classical stochastic problem, with the characterictic property of relating time-Wigner averages directly to classical averages.

The question then is if, rather than taking the usual route from q-number to c-number equations based on phase-space techniques , the positive-P and the positive-W representations may be derived by directly linking q-number Heisenberg equations to c-number Langevin equations. In this paper, we show that such a derivation is indeed possible, by employing methods  based on the techniques of quantum field theory. Our considerations closely follow the way in which Feynman diagram techniques were derived in textbooks dating back to the fifties and sixties . We derive Keldysh diagram series  for the Kerr oscillator, then recast them as a Wyld-type series , otherwise termed causal series  (see also Ref. ). (More precisely speaking, we derive generating expressions  for these series, avoiding diagram notation as such.) Causal series emerge as solutions to c-number stochastic problems . On the other hand, given a causal series, it is easy to write a stochastic differential (or difference) equation for which this series is a solution . This yields strikingly simple and powerful techniques for obtaining stochastic representations of quantum problems. (Cf., e.g., section V B 1 where positive-P equations are derived for an optical parametric oscillator, the essence of the derivation being a two-line calculation.) Not the least important property of these techniques is that they can be formulated using simple recipes and then used without any reference to the advanced methods employed in their derivation.

The paper is structured as follows. In section II, we reiterate the standard techniques of quantum field theory as applied to a nonlinear oscillator. We discuss in detail causal regularisation  of the propagator, which is necessary to order to make our relations unambiguous. A functional (Hori’s) form of Wick’s theorem  is introduced for the normal ordering and then generalised to the case of symmetric ordering. The result of the section is a pair of closed perturbative relations. One of them applies when the inital state of the field is characterised via a P-distribution; the other applies if this state is characterised via a W-distribution . In section III, we investigate the causal structure of the perturbative relations derived in section II. The concept of causality is introduced via the retarded Green’s function of the free Schrödinger equation. Following it, we are able to define the input and output of a quantum system. We then show that, physically, the input and output thus introduced correspond to a generalisation of Kubo’s linear reaction approach  to a full nonlinear quantum-stochastic response problem. In section IV, we develop techniques (based on the Hubbard-Stratonovich transformation ) which allow for a constructive mapping of the quantum nonlinear response problem onto a classical nonlinear stochastic response problem. For the oscillator, this results in the positive-P and positive-W representations which we are seeking. We complete the section with results of computer simulations of the Kerr oscillator using these representations. Finally, in section V we reformulate our results recipe-style, the way they should be applied in calculations. Despite all the tediousness of their derivation, we end up with two simple relations, which allow one to derive positive-P and positive-W representations for an arbitrary nonlinear quantum system. We illustrate their use by applying them to a number of systems common in quantum optics.

## II Quantum field theory of the Kerr oscillator

### A The model

The techniques we introduce in this paper are applicable to any Hamiltonians which have a polynomial form in the field operators. We also assume that the Hamiltonian can be divided into a quadratic part, called the free Hamiltonian, and the remainder, termed the interaction Hamiltonian. This allows us to introduce, in the usual manner, Schrödinger, Heisenberg, and interaction-picture field operators. However, this requirement may always be satisfied by subtracting a suitable quadratic term from the Hamiltonian, and declaring the remainder as being the “interaction” part (cf. the way a Higgs-type phase transition in an anharmonic oscillator was treated in Ref. ).

This generality notwithstanding, our techniques may be effectively demonstrated for a 1D oscillator (as is usual in quantum field theory). We therefore consider a nonlinear quantum oscillator with the Hamiltonian,

$``$ $`=`$ $`_0+_{\text{int}}=\omega \widehat{a}_\text{S}^{}\widehat{a}_\text{S}+{\displaystyle \frac{\kappa }{4}}\widehat{a}_\text{S}^2\widehat{a}_\text{S}^2,`$ (1)
using units such that $`\mathrm{}=1`$. In (1), $`\widehat{a}_\text{S}^{}`$ and $`\widehat{a}_\text{S}`$ are the usual pair of creation and annihilation operators with commutator $`[\widehat{a}_\text{S},\widehat{a}_\text{S}^{}]=1`$. They play the role of Schrödinger-picture field operators for this illustrative system. The field operators in the interaction picture are simply

$`\widehat{a}(t)=\text{e}^{i\omega t}\widehat{a}_\text{S},\widehat{a}^{}(t)=\text{e}^{i\omega t}\widehat{a}_\text{S}^{},`$ (2)
while Heisenberg picture operators will be denoted in a slanted font as $`\widehat{a}^{}(t)`$ and $`\widehat{a}(t)`$.

### B Time orderings of operators

Time ordering of operators puts operators from right to left in the order of increasing time arguments, e.g.,

$`T_+\widehat{a}(t^{})\widehat{a}^{}(t)=\{\begin{array}{cc}\widehat{a}^{}(t)\widehat{a}(t^{}),& tt^{},\\ \widehat{a}(t^{})\widehat{a}^{}(t),& t<t^{}.\end{array}`$ (5)
For equal times, the time ordering is specified as normal ordering (which places all creation operators on the left of annihilation operators). Further, the reverse time ordering, $`T_{}`$, places operators in the order of decreasing time arguments. Formally, it may be defined as a conjugate of $`T_+`$:

$`T_{}\widehat{P}=\left[T_+\left(\widehat{P}^{}\right)\right]^{},`$ (6)
where $`\widehat{P}`$ is a product of field operators. Then, e.g.,

$`T_{}\widehat{a}(t^{})\widehat{a}^{}(t)=\{\begin{array}{cc}\widehat{a}^{}(t)\widehat{a}(t^{}),& tt^{},\\ \widehat{a}(t^{})\widehat{a}^{}(t),& t>t^{}.\end{array}`$ (9)
For equal times, $`T_{}`$ also becomes normal ordering. Double time ordering is the combination of the $`T_+`$ and $`T_{}`$-orderings,

$`T_{}\widehat{P}_{}T_+\widehat{P}_+,`$ (10)
where $`\widehat{P}_{}`$ and $`\widehat{P}_+`$ are operator products (by definition, $`T_{}`$ acts on $`\widehat{P}_{}`$ and $`T_+`$ acts on $`\widehat{P}_+`$).

### C Non-stationary perturbation approach

#### 1 Perturbation expressions for time-ordered operator products

Heisenberg operators are related to those in the interaction picture via the evolution operator,

$`\widehat{a}(t)`$ $`=`$ $`𝒰^{}(t,\mathrm{})\widehat{a}(t)𝒰(t,\mathrm{}),`$ (11)
$`\widehat{a}^{}(t)`$ $`=`$ $`𝒰^{}(t,\mathrm{})\widehat{a}^{}(t)𝒰(t,\mathrm{}),`$ (12)
which is a solution to the Schrödinger equation,

$`i{\displaystyle \frac{𝒰(t,t_0)}{t}}=_{\text{int}}(t)𝒰(t,t_0),𝒰(t_0,t_0)=𝟙,`$ (13)
where

$`_{\text{int}}(t)={\displaystyle \frac{\kappa }{4}}\widehat{a}^2(t)\widehat{a}^2(t).`$ (14)
The evolution operator may be written as a time-ordered operator expression (T-exponent):

$`𝒰(t,t_0)=T_+\mathrm{exp}\left[i{\displaystyle _{t_0}^t}𝑑t^{}_{\text{int}}(t^{})\right].`$ (15)
It is unitary, and has the group property,

$`𝒰(t,t^{})𝒰(t^{},t^{\prime \prime })=𝒰(t,t^{\prime \prime }),tt^{}t^{\prime \prime }.`$ (16)
In particular it follows that

$`𝒰(t,t^{})=𝒰(t,\mathrm{})𝒰^{}(t^{},\mathrm{}).`$ (17)
Then, assuming that times are ordered in a certain way, $`t>t^{}`$, we have:

$`T_+\widehat{a}^{}(t)\widehat{a}(t^{})=𝒰^{}(\mathrm{},\mathrm{})𝒰(\mathrm{},t)\widehat{a}^{}(t)𝒰(t,t^{})\widehat{a}(t^{})𝒰(t^{},\mathrm{}).`$ (18)
It should be noted here that the combination

$`𝒰(\mathrm{},t)\widehat{a}^{}(t)𝒰(t,t^{})\widehat{a}(t^{})𝒰(t^{},\mathrm{})`$ (19)
is time-ordered. Finally, introducing the S-matrix,

$`𝒮=𝒰(\mathrm{},\mathrm{}),`$ (20)
we have,

$`T_+\widehat{a}^{}(t)\widehat{a}(t^{})`$ $`=`$ $`𝒮^{}T_+\widehat{a}(t^{})\widehat{a}^{}(t)\mathrm{exp}\left[{\displaystyle \frac{i\kappa }{4}}{\displaystyle _{\mathrm{}}^+\mathrm{}}𝑑\tau \widehat{a}^2(\tau )\widehat{a}^2(\tau )\right].`$ (21)
There is a certain subtlety in the way in which relations such as (15) and (21) should be understood. Quantities under ordering should be regarded as functions of $`\widehat{a}(t)`$ and $`\widehat{a}^{}(t)`$, whereas (2) may be used only after the time ordering has been completed. So calculating the time integral in (21) should follow the ordering rather than precede it. Were one to take the integral before the time ordering, it would prevent $`T_+`$ from acting on the operators comprising the interaction Hamiltonian.

#### 2 Functional perturbation techniques

For our purposes we need relations applying to the whole assemblage of time-ordered operator products. An arbitrary time-ordered operator product may be found by differentiating the following time-ordered operator exponent,

$`T_+\mathrm{exp}\left\{{\displaystyle _{\mathrm{}}^+\mathrm{}}𝑑t\left[\zeta (t)\widehat{a}^{}(t)+\zeta ^{}(t)\widehat{a}(t)\right]\right\}T_+\mathrm{exp}\left(\zeta \widehat{a}^{}+\zeta ^{}\widehat{a}\right),`$ (22)
where $`\zeta (t),\zeta ^{}(t)`$ is a pair of c-number functions. Then, e.g.,

$`T_+\widehat{a}(t^{})\widehat{a}^{}(t)={\displaystyle \frac{\delta ^2}{\delta \zeta (t)\delta \zeta ^{}(t^{})}}T_+\mathrm{exp}\left(\zeta \widehat{a}^{}+\zeta ^{}\widehat{a}\right)|_{\zeta =\zeta ^{}=0},`$ (23)
and so on. Equation (21) is readily generalised, resulting in

$`T_+\mathrm{exp}\left(\zeta \widehat{a}^{}+\zeta ^{}\widehat{a}\right)`$ $`=`$ $`𝒮^{}T_+\mathrm{exp}\left\{{\displaystyle _{\mathrm{}}^+\mathrm{}}𝑑t\left[\zeta (t)\widehat{a}^{}(t)+\zeta ^{}(t)\widehat{a}(t){\displaystyle \frac{i\kappa }{4}}\widehat{a}^2(t)\widehat{a}^2(t)\right]\right\}`$ (24)
$``$ $`𝒮^{}T_+\mathrm{exp}\left(\zeta \widehat{a}^{}+\zeta ^{}\widehat{a}{\displaystyle \frac{i\kappa }{4}}\widehat{a}^2\widehat{a}^2\right).`$ (25)
Relations for the inverse and double time orderings then read

$`T_{}\mathrm{exp}\left(\zeta \widehat{a}^{}+\zeta ^{}\widehat{a}\right)=T_{}\mathrm{exp}\left(\zeta \widehat{a}^{}+\zeta ^{}\widehat{a}+{\displaystyle \frac{i\kappa }{4}}\widehat{a}^2\widehat{a}^2\right)𝒮,`$ (26)
$`T_{}\mathrm{exp}\left(\zeta _{}\widehat{a}^{}+\zeta _{}^{}\widehat{a}\right)T_+\mathrm{exp}\left(\zeta _+\widehat{a}^{}+\zeta _+^{}\widehat{a}\right)`$ (27)
$`=T_{}\mathrm{exp}\left(\zeta _{}\widehat{a}^{}+\zeta _{}^{}\widehat{a}+{\displaystyle \frac{i\kappa }{4}}\widehat{a}^2\widehat{a}^2\right)T_+\mathrm{exp}\left(\zeta _+\widehat{a}^{}+\zeta _+^{}\widehat{a}{\displaystyle \frac{i\kappa }{4}}\widehat{a}^2\widehat{a}^2\right).`$ (28)
The relation for the double time ordering (naturally) involves four arbitrary c-number functions, $`\zeta _\pm (t),\zeta _\pm ^{}(t)`$. Note that, in this relation, the factors $`𝒮`$ and $`𝒮^{}`$ outside the orderings have cancelled each other. This results in a genuine double-time-ordered structure on the RHS of (28).

In what follows, we will make wide use of the brief notation as in Eqs. (22), (25) and (28), which will always be signalled by the absence of time arguments of the fields. Apart from these cases, we will not omit time arguments where they apply, as in Eq. (14), so as to remove any ambuguity. (Note that this is also the reason why the Scrödinger-picture field operators were denoted $`\widehat{a}_\text{S},\widehat{a}_\text{S}^{}`$, not $`\widehat{a},\widehat{a}^{}`$.)

### D Wick’s theorems

#### 1 Hori’s form of Wick’s theorem proper

A common way of dealing with time-ordered expressions is by converting them to normally ordered form. This is done following Wick’s theorem, which states that “a time ordered product of interaction-picture field operators equals the sum of all possible normally ordered operator products, obtained by replacing pairs of operators in the initial product by corresponding contractions (including the term without contractions)”. For the oscillator, we have only one nonzero contraction, namely,

$`T_+\widehat{a}(t^{})\widehat{a}^{}(t):\widehat{a}^{}(t)\widehat{a}(t^{}):`$ $`=`$ $`0\left|T_+\widehat{a}^{}(t)\widehat{a}(t^{})\right|0`$ (29)
$`=`$ $`\theta (t^{}t)\text{e}^{i\omega (t^{}t)}iG(t^{}t).`$ (30)
In the above, $`G(t)`$ is a retarded Green’s function of the free Schrödinger equation:

$`\left(i{\displaystyle \frac{}{t}}\omega \right)G(t)=\delta (t).`$ (31)
It then follows that

$`T_+\widehat{a}^{}(t)\widehat{a}(t^{})`$ $`=`$ $`\widehat{a}^{}(t)\widehat{a}(t^{})+iG(t^{}t),`$ (32)
$`T_+\widehat{a}^{}(t)\widehat{a}(t^{})\widehat{a}(t^{\prime \prime })`$ $`=`$ $`\widehat{a}^{}(t)\widehat{a}(t^{})\widehat{a}(t^{\prime \prime })+iG(t^{}t)\widehat{a}(t^{\prime \prime })+iG(t^{\prime \prime }t)\widehat{a}(t^{}),`$ (33)
$`T_+\widehat{a}^{}(t)\widehat{a}^{}(t^{})\widehat{a}(t^{\prime \prime })\widehat{a}(t^{\prime \prime \prime })`$ $`=`$ $`\widehat{a}^{}(t)\widehat{a}^{}(t^{})\widehat{a}(t^{\prime \prime })\widehat{a}(t^{\prime \prime \prime })`$ (37)
$`+iG(t^{\prime \prime }t)\widehat{a}^{}(t^{})\widehat{a}(t^{\prime \prime \prime })+iG(t^{\prime \prime \prime }t)\widehat{a}^{}(t^{})\widehat{a}(t^{\prime \prime })+`$
$`+iG(t^{\prime \prime }t^{})\widehat{a}^{}(t)\widehat{a}(t^{\prime \prime \prime })+iG(t^{\prime \prime \prime }t^{})\widehat{a}^{}(t)\widehat{a}(t^{\prime \prime })`$
$`G(t^{\prime \prime }t)G(t^{\prime \prime \prime }t^{})G(t^{\prime \prime }t^{})G(t^{\prime \prime \prime }t),`$
and so on.

As was noticed by Hori , the pattern of products with contractions is exactly that produced by a quadratic form of functional derivatives,

$`\mathrm{\Delta }={\displaystyle 𝑑t𝑑t^{}iG(t^{}t)\frac{\delta ^2}{\delta a(t^{})\delta a^{}(t)}},`$ (38)
acting repeatedly on products of c-number functions $`a(t),a^{}(t)`$. Applying $`\mathrm{\Delta }`$ once produces terms with a single contraction:

$`\mathrm{\Delta }a^{}(t)a(t^{})`$ $`=`$ $`iG(t^{}t),`$ (39)
$`\mathrm{\Delta }a^{}(t)a(t^{})a(t^{\prime \prime })`$ $`=`$ $`iG(t^{}t)a(t^{\prime \prime })+iG(t^{\prime \prime }t)a(t^{}),`$ (40)
$`\mathrm{\Delta }a^{}(t)a^{}(t^{})a(t^{\prime \prime })a(t^{\prime \prime \prime })`$ $`=`$ $`iG(t^{\prime \prime }t)a^{}(t^{})a(t^{\prime \prime \prime })+iG(t^{\prime \prime \prime }t)a^{}(t^{})a(t^{\prime \prime })+`$ (42)
$`iG(t^{\prime \prime }t^{})a^{}(t)a(t^{\prime \prime \prime })+iG(t^{\prime \prime \prime }t^{})a^{}(t)a(t^{\prime \prime }),`$
etc. Terms with $`n`$ contractions are produced by $`\mathrm{\Delta }^n/n!`$. As an example, consider the case where $`n=2`$:

$`{\displaystyle \frac{\mathrm{\Delta }^2}{2}}a^{}(t)a^{}(t^{})a(t^{\prime \prime })a(t^{\prime \prime \prime })`$ $`=`$ $`G(t^{\prime \prime }t)G(t^{\prime \prime \prime }t^{})G(t^{\prime \prime }t^{})G(t^{\prime \prime \prime }t).`$ (43)
The whole assemblage of terms required by Wick’s theorem (including those without contractions) is thus found by applying a differential operator,

$`1+\mathrm{\Delta }+{\displaystyle \frac{\mathrm{\Delta }^2}{2}}+\mathrm{}=\text{e}^\mathrm{\Delta }.`$ (44)
This allows one to write Wick’s theorem in a compact form as (with $`\widehat{P}`$ being an arbitrary operator product)

$`T_+\widehat{P}=:\left[\text{e}^\mathrm{\Delta }\left(\widehat{P}|_{\widehat{a}(t)a(t),\widehat{a}^{}(t)a^{}(t)}\right)\right]|_{a(t)\widehat{a}(t),a^{}(t)\widehat{a}^{}(t)}:.`$ (45)
#### 2 Causal regularisation

Wick’s theorem requires that no contractions should occur between operators with equal time arguments. (This is simply because the time ordering was specified for equal times as normal ordering.) However, applying (45) results in ambiguities, e.g.,

$`\mathrm{\Delta }a^{}(t)a(t)=iG(0)=i\theta (0)\text{ (undefined)}.`$ (46)
A convenient way around this problem is a causal regularisation of $`G(t)`$. To this end, we shall always assume that $`G(t)`$ is somehow smoothed while still preserving its causal nature: $`G(t)=0`$ for $`t<0`$. For a continuous function this also means that $`G(0)=0`$. With this specification (and implying a final limit of unsmoothed $`G`$), Eq. (45) becomes fully equivalent to Wick’s theorem.

#### 3 Wick’s theorem for double time ordering

It may be checked that proof of Wick’s theorem  is based only on the linear ordering of the time axis; consequently Wick’s theorem may be generalised to operators defined formally on any linearly ordered set. This clearly applies to the double time ordering, Eq. (10), which may alternatively be introduced as an ordering on the so-called C-contour . The C-contour (see Fig. 1) first travels from $`t=\mathrm{}`$ to $`t=+\mathrm{}`$ (direct branch) and then back to $`t=\mathrm{}`$ (reverse branch). Then,

$`T_{}\widehat{P}_{}T_+\widehat{P}_+=T_C\widehat{P}_{}|_{\widehat{a}\widehat{a}_{},\widehat{a}^{}\widehat{a}_{}^{}}\widehat{P}_+|_{\widehat{a}\widehat{a}_+,\widehat{a}^{}\widehat{a}_+^{}},`$ (47)
with the subscripts ‘$`+`$’ and ‘$``$’ specifying operators as being defined on the direct and reverse branches of the C-contour. (In other words, operators acquire an additional argument; a C-contour index. This argument is useful purely for labelling purposes and should be disregarded after the ordering has been performed.) For the $`T_C`$ ordering, operator contraction becomes a matrix in respect of the C-countour indices, ($`\alpha ,\beta =+,`$)

$`iG_{\alpha \beta }(t^{}t)=0\left|T_C\widehat{a}_\alpha (t^{})\widehat{a}_\beta ^{}(t)\right|0,`$ (48)
and the functional (Hori’s) form of Wick’s theorem generalised to the double time ordering is found to be:

$`T_{}\widehat{P}_{}T_+\widehat{P}_+=:[\text{e}^{\mathrm{\Delta }_C}(\widehat{P}_{}|_{\widehat{a}(t)a_{}(t),\widehat{a}^{}(t)a_{}^{}(t)}`$ (49)
$`\times \widehat{P}_+|_{\widehat{a}(t)a_+(t),\widehat{a}^{}(t)a_+^{}(t)})]|_{a_{}(t),a_+(t)\widehat{a}(t),a_{}^{}(t),a_+^{}(t)\widehat{a}^{}(t)}:.`$ (50)
Here, $`\widehat{P}_{}`$ and $`\widehat{P}_+`$ are arbitrary operator products, $`a_\pm (t),a_\pm ^{}(t)`$ are four independent c-number functions, and $`\mathrm{\Delta }_C`$ is a quadratic form of functional derivatives,

$`\mathrm{\Delta }_C={\displaystyle }dtdt^{}[iG_{++}(t^{}t){\displaystyle \frac{\delta ^2}{\delta a_+(t^{})\delta a_+^{}(t)}}`$ (51)
$`+iG_{}(t^{}t){\displaystyle \frac{\delta ^2}{\delta a_{}(t^{})\delta a_{}^{}(t)}}+iG_+(t^{}t){\displaystyle \frac{\delta ^2}{\delta a_{}(t^{})\delta a_+^{}(t)}}],`$ (52)
where the kernels are the three nonzero components of $`G_{\alpha \beta }`$ (shown schematically in Fig. 1). In turn, these are conveniently expressed in terms of $`G(t)`$:

$`G_{++}(t)=G(t),G_{}(t)=G^{}(t),G_+(t)=G(t)G^{}(t).`$ (53)
Note that Eqs. (50) and (53) are formulated so as to eliminate the concept of the C-contour from our considerations from here on.

Relations (53) require a word of caution. For equal times, both the $`T_+`$ and the $`T_{}`$ orderings are specified as normal ordering, hence the no-contractions-between-the-same-time-operators caveat of Wick’s theorem applies equally to operators under the $`T_{}`$ ordering. In (50), it is enforced through the equation relating $`G_{}(t)`$ to $`G(t)`$ and the (implied) causal regularisation of $`G(t)`$. At the same time, through the equation relating $`G_+(t)`$ to $`G(t)`$, the regularisation also modifies $`G_+(t)`$, “burning a hole” in it in the vicinity of $`t=0`$. Both the causal regularisation and relations (53) play central roles in our considerations, so we have to make sure that the stated modification of $`G_+(t)`$ does not lead to incorrect results. The fact that weird results may indeed follow demonstrates, e.g., a “proof” that $`\widehat{a}_\text{S}`$ and $`\widehat{a}_\text{S}^{}`$ commute. With regularisation, $`G_+(0)=0`$, and we “obtain”,

$`\widehat{a}_\text{S}\widehat{a}_\text{S}^{}=T_C\widehat{a}_{}(t)\widehat{a}_+^{}(t)=\widehat{a}^{}(t)\widehat{a}(t)+iG_+(0)=\widehat{a}_\text{S}^{}\widehat{a}_\text{S}.`$ (54)
The flaw in this “proof” is that all quantities we deal with should be regarded as generalised functions (distributions) and not pointwise functions. This is especially true under regularisation. “Holes” in continuous functions which emerge due to regularisation should be simply ignored (smeared out). We would only expect problems associated with the “holes” if they were overlapping with sufficiently strong singularities, $`\delta `$-functions or worse, whereas the worst type of singularity that we may expect to occur is a step-function. The “hole” in $`G_+(t)`$ should thus be of no consequence.

#### 4 Reordering time-ordered operators symmetrically

Wick’s theorem also holds if the desired result is symmetric, or Wigner, operator ordering. (For a definition of the symmetric ordering see, e.g. .) The easiest way to demonstrate this is to derive an analog of Eq. (50). A relation between a normal and a symmetric representation of a particular Scrödinger operator, most suitable for our purposes, is given in Ref. . It can be readily adapted to operators in the interaction picture resulting in

$`:\widehat{P}:=\text{W}\left\{\left[\mathrm{exp}\left(i{\displaystyle \frac{\delta }{\delta a}}G_W{\displaystyle \frac{\delta }{\delta a^{}}}\right)\left(\widehat{P}|_{\widehat{a}(t)a(t),\widehat{a}^{}(t)a^{}(t)}\right)\right]|_{a(t)\widehat{a}(t),a^{}(t)\widehat{a}^{}(t)}\right\}.`$ (55)
Here, $`\widehat{P}`$ is an arbitrary operator product, W denotes Wigner ordering, and

$`{\displaystyle \frac{\delta }{\delta a}}G_W{\displaystyle \frac{\delta }{\delta a^{}}}={\displaystyle 𝑑t𝑑t^{}G_W(t^{}t)\frac{\delta ^2}{\delta a(t^{})\delta a^{}(t)}},`$ (56)
where the kernel $`G_W`$ is defined as

$`iG_W(tt^{})=\widehat{a}^{}(t)\widehat{a}(t^{})\text{W}\widehat{a}(t^{})a^{}(t)=0\left|\text{W}\widehat{a}(t^{})a^{}(t)\right|0={\displaystyle \frac{\text{e}^{i\omega (t^{}t)}}{2}}.`$ (57)
Note that the type of brief notation introduced by (56) accounts for asymmetry of the kernel $`G_W`$. We now make use of the following fancy form of the rule of product differentiation,

$`\mathrm{\Phi }\left({\displaystyle \frac{\delta }{\delta \phi }}\right)\mathrm{\Phi }_1\left(\phi \right)\mathrm{\Phi }_2\left(\phi \right)=\mathrm{\Phi }\left({\displaystyle \frac{\delta }{\delta \phi _1}}+{\displaystyle \frac{\delta }{\delta \phi _2}}\right)\mathrm{\Phi }_1\left(\phi _1\right)\mathrm{\Phi }_2\left(\phi _2\right)|_{\phi _1,\phi _2=\phi }.`$ (58)
where $`\phi (t),\phi _1(t),\phi _2(t)`$ are c-number functions and $`\mathrm{\Phi }(\phi ),\mathrm{\Phi }_1(\phi ),\mathrm{\Phi }_2(\phi )`$ are functionals of such functions. Combining (50) and (55) and employing (58), we arrive at a relation,

$`T_{}\widehat{P}_{}T_+\widehat{P}_+=\text{W}\{[\text{e}^{\stackrel{~}{\mathrm{\Delta }}_C^W}(\widehat{P}_{}|_{\widehat{a}(t)a_{}(t),\widehat{a}^{}(t)a_{}^{}(t)}`$ (59)
$`\times \widehat{P}_+|_{\widehat{a}(t)a_+(t),\widehat{a}^{}(t)a_+^{}(t)})]|_{a_{}(t),a_+(t)\widehat{a}(t),a_{}^{}(t),a_+^{}(t)\widehat{a}^{}(t)}\},`$ (60)
where

$`\stackrel{~}{\mathrm{\Delta }}_C^W=i{\displaystyle \underset{\alpha ,\beta =+,}{}}{\displaystyle \frac{\delta }{\delta a_\alpha }}\left(G_{\alpha \beta }+G_W\right){\displaystyle \frac{\delta }{\delta a_\beta ^{}}}.`$ (61)
For nonzero time arguments, kernels comprising (61) may be expressed by the retarded Green’s function, $`G(t)`$. We note that

$`G_{\alpha \beta }(t)+G_W(t)=G_{\alpha \beta }(t){\displaystyle \frac{1}{2}}[G(t)G^{}(t)]G_{\alpha \beta }^W(t).`$ (62)
Then, employing (53), we find,

$`G_{++}^W(t)=G_{}^W(t)={\displaystyle \frac{1}{2}}\left[G(t)+G^{}(t)\right],G_+^W(t)=G_+^W(t)={\displaystyle \frac{1}{2}}\left[G(t)G^{}(t)\right].`$ (63)
Extending these relations to all times, we define,

$`\mathrm{\Delta }_C^W=i{\displaystyle \underset{\alpha ,\beta =+,}{}}{\displaystyle \frac{\delta }{\delta a_\alpha }}G_{\alpha \beta }^W{\displaystyle \frac{\delta }{\delta a_\beta ^{}}}.`$ (64)
The question then is if $`\stackrel{~}{\mathrm{\Delta }}_C^W`$ in (60) may be replaced by $`\mathrm{\Delta }_C^W`$. The same arguments as above show that the “holes” in $`G_+^W(t)`$ and $`G_+^W(t)`$ are of no concern. This is not the case for $`G_{++}^W(t)`$ and $`G_{}^W(t)`$. Indeed, recall that the $`T_+`$ and $`T_{}`$-ordering are both specified for equal times as normal ordering. $`\stackrel{~}{\mathrm{\Delta }}_C^W`$ takes care of such same-time operator groups reordering them symmetrically, whereas $`\mathrm{\Delta }_C^W`$ misses them. That is, $`\mathrm{\Delta }_C^W`$ enforces the no-contractions-between-the-same-time-operators caveat of Wick’s theorem, whereas this caveat does not apply if a double-time-ordered operator product is reordered symmetrically. We should thus either amend Wick’s theorem introducing same-time contractions, or, which is more convenient and suitable for our purposes, redefine the time orderings. We then obtain

$`T_{}^W\widehat{P}_{}T_+^W\widehat{P}_+=\text{W}\{[\text{e}^{\mathrm{\Delta }_C^W}(\widehat{P}_{}|_{\widehat{a}(t)a_{}(t),\widehat{a}^{}(t)a_{}^{}(t)}`$ (65)
$`\times \widehat{P}_+|_{\widehat{a}(t)a_+(t),\widehat{a}^{}(t)a_+^{}(t)})]|_{a_{}(t),a_+(t)\widehat{a}(t),a_{}^{}(t),a_+^{}(t)\widehat{a}^{}(t)}\},`$ (66)
were $`T_{}^W`$ and $`T_+^W`$ differ from $`T_{}`$ and $`T_+`$ in that same-time operators are ordered symmetrically rather than normally.

### E Closed perturbative relations for quantum field averages

#### 1 The method of normal ordering

For practical purposes, the quantities of interest are operator averages rather than the operators themselves. We therefore consider a characteristic functional of averages of double-time-ordered operator products:

$`\mathrm{\Xi }(\zeta _{},\zeta _+,\zeta _{}^{},\zeta _+^{})=T_{}\mathrm{exp}\left(\zeta _{}\widehat{a}^{}+\zeta _{}^{}\widehat{a}\right)T_+\mathrm{exp}\left(\zeta _+\widehat{a}^{}+\zeta _+^{}\widehat{a}\right).`$ (67)
Angle brackets here define an averaging over the Heisenberg $`\rho `$-matrix of the quantum field, (or over the field’s initial state, which is the same thing):

$`\mathrm{}=\text{Tr}\widehat{\rho }\left(\mathrm{}\right).`$ (68)
Although defined in q-number terms, in itself functional (67) is a c-number object. Applying the closed perturbative relation, Eq. (28), and then Wick’s theorem, Eq. (50), allows one to express it as an average of a normally ordered operator expression. To eliminate q-numbers completely, we shall characterise the initial state of the field by the corresponding P-distribution,

$`P(\alpha )`$ $`=`$ $`{\displaystyle \frac{1}{\pi ^2}}{\displaystyle d^2\eta \text{e}^{\eta \left(\widehat{a}^{}\alpha ^{}\right)\eta ^{}\left(\widehat{a}\alpha \right)|\eta |^2/2}},`$ (69)
$`\widehat{\rho }`$ $`=`$ $`{\displaystyle d^2\alpha P(\alpha )|\alpha \alpha |},`$ (70)
where $`|\alpha `$ is a coherent state with the amplitude $`\alpha `$:

$`\widehat{a}|\alpha =\alpha |\alpha .`$ (71)
For any normally ordered operator expression, $`:\widehat{X}:`$, we then have

$`:\widehat{X}:={\displaystyle }d^2\alpha P(\alpha )\alpha |:\widehat{X}:|\alpha ={\displaystyle }d^2\alpha P(\alpha )\widehat{X}|_{\widehat{a}(t)\alpha (t),\widehat{a}^{}(t)\alpha ^{}(t)},`$ (72)
where $`\alpha (t)=\alpha \text{e}^{i\omega t}`$ is the coherent amplitude of the interaction-picture field operator, $`\widehat{a}(t)|\alpha =\alpha (t)|\alpha `$. Combining Eqs. (28), (50) and (72) then yields:

$`\mathrm{\Xi }(\zeta _{},\zeta _+,\zeta _{}^{},\zeta _+^{})`$ $`=`$ $`{\displaystyle d^2\alpha P(\alpha )\text{e}^{\mathrm{\Delta }_C}\text{e}^{\zeta _{}a_{}^{}+\zeta _{}^{}a_{}+\zeta _+a_+^{}+\zeta _+^{}a_+}}`$ (74)
$`\times \text{e}^{\frac{i\kappa }{4}\left(a_{}^2a_{}^2a_+^2a_+^2\right)}|_{a_{}(t)=a_+(t)=\alpha (t),a_{}^{}(t)=a_+^{}(t)=\alpha ^{}(t)}.`$
#### 2 The method of Wigner ordering

In order to derive an analog of Eq. (74) for symmetric ordering, consider, to start with, relation (15) for the evolution operator. By derivation , it is a sum of $`T`$-ordered terms, found by expanding the exponent in a power series:

$`𝒰(t,t_0)`$ $`=`$ $`T\mathrm{exp}\left[i{\displaystyle _{t_0}^t}𝑑t^{}_{\text{int}}(t^{})\right]`$ (75)
$`=`$ $`𝟙𝕚{\displaystyle _{𝕥_\mathrm{𝟘}}^𝕥}𝕕𝕥^{}_{\text{int}}(𝕥^{})+{\displaystyle \frac{(𝕚)^\mathrm{𝟚}}{\mathrm{𝟚}}}{\displaystyle _{𝕥_\mathrm{𝟘}}^𝕥}𝕕𝕥^{}𝕕𝕥^{\prime \prime }𝕋_{\text{int}}(𝕥^{})_{\text{int}}(𝕥^{\prime \prime })+\mathrm{}.`$ (76)
The $`T`$-ordering here applies, strictly speaking, to the Hamiltonian operators regarded as entities rather than to the field operators. Redefining the time ordering for the field operators via Eq. (5) changes the integrands, because, e.g.,

$`T_{\text{int}}(t)_{\text{int}}(t)=\left[_{\text{int}}(t)\right]^2T_+_{\text{int}}(t)_{\text{int}}(t)=:\left[_{\text{int}}(t)\right]^2:.`$ (77)
This change, however, (i) is finite and (ii) affects only the submanifold of measure zero of the integration manifold. Hence it has no effect on the result of integration (this and further arguments below may become flawed for continuous-space problems where expressions like $`\left[_{\text{int}}(t)\right]^2:\left[_{\text{int}}(t)\right]^2:`$ often contain infinities).

The key property of the $`T_+`$-ordering is that it does not affect the Hamiltonian operators themselves,

$`T_+\left\{_{\text{int}}(t)\right\}=_{\text{int}}(t),`$ (78)
which in turn is due to the normal form of $`_{\text{int}}(t)`$. It is then immediately clear that (15) may be rewritten in terms of $`T_+^W`$ by using a symmetric form of the interaction:

$`{\displaystyle \frac{\kappa }{4}}\widehat{a}^2(t)\widehat{a}^2(t)={\displaystyle \frac{\kappa }{4}}\text{W}\left\{\widehat{a}^2(t)\widehat{a}^2(t)2\widehat{a}^{}(t)\widehat{a}(t)+{\displaystyle \frac{1}{2}}\right\}_{\text{int}}^W(t),`$ (79)
so that

$`T_+^W\left\{_{\text{int}}^W(t)\right\}=_{\text{int}}^W(t),`$ (80)
When deriving (79) we have used the property that

$`\text{W}\left\{\widehat{a}^2\widehat{a}^2\right\}`$ $`=`$ $`{\displaystyle \frac{1}{6}}\left(\widehat{a}^2\widehat{a}^2+\widehat{a}^{}\widehat{a}\widehat{a}^{}\widehat{a}+\widehat{a}^{}\widehat{a}^2\widehat{a}^{}+\widehat{a}\widehat{a}^2\widehat{a}+\widehat{a}\widehat{a}^{}\widehat{a}\widehat{a}^{}+\widehat{a}^2\widehat{a}^2\right),`$ (81)
$`\text{W}\left\{\widehat{a}^{}\widehat{a}\right\}`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left(\widehat{a}^{}\widehat{a}+\widehat{a}\widehat{a}^{}\right).`$ (82)
For the evolution operator we then find,

$`𝒰(t,t_0)`$ $`=`$ $`T_+^W\mathrm{exp}\left[i{\displaystyle _{t_0}^t}𝑑t^{}_{\text{int}}^W(t^{})\right]`$ (83)
$`=`$ $`T_+^W\mathrm{exp}\left\{{\displaystyle \frac{i\kappa }{4}}{\displaystyle _{t_0}^t}𝑑t^{}\left[\widehat{a}^2(t^{})\widehat{a}^2(t^{})2\widehat{a}^{}(t^{})\widehat{a}(t^{})+{\displaystyle \frac{1}{2}}\right]\right\}.`$ (84)
Relation (28) also remains valid after replacing $`T_+T_+^W`$ and $`_{\text{int}}(t)_{\text{int}}^W(t)`$. For any symmetrically ordered operator expression, $`\text{W}\left\{\widehat{X}\right\}`$, we have

$`\text{W}\left\{\widehat{X}\right\}={\displaystyle d^2\alpha W(\alpha )\widehat{X}}|_{\widehat{a}(t)\alpha (t),\widehat{a}^{}(t)\alpha ^{}(t)},`$ (85)
where $`W(\alpha )`$ is the Wigner distribution characterising the initial state,

$`W(\alpha )={\displaystyle \frac{1}{\pi ^2}}{\displaystyle d^2\eta \text{e}^{\eta \left(\widehat{a}^{}\alpha ^{}\right)\eta ^{}\left(\widehat{a}\alpha \right)}}.`$ (86)
Finally, we arrive at:

$`\mathrm{\Xi }^W(\zeta _{},\zeta _+,\zeta _{}^{},\zeta _+^{})`$ $`=`$ $`{\displaystyle d^2\alpha W(\alpha )\text{e}^{\mathrm{\Delta }_C^W}\text{e}^{\zeta _{}a_{}^{}+\zeta _{}^{}a_{}+\zeta _+a_+^{}+\zeta _+^{}a_+}}`$ (88)
$`\times \text{e}^{\frac{i\kappa }{4}\left(a_{}^2a_{}^2a_+^2a_+^2+2a_+^{}a_+2a_{}^{}a_{}\right)}|_{a_{}(t)=a_+(t)=\alpha (t),a_{}^{}(t)=a_+^{}(t)=\alpha ^{}(t)},`$
where

$`\mathrm{\Xi }^W(\zeta _{},\zeta _+,\zeta _{}^{},\zeta _+^{})=T_{}^W\mathrm{exp}\left(\zeta _{}\widehat{a}^{}+\zeta _{}^{}\widehat{a}\right)T_+^W\mathrm{exp}\left(\zeta _+\widehat{a}^{}+\zeta _+^{}\widehat{a}\right).`$ (89)
Note that averages generated by functionals (67) and (89) may indeed differ. For instance, with $`\widehat{n}(t)=\widehat{a}^{}(t)\widehat{a}(t)`$,

$`T_{}\widehat{n}(t)T_+\widehat{n}(t)`$ $`=`$ $`\left[\widehat{n}(t)\right]^2,`$ (90)
$`T_{}^W\widehat{n}(t)T_+^W\widehat{n}(t)`$ $`=`$ $`\left[\widehat{n}(t)+{\displaystyle \frac{1}{2}}\right]^2.`$ (91)
To obtain $`\left[\widehat{n}(t)\right]^2`$ directly from (89), we may use, e.g., the property

$`\left[\widehat{n}(t)\right]^2=\underset{\delta t0}{lim}T_{}^W\widehat{a}^{}(t)\widehat{a}(t+\delta t)T_+^W\widehat{a}^{}(t+\delta t)\widehat{a}(t).`$ (92)
This limit should be preceded by the one associated with causal regularisation.

## III Causal variables

Our next goal is to investigate the causal structure of Eqs. (74) and (88). The concept of causality is introduced via the retarded Green’s function, $`G(t)`$. Following this, we are able to define the input and output of a quantum system. We then show that, physically, the input and output thus introduced correspond to a generalisation of Kubo’s linear reaction approach  to a full nonlinear quantum-stochastic response problem.

### A The method of normal ordering

Consider in more detail the differential quadratic form in (74). Making use of relations (53), and utilising the notation introduced by Eq. (56), we find

$`\mathrm{\Delta }_C=i\left({\displaystyle \frac{\delta }{\delta a_+}}+{\displaystyle \frac{\delta }{\delta a_{}}}\right)G{\displaystyle \frac{\delta }{\delta a_+^{}}}i\left({\displaystyle \frac{\delta }{\delta a_+^{}}}+{\displaystyle \frac{\delta }{\delta a_{}^{}}}\right)G^{}{\displaystyle \frac{\delta }{\delta a_{}}}`$ (93)
We now change the functional variables, $`a_\pm (t),a_\pm ^{}(t)a(t),a^{}(t),\xi (t),\xi ^{}(t)`$, in order to obtain

$`\mathrm{\Delta }_C={\displaystyle \frac{\delta }{\delta a}}G{\displaystyle \frac{\delta }{\delta \xi ^{}}}+{\displaystyle \frac{\delta }{\delta a^{}}}G^{}{\displaystyle \frac{\delta }{\delta \xi }}.`$ (94)
That is,

$`{\displaystyle \frac{\delta }{\delta \xi ^{}(t)}}`$ $`=`$ $`i{\displaystyle \frac{\delta }{\delta a_+^{}(t)}},`$ (95)
$`{\displaystyle \frac{\delta }{\delta \xi (t)}}`$ $`=`$ $`i{\displaystyle \frac{\delta }{\delta a_{}(t)}},`$ (96)
$`{\displaystyle \frac{\delta }{\delta a(t)}}`$ $`=`$ $`{\displaystyle \frac{\delta }{\delta a_+(t)}}+{\displaystyle \frac{\delta }{\delta a_{}(t)}},`$ (97)
$`{\displaystyle \frac{\delta }{\delta a^{}(t)}}`$ $`=`$ $`{\displaystyle \frac{\delta }{\delta a_+^{}(t)}}+{\displaystyle \frac{\delta }{\delta a_{}^{}(t)}}.`$ (98)
These relations determine the new variables up to given functions which we chose to be zero:

$`a_+(t)`$ $`=`$ $`a(t),`$ (99)
$`a_{}(t)`$ $`=`$ $`a(t)i\xi (t),`$ (100)
$`a_+^{}(t)`$ $`=`$ $`a^{}(t)+i\xi ^{}(t),`$ (101)
$`a_{}^{}(t)`$ $`=`$ $`a^{}(t).`$ (102)
Consider now the second exponent in (74):

$`\zeta _{}a_{}^{}+\zeta _{}^{}a_{}+\zeta _+a_+^{}+\zeta _+^{}a_+=\left(\zeta _{}+\zeta _+\right)a^{}+\left(\zeta _{}^{}+\zeta _+^{}\right)ai\zeta _{}^{}\xi +i\zeta _+\xi ^{}.`$ (103)
This clearly suggests another substitution, this time in the functional $`\mathrm{\Xi }`$ itself:

$`\mathrm{\Xi }(\zeta _{},\zeta _+,\zeta _{}^{},\zeta _+^{})\mathrm{\Phi }(\zeta ,\zeta ^{},\sigma ,\sigma ^{}),`$ (104)
where

$`\zeta _{}(t)+\zeta _+(t)`$ $`=`$ $`\zeta (t),`$ (105)
$`\zeta _{}^{}(t)+\zeta _+^{}(t)`$ $`=`$ $`\zeta ^{}(t),`$ (106)
$`i\zeta _{}^{}(t)`$ $`=`$ $`\sigma ^{}(t),`$ (107)
$`+i\zeta _+(t)`$ $`=`$ $`\sigma (t),`$ (108)
and

$`\zeta _+(t)`$ $`=`$ $`i\sigma (t),`$ (109)
$`\zeta _{}(t)`$ $`=`$ $`\zeta (t)+i\sigma (t),`$ (110)
$`\zeta _+^{}(t)`$ $`=`$ $`\zeta ^{}(t)i\sigma ^{}(t),`$ (111)
$`\zeta _{}^{}(t)`$ $`=`$ $`i\sigma ^{}(t).`$ (112)
In causal variables thus introduced, Eq. (74) becomes:

$`\mathrm{\Phi }(\zeta ,\zeta ^{},\sigma ,\sigma ^{})`$ $`=`$ $`{\displaystyle d^2\alpha P(\alpha )\mathrm{exp}\left(\frac{\delta }{\delta a}G\frac{\delta }{\delta \xi ^{}}+\frac{\delta }{\delta a^{}}G^{}\frac{\delta }{\delta \xi }\right)\mathrm{exp}\left(\zeta a^{}+\zeta ^{}a+\xi \sigma ^{}+\xi ^{}\sigma \right)}`$ (114)
$`\times \mathrm{exp}S_{\text{int}}(\xi ,\xi ^{},a,a^{})|_{a(t)=\alpha (t),a^{}(t)=\alpha ^{}(t),\xi (t)=\xi ^{}(t)=0},`$
where

$`S_{\text{int}}(\xi ,\xi ^{},a,a^{})={\displaystyle \frac{\kappa }{2}}\left(\xi a^2a+\xi ^{}a^2a^{}\right)+{\displaystyle \frac{i\kappa }{4}}\left(\xi ^2a^2\xi ^2a^2\right).`$ (115)
It should be stressed that equation (114) is universal, while all details of the problem enter through $`S_{\text{int}}`$. In general $`S_{\text{int}}`$ is found as

$`S_{\text{int}}(\xi ,\xi ^{},a,a^{})=i{\displaystyle 𝑑t\left[h(a^{}(t),a(t)i\xi (t))h(a^{}(t)+i\xi ^{}(t),a(t))\right]},`$ (116)
where $`h`$ is the normal representation of the interaction Hamiltonian,

$`_{\text{int}}=:h(\widehat{a}^{},\widehat{a}):.`$ (117)
Generalisation of (114) to multimode problems (see Sec. V) is also straightforward.

### B The method of Wigner ordering

Causal variables for the case of symmetric ordering are defined in a similar way. We require that

$`\mathrm{\Delta }_C^W={\displaystyle \frac{\delta }{\delta a}}G{\displaystyle \frac{\delta }{\delta \xi ^{}}}+{\displaystyle \frac{\delta }{\delta a^{}}}G^{}{\displaystyle \frac{\delta }{\delta \xi }},`$ (118)
and then proceed as above to yield two sets of substitutions:

$`a_+(t)`$ $`=`$ $`a(t)+{\displaystyle \frac{i}{2}}\xi (t),`$ (119)
$`a_{}(t)`$ $`=`$ $`a(t){\displaystyle \frac{i}{2}}\xi (t),`$ (120)
$`a_+^{}(t)`$ $`=`$ $`a^{}(t)+{\displaystyle \frac{i}{2}}\xi ^{}(t),`$ (121)
$`a_{}^{}(t)`$ $`=`$ $`a^{}(t){\displaystyle \frac{i}{2}}\xi ^{}(t),`$ (122)
and

$`\zeta _+(t)`$ $`=`$ $`{\displaystyle \frac{1}{2}}\zeta (t)i\sigma (t),`$ (123)
$`\zeta _{}(t)`$ $`=`$ $`{\displaystyle \frac{1}{2}}\zeta (t)+i\sigma (t),`$ (124)
$`\zeta _+^{}(t)`$ $`=`$ $`{\displaystyle \frac{1}{2}}\zeta ^{}(t)i\sigma ^{}(t),`$ (125)
$`\zeta _{}^{}(t)`$ $`=`$ $`{\displaystyle \frac{1}{2}}\zeta ^{}(t)+i\sigma ^{}(t).`$ (126)
Further, on rewriting functional (89) in the variables (123)–(126),

$`\mathrm{\Xi }^W(\zeta _{},\zeta _+,\zeta _{}^{},\zeta _+^{})\mathrm{\Phi }^W(\zeta ,\zeta ^{},\sigma ,\sigma ^{}),`$ (127)
we find that

$`\mathrm{\Phi }^W(\zeta ,\zeta ^{},\sigma ,\sigma ^{})`$ $`=`$ $`{\displaystyle d^2\alpha W(\alpha )\mathrm{exp}\left(\frac{\delta }{\delta a}G\frac{\delta }{\delta \xi ^{}}+\frac{\delta }{\delta a^{}}G^{}\frac{\delta }{\delta \xi }\right)\mathrm{exp}\left(\zeta a^{}+\zeta ^{}a+\xi \sigma ^{}+\xi ^{}\sigma \right)}`$ (129)
$`\times \mathrm{exp}S_{\text{int}}^W(\xi ,\xi ^{},a,a^{})|_{a(t)=\alpha (t),a^{}(t)=\alpha ^{}(t),\xi (t)=\xi ^{}(t)=0}.`$
Here,

$`S_{\text{int}}^W(\xi ,\xi ^{},a,a^{})`$ $`=`$ $`i{\displaystyle }dt[h^W(a^{}(t){\displaystyle \frac{i}{2}}\xi ^{}(t),a(t){\displaystyle \frac{i}{2}}\xi (t))`$ (131)
$`h^W(a^{}(t)+{\displaystyle \frac{i}{2}}\xi ^{}(t),a(t)+{\displaystyle \frac{i}{2}}\xi (t))],`$
where $`h^W`$ is the symmetric representation of the interaction Hamiltonian,

$`_{\text{int}}=\text{W}\left\{h^W(\widehat{a}^{},\widehat{a})\right\}.`$ (132)
For the Kerr nonlinearity,

$`S_{\text{int}}^W(\xi ,\xi ^{},a,a^{})={\displaystyle \frac{\kappa }{2}}\left[\left(\xi a^{}+\xi ^{}a\right)\left(a^{}a2\right)\right]{\displaystyle \frac{\kappa }{8}}\left(\xi ^2\xi ^{}a^{}+\xi ^2\xi a\right).`$ (133)
It should not be overlooked that the definition of causal variables is ordering-specific. In fact, the notation $`\zeta ,\zeta ^{},\sigma ,\sigma ^{}`$ is used for two different sets of causal variables: (109)–(112) in the case of normal ordering, and (123)–(126) in the case of symmetric ordering. Since it is always clear by the context which case is being discussed, no confusion should occur.

### C Quantum nonlinear-reaction problem

To gain more insight into the causal variables, we now introduce a variable pump term into the interaction Hamiltonian, which then reads:

$`\stackrel{~}{H}_{\text{int}}(t)={\displaystyle \frac{\kappa }{4}}\widehat{a}^2(t)\widehat{a}^2(t)+s(t)\widehat{a}^{}(t)+s^{}(t)\widehat{a}(t).`$ (134)
The external source, $`s(t)`$, is a given c-number function. We mark by tilde all quantitities defined in the presence of the source. With the source,

$`\stackrel{~}{S}_{\text{int}}=S_{\text{int}}+\xi s^{}+\xi ^{}s.`$ (135)
Moving $`\xi s^{}+\xi ^{}s`$ to the second exponent in (114) results in an identity,

$`\stackrel{~}{\mathrm{\Phi }}(\zeta ,\zeta ^{},\sigma ,\sigma ^{})=\mathrm{\Phi }(\zeta ,\zeta ^{},\sigma +s,\sigma ^{}+s^{}),`$ (136)
which also holds for the $`\stackrel{~}{\mathrm{\Phi }}^W/\mathrm{\Phi }^W`$ pair.

This way, the physical information contained in the dependence of both $`\mathrm{\Phi }(\zeta ,\zeta ^{},\sigma ,\sigma ^{})`$ and $`\mathrm{\Phi }^W(\zeta ,\zeta ^{},\sigma ,\sigma ^{})`$ on $`\sigma `$ and $`\sigma ^{}`$ is the system’s reaction to an external perturbations, while the variables $`\sigma ,\sigma ^{}`$ define an input of the system. The fact that this input is also determined dynamically makes it independent of ordering. The variables $`\zeta ,\zeta ^{}`$ define an output of the system. Unlike the input, the output turns out to be ordering-specific. Consider firstly the case of normal ordering; assuming the source to be arbitrary, the full physical information about the system is obtainable from

$`\stackrel{~}{\mathrm{\Phi }}(\zeta ,\zeta ^{},0,0)=\mathrm{\Phi }(\zeta ,\zeta ^{},s,s^{}).`$ (137)
In turn, making use of (67) we get:

$`\stackrel{~}{\mathrm{\Phi }}(\zeta ,\zeta ^{},0,0)=T_{}\mathrm{exp}\left(\zeta ^{}\widehat{\stackrel{~}{a}}\right)T_+\mathrm{exp}\left(\zeta \widehat{\stackrel{~}{a}}^{}\right).`$ (138)
This is nothing but a characteristic functional of Glauber’s renowned time-normal averages of the Heisenberg field operator in the presence of the source. The source terms in the Hamiltonian are also quite recognisable; they appear in Kubo’s linear reaction theory . Introducing causal variables is thus equivalent to a nonlinear-reaction reformulation of a quantum system.

Unlike the linear reaction theory, the nonlinear reaction theory depends explicitly on which quantities are to be “measured”. Causal variables for the normal ordering, Eqs. (109)–(112), correspond to “measuring” time-normal averages of the field operators. Causal variables for the symmetric ordering, Eqs. (123)–(126), introduce another set of “measured” quantities. For example, relation (137) holds equally for $`\stackrel{~}{\mathrm{\Phi }}^W`$ and $`\mathrm{\Phi }^W`$ and in place of (138), one finds

$`\stackrel{~}{\mathrm{\Phi }}^W(\zeta ,\zeta ^{},0,0)=T_{}^W\mathrm{exp}\left({\displaystyle \frac{1}{2}}\zeta ^{}\widehat{\stackrel{~}{a}}+{\displaystyle \frac{1}{2}}\zeta \widehat{\stackrel{~}{a}}^{}\right)T_+^W\mathrm{exp}\left({\displaystyle \frac{1}{2}}\zeta ^{}\widehat{\stackrel{~}{a}}+{\displaystyle \frac{1}{2}}\zeta \widehat{\stackrel{~}{a}}^{}\right).`$ (139)
Disregarding the time orderings turns (139) into a characteristic functional of symmetrically ordered operator averages. We shall therefore term the operator ordering introduced by (139) a time-Wigner ($`T_W`$) ordering. That is, by definition,

$`\stackrel{~}{\mathrm{\Phi }}^W(\zeta ,\zeta ^{},0,0)`$ $`=`$ $`T_W\mathrm{exp}\left(\zeta ^{}\widehat{\stackrel{~}{a}}+\zeta \widehat{\stackrel{~}{a}}^{}\right),`$ (140)
where

$`T_W\mathrm{exp}\left(\zeta ^{}\widehat{a}+\zeta \widehat{a}^{}\right)`$ $`=`$ $`T_{}^W\mathrm{exp}\left({\displaystyle \frac{1}{2}}\zeta ^{}\widehat{a}+{\displaystyle \frac{1}{2}}\zeta \widehat{a}^{}\right)T_+^W\mathrm{exp}\left({\displaystyle \frac{1}{2}}\zeta ^{}\widehat{a}+{\displaystyle \frac{1}{2}}\zeta \widehat{a}^{}\right).`$ (141)
We have dropped the tildes here to emphasise that, in itself, the definition of the time-Wigner ordering does not depend on the presence or absence of the sources (nor on any other details of the dynamics).

Under time-normal ordering, the “most recent” creation (annihilation) operator becomes the leftmost (rightmost) in the product. The $`T_W`$-ordering acts in a similar manner, only symmetrised in respect of the creation and annihilation operators:

$`T_W\left\{\widehat{x}(t)\widehat{y}(t^{})\mathrm{}\widehat{z}(t^{\prime \prime })\right\}`$ $`=`$ $`{\displaystyle \frac{1}{2}}[T_W\left\{\widehat{y}(t^{})\mathrm{}\widehat{z}(t^{\prime \prime })\right\}\widehat{x}(t)`$ (143)
$`+\widehat{x}(t)T_W\left\{\widehat{y}(t^{})\mathrm{}\widehat{z}(t^{\prime \prime })\right\}].`$
Here, $`\widehat{x},\widehat{y},\mathrm{},\widehat{z}`$ stand for field operators (i.e., $`\widehat{a}`$ or $`\widehat{a}^{}`$), and $`t`$ should exceed all other time arguments in the product, $`t>t^{},\mathrm{},t^{\prime \prime }`$. Equation (143) is a recurrence relation, defining $`T_W`$ for the case of time arguments which are all different. For products of two operators, time-Wigner ordering is merely symmetric ordering,

$`T_W\widehat{x}(t)\widehat{y}(t^{})={\displaystyle \frac{1}{2}}\left[\widehat{x}(t)\widehat{y}(t^{})+\widehat{y}(t^{})\widehat{x}(t)\right]={\displaystyle \frac{1}{2}}[\widehat{x}(t),\widehat{y}(t^{})]_+,`$ (144)
where $`[\mathrm{}]_+`$ denotes an anticommutator. For higher-order products, time ordering becomes essential, so that, for three operators, assuming that $`t>t^{},t^{\prime \prime }`$,

$`T_W\widehat{x}(t)\widehat{y}(t^{})\widehat{z}(t^{\prime \prime })={\displaystyle \frac{1}{4}}[\widehat{x}(t),[\widehat{y}(t^{}),\widehat{z}(t^{\prime \prime })]_+]_+.`$ (145)
It is easy to see that the following rule holds for an arbitrary $`T_W`$-ordered product: write all possible permutations of the operators in a product, then retain only those where time arguments are arranged in a C-contour sequence (i.e., times increase then decrease).

Consider now the behaviour of time-Wigner averages for equal time arguments. Assuming $`t^{}>t`$, from (145) we find,

$`T_W\widehat{x}(t)\widehat{y}(t^{})\widehat{z}(t^{}+0)T_W\widehat{x}(t)\widehat{y}(t^{})\widehat{z}(t^{}0)={\displaystyle \frac{1}{4}}[\widehat{x}(t),[\widehat{y}(t^{}),\widehat{z}(t^{})]]_+,`$ (146)
where $`[\mathrm{}]`$ is a commutator, and $`T_W\widehat{x}(t)\widehat{y}(t^{})\widehat{z}(t^{}\pm 0)lim_{\delta t0}T_W\widehat{x}(t)\widehat{y}(t^{})\widehat{z}(t^{}\pm \delta t)`$ (i.e. the limit applies to the $`T_W`$-ordered product as a whole). This means that, unlike time-normal products, time-Wigner products are not continuous at coinciding times, nor do the limits of $`T_W`$-ordered products at equal times coincide with symmetric products. For instance,

$`T_W\widehat{a}^{}(t+0)\widehat{a}^2(t)`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left[\widehat{a}^{}(t)\widehat{a}^2(t)+\widehat{a}^2(t)\widehat{a}^{}(t)\right]`$ (147)
$`\text{W}\widehat{a}^{}(t)\widehat{a}^2(t)`$ $`=`$ $`{\displaystyle \frac{1}{3}}\left[\widehat{a}^{}(t)\widehat{a}^2(t)+\widehat{a}(t)\widehat{a}^{}(t)\widehat{a}(t)+\widehat{a}^2(t)\widehat{a}^{}(t)\right].`$ (148)
Symmetric averages may only be recovered as combinations of such limits from different directions, e.g.,

$`\text{W}\widehat{a}^{}(t)\widehat{a}^2(t)={\displaystyle \frac{1}{3}}T_W\widehat{a}^{}(t+0)\widehat{a}^2(t)+{\displaystyle \frac{2}{3}}T_W\widehat{a}^{}(t)\widehat{a}(t)\widehat{a}(t+0).`$ (149)
## IV Calculating quantum averages as classical stochastic averages

In this section, we show that there exist classical stochastic problems, such that Eqs. (74) and (88) may be interpreted as stochastic averages. In particular, Eq. (74) yields a generalisation of the positive-P representation (well known in quantum optics) to a quantum nonlinear response problem.

### A Classical stochastic response problem

Relations (114) and (129) give formal solutions to quantum response problems (formulated, respectively, in terms of time-normal and time-Wigner averages). It is instructive to compare these relations to a solution to a classical stochastic response problem. To this end, consider a c-number stochastic field $`a(t)`$, which obeys an integral equation,

$`a(t)=\alpha (t)+{\displaystyle 𝑑t^{}G(tt^{})\sigma _{\text{tot}}(t^{})},`$ (150)
where $`\alpha (t)`$ is the in-field and $`\sigma _{\text{tot}}(t)`$ is the field source. For the purposes of this paragraph, the kernel $`G(t)`$ is assumed to be regular and retarded, $`G(t)=0,t0`$, and otherwise arbitrary. We assume that the in-field, $`\alpha (t)`$, is also arbitrary. (This allows one to regularise $`G(t)`$ without changing $`\alpha (t)`$.) The full field source $`\sigma _{\text{tot}}(t)`$ consists of two parts,

$`\sigma _{\text{tot}}(t)=\sigma (t)+\sigma ^{}(t).`$ (151)
The external source, $`\sigma (t)`$, is regarded as given, while the random source, $`\sigma ^{}(t)`$, depends on the field. That is, the random source describes the field’s effective self-action (which physically originates, e.g., in interaction with a medium). As a random quantity, $`\sigma ^{}(t)`$ is fully characterised by a probability distribution, conditional on the field at the same time $`t`$:

$`\mathrm{\Pi }\left(\sigma ^{}(t)|a(t)\right).`$ (152)
Resolving formally the self-action problem results in a probability distribution over the random source, $`\sigma ^{}(t)`$, conditional on the in-field, $`\alpha (t)`$, and the external source, $`\sigma (t)`$. This is found by substituting Eq. (150) for $`a(t)`$,

$`\mathrm{\Pi }\left(\sigma ^{}(t)|\left[\alpha +G(\sigma +\sigma ^{})\right](t)\right).`$ (153)
Importantly, this expression does not contain a vicious cycle because of the assumed regular-and-retarded nature of $`G(t)`$: $`\sigma ^{}(t)`$ depends on $`\sigma ^{}(t^{})`$ only for $`t^{}<t`$.

Consider now statistical properties of the field. With the self-action resolved, these are also conditional on the in-field, $`\alpha (t)`$, and the external source, $`\sigma (t)`$. For the characteristic functional of multi-time stochastic field averages we find, (with $`\zeta (t)`$ being an arbitrary function)

$`\mathrm{\Sigma }\left(\zeta |\alpha ,\sigma \right)`$ $`=`$ $`\overline{\mathrm{exp}\left({\displaystyle 𝑑t\zeta (t)a(t)}\right)}=\overline{\text{e}^{\zeta \left[\alpha +G(\sigma +\sigma ^{})\right]}}`$ (154)
$`=`$ $`{\displaystyle \text{D}^{\mathrm{}}\sigma ^{}\text{e}^{\zeta \left[\alpha +G(\sigma +\sigma ^{})\right]}\mathrm{\Pi }\left(\sigma ^{}|\alpha +G[\sigma +\sigma ^{}]\right)}.`$ (155)
The upper bar here denotes an averaging over the statistics of $`\sigma ^{}`$, which is afterwards explicitly rewritten as a trajectorial (functional) integral. (Here and in what follows, we will make extensive use of a brief notation so as to prevent our relations from growing bulky.) Then, firstly, we pull $`G(\sigma +\sigma ^{})`$ out of $`\mathrm{\Pi }\left(\sigma ^{}|\alpha +G[\sigma +\sigma ^{}]\right)`$ by applying a shift operator:

$`\mathrm{\Pi }\left(\sigma ^{}|\alpha +G[\sigma +\sigma ^{}]\right)=\text{e}^{\frac{\delta }{\delta \alpha }G(\sigma +\sigma ^{})}\mathrm{\Pi }\left(\sigma ^{}|\alpha \right).`$ (156)
Secondly, we pull all factors except $`\mathrm{\Pi }`$ out of the functional integral, resulting in: (with $`\xi (t)`$ being another arbitrary function)

$`\mathrm{\Sigma }\left(\zeta |\alpha ,\sigma \right)`$ $`=`$ $`{\displaystyle \text{D}^{\mathrm{}}\sigma ^{}\text{e}^{\zeta \alpha +\left(\zeta +\frac{\delta }{\delta \alpha }\right)G\left(\sigma +\sigma ^{}\right)}\mathrm{\Pi }\left(\sigma ^{}|\alpha \right)}`$ (157)
$`=`$ $`\text{e}^{\zeta \alpha +\left(\zeta +\frac{\delta }{\delta \alpha }\right)G\left(\sigma +\frac{\delta }{\delta \xi }\right)}{\displaystyle \text{D}^{\mathrm{}}\sigma ^{}\text{e}^{\xi \sigma ^{}}\mathrm{\Pi }\left(\sigma ^{}|\alpha \right)}|_{\xi =0}`$ (158)
$`=`$ $`\text{e}^{\zeta \alpha +\left(\zeta +\frac{\delta }{\delta \alpha }\right)G\left(\sigma +\frac{\delta }{\delta \xi }\right)}\text{e}^{S\left(\zeta |\alpha \right)}|_{\xi =0}.`$ (159)
We have introduced a characteristic functional of cumulants of the random source conditional on the full field,

$`S\left(\zeta |a\right)=\mathrm{ln}{\displaystyle \text{D}^{\mathrm{}}\sigma ^{}\text{e}^{\xi \sigma ^{}}\mathrm{\Pi }\left(\sigma ^{}|a\right)}.`$ (160)
After further algebra, which is much assisted by Eq. (58), we arrive at

$`\mathrm{\Sigma }\left(\zeta |\alpha ,\sigma \right)=\text{e}^{\frac{\delta }{\delta a}G\frac{\delta }{\delta \xi }}\text{e}^{\zeta a+\xi \sigma }\text{e}^{S\left(\zeta |a\right)}|_{\xi =0,a=\alpha }.`$ (161)
Disregarding the averaging over the pseudo-distributions, Eqs. (114) and (129) look very much like generalisations of Eq. (161) to a pair of random fields. It is not however immediately obvious that $`S_{\text{int}}`$ and $`S_{\text{int}}^W`$ may be interpreted as characteristic functionals of cumulants of some classical noises. In fact, this turns out to be unconditionally the case for $`S_{\text{int}}`$, whereas for $`S_{\text{int}}^W`$ this interpretation holds only with discretisation of the time axis.

### B Hubbard-Stratonovich transformation: introducing noise sources constructively

#### 1 Generalising the positive-P representation to a quantum problem of nonlinear reaction

Consider a standardised real delta-correlated Gaussian noise, $`\chi (t)`$,

$`\overline{\chi (t)\chi (t^{})}=\delta (tt^{}).`$ (162)
The following relation (Hubbard-Stratonovich Transformation—HST) holds for an arbitrary function, $`x(t)`$,

$`\overline{\text{e}^{\chi x}}=\text{e}^{x^2/2},`$ (163)
so we find,

$`\mathrm{exp}S_{\text{int}}(\xi ,\xi ^{},a,a^{})=\overline{\mathrm{exp}\left[\xi \left({\displaystyle \frac{\kappa }{2}}a^2a+\sqrt{{\displaystyle \frac{i\kappa }{2}}}\chi ^{}a^{}\right)+\xi ^{}\left({\displaystyle \frac{\kappa }{2}}a^2a^{}+\sqrt{{\displaystyle \frac{i\kappa }{2}}}\chi a\right)\right]},`$ (164)
where $`\chi ^{}(t)`$ is another standardised real delta-correlated Gaussian noise, uncorrelated with $`\chi (t)`$, and the upper bar denotes averaging over the statistics of the noises. The exponent on the RHS of (164) is linear in $`\xi (t)`$ and $`\xi ^{}(t)`$. This means that for given $`\chi (t)`$ and $`\chi ^{}(t)`$, the equations for $`a(t)`$ and $`a^{}(t)`$ become regular, with the random sources being, respectively,

$`\sigma ^{}(t)`$ $`=`$ $`{\displaystyle \frac{\kappa }{2}}a^2(t)a^{}(t)+\sqrt{{\displaystyle \frac{i\kappa }{2}}}\chi (t)a(t),`$ (165)
$`\sigma ^{}(t)`$ $`=`$ $`{\displaystyle \frac{\kappa }{2}}a^2(t)a(t)+\sqrt{{\displaystyle \frac{i\kappa }{2}}}\chi ^{}(t)a^{}(t).`$ (166)
In other words, $`a(t)`$ and $`a^{}(t)`$ obey a pair of coupled stochastic integral equations,

$`a(t)`$ $`=`$ $`\alpha (t)+{\displaystyle 𝑑t^{}G(tt^{})\left[\sigma (t)+\frac{\kappa }{2}a^2(t^{})a^{}(t^{})+\sqrt{\frac{i\kappa }{2}}\chi (t^{})a(t^{})\right]},`$ (167)
$`a^{}(t)`$ $`=`$ $`\alpha ^{}(t)+{\displaystyle 𝑑t^{}G^{}(tt^{})\left[\sigma ^{}(t)+\frac{\kappa }{2}a^2(t^{})a(t^{})+\sqrt{\frac{i\kappa }{2}}\chi ^{}(t^{})a^{}(t^{})\right]}.`$ (168)
On removing the causal regularisation of $`G(t)`$, one recovers a pair of Itô stochastic differential equations,

$`i{\displaystyle \frac{da(t)}{dt}}`$ $`=`$ $`\sigma (t)+{\displaystyle \frac{\kappa }{2}}a^2(t^{})a^{}(t^{})+\sqrt{{\displaystyle \frac{i\kappa }{2}}}\chi (t^{})a(t^{}),`$ (169)
$`i{\displaystyle \frac{da^{}(t)}{dt}}`$ $`=`$ $`\sigma ^{}(t)+{\displaystyle \frac{\kappa }{2}}a^2(t^{})a(t^{})+\sqrt{{\displaystyle \frac{i\kappa }{2}}}\chi ^{}(t^{})a^{}(t^{}).`$ (170)
The fact that Itô calculus should be chosen is due to the causal regularisation of $`G(t)`$, which makes sources at time $`t`$ independent of fields at the same time, which is the characteristic property of Itô calculus.

Without the external sources (i.e., with $`\sigma =\sigma ^{}=0`$), equations (169), (170) are well known in quantum optics under the name of the positive-P representation. In that form, they allow one to calculate time-normal averages of the field operators, describing radiation properties of the system. With sources included, they become applicable to a much wider assemblage of operator averages, covering the full quantum-stochastic nonlinear reaction problem.

#### 2 Positive-W representation

At a first glance, factorising noise terms in $`S_{\text{int}}^W`$ does not pose a problem either. Introducing a standardised complex delta-correlated Gaussian noise, $`\eta (t)`$, such that

$`\overline{\eta (t)\eta (t^{})}=0,\overline{\eta (t)\eta ^{}(t^{})}=\delta (tt^{}),`$ (171)
allows one to factorise arbitrary products:

$`\text{e}^{xy}=\overline{\text{e}^{x\eta +y\eta ^{}}}.`$ (172)
We shall write the real and complex HST’s, Eqs. (163) and (172), as:

$`{\displaystyle \frac{x^2}{2}}\stackrel{\chi }{}x\chi ,xy\stackrel{\eta }{}x\eta +y\eta ^{}.`$ (173)
These relations imply the definition of the corresponding noises as standardised Gaussian noises (respectively, real and complex ones). Then, e.g., (cf. Eq. (133))

$`{\displaystyle \frac{\kappa }{8}}\xi ^2\xi a\stackrel{\eta }{}\eta {\displaystyle \frac{\xi ^2}{2}}\eta ^{}{\displaystyle \frac{\kappa }{4}}\xi a\stackrel{\chi }{}\chi \sqrt{\eta }\xi ^{}\eta ^{}{\displaystyle \frac{\kappa }{4}}\xi a.`$ (174)
Consider now the random quantity, $`\nu (t)=\chi (t)\sqrt{\eta (t)}`$. $`\chi (t)`$ and $`\eta (t)`$ are independent, and we have,

$`\overline{\nu (t)\nu ^{}(t^{})}=\overline{\chi (t)\chi (t^{})}\overline{\sqrt{\eta (t)\eta ^{}(t^{})}}=\delta (tt^{})\sqrt{\pi \delta (tt^{})}=\delta (tt^{})\sqrt{\pi \delta (0)}.`$ (175)
That is, we have encounted an infinity (divergence). The underlying reason is certainly that both $`\chi (t)`$ and $`\eta (t)`$ are highly singular functions and the “quantity” $`\nu (t)`$ is simply not defined.

For all practical purposes (such as computer simulations) stochastic differential equations will be replaced by finite-difference equations. We therefore assume an equidistant discretisation of the time axis, with $`\mathrm{\Delta }t`$ being the step of the time grid. The integration and the $`\delta `$-function are understood as usual,

$`{\displaystyle 𝑑t}\mathrm{\Delta }t{\displaystyle \underset{t}{}},\delta (tt^{})\mathrm{\Delta }t^1\delta _{tt^{}},`$ (176)
where $`\delta _{tt^{}}`$ is the Kronecker symbol on the time grid (i.e. $`\delta _{tt^{}}=1`$ if $`t=t^{}`$ and 0 otherwise). We renormalise the standardised Gaussian sources,

$`\eta (t)=\mathrm{\Delta }t^{1/2}\overline{\eta }(t),\chi (t)=\mathrm{\Delta }t^{1/2}\overline{\chi }(t),`$ (177)
making them $`\delta _{tt^{}}`$ (Kronecker symbol) correlated, not $`\delta (tt^{})`$ ($`\delta `$-function) correlated:

$`\overline{\overline{\eta }(t)\overline{\eta }^{}(t^{})}=\delta _{tt^{}},\overline{\overline{\chi }(t)\overline{\chi }(t^{})}=\delta _{tt^{}}.`$ (178)
The short upper bar always marks the Kronecker-correlated noises (and should not be confused with the long one denoting the averaging; the latter always extends over the whole expression). The HSTs then become,

$`{\displaystyle \frac{x^2}{2}}\stackrel{\overline{\chi }}{}{\displaystyle \frac{x\overline{\chi }}{\sqrt{\mathrm{\Delta }t}}},xy\stackrel{\overline{\eta }}{}{\displaystyle \frac{x\overline{\eta }+y\overline{\eta }^{}}{\sqrt{\mathrm{\Delta }t}}}.`$ (179)
There are at least two ways in which the third order terms may be factorised. One of them, generalising (174), employs two complex and two real HST’s:

$`{\displaystyle \frac{\kappa }{8}}\xi ^2\xi a\stackrel{\overline{\eta }}{}p\mathrm{\Delta }t\overline{\eta }{\displaystyle \frac{\xi ^2}{2}}+q\overline{\eta }^{}\xi \stackrel{\overline{\chi }}{}\overline{\chi }\sqrt{p\overline{\eta }}\xi ^{}+q\overline{\eta }^{}\xi ,`$ (180)
$`{\displaystyle \frac{\kappa }{8}}\xi ^2\xi ^{}a^{}\stackrel{\overline{\eta }^{}}{}p^{}\mathrm{\Delta }t\overline{\eta }^{}{\displaystyle \frac{\xi ^2}{2}}+q^{}\overline{\eta }^{}\xi ^{}\stackrel{\overline{\chi }^{}}{}\overline{\chi }^{}\sqrt{p^{}\overline{\eta }^{}}\xi +q^{}\overline{\eta }^{}\xi ^{},`$ (181)
where the functions $`p(t),q(t),p^{}(t),q^{}(t)`$ obey the relations

$`p(t)q(t)={\displaystyle \frac{\kappa a(t)}{4\mathrm{\Delta }t^2}},p^{}(t)q^{}(t)={\displaystyle \frac{\kappa a^{}(t)}{4\mathrm{\Delta }t^2}},`$ (182)
and are otherwise arbitrary (they can even be random, given they are not correlated with the Gaussian noises, $`\overline{\eta }(t),\overline{\chi }(t),\overline{\eta }^{}(t),\overline{\chi }^{}(t)`$). This results in the following coupled finite-difference equations for $`a(t)`$ and $`a^{}(t)`$:

$`i\left[a(t+\mathrm{\Delta }t)a(t)\right]`$ $`=`$ $`\left\{{\displaystyle \frac{\kappa }{2}}\left[a^{}(t)a(t)2\right]a(t)+\mu (t)\right\}\mathrm{\Delta }t,`$ (183)
$`i\left[a^{}(t+\mathrm{\Delta }t)a^{}(t)\right]`$ $`=`$ $`\left\{{\displaystyle \frac{\kappa }{2}}\left[a^{}(t)a(t)2\right]a^{}(t)+\mu ^{}(t)\right\}\mathrm{\Delta }t,`$ (184)
where

$`\mu (t)`$ $`=`$ $`\overline{\chi }(t)\sqrt{p(t)\overline{\eta }(t)}+q^{}(t)\overline{\eta }^{}(t),`$ (185)
$`\mu ^{}(t)`$ $`=`$ $`\overline{\chi }^{}(t)\sqrt{p^{}(t)\overline{\eta }^{}(t)}+q(t)\overline{\eta }^{}(t).`$ (186)
If we require the average random excursion squared at each time step,

$`\mathrm{\Delta }t^2\left[\overline{\left|\mu (t)\right|^2+\left|\mu ^{}(t)\right|^2}\right]=\mathrm{\Delta }t^2\left\{\sqrt{\pi }\left[\left|p(t)\right|+\left|p^{}(t)\right|\right]+\left|q(t)\right|^2+\left|q^{}(t)\right|^2\right\},`$ (187)
to be minimal, we find,

$`\left|q(t)\right|^3={\displaystyle \frac{\sqrt{\pi }\kappa \left|a(t)\right|}{8\mathrm{\Delta }t^2}},\left|q^{}(t)\right|^3={\displaystyle \frac{\sqrt{\pi }\kappa \left|a^{}(t)\right|}{8\mathrm{\Delta }t^2}}.`$ (188)
Apart from Eq. (182)), phases of the $`p`$’s and $`q`$’s are of no physical consequence. We choose:

$`q(t)`$ $`=`$ $`\left[{\displaystyle \frac{\sqrt{\pi }\kappa a(t)}{8\mathrm{\Delta }t^2}}\right]^{1/3},p(t)={\displaystyle \frac{2}{\sqrt{\pi }}}\left[{\displaystyle \frac{\sqrt{\pi }\kappa a(t)}{8\mathrm{\Delta }t^2}}\right]^{2/3},`$ (189)
$`q^{}(t)`$ $`=`$ $`\left[{\displaystyle \frac{\sqrt{\pi }\kappa a^{}(t)}{8\mathrm{\Delta }t^2}}\right]^{1/3},p^{}(t)={\displaystyle \frac{2}{\sqrt{\pi }}}\left[{\displaystyle \frac{\sqrt{\pi }\kappa a^{}(t)}{8\mathrm{\Delta }t^2}}\right]^{2/3}.`$ (190)
For (187), we then find,

$`\mathrm{\Delta }t^2\left[\overline{\left|\mu (t)\right|^2+\left|\mu ^{}(t)\right|^2}\right]\mathrm{\Delta }t^{2/3}.`$ (191)
This should be compared with the scaling characteristic of a Wiener process, $`\mathrm{\Delta }t`$. The “shortage” of $`\mathrm{\Delta }t^{1/3}`$ results in an overall random excursion squared over time $`T`$ scaling as, $`T/\mathrm{\Delta }t^{1/3}`$. This is indeed better that the scaling, $`T/\mathrm{\Delta }t^{1/2}`$, characteristic of the “naïve” factorisation, Eq. (174), yet the improvement is only quantitative, the process still diverging in the continuous time limit. From a practical perspective, this means that sampling noise increases with decreasing time step.

Another way of factorising the third order terms employs three complex HSTs:

$`{\displaystyle \frac{\kappa }{8}}\xi ^2\xi a`$ $`\stackrel{\overline{\eta }}{}`$ $`{\displaystyle \frac{p\mathrm{\Delta }t}{2}}\overline{\eta }\xi ^{}\xi +q\xi ^{}\overline{\eta }^{},`$ (192)
$`{\displaystyle \frac{\kappa }{8}}\xi ^2\xi ^{}a^{}`$ $`\stackrel{\overline{\eta }^{}}{}`$ $`{\displaystyle \frac{p^{}\mathrm{\Delta }t}{2}}\overline{\eta }^{}\xi ^{}\xi +q^{}\xi \overline{\eta }^{},`$ (193)
$`{\displaystyle \frac{p\overline{\eta }+p^{}\overline{\eta }^{}}{2}}\xi ^{}\xi \mathrm{\Delta }t`$ $`\stackrel{\overline{\eta }^{}}{}`$ $`\sqrt{p\overline{\eta }+p^{}\overline{\eta }^{}}\left(r\overline{\eta }^{}\xi ^{}+r^{}\overline{\eta }^{}\xi \right).`$ (194)
As above, the $`p`$’s and $`q`$’s here are arbitrary functions obeying Eq. (182), and $`r(t)r^{}(t)=1/2`$. That is, we recover the generic equations (180) and (181), this time with

$`\mu (t)`$ $`=`$ $`r\overline{\eta }^{}\sqrt{p\overline{\eta }+p^{}\overline{\eta }^{}}+q\overline{\eta }^{},`$ (195)
$`\mu ^{}(t)`$ $`=`$ $`r^{}\overline{\eta }^{}\sqrt{p\overline{\eta }+p^{}\overline{\eta }^{}}+q^{}\overline{\eta }^{}.`$ (196)
(Note that the terms proportional to the $`q`$s have exchanged places.) Minimising the random excursion squared, Eq. (187), we find the optimised values of free parameters to be

$`r(t)`$ $`=`$ $`r^{}(t)={\displaystyle \frac{1}{\sqrt{2}}},`$ (197)
$`q(t)`$ $`=`$ $`\left[{\displaystyle \frac{\pi \kappa ^2a^3(t)}{64\mathrm{\Delta }t^4\left(\left|a(t)\right|+\left|a^{}(t)\right|\right)}}\right]^{1/6},`$ (198)
$`q^{}(t)`$ $`=`$ $`\left[{\displaystyle \frac{\pi \kappa ^2a^3(t)}{64\mathrm{\Delta }t^4\left(\left|a(t)\right|+\left|a^{}(t)\right|\right)}}\right]^{1/6},`$ (199)
with the $`p`$s being recovered from Eq. (182). Evidently both ways of factorising the third order terms result in the same scaling of the noises.

### C Numerical experiment

We have calculated the time dependence of the quadrature amplitude,

$`X(t)=\widehat{a}(t)+\widehat{a}^{}(t)=\overline{a(t)+a^{}(t)},`$ (200)
assuming that the initial state of the oscillator is a coherent state $`|\alpha `$. An analytical expression for this quantity may be found in . In Fig. 2, we plot the relative error in $`X(t)`$,

$`\delta X={\displaystyle \frac{X_{\text{simulated}}(t)X_{\text{exact}}(t)}{X_{\text{exact}}(t)}},`$ (201)
where $`X_{\text{simulated}}(t)`$ is found via stochastic simulations of the “$`+`$W” equations, (183)–(184), for $`\kappa =4`$ and $`\alpha =1`$. To have an idea how important the third order noises are, we also calculated $`X(t)`$ using the truncated Wigner representation . The latter is obtained by dropping the noise sources in Eqs. (183)–(184). The resulting equations are readily solved analytically; the averaging over the initial W-distribution does not pose a problem. The relative error in $`X(t)`$ found via the truncated Wigner (“$``$W”) representation is also shown in Fig. 2. As is evident from inspection of Fig. 2, the “$`+`$W” result is a manifest improvement over the “$``$W” one. The bad news is, however, that the sampling error remains rather large, even after averaging over more than half a billion trajectories. A further piece of bad news is that simulations fell victim to numerical instability shortly after the maximal time shown in Fig. 2. These problems notwithstanding, the very fact that third order noises my be subject to stochastic simulations has been clearly established.

## V The cookbook

### A The general recipe

We now summarise our results, recipe-style, the way they should be applied to practical calculations. For an $`n`$-mode system, Schrödinger-picture annihilation and creation operators are vectors in respect of the mode index,

$`\widehat{𝒂}_\text{S}=\left\{\widehat{a}_{\text{S}k}\right\},\widehat{𝒂}_\text{S}^{}=\left\{\widehat{a}_{\text{S}k}^{}\right\},k=1,\mathrm{},n,`$ (202)
as are the interaction-picture field operators,

$`\widehat{𝒂}(t)=\left\{\widehat{a}_k(t)\right\},\widehat{𝒂}^{}(t)=\left\{\widehat{a}_k^{}(t)\right\},k=1,\mathrm{},n.`$ (203)
The system Hamiltonian consists, as usual, of the free and interaction Hamiltonians, cf. Eq. (1). The free Hamiltonian is an arbitrary matrix, $`H=\left\{H_{kk^{}}\right\}`$, $`k,k^{}=1,\mathrm{},n`$, in the mode indices. The free Schrödinger equation thus reads

$`i{\displaystyle \frac{d\widehat{𝒂}(t)}{dt}}=H\widehat{𝒂}(t).`$ (204)
Quantum-classical mappings yield the generic system of $`2\times n`$ equations,

$`i{\displaystyle \frac{d𝒂(t)}{dt}}`$ $`=`$ $`H𝒂(t)+𝝈(t)+𝝈^{}(t),`$ (205)
$`i{\displaystyle \frac{d𝒂^{}(t)}{dt}}`$ $`=`$ $`H𝒂^{}(t)+𝝈^{}(t)+𝝈^{}(t),`$ (206)
for $`2\times n`$ c-number random fields,

$`𝒂(t)=\left\{a_k(t)\right\},𝒂^{}(t)=\left\{a_k^{}(t)\right\},k=1,\mathrm{},n.`$ (207)
The vectors of random sources, $`𝝈^{}(t),𝝈^{}(t)`$, depend on the interaction and on the actual type of operator ordering underlying the mapping.

Consider firstly the case of time-normal ordering. Without the given sources (i.e., with $`𝝈(t)=𝝈^{}(t)=0`$), stochastic averages of the random fields coincide with the quantum averages of the time-normally-ordered products of the Heisenberg field operators, $`\widehat{a}_k(t),\widehat{a}_k^{}(t)`$, $`k=1,\mathrm{},n`$. With the given sources, Eqs. (205)–(206) constitute a constructive nonlinear-reaction formulation of the quantum system. The random sources are deduced from the normal form of the interaction Hamiltonian,

$`\widehat{}_{\text{int}}=:h(\widehat{𝒂}_\text{S},\widehat{𝒂}_\text{S}^{}):.`$ (208)
Namely, one should calculate the functional, (cf. Eq. (115))

$`S_{\text{int}}(𝝃,𝝃^{},𝒂,𝒂^{})`$ $`=`$ $`i{\displaystyle 𝑑t\left[h(𝒂^{}(t),𝒂(t)i𝝃(t))h(𝒂^{}(t)+i𝝃^{}(t),𝒂(t))\right]}.`$ (209)
Since $`\widehat{}_{\text{int}}`$ is Hermitian, this can be further simplified resulting in

$`S_{\text{int}}(𝝃,𝝃^{},𝒂,𝒂^{})`$ $`=`$ $`i{\displaystyle 𝑑th(𝒂^{}(t),𝒂(t)i𝝃(t))}+\text{conj.},`$ (210)
where conjugation acts as a formal Hermitian transformation (i.e., it interchanges quantities with and without dagger, $`𝒂(t)𝒂^{}(t)`$, $`𝝃(t)𝝃^{}(t)`$, and complex-conjugates other c-numbers). One should then factorise powers and products of $`\xi `$s, employing a suitable set of Hubbard-Stratonovich transformations (see Eqs. (173)), until an expression emerges which is already linear in respect of the $`\xi `$s:

$`S_{\text{int}}(𝝃,𝝃^{},𝒂,𝒂^{})\mathrm{}\stackrel{\mathrm{}}{}\mathrm{}\stackrel{\mathrm{}}{}\mathrm{}{\displaystyle 𝑑t\underset{k=1}{\overset{n}{}}\left[\xi _k(t)\sigma _k^{}(t)+\xi _k^{}(t)\sigma _k^{}(t)\right]}.`$ (211)
(Note that terms without $`\xi `$s in $`S_{\text{int}}`$ always cancel.) The $`\sigma `$s thus obtained are exactly those which appear in Eqs. (205)–(206). If $`S_{\text{int}}`$ is at most quadratic in $`\xi `$’s, Eqs. (205)–(206) are a system of genuine Itô stochastic differential equations. Otherwise they can only be interpreted as difference equations over a finite time step $`\mathrm{\Delta }t`$.

In the case of the time-Wigner ordering, this recipe applies with the replacement $`S_{\text{int}}S_{\text{int}}^W`$, where

$`S_{\text{int}}^W(𝝃,𝝃^{},𝒂,𝒂^{})`$ $`=`$ $`i{\displaystyle }dt[h^W(𝒂^{}(t){\displaystyle \frac{i}{2}}𝝃^{}(t),𝒂(t){\displaystyle \frac{i}{2}}𝝃(t))`$ (213)
$`h^W(𝒂^{}(t)+{\displaystyle \frac{i}{2}}𝝃^{}(t),𝒂(t)+{\displaystyle \frac{i}{2}}𝝃(t))]`$
$`=`$ $`i{\displaystyle 𝑑th^W(𝒂^{}(t)\frac{i}{2}𝝃^{}(t),𝒂(t)\frac{i}{2}𝝃(t))}+\text{conj.},`$ (214)
and $`h^W`$ is the symmetric form of the interaction Hamiltonian,

$`\widehat{}_{\text{int}}=\text{W}\left\{h^W(\widehat{𝒂}_\text{S},\widehat{𝒂}_\text{S}^{})\right\}.`$ (215)
We are however not aware of any nonlinear system where $`S_{\text{int}}^W`$ would turn out to be quadratic in the $`\xi `$’s.

### B Degenerate OPO

#### 1 Generalised positive-P representation

We shall now illustrate the general recipe by deriving the positive-P and positive-W representations for a degenerate optical parametric oscillator (OPO). The degenerate OPO consists of two coupled oscillators described by the Hamiltonian,

$`=\omega \widehat{a}_{\text{S}1}\widehat{a}_{\text{S}1}^{}+2\omega \widehat{a}_{\text{S}2}\widehat{a}_{\text{S}2}^{}+{\displaystyle \frac{i\kappa }{2}}\left[\widehat{a}_{\text{S}1}^2\widehat{a}_{\text{S}2}\widehat{a}_{\text{S}1}^2\widehat{a}_{\text{S}2}^{}\right].`$ (216)
Then, in brief notation,

$`S_{\text{int}}(𝝃,𝝃^{},𝒂,𝒂^{})=`$ $`i{\displaystyle \frac{i\kappa }{2}}\left[a_1^2\left(a_2i\xi _2\right)\left(a_1i\xi _1\right)^2a_2^{}\right]+\text{conj.}`$ (217)
$`=`$ $`i\kappa \xi _1a_1a_2^{}+i{\displaystyle \frac{\kappa }{2}}\xi _2a_1^2{\displaystyle \frac{\kappa }{2}}\xi _1^2a_2^{}+\text{conj.}`$ (218)
$`\stackrel{\chi ,\chi ^{}}{}`$ $`+i\xi _1^{}\left(\kappa a_1^{}a_2+\chi \sqrt{\kappa a_2}\right)`$ (222)
$`i\xi _1\left(\kappa a_1a_2^{}+\chi ^{}\sqrt{\kappa a_2^{}}\right)`$
$`i\xi _2^{}{\displaystyle \frac{\kappa }{2}}a_1^2`$
$`+i\xi _2{\displaystyle \frac{\kappa }{2}}a_1^2,`$
where $`\chi (t),\chi ^{}(t)`$ are a pair of independent real $`\delta `$-correlated Gaussian noises. The positive-P representation for the OPO, generalised to the response problem, then reads,

$`{\displaystyle \frac{da_1(t)}{dt}}`$ $`=`$ $`i\sigma _1(t)+\kappa a_1^{}(t)a_2(t)+\chi (t)\sqrt{\kappa a_2(t)},`$ (223)
$`{\displaystyle \frac{da_1^{}(t)}{dt}}`$ $`=`$ $`+i\sigma _1^{}(t)+\kappa a_1(t)a_2^{}(t)+\chi ^{}(t)\sqrt{\kappa a_2^{}(t)},`$ (224)
$`{\displaystyle \frac{da_2(t)}{dt}}`$ $`=`$ $`i\sigma _2(t){\displaystyle \frac{\kappa }{2}}a_1^2(t),`$ (225)
$`{\displaystyle \frac{da_2^{}(t)}{dt}}`$ $`=`$ $`+i\sigma _2^{}(t){\displaystyle \frac{\kappa }{2}}a_1^2(t).`$ (226)
This derivation is strikingly simple and straightforward, and compares very favorably to the common derivation based on phase-space techniques  (not to mention that Eqs. (223)–(226) allow for a much deeper insight).

#### 2 Positive-W representation

Similarly,

$`S_{\text{int}}^W(𝝃,𝝃^{},𝒂,𝒂^{})=`$ $`i{\displaystyle \frac{i\kappa }{2}}\left[\left(a_1^{}i\xi _1^{}/2\right)^2\left(a_2i\xi _2/2\right)\left(a_1i\xi _1/2\right)^2\left(a_2^{}i\xi _2^{}/2\right)\right]+\text{conj.}`$ (227)
$`=`$ $`i\kappa \xi _1a_1a_2^{}+i{\displaystyle \frac{\kappa }{2}}\xi _2a_1^2+{\displaystyle \frac{i\kappa }{8}}\xi _1^2\xi _2^{}+\text{conj.}`$ (228)
$`\stackrel{\eta ,\eta ^{}}{}`$ $`i\kappa \xi _1a_1a_2^{}+i{\displaystyle \frac{\kappa }{2}}\xi _2a_1^2{\displaystyle \frac{p^{}}{2}}\eta ^{}\xi _1^2+iq\eta \xi _2^{}+\text{conj.}`$ (229)
$`\stackrel{\chi ,\chi ^{}}{}`$ $`+i\xi _1^{}\left(\kappa a_1^{}a_2+\chi \sqrt{p\eta ^{}}\right)`$ (233)
$`i\xi _1\left(\kappa a_1a_2^{}+\chi ^{}\sqrt{p^{}\eta ^{}}\right)`$
$`+i\xi _2^{}\left({\displaystyle \frac{\kappa }{2}}a_1^2+q\eta \right)`$
$`i\xi _2\left({\displaystyle \frac{\kappa }{2}}a_1^2+q^{}\eta ^{}\right),`$
resulting in the positive-W representation for the OPO:

$`{\displaystyle \frac{da_1(t)}{dt}}`$ $`=`$ $`i\sigma _1(t)+\kappa a_1^{}(t)a_2(t)+\chi (t)\sqrt{p(t)\eta ^{}(t)},`$ (234)
$`{\displaystyle \frac{da_1^{}(t)}{dt}}`$ $`=`$ $`+i\sigma _1^{}(t)+\kappa a_1(t)a_2^{}(t)+\chi ^{}(t)\sqrt{p^{}(t)\eta ^{}(t)},`$ (235)
$`{\displaystyle \frac{da_2(t)}{dt}}`$ $`=`$ $`i\sigma _2(t){\displaystyle \frac{\kappa }{2}}a_1^2(t)+q(t)\eta (t),`$ (236)
$`{\displaystyle \frac{da_2^{}(t)}{dt}}`$ $`=`$ $`+i\sigma _2^{}(t){\displaystyle \frac{\kappa }{2}}a_1^2(t)+q^{}(t)\eta ^{}(t).`$ (237)
Here, $`\eta (t),\eta ^{}(t)`$ and $`\chi (t),\chi ^{}(t)`$ are pairs of standardised Gaussian noises (respectively, complex and real ones), and the functions $`p(t),p^{}(t),q(t),q^{}(t)`$ obey the relations,

$`p(t)q^{}(t)=p^{}(t)q(t)={\displaystyle \frac{\kappa }{4}}.`$ (238)
Equations (234)–(237) are difference equations over a finite time step $`dt=\mathrm{\Delta }t`$, with the noises normalised such that

$`\overline{\chi (t)\chi (t^{})}=\overline{\chi ^{}(t)\chi ^{}(t^{})}=\overline{\eta (t)\eta ^{}(t^{})}=\overline{\eta ^{}(t)\eta ^{}(t^{})}=\delta (tt^{})={\displaystyle \frac{\delta _{tt^{}}}{\sqrt{dt}}}.`$ (239)
The $`p`$’s and $`q`$’s will be chosen so as to minimise the weighted sum of random increments squared,

$`dt^2\left\{\overline{\left|\chi (t)\sqrt{p(t)\eta ^{}(t)}\right|^2}+\overline{\left|\chi ^{}(t)\sqrt{p^{}(t)\eta ^{}(t)}\right|^2}+w\left[\overline{\left|q(t)\eta (t)\right|^2}+\overline{\left|q^{}(t)\eta ^{}(t)\right|^2}\right]\right\}`$ (240)
$`=p(t)\sqrt{\pi dt}+{\displaystyle \frac{w\kappa ^2dt}{16p^2(t)}}+\text{conj.},`$ (241)
where we have chosen the $`p`$’s to be real and positive. The weighting factor $`w>0`$ allows one to redistribute the noise between the fundamental and the subharmonic, and may be fine-tuned by trial and error when running stochastic simulations. Minimising (241) yields:

$`p(t)=p^{}(t)={\displaystyle \frac{1}{2}}\sqrt[3]{w\kappa ^2\sqrt{dt/\pi }}.`$ (242)
The $`q`$’s are then recovered from (238). Accounting for (242), all noise increments in the positive-W equations scale as, $`dt^{1/3}`$.