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

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


Introduction
Phase-space representations of quantum mechanics such as the Wigner, P, positive-P, Q and related approaches are a powerful tool for the study and understanding of quantum mechanics [14,50,62,93,118]. Their use in recent times has been directed particularly as a tool Piotr Deuar: deuar@ifpan.edu.pl in quantum information science (see [58] and [129] for recent reviews) and for the simulation of large-scale quantum dynamics. Negativity of Wigner, P, and other phase-space quasi-distributions is a major criterion for quantumness and closely related to contextuality and nonlocality in quantum mechanics [58,65]. The inability to interpret the Glauber-Sudarshan P [67,132] and Wigner distributions in terms of a classical probability density is the fundamental benchmark for quantum light [88,128]. Quasi-probability representations arise naturally when looking for hidden variable descriptions and an ontological model of quantum theory [16,58,93,127]. Wigner and P-distribution negativity are considered a resource for quantum computation [8,58] and closely related to magic state distillation and quantum advantage [138]; the Glauber-Sudarshan P distribution defines a hierarchy of nonclassicalities and nonclassicality witnesses [88,139] and can be experimentally reconstructed [3,83]. Moreover, P, Q and Wigner distributions and have been used for simulations of models important for quantum information topics such as open spin-qubit systems [136] and boson sampling [53,112,113].
For simulation of quantum mechanics, the prime advantage of the phase-space approach is that its computational cost typically grows only linearly with system size even in interacting systems of many particles. Therefore it provides a route to the non-perturbative treatment of the quantum dynamics of large systems. It has been used for solving and simulating a multitude of problems in various physical fields: e.g. quantum optics [21,28,30,32,52,54,84,91,110], ultracold atoms [22,39,40,49,74,82,94,95,101,102,108,125,131,134,144], fermionic systems [6,7,27,29], spin systems [11,106,107], nuclear physics [137], dissipative systems in condensed matter [26,45,143], or cosmology [111]. The vast majority of such calculations so far have considered only equal time correlations and observables.
Multi-time correlations are also important for answering many physical questions which cannot necessarily be dealt with by monitoring the time dependence of equal time observables [13,62,92,118,139]. For example, the determination of lifetimes in an equi-librium or stationary state, response functions, out-oftime-order correlations (OTOCs), or finding the time resolution required to observe a transient phenomenon. However their calculation in the phase-space framework has been restricted because few general results have been available. What is known is largely restricted to time-and normal-ordered correlations in the Glauber-Sudarshan P representation [5,62], time-symmetric ones in the truncated Wigner representation [13,115], or the corresponding linear response corrections for closed systems. The Glauber-Sudarshan P representation results provide a particularly intuitive framework, simply replacing Heisenberg operators a † (t) and a(t) with timedependent phase-space variables in normal-ordered observables such as a † (t 1 ) a(t 2 ). They also apply to open systems [2,62]. Such ordering corresponds to general photon counting measurements [67,68,81]. Still, while scaling with system size is excellent, both above approaches usually end up being only approximate because one is forced to omit part of the full quantum mechanics in their numerical implementations [62]. Hence, a broader application of full quantum phase-space methods to multi-time observables would be advantageous. This is especially so, since recent years have shown a lot of interest in such quantities, be it in the study of nonclassicality in the time dimension [9,89,104,119], time crystals [55,135,142,145], quantum technologies [78,96,103] where time correlations and time resolution are crucial, methods development [86,115], and the OTOCs that have recently found to be important for the study of quantum scrambling [15,64,133], quantum chaos [98,124] and many-body localisation [56,70]. Therefore, this paper sets out to extend the infrastructure available for multi-time observables with phasespace methods.
The positive-P representation [50] does not suffer from the limitations of the Glauber-Sudarshan P (seen e.g. in Sec. 7.3) due to its anchoring in doubled phase space. This allows it to incorporate the full quantum mechanics for most Hamiltonians, including all twobody interactions. Importantly, as we shall see, the same applies for its differently ordered analogues like the doubled-Wigner and doubled-Q. While the positive-P is known to be limited to short times for closed systems due to a noise amplification instability [38,66], this abates in dissipative systems [45,66]. Recent work has shown that a very broad range of driven dissipative Bose-Hubbard models can be simulated with the positive-P into and beyond the stationary state [45]. This covers systems of current interest such as micropillars, transmon qubits, and includes strong quantum effects such two-photon interference in the unconventional photon blockade [10,96], making its extension quite timely. Therefore, after basic background in Sec. 2, the first order of business in this paper in Sec. 3 is to extend the long known Glauber-Sudarshan P results on time-and normal-ordered multi-time correlations to the positive-P representation. This is sort of an obvious extension, but was missing from the literature, and allows complete quantum calculations. It also lets us explain the ins and outs and set the stage for the less obvious and broader extensions that follow: The first one of those in Sec. 4.1 considers general representations and reformulates the case of Q representations in an operational form that can be extended to doubled phase-space. This gives stochastic access to anti-normal ordered observables such as a(t 1 ) a † (t 2 ). It is then shown in Sec. 5 how the Gaussian convolution relationship between Wigner, P, and Q distributions allows one to easily evaluate anti-normal ordered correlations in Wigner and P. This is extended also to a wide range of mixed ordered observables such as a † (t 1 ) a(t 2 ) a † (t 3 ) a(t 4 ) in Sec. 5.4, and further extended to cover doubled phase-space representations in Sec. 6. In Sec. 7, the use and accuracy of this approach is demonstrated on the numerical example of the unconventional photon blockade, a system that may have application to the creation of single-photon sources. Finally, Appendix A gives details of a robust algorithm for integrating the resulting phase-space stochastic equations, a matter that has been somewhat neglected in the literature to date.

Positive-P representation
Consider an M -mode system (modes 1, . . . , j, . . . , M ) with bosonic annihilation operators a j for the jth mode. The density operator of the system is written in terms of bosonic coherent states of complex amplitude α j defined relative to the vacuum: using [50]: The kernel Λ involves "ket" |α j j and "bra" β * j | j states described by separate and independent variables. Let us define the container variable as a shorthand to hold the full configuration information. Individual complex variables λ µ can be indexed by µ = 1, . . . , 2M . The positive-P distribution P + is guaranteed real and non-negative [50]. This is essential for the utility of the method, which lies in its ability to represent full quantum mechanics using stochastic trajectories of the samples λ of the distribution P + . Notably, the size of these samples scales merely linearly with M , allowing for the simulation of the full quantum dynamics of systems with even millions of modes [39,82]. The distribution (2) can be compared to the simpler and more widely known Glauber-Sudarshan P representation [67,132], which uses a single set of coherent state amplitudes and ties the "bra" and "ket" amplitudes to be equal. It is written where the multimode coherent amplitude is α = {α 1 , . . . , α M }. This representation gives a well behaved representation of only a subset of possible quantum states [20,128], though a very useful one that includes Gaussian density operators [4]. For full computational utility of the positive-P method, all quantum mechanical actions should be able to be written in terms of the samples λ. Central to this are the following identities: A general master equation for the density matrix written in Lindblad form with reservoir operators R n is This can be transformed by standard methods [62], using identities (5), to a partial differential equation (PDE) for P + of the following general form: with maximum differential order n max and coefficients F that depend on the details of (6). Most cases have n max 2, in which case this is a Fokker-Planck equation (FPE): An exception occurs if irreducible higher order partial derivatives appear, e.g. due to explicit three-body interactions in the Hamiltonian. Most first principles models use only two-particle interactions, though. Standard methods convert an FPE like (8) to the following Ito stochastic equations of the samples [62]: where i.e. D = BB T in matrix notation, and ξ σ (t) are independent real white noises with variances The notation · stoch indicates a stochastic average over the samples. The equations (9) give us quantum mechanical evolution in terms of the samples, which becomes ever more exact as the number of samples grows.
Quantum expectation values at a given time are evaluated via The equivalence can be written more explicitly in terms of S individual samples λ (u) labelled by u = 1, . . . , S as: The mode labels in the above can be in any combination, provided all annihilation operators are to the right of all creation operators (which is called "normal ordering"). Any single-time operator can be expressed as a sum of normally ordered terms like (12).

Multi-time averages in open systems
The situation is much more complicated when the operators in the expectation value are not evaluated all at the same time. The root of the difficulty is that multi-time commutation relations usually depend nontrivially on the full system dynamics, and a reduction of arbitrary multi-time operators to a normal-ordered form is not generally possible.
To describe what is meant at the operator level it is first helpful to introduce the two-sided evolution opera-torV (t 1 , t 2 ) such that evolution by the master equation (6) can be summarised as ρ(t 2 ) =V (t 2 , t 1 ) ρ(t 1 ). (14) Notably, the evolution operator has the semigroup property [130]V (t 3 , t 2 )V (t 2 , t 1 ) =V (t 3 , t 1 ). We will use the convention that two-sided operators, indicated by a breve˘, act on everything to their right. Hencȇ Let us also define the time-ordered form via where the A j (t) and B j (t) are Heisenberg picture operators, and the times obey There are N operators A with times labelled t 1 , . . . , t N increasing to the right, inward, and M operators B with times labelled s 1 , . . . , s M increasing to the left, also inward. The location of the inner "meeting point" is arbitrary, and the number of A or B operators can also be zero. No particular constraints are imposed on the A and B operators, except that they should be single time quantities. These operators refer to measurements made at the respective times τ r , when the density operator was ρ(τ r ). Between measurements the state evolves according to the master equation (6) (i.e. (14)). It has been shown [62] that multi-time correlations that correspond to sequences of measurements can always be written in the above form (16). Therefore this time ordering is not an arbitrary one, and not particularly restrictive in itself. It has also been shown that a general time-ordered correlation function (16) obeying (17) can be written in a form that uses theV . To do so it is necessary to be careful about operator ordering. Following [62], let us order all the times t p and s q in the correlation in sequence from earliest to latest. Let us then rename them τ r so that with R = N + M. We also define the corresponding two-sided operatorsF ȓ Then, · · ·F 2V (τ 2 , τ 1 )F 1 ρ(τ 1 )) . (20) 2.3 Multi-time correlations in the Glauber-Sudarshan P representation provided the times are ordered according to (17). The requirement that the operators in (21) be normally ordered is an additional constraint on top of time-ordering (17), but one that leads to all operators sorted as they occur in photo-counting theory [81]. It covers a very large subset of the potentially physically interesting correlations. The Gardiner approach is more amenable to extension to doubled phase space and will be used in what follows. This ordering can be contrasted with the time-symmetric ordering for which straightforward truncated Wigner correspondences for closed system evolution were found in [13,115,118]. Examples of time-symmetric ordered quantities are 1 2 a † (t 2 ) a(t 1 ) + a(t 1 ) a † (t 2 ) and 3 Time ordered moments in the positive-P representation

Normally ordered observables
To derive an expression like (21) for the positive-P representation, we will follow Gardiner's approach [62] that was previously used for the Glauber-Sudarshan P. Smaller steps will, however, be taken here to draw attention to a few subtleties that will be necessary later.

First order correlation function
First consider the correlation Comparing to (20), we can identify A 1 = a † j , B 1 = a k , and the times t 1 = t , s 1 = t. For this low order correlation, the time ordering (17) sets no additional conditions. There are two possibilities for τ 1 , depending on whether time t or t is later. Consider first t t, so that τ 1 = t, τ 2 = t . Using (20) and (19) we have that The 2nd line follows from the cyclic property of traces.
The {·} is kept for now for clarity. We can see that the evolution operator from t to t acts to the right on all the quantities at time t, while the later-time operator a † j (t ) acts only on the evolved quantities to its right. This makes intuitive physical sense. In the second case of t < t, we have τ 1 = t , τ 2 = t, and get that This again has the intuitive form of the operatorV acting to the right on all the earlier-time quantities. Take the first case with t t. Expressing the density matrix in (24) in the positive-P representation (2a), The 2nd line follows from application of (5a). We cannot do the same for a † j (t ) yet, because the kernel Λ finds itself inside the prior action of theV operator.
To deal with this, consider now the action of the evolution operatorV on distributions P + . i.e. the action of the PDE (7). If we define the conditional distribution P(λ, t |λ, t) as the solution of this PDE at time t t starting from the initial condition δ 4M (λ − λ), i.e. the "propagator", then it can be used to formally writȇ This now contains no more two-sided operators. Through this convolution, P + is expressed in λ variables which accompany the earliest time t to aid for later interpretation as part of a joint probability, while the kernel Λ is expressed in the variables λ that accompany later times, ready for application of the next operator identity. Notice that there are no particularly stringent assumptions about P + for (27) itself to apply. For example, (27) applies equally well if one replaces P + with some complex distribution function P . This point will soon be useful. However, there was an assumption that the propagator P is well behaved. This is certainly true if the PDE was of a Fokker-Planck form (8), and therefore is always justified if we have an exact mapping of the master equation (6) to stochastic equations (9). However, in some other cases of third/higher order terms in the PDE, it might not. We will not be concerned with such cases here. Now in (26),V is acting on a distribution P (λ, t) = α k P + (λ, t). Using (27) we get The two-way operator acting on the right that wasV , has now been gotten rid of, by virtue of being incorporated in the propagator P. Therefore, the remaining operator a † j can now be shifted to the right due to the cyclic property of the trace, and then processed via (5d) as so: The last line follows from Tr Λ = 1, which is pre-set by the definition (2b). The quantity PP + is the just the joint probability of having configuration λ at time t and configuration λ at time t . (Provided P is well behaved, positive, real, as mentioned before, which is the case for any model fully described by an FPE (8)). Therefore, At this stage we can identify probability with stochastic realisations. The noises ξ µ (t) introduced during time evolution are independent from each other, independent at each time step, independent for each sample's trajectory. They are also independent of any other random variables used to produce the initial ensemble {λ (1) , . . . , λ (S) }(t) that samples P + (λ, t). Therefore, the evolved configuration λ (u) (t ) at time t depends only on mutually independent random variables that consist of λ (u) (t) and the noise history of the uth trajectory. As a result, the combination of initial configuration λ (u) (t) and the evolved configuration λ (u) (t ), together form an unbiased sample of the joint distribution P (λ, t ; λ, t). We arrive then at the final result that The procedure for the case t < t gives a result analogous to (31): which once again leads to (32). It is especially important to note that the quantity to be averaged comes from the time evolution of individual sample trajectories. Explicitly: This allows for very efficient calculations.

Higher order correlations
Other time and normal ordered correlations follow a similar pattern. For example, when t t t, In this case, following a similar procedure to before, one can act alternately with (5a) on the kernel to extract a factor of α, and (27) to convert the evolution operators to propagators. One finds Since the times are ordered t t t, and the conditional probabilities follow from the evolution of the FPE, λ are parent variables of the λ and so on, and the product of conditional probabilities is just the joint probability. Hence (36) is (37) and Working similarly, using just the 1st and 4th identities in (5), one readily but somewhat cumbersomely finds that the stochastic estimator for the general time-andnormal-ordered correlation function is Here, of course (17) must hold, and the stochastic averaging is over products constructed using values from the evolution of a single sample, as in (34). It confirms the suspicion and intuition that the behaviour of the positive-P representation in this regard should be similar to the earlier expression (21), for the Glauber-Sudarshan P.

Other ordering in the positive-P
Now to see the limitations of this scheme, consider the anti-normally (but time-ordered) ordered correlation with t > t. The first point to make is that we cannot rearrange this to a normal-ordered form like a † k (t) a j (t ) + δ jk and then use (39) because generally a j (t ), a † k (t) = δ jk when t = t . Instead it is some time-dependent operator. Now applying (20), (40) can be written Upon expansion, we will need to act on Λ using the 2nd identity in (5), (5b), to convert a † k to variable form. Thus This is not of a form amenable to (27). We can soldier on applying integration by parts and assuming negligible boundary terms to obtain The matter of whether boundary terms can be discarded has been studied in depth [23,34,37,66,85,126]. The summary is that one can determine operationally in a stochastic simulation whether boundary terms are negligible or not. If deemed negligible, then integration by parts is justified. In (43) we can now identify a distribution P (λ, t) = [β k − ∂ ∂α k ]P + (λ, t) to act on with (27). Doing so gives: The second operator a j (t ) can now act on Λ from the left. One obtains While this is formally acceptable (given those negligible boundary terms), and could be used for some analytic work in small systems such as demonstrated in [62], unfortunately the derivative of P + is not amenable to interpretation in terms of stochastic samples. At least not the direct samples we are investigating here. It may be partially treatable using the quantum jump and response theory approach previously applied to truncated Wigner [13,115,118], which is a topic for another time. However -in Sections 5 and 6, a different direct way to evaluate anti-normal ordered observables such as A = a j (t ) a † k (t) will be demonstrated.

Other phase-space representations
In the derivations of Sec. 3.1, one can see that the crucial aspect for obtaining a stochastically useful expression is to use only those identities which do not contain derivatives 1 . This suggests that convenient expressions for multi-time correlations similar to (39) will be obtainable whenever the operators in the correlation can be converted to phase-space variables without resorting to identities with derivatives. However, it happens that such identities without derivatives are not particularly abundant in other phase-space representations. For example, the Wigner representation [105,118,141] has derivatives in all identities, as does its dimension-doubled analogue [71,116].
1 Strictly speaking, a small exception to this appears if the derivative appears only at the final time when the only remaining operator is Λ since Tr[ ∂ ∂λµ Λ] = ∂ ∂λµ Tr[ Λ] = 0 removes any awkward terms. Such a case can, however, also be treated by an identity without any derivatives after using the cyclic property of the trace at the right step.
Notably, while the trace of the symmetric form in the Wigner representation corresponds to |α| 2 with no corrections: Tr 1 2 a † a + a a † Λ Wig = |α| 2 , this does not remove derivatives in the corresponding identity 2 . The best that appears to be achievable in this way is Hence, the path-integral and time-symmetric approach [13] is more suited to the Wigner representation. Also the phase-space representations developed for spin systems [11,100,107], contain derivatives for all identities. One notable exception is the Q representation, which admits derivative-free identities similar to (5) for antinormally ordered operators, and so is a good candidate for convenient phase-space expressions. Time-ordered anti-normal operators occur for example in the theory of photon detectors that operate via emission rather than absorption of photons [99]. Formal integral expressions for this kind of correlation were also provided in [5]. Below, the operational stochastic expressions are found using the Gardiner approach.

The case of the Q representation
The Husimi Q representation [75] is defined as and is positive for any ρ. Due to the Q distribution being defined in this explicit way, rather than the implicit form (2), observable expressions in the Q distribution have traditionally been found by simply expanding the trace: and applying the eigenvalue equation for coherent states For anti-normal ordered correlations at equal times, the cyclic property of traces gives However, this traditional approach fails with timeordered correlations. Whatever way one orders the operators, working this way on expression (20) will lead to a form 1 in which it is the latest time operatorF R (τ R ) that acts on the outer coherent states. One can convert F R to phase-space variables using (49), and move the state vectors |α or α| closer to the density matrix and the form (47). However, in the next putative step, there is no clear way to convert the evolution opera-torV (τ R , τ R−1 ) to phase space form. Moreover, when there is phase-space diffusion, the propagator P is wellbehaved only in the forward time direction, so there is no way to act with P on the outer state vectors and variables which correspond to later times than the inner ones. Therefore, we proceed in a non-traditional way, using an implicit form similar to what was done for the positive-P distribution in Sec. 3.1. The Q representation is the s → −1 limiting case of the family of s-ordered representations W s (α) studied by Cahill and Glauber [19,20] (the Glauber-Sudarshan P and Wigner correspond to s = 1 and s = 0, respectively). All these distributions can be written using coherent displacement operators in an implicit form similar to (2): with the base operator [20] T One has Tr Λ s = 1 and the operator identities which can be verified by equating the LHS and RHS when T (0, −s) is expanded in number states. Therefore, in the limit s → −1 corresponding to the Q representation with Λ −1 = Λ Q , We see then that multi-time correlations in which only the orderings a † j Λ Q and Λ Q a j appear could have the capacity to correspond to simple stochastic expressions. This implies anti-normal ordered moments since by the cyclic property of traces The implicit form (53) for ρ allows us to proceed in a similar fashion to Sec. 3.1. Consider then the antinormal and time-ordered correlation with times obeying (17). A potentially awkward issue is that the kernel Λ s is not very well bounded in the s → −1 limit, with β| Λ s (α)|β = 2 ε exp −2|α − β| 2 /ε where s = ε−1. To gauge whether this is a problem, we will work using infinitesimal ε > 0 in which case Using the form (20) on (57) with now B = a † and A = a, one obtains The part in the inner { } brackets will be either Both cases allow us to use derivative-free identities (56c) or (56b), to replace (61) with P s (α, τ 1 ) Λ Q where Now assuming an acceptable propagator P s (α, τ 2 |α, τ 1 ) exists for the evolution of the s-ordered distribution W s , use of (27) leads to in which the objects related to the first time interval have been fully converted to phase-space quantities. Proceeding in this fashion with increasing time for the remaining a and a † operators, one arrives at Variables α, α etc were relabelled to α(τ R ), α(τ R−1 ).
In the limit ε → 0 that we are considering, W s (α) → Q(α), which is real non-negative. To be consistent, the propagator P −1 must also be real nonnegative. When the master equation (6) is faithfully reproduced by a Fokker-Planck equation for Q(α), this will be the case. Then, the P s . . . P s W s factors can be interpreted similarly to (30) as the joint probability of samples α(τ 1 ), α(τ 2 ), . . . at successive times. With that, we arrive at the hoped for result that a time ordered (as (17)) and anti -normally ordered correlation is evaluated as Stochastic averaging is over products constructed using values from the evolution of a single sample, like in (34).

Conversion between samples of s-ordered distributions
The s-ordered distributions W s introduced in (53) are mutually related by [20] W s (α ) = 2 in the sense that W s and W s0 represent the same quantum density matrix ρ. When s 0 > s, this is a Gaussian convolution of the more normally-ordered distribution W s0 , and reflects the well known property that Q distributions (s = −1) are broader and more smoothed than Wigner (s = 0), which are in turn broader than Glauber-Sudarshan P distributions (s = 1) for the same state. Importantly for us here, this means that if we have samples α of a more normally ordered distribution, we can easily also obtain samples α of the less normally ordered distributions simply by adding Gaussian noise. The prescription is for each mode j, with each ζ j a complex random variable of variance 1: In particular, converting samples of a P distribution to samples of Q requires Converting samples of a Wigner distribution (if it is nonnegative to begin with) to samples of Q can be done with

Evaluation of anti-normal ordered moments starting from P and Wigner representations
Therefore, if one has samples α(t 0 ) of a P or Wigner representation up to a time t 0 (e.g. from a prior stochastic evolution), the prescription (67) can be used to convert them to samples α (t 0 ) of the Q representation at that time. The noises ζ j are generated just once at this time. Subsequent evolution according to the Q stochastic equations then leads to Q samples α (t) at later times t > t 0 . These can be directly used in expression (65) to evaluate anti-normal ordered multi-time observables. This is potentially a little less convenient that remaining in one distribution throughout because the conversion time t 0 has to be chosen before starting a simulation, and the anti-normal ordered operators cannot extend to times before t 0 . It does preserve the principal advantages of phase-space simulation, though: intuitive and computationally tractable expressions for observables, stochasticity, gentle scaling with system size. The evolution equations (9) are usually of similar form in all s-ordered representations, apart from simplification at special s values.
One cannot convert the other way with this procedure (e.g. from Wigner to P), so normally ordered multi-time correlations cannot be extracted this way from samples of Wigner or Q distributions.

Mixed-order moments
The above discussion suggests a way that some mixedorder multi-time correlations that contain both normal and anti-normal factors could be evaluated. Suppose early time operators (t t 0 ) are normally ordered (can be evaluated using the P distribution), while later time operators (t > t 0 ) are anti-normally ordered (can be evaluated using the Q representation). A switch according to (69) could be made at t 0 after evaluating any normally-ordered factors, and the Q distribution evolved and used at later times to obtain the remaining inner factors with t > t 0 . Let us check in detail whether this is feasible.
we have Now using (27), and taking care with variable labels: After the variable change caused by (27), now Notice also the −1 label on the latest propagator in (77), indicating that it is according to the Q representation equations in the time interval (t 2 , s 2 ]. Finally, applying (56b), one arrives at: This encodes the following sequence of operations: 1. Start with initial samples α at t 1 .
2. Propagate P representation equations to s 2 obtaining samples α. 3. Propagate P representation equations to t 2 obtaining samples α. 4. Add Gaussian noise as per (78), to get samples α . 5. Propagate Q representation equations to s 1 obtaining samples α . Along the way, samples are collected to use in the final stochastic average. The P s factors in (79) together with the Gaussian factor in the top line form the joint probability P (α , s 1 ; α , t 2 ; α, s 2 ; α, t 1 ). Therefore, the final estimator for the observable A is Primed variables α are samples of the Q distribution, while un-primed ones α are samples of the P distribution. This all matches intuitively with the evolution and the expectation that anti-normally ordered elements will use Q distribution samples and normally-ordered elements samples of the P distribution.

Most general ordering case
The widest generalisation of this procedure to other (time-ordered) cases is as follows: If there is an early time set of normally ordered operators, on either side of the correlation, it can be dealt with by sampling the P distribution according to the replacements Once this avenue becomes exhausted, one adds noise via (69) to convert α's to Q distribution samples α . If the remaining later time inner factors are anti-normally ordered, they can then be dealt with using the replacement The above covers both fully normal and fully antinormal ordered products as special cases. On the other hand, the case of an early time anti-normally ordered block and a later time normally ordered block containing several times is not amenable to this approach because one cannot stochastically convert Q samples to P samples. A large number of cases can also be reduced to a form amenable to this procedure by the use of the commutator [ a, a † ] = 1 for equal time factors. Also, simulations starting in a Wigner representation can be used for evaluation of outer symmetric ordered parts, and then a switch can be made to the Q representation via (70) to evaluate any remaining inner anti-normal ordered parts. Table 1 counts the number of distinct 1, 2, and 3 annihilation/creation operator products that can be evaluated using the various methodology described above. All one and two-factor combinations can be evaluated (though 4/12 of the latter require use of the Q representation). For third order correlations, almost all time-ordered cases can be evaluated (74 out of 80). The majority require use of the Q representation or a shift from a P to Q representation as described in Sec. 5.4 to work. There are 6 exceptions that cannot be evaluated:

Correlation function coverage
These all share the feature that the earliest factor already requires the Q representation, while the next operator in time requires the P representation to which one cannot return. There are also a number of correlations that are not time ordered, for which the earliest time t 1 is on the middle operator, and these are not possible to evaluate according to the schemes presented here.
The greatest interest in multi-time correlations usually concerns those involving two times (say t = 0 and t = τ > 0). Examples are counting correlations like a † (0) a † (τ ) a(τ ) a(0) , and pair correlations such as a † (0) a † (0) a(τ ) a(τ ) . The case count for these is summarised by  time ordered four-operator products of this form can be evaluated, including atypical combinations such as a(0) a † (τ ) a † (τ ) a(0) , but very many require the Q representation. (72 out of the 160 accessible ones require the use of the mid-simulation switching to the Q representation described in Sec. 5.4). The only correlations that are inaccessible are the non time ordered ones such as e.g. a † (τ ) a † (0) a(0) a(τ ) .
The change of distribution can introduce additional restrictions. In particular, Wigner and Q distributions have a tendency to involve higher order derivatives in the PDE (7) than P distributions, so that a standard diffusion process no longer captures the full quantum dynamics. An example are two-photon losses with operators R = a 2 for which P distributions produce only 2nd order derivatives but Q distributions 4th order terms, and Wigner distributions 3rd order ones. This precludes fully accurate calculation of correlation functions such as a † (0) a(t) a(t) a † (0) = α * (0)α (t)α (t)α * (0) stoch requiring Q evolution after the changeover, but not cases where the change to a Q distribution is only needed at the final time such as a † (0) a(t) a † (t) a(0) = α * (0)α (t)α * (t)α(0) stoch . Stochastic techniques for dealing with higher order PDE terms have been investigated [79,80,109,117] particularly in [48] for doubled phase space, though attempts to date have shown strong time and stability limitations.
, where A, B, C, D can be either of a or a † (same mode), and the time arguments can take up to two distinct times t = 0 and t = τ > 0.

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

Doubled phase space s-ordered representations
The s-ordered representations were first generalised to doubled phase space by de Oliveira [31]. Later studies [35,71,116] used a different normalisation that is closer to the original single-phase space formulation (52) of Cahill and Glauber [19]. It will also be used here. The explicit expansion of the density matrix was given in [35] as 3 : where the displacement-like operator d is obtained by the replacement α * → β in (52). To distinguish from single phase space, the superscript + is used 3 In the supplemental material therein.
where necessary. The usual properties Tr Λ + s = 1 and d −1 (α, β) = d(−α, −β) continue to apply. The operator identities are now The s → 1 limit of all the above gives the positive-P representation, s = 0 the doubled-Wigner representation of [71], and the limit s → −1 a doubled phase space analogue to the Q representation ("doubled-Q"). Eqs. (85) are equivalent to the correspondences found in [31]. An advantage of the doubled-Q and doubled Wigner representations relative to their single-phase space analogues is that all 2nd order derivative terms in the PDE can be made positive-definite and converted fully to stochastic equations via the same analytic kernel trick as for the positive-P representation. For example, spontaneous emission in a Schwinger boson system need not be amputated like in [73] for the Wigner representation. The kernel transform between different orderings in doubled phase-space is found to be = 2 which is easily verified with the help of from [19], the easy to show identities and Gaussian integrals. The distribution transform = 2 is readily found by equating two (83a) expansions of ρ which have different s, and applying (86) to the one with higher s (= s 0 ).

Non normally-ordered correlations in the positive-P representation
The ideas from Sec. 5 can be used to treat non normallyordered correlations in the positive-P representation, which -unlike the Glauber-Sudarshan P -is applicable for all quantum states and systems. With the tools for the doubled s-ordered phase space described in Sec. 6.1, derivation of the expressions for mixed-time expectation values follows the same path as in Secs. 5.1 and 5.3, but some care needs to be taken to incorporate the doubled phase-space. We obtain that: (1) To shift samples from a more to a less normally ordered doubled-phase space representation, one adds the following noise: Notably -the same noise for α and β * , which was not obvious a priori. The prefactor is 1 for positive-P to "doubled"-Q, and 1/ √ 2 for positive-P to doubled-Wigner.
(2) Anti-normal ordered mixed-time correlations are evaluated as where α and β are doubled-Q representation samples, possibly created via (90) (with s 0 − s = 2) from initial positive-P samples, and later evolved via the appropriate doubled-Q representation evolution equations.
(3) For mixed ordering, the procedure in Sec. 5.4 follows with the same structure, except that the correspondences are in the positive-P and in the doubled-Q. For example, the expression for the correlation A (80) using samples starting in the positive-P is (4) The correlation tallies of Sec. 5.5 apply without change to the doubled phase space representations.
with standard annihilation operators a j for modes j using units = m = 1. This describes two sites with (real) tunnelling J, local on-site interaction constant U , detuning ∆, and a coherent drive F (real) only on the first site. There is also a decay process with rate γ that is treated by describing the evolution of the system with the master equation To test the performance a bit more beyond the standard model, we have also addded a thermal bath with mean occupation N . Such a system can be realised using e.g. micropillars [69,123] or transmon qubits [72,122]. The nontrivial feature here is a two-boson destructive interference effect between photons injected by the drive, and other photons that have been previously injected, tunnelled to site 2 and then back, returning with a relative phase of π [10]. The result of this is that in the steady state only single photons can be present at the driven site 1, providing possibly an avenue to create single-photon sources [96]. The lack of double occupation is evidenced in a single-time two body correlation function g 11 = a †2 1 a 2 1 / a † 1 a 1 2 , which is very close to zero. However, for practical application, it is of particular interest to find out how large a time mismatch between measured photons can be accommodated without significantly increasing g (2) 11 from this low level. If it is too short, then the single-photon source will have a too short active time for practical applications. Therefore the correlation function of particular interest is in the stationary state, with delay time τ . The mean occupation is n 1 = a † 1 a 1 . The Ito stochastic evolution equations in the positive-P (s = 1) and doubled-Q (s = −1) representations are These were derived in [45] for the positive-P. Here the matrix elements of the one-particle Hamiltonian are H sp jj = −∆, H sp 12 = H 21 = −J, the drives are F 1 = F , F 2 = 0, while ξ j and ξ j are delta-correlated independent white real noises with variance and η j are similarly correlated independent complex noises: Appendix A.3 gives some detail on generation of the noise. In some cases, calculations with these stochastic equations can already be faster than brute force calculations directly with the density matrix ρ in a suitably truncated number state basis. For large systems, they are the only tractable way to access full quantum mechanics. Appendix A gives details of an integration algorithm that is robust to the multiplicative noise appearing in (98), and was used for the simulations reported here and elsewhere [35, 36, 41-45, 82, 107, 120, 134].

Two-time photon-photon correlations
Consider first the strong antibunching parameters studied in [10] and [45]: U = 0.0856, J = 3, ∆ = −0.275, γ = 1, F = 0.01. We will study the correlations primarily in the stationary state. To this end, a positive-P simulation is initialised in the vacuum α = β = 0, and evolved up to t = 30. The stationary state is attained after t 15, and we calculate multi time correlation functions from times t 0 = 20 to t 0 + τ using the available samples. (S = 2 16 in all cases). Uncertainty in predictions is calculated using sub-ensemble averaging, as explained in Appendix A.4. These errors are shown as triple lines in all the plots. The two-time photon-photon correlations between one photon measured at site 1 at time t 0 and the other at site j after a delay time of τ are where n j (t) = a † (t) a(t) = Re α j (t)β j (t) stoch is the mean occupation of site j. In the positive-P representation the correlation (101) is calculated via We can take the real parts above, because the imaginary parts must converge to zero in the S → ∞ limit. This behaviour is verified in the simulations.  Fig. 1 (top) shows the multi-time local photon correlation (101) at site 1 (i.e. j = 1) for the clean case (N = 0), and verifies that the result perfectly matches the exact brute force solution. The desired anticorrelation dip is seen around τ = 0, along with characteristic oscillations out to delay times of about τ = 5. The bottom panel shows how the anti-correlation dip degrades when the system is linked to a particle reservoir (growing N ). Notably the dip timescale does not change appreciably as the minimum correlation rises. A 10% remnant correlation which might still be acceptable for applications is found when N = 10 −8 = 0.026n 1 . This sets a limit on how much background photon reservoir is acceptable.

Breakdown of the Glauber-Sudarshan P
This system also shows the breakdown of the singlephase-space P representation very clearly. The evolution equation is (98a) with s = 1, but it is a correct representation of the quantum FPE (8) only when the diffusion matrix D µν in the FPE has all nonnegative eigenvalues [62]. Here, the diffusion matrix for real, imaginary parts of α j = α jr +iα ji has elements D jr,jr = U 2 Im(α 2 j )+ γN 2 , . with eigenvalues λ j,± = γN 2 ± U 2 |α j | 2 . These only become non-negative once γN > U |α j | 2 , i.e. γN U n j . The question then is: for what parameters does the evolution remain well described? When N = 0 The antibunched mode has mean occupation n 1 = 3.87 × 10 −7 , naively suggesting N ∼ 5×10 −8 = 1.5U n 1 /γ to already be a value for which the description is good. However, we can see in Fig. 2 that it does not give correct results at all. This is because the problem lies in the nr. 2 mode with n 2 = 1.07×10 −5 . Taking a far larger N = 2×10 −5 (in which case γN ∼ 10U n 1 , 7U n 2 ) gives almost correct results with the Glauber-Sudarshan P (though still not fully), but of course antibunching in mode 1 is long gone for such a relatively high thermal noise level. This is an indication that skimping on full quantum effects by trying approximate semiclassical methods is not a good strategy for this kind of system.

Differently ordered correlations
To test how the prescriptions developed in Secs. 4.1 and 5 work, we use a different driving, F = 3, which generates larger occupations (in the stationary state n 1 ≈ 0.043 and n 2 ≈ 0.98) and as a result more interesting anti-normal and mixed-order correlations. The doubled-Q simulations are also too noisy to get useful predictions for the parameters in Fig. 1 because Q distributions have a width of O(1) even in vacuum. This is a certain limitation to the Q distribution approach. The difference in sampling accuracy can be nicely seen in Fig. 3 which shows the mean and uncertainty of n 1 during the simulation used to generate Figs. 4-5. The samples are switched from positive-P to doubled-Q at t = t 0 = 20.  Arguably the primary practical interest lies in cases with four factors and two times -corresponding roughly to either detection of single particles at two times, or creation and destruction of pairs. Fig. 4 shows the simulation of some anti-normally ordered correlations of this kind (not normalised), evaluated using the doubled-Q representation: The latter G b is an anomalous correlation that has both real and imaginary components. Similarly, Fig. 5 shows simulations of a number of mixed-order correlations which absolutely require a swapping from positive-P to doubled-Q representation using the procedure of Sec. 5.4: In all cases, primed variables are evaluated in the doubled-Q representation, un-primed in the positive-P. Notably, in both figures the stochastic simulations perfectly and very accurately agree with brute force calculations using the density matrix. This is despite the regime being one which is very poorly treated by approximate semiclassical methods (occupations are O(1) or smaller). This is strong evidence that the intuitive expressions and approach laid out in the previous sections is appropriate.

Correlation dynamics in the stationary state
As a larger-sized example, consider the same Hamiltonian and master equation as (95) and (96) (take N = 0), but with a longer chain of 32 sites, which is already beyond or at the limit of the capabilities of alternative methods such as brute force or corner space renormalisation [24,59]. The tunnelling term in (95) now becomes −J 31 j=1 a † j a j+1 + a † j+1 a j , and a driving of F = 3 remains only at the j = 1 site. While the stationary state does not show any time-dependent change in expectation values, we should expect to still see signatures of transport in its multi-time correlations as a function of distance and delay time if the method is good.   6 shows such a calculation in the 32 site system, plotting the normalised photon-photon correlation g 1,j (τ ) with spatial separation j − 1 and delay time τ , as defined in (101). The correlations are evaluated from t 0 = 20. Despite this being the stationary state, a very clear anti-correlation signal can be seen moving steadily with delay time. Its speed is approximately 2J, twice the tunnelling rate. Correlation waves often travel at twice the characteristic speed for single particles [25,97,114], so this is not unexpected. However, the speed is clearly not related to the superfluid speed of sound, here ∼ U n j , which is far lower.

Summary
The known framework for evaluation of multi-time observables from the work of Gardiner in the P representation [62] has been extended here to include firstly: the much more widely applicable positive-P representation (39); Secondly the Q (65), doubled-Q and other single and double phase space s-ordered representations; Thirdly: other orderings such as anti-normal and mixedordered observable products (Sec. 5, and especially (80) and the algorithm of Sec. 5.4). These results allow the evaluation of a very wide range of quantum multi-time observables in bosonic systems (classified in Tables. 1-2) in a way that contains the full quantum mechanics, and is scalable to large systems. Systems with localised interactions require computational effort proportional to M or M log M with the number of modes or sites, M .
While out-of-time correlations (OTOCs) are not directly accessible through the mechanism outlined here, time reversal schemes like those used in experiment [63] or theory [15,46,133] could be attempted by changing the signs of constants in the equation of motion. When combined with the techniques developed here, this would give access to a wide range of information about quantum chaos [15,60,98], many body localisation [56,70], and quantum phase transitions [124] in larger systems than were accessible to date [15,63].
Along the way, a number of additional results were obtained regarding conversion formulae between different orderings in the doubled-phase space representations ( (86) and (89)) and their stochastic samples (90). A clear case of the breakdown of the Glauber-Sudarshan P representation was seen in Fig. 2. Also, in Sec. 7 some results regarding the correlation functions (Figs. 4 and 5), susceptibility to background thermal density (Fig. 1) and signal speed (Fig. 6) in the unconventional photon blockade system were obtained.
Finally, Appendix A details a convenient algorithm for simulation of phase-space stochastic equations, an item that has been hard to find in the literature in the past.

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

A Numerical techniques for stochastic simulation of phase-space trajectories
This appendix details how the stochastic equations for the examples in Sec. 7 were integrated, and describes a robust and general algorithm that is particularly applicable to stochastic trajectories generated by phasespace descriptions of quantum systems. The bulk of the method is based on the semi-implicit midstep algorithm described in [51,140], adding some modified propagators as per [44] which allow for more efficient treatment of exponential and dominant deterministic processes. Some practical elements that are usually skipped over in the literature are also pointed out. Example textbooks for the broad topic of stochastic integration include [76,87].

A.1 Integration algorithm
Let us introduce the following notation. The Ito stochastic equations for a set of variables v with elements v j indexed by j will be written where X j are noise terms with zero mean X j stoch = 0, whereas A j are deterministic "drift terms" that contain no explicit noise contribution, apart from the statistical spread of the variables v j themselves. At the level of the Ito equations (105) we also require that the noise is not correlated in time i.e. proportional to Wiener increments. For example, in (98): For various reasons, the standard mathematical environments and numerical packages tend not to perform well when integrating stochastic equations 4 , especially ones that contain variable-dependent diffusion coefficients such as in (106) or the prototypical Kubo oscillator [90] with multiplicative noise. Phase space methods that treat the full quantum dynamics also typically contain such terms. Once one leaves the simplest Ito-Euler algorithm which requires very small step sizes, two of the recurring issues are that advanced integration algorithms usually assume sufficiently continuous derivatives (completely violated by white noise) or try to adapt the timestep in real time (which can introduce systematic errors when noise is involved). A third difficulty is the need for autocorrelation corrections (seen Sec. A.2) which become much more difficult to derive as the complexity of the algorithm grows.
The semi-implicit midpoint approach [47,126] has been shown to be robust against spurious numerical instability as carefully analysed in [51,140], but remains simple enough for any autocorrelation corrections to be readily calculated.
Suppose we are to advance time by ∆t, starting from variables v  (107) which just uses the values at the beginning of the step for evaluating the derivative, as per the definition of Ito stochastic calculus.
The semi-implicit midpoint algorithm can be viewed as a container procedure that calculates the final variables v j via a series of iterations [140] v (1) The numbered iterations always start from the initial point, but iteratively find a self-consistent estimate of the midpoint value of the derivatives. The last iteration (109) uses this midpoint value and midstep time to make the final advancement of the variables. Usually, a single midpoint iteration (i.e. n = 1) is sufficient. When using the basic Ito-Euler propagator (107) at time argument t + ∆t/2 in (109), the final step is accurate to O(∆t 2 ) in the deterministic parts A j , and the step's mean to order O(∆t) in the noise terms X j [51], provided the autocorrelation correction is incorporated into v step j , as explained in Sec. A.2. This is so-called "weak" stochastic convergence to O(∆t) in the terminology of [51] because the accuracy for single trajectories is O( √ ∆t). "Strong" stochastic convergence to O(∆t) for single trajectories requires a more involved and more inconvenient autocorrelation correction than that presented in Sec. A.2. Nevertheless, the "weak" level of accuracy in the noise term is actually quite high already. For example, even 4th or 9th order Runge-Kutta implementations made available are generally only accurate to O( √ ∆t) in this regard [33].
Moreover, in practice, one can often easily but substantially improve on the Euler step (107). For example, if ∆, γ, or U are large in (98), the evolution is primarily exponential. Another typical case is in continuum models when the kinetic energy or external potential energy contributions are dominant and trivial to integrate [44]. A pragmatic and flexible propagator choice is to separate the evolution into parts that will be treated exponentially (E), strictly linearly (L) and the rest (R): and solve the equation for the E and R parts (111) as if the coefficients were constant [44]. This gives the full propagator v step Then, the leading processes are often integrated exactly, with only higher order corrections needed to be worked on by the midpoint algorithm. This is particularly efficient when the stochastic terms are a perturbation on the dominant exponential deterministic evolution. It is usually optimal to place all non-exponential parts into the "R" part, but some special cases can require avoidance of the nonlinear solution [44]. Both A drift and X noise parts can enter the coefficients as one chooses.
In the case of the calculations in this paper using (98), and D L j = 0 were used for the propagator (112) and combined with a single iteration (n = 1) of the midpoint algorithm (108). Note the presence of the autocorrelation corrections as explained in Sec. A.2. Importantly -the form of the corrections (114) above arises when the noises ξ, η etc. appearing in the X j are calculated just once at the beginning of the entire procedure for making the full ∆t step. Other noise input can change the expressions (114).
The actual timestep ∆t that must be used is constrained by the need for the coefficients D j to change little over a single timestep ∆t. A good practical rule of thumb is to require for each term D (ν) in D = ν D (ν) . One caveat is that one should take care not to estimate D( v (0) , t) on the run using noise or even variable values from individual trajectories. Doing so can, and often does, introduce autocorrelations between noise values and time step lengths, leading to systematic errors. For treating noise terms, it is useful to use the RMS estimate X j ≈ |X j | 2 ∝ √ ∆t based on expected mean variable values v j for timestep estimates. Finally, for completeness and general applicability beyond the examples of Sec. 7, it can be advantageous to treat widely differing processes using a split-step approach [140]. In particular, the treatment of kinetic energy and tunnelling is much bettered for many-mode systems. A split step allows to take advantage of the fact that kinetic energy and tunnelling are diagonal in kspace and can be done exactly there. k-space is reached using a discrete Fast Fourier Transform that has computational cost of only M log M for an M-site system [61]. The split-step is applied in three parts: (i) First evolve terms diagonal in k-space for ∆t/2; (ii) then change variables to x-space and evolve the remaining terms there for ∆t; (iii) finally return to k-space for another k-space evolution of ∆t/2.
For cases like (98) in which the main mean-field evolution has a Gross-Pitaevskii form, the split step method is symplectic (conserves energy) and has been shown to have an accuracy of O(∆t) 2 for the long time solution [77] provided that the latest available copy of the field is input as v (0) after each Fourier transform [57]. Within each k-or x-space split step, the midpoint method (108) can employed if the coefficients D of the substep are variable-dependent, or otherwise if the coefficients are constant just a plain propagator substep (112) is done. Thus, one has the split-step algorithm as a container for midpoint iterations, which themselves are containers for the base propagators. Long-range interactions (eg. dipole-dipole) can also often be treated efficiently through Fourier transforms [144].

A.2 Timestep autocorrelation correction
Realisations of Ito stochastic equations (105) must at the least produce the required average and noise variance to lowest order. That is, writing In orders of ∆t, the drift is A j ∆t ∼ O(∆t) and the noise X j ∆t ∼ O( √ ∆t). The time-dependence of coefficients A and X does not enter at this order of ∆t. A simple Ito- ) ∆t meets both conditions (117).
However, the midstep algorithm with an Ito-Euler propagator (107) can introduce an extra correlation of order ∆t by virtue of using the same noise for both the halfstep of ∆t/2 and then for the final full step. The first halfstep advance using (107) produces Note that the X j themselves are O( 1 √ ∆t ), which is important below. The next step in (108) involves coefficients A j ( v (1) ) and X j ( v (1) ). Taylor expanding around v (0) , and analogously for A j . Using (118) and discarding terms of subleading order O( √ ∆t) in the new coefficients one finds Notice the appearance of a term in X j ( v (1) ) of the same order as the original coefficients X j ( v (0) ). Assuming first just one n = 1 iteration, the final step (109) becomes It turns out that iterations of n > 1 do not affect the expression (121) at O(∆t). The new term in (121) breaks the equivalence of the algorithm (117) because it is of the same order as the drift A j ∆t.
To restore equivalence, one can add a correction C j to the drift used in the algorithm as per in either the "strongly convergent" form (123) or the "weakly convergent" averaged form C j = C strong j stoch . For the case of Ito-Euler propagator in the midstep algorithm, The strong forms (123) have different values for each stochastic realization [51] and obtain O(∆t) timestep accuracy for the noise terms, while the weak forms use a pre-calculated mean (124). The weak variant is more commonly used in practice and sufficient to obtain overall stable and accurate integration. It is what is applied in (114) and the calculations of Sec. 7. The weak autocorrelation correction C j has often been dubbed a "Stratonovich correction" because the form (124) is identical to the correction used to move between Ito and Stratonovich stochastic calculus. However, this is a misnomer because the match is merely a coincidence that occurs when an Ito-Euler propagator is used with the midpoint algorithm. For the case of the partly exponential propagator (112), the same procedure as above leads to the following, different, form of the autocorrelation corrections. When n 1 in the midstep algorithm: or when one uses no midstepping (n = 0): As it happens, in the case of the equations (98), the corrections are the same whether using an Ito-Euler propagator ( (114)) or the partially exponential one (125)- (126), but that is not always the case.

A.3 Noise implementation
Since the underlying stochastic equations are defined in the infinitesimal limit ∆t → 0, in principle any implementation of the noise that has zero mean and satisfies the variance conditions (99) and (100) in the ∆t → 0 limit is correct. By the central limit theorem, after several timesteps the effective distribution of the sum of the noises will always converge to a Gaussian one. In practice, however, it is usually desirable to use explicitly Gaussian distributed noises in the implementation. Doing so already accurately depicts the limiting distribution for each ∆t step and avoids a potential need to reduce timestep further to capture the right noise distribution. For real noises ξ, one generates fresh Gaussian random variables of variance 1/∆t for each j at the beginning of each timestep. Generating new ones for each time step ensures the δ(t − t ) limit as ∆t → 0. The Box-Muller algorithm [17] is a simple way to obtain two such independent Gaussian variables from two uniformly distributed noises, whereas a variety of more efficient though less transparent algorithms are also known [12,18]. In our case Box-Muller was used, while the underlying uniformly distributed noise was generated using the SFMT fast Mersenne twister method [121] 5 . Complex noise η is simply constructed as η = (ξ r + iξ i )/ √ 2 using two real Gaussian noises ξ i,r .
When using efficient random number generation such as [121] it turns out that the creation of Gaussian noises from uniform ones is the most computationally costly step. The Brent algorithm [18] alleviates this appreciably, while a simple alternative route applicable at least for stochastic simulations -due to their central limit properties -is to produce binomially distributed noise instead. One adds n B uniform noises r n on [0, 1] as per ξ ≈ 12 n B ∆t n B n=1 (r n − 1 2 ). In practice n B = 4 or n B = 3 is already sufficient.

A.4 Error estimation
Statistical uncertainty in quantities Q calculated from the simulations can be robustly estimated via subensemble averaging [140] even when the distribution of samples is unknown. The full ensemble of, say S, realisations is divided into a smaller number u of subensembles. Auxiliary subensemble means Q (i=1,...,u) are calculated for each sub-ensemble individually. Then by the central limit theorem, the 1σ uncertainty in the full ensemble prediction is In practice u √ S is useful. In the simulations reported in this paper, u = 32 was used.
A separate issue is testing the time discretization accuracy in stochastic equations. In a simple implementation where noise is generated sequentially, changing timestep changes also the noise history, making comparison of single trajectories with different ∆t useless for determination of accuracy. Yet, comparing the ensemble means for different timesteps can be onerous in large systems. A solution is to compare two runs of a single realisation with identical underlying noise sources [51]: One with time step ∆t and noises ξ (i) j generated sequentially for time steps numbered i = 1, . . . . The other with integration timestep 2∆t but noises generated by adding pairs of noises from the first simulation: ξ . This allows separation of the time discretization error from the stochastic randomness.