Fully quantum scalable description of driven dissipative lattice models

Methods for modeling large driven dissipative quantum systems are becoming increasingly urgent due to recent experimental progress in a number of photonic platforms. We demonstrate the positive-P method to be ideal for this purpose across a wide range of parameters, focusing on the archetypal driven dissipative Bose-Hubbard model. Notably, these parameters include intermediate regimes where interactions and dissipation are comparable, and especially cases with low occupations for which common semiclassical approximations can break down. The presence of dissipation can alleviate instabilities in the method that are known to occur for closed systems, allowing the simulation of dynamics up to and including the steady state. Throughout the parameter space of the model, we determine the magnitude of dissipation that is sufficient to make the method useful and stable, finding its region of applicability to be complementary to that of truncated Wigner. We then demonstrate its use in a number of examples with nontrivial quantum correlations, including a demonstration of solving the urgent open problem of large and highly non-uniform systems with even tens of thousands of sites.

Unbiased quantum methods, including corner-space renormalization [32] and quantum trajectories [33,34] can successfully treat small systems, but suffer from the usual runaway complexity problems once larger numbers of modes or sites are present. This issue is exacerbated even further for open systems since density matrices are needed, where the number of variables scales as (e M ) 2 with the configuration size M rather than "only" e M for pure states. Matrix product states and related techniques [35,36] offer one way around this for closed systems, but their extension to include drive and dissipation is difficult [37].
In contrast, techniques known as phase-space methods, in which quantum expectation values are calculated from averages over stochastic trajectories in phase-space, are readily scalable to quantum problems with large numbers of sites or modes, and are naturally adapted to open systems due to already being based on a density matrix formalism. Their performance does not depend much on * deuar@ifpan.edu.pl dimensionality. Indeed, the use of the approximate truncated Wigner method has become common for studying semiclassical phenomena in ultra-cold atoms and microcavity polaritons [38][39][40][41][42][43][44][45][46][47][48][49][50][51][52]. However, as lattice experiments increasingly aim to delve further into the quantum regime in these media, other techniques are needed to study quantum effects beyond the reach of the truncated Wigner approximation.
An alternative phase-space method, the positive-P approach [53] allows for the full quantum mechanics of systems with up to two-body interactions to be simulated in an unbiased way without approximations. It has already found significant application in quantum optics [54], and in ultracold atoms [55], where it has been successfully applied to cases with hundreds or even millions of sites [56,57]. For closed systems, the trade-off has always been that while results for short evolution times are accessible, a nonlinear amplification of the trajectory spread eventually appears at sufficiently long times to obscure predictions below a rising noise floor [58,59]. However, it is already known that dissipation is beneficial to the stability of the method, and simulations can stabilize fully if it is sufficiently large [58].
It is with this in mind, and with the increasing relevance of the physics of open quantum systems to a number of experimental platforms, that we propose positive-P as an ideal method for simulating such systems in intermediate regimes, relevant to current experiments, where driving, dissipation and quantum correlations are all relevant effects. To demonstrate this, we focus on the archetypal driven dissipative Bose-Hubbard model, which is directly applicable to a number of the different experimental realizations [1,60]. We firstly endeavor to thoroughly characterize the regimes of applicability of positive-P in the parameter space of the driven dissipative Bose-Hubbard model, before also demonstrat-ing a number of specific examples of nontrivial effects accessible to the method, some of which may be difficult to solve accurately by other means due to the very large or highly non-uniform systems considered. The success of the positive-P method demonstrated here for the driven dissipative Bose-Hubbard model also implies that the stabilizing effect of dissipation on the method should likely allow it to be useful for simulating a number of related models of open quantum systems in future. We also demonstrate that the regions of applicability of positive-P and truncated Wigner happen to be complementary to each other, with the truncated Wigner approximation being fairly accurate for large occupations (i.e. strong drive) and positive-P being stable for strong dissipation. Between them they provide a viable phase-space method for almost all regimes where external drive and/or dissipation are significant effects.
The paper is organized as follows: In Sec. II we describe the driven dissipative Bose-Hubbard model, and then in Sec. III present its mapping to the positive-P representation (5). Sec. IV studies the single-site case and determines the level of damping (IV A) needed for successful simulation, while benchmarking against known exact solutions. We then investigate use cases in multimode models (Sec. V), including Lieb lattices with dark sites and large non-uniformly driven 2d square lattices, demonstrating scalability to huge systems (Fig. 11). Extension to nonzero temperature is given in Sec. VI before concluding in Sec. VII. An illustration indicating the key messages of this paper is presented in Fig. 1.

II. MODEL
The Bose-Hubbard model is the standard go-to description for bosonic driven-dissipative lattice systems. In dimensionless units the Hamiltonian can be written in the tight-binding form: Here, the local part of the Hamiltonian at site j is where a j is the bosonic annihilation operator at site j, U j 0 the local two-body interaction, F j the strength of coherent driving (can be complex), and −∆ j is the local energy bias. For example, for polaritons in micropillars [18,19,[21][22][23][24][25][26][27] with pumping frequency ω p and natural mode frequency ω j , the ∆ j = ω p − ω j plays the role of an effective chemical potential [34]. Returning to (1), J ij = J * ji is the tunneling amplitude for a transfer i → j between connected sites. For definiteness, in this notation, each connection occurs only once in the sum, so that e.g. a system consisting of just two connected sites has the tunneling terms −J 12 a † 2 a 1 − J * 12 a † 1 a 2 . Complicated connections and lattices can also be trivially incorporated into the model via the general form in (1). While in this work we consider examples with nearest neighbour connections in one or two dimensions, there is no reason in principle that these methods should be any less effective for higher dimensions, all-to-all connections, or longer range tunneling that could be represented by arbitrary J ij . The local single-particle dissipation rate is γ j . The system is then described via the density matrix ρ and evolves according to the master equation (3) This assumes dissipation into empty modes. The case of non-empty reservoir modes is described in Sec. VI.
For a single mode (site) with parameters F , ∆, U , γ, the observables of most interest are the mode occupation N = a † a , mean amplitude a , and normalised twobody correlation g 2 = a † a † a a / a † a 2 . Bunching is indicated by g 2 > 1 and antibunching by g 2 < 1. Strongly antibunched modes can in principle be good quantum sources of single photons. The steady state solution of the single mode has been calculated analytically by Drummond and Walls [61]. Several regimes can be identified based on which process is dominant on the observables N and g 2 : 1. A strongly driven regime when |F | U and |F | γ with coherent high occupation in the stationary state N ≈ (|F |/U ) 2/3 , g 2 ∼ 1.
Coupling different sites will inevitably mix the different regimes, leading to novel quantum phenomena [34,60,62,63], including more exotic physics with hysteresis and large collective fluctuations [34,64,65]. Needless to say, no exact solution of the steady state of the many-site problem is currently available, even in one dimension. Models with space-dependent parameters are certainly possible and often demonstrated experimentally (e.g. micropillars allow for the fabrication of systems with parameters that are very flexible from site to site [66]), but have been much less studied and simulated. Time dependence is also possible -most readily for F (t). Positive-P works especially well for low occupations and/or strong dissipation, while the truncated Wigner approximation is accurate for large occupations (see Fig. 6).
(f) In positive-P, normally ordered quantum observables are calculated by averaging the corresponding stochastic phase space variables over realizations. This correspondence is exact in the limit of large numbers of realizations.

III. POSITIVE-P REPRESENTATION
The application of the positive-P representation [53] to the model (1)-(3) generally follows the standard procedure applied to the related ultracold Bose gas systems without drive and dissipation [59,67]. One expresses the density matrix of an M mode/site system as in terms of local coherent state kernels Λ j at each site j, with Tr[ Λ j ] = 1. The |α j j and | α j j are local coherent states |α j j = exp α j a † j |vac . The bold notation α α α indicates a vector of all α j values. As a result of the properties of Λ, the distribution P can be made positive real for any density matrix, hence it is a true probability distribution of the configurations v = {α α α, α α α * } [53]. For this to be possible, however, the α j "bra" duals to the "ket" amplitudes α j must be independent, leading to an off-diagonal kernel operator Λ j . There is a full equivalence between the density matrix ρ and the distribution P ( v). Moreover, a set of S samples of the configuration v, distributed according to P , is also equivalent to the full density matrix in the limit S → ∞. Therefore, a set of such samples can in principle be used to approximate full quantum mechanics with increasing and unbiased precision as S grows.
We can then use the properties of the projector Λ to convert the master equation (3) into a Fokker-Planck equation for the evolution of the distribution P (see Appendix A for details), which in turn leads to stochastic differential equations for trajectories of the phase space variables v. The resulting (Itô) stochastic equations for the samples of v are where the final sum is over all sites k connected to j. The real random variables ξ j (t) and ξ j (t) are independent white noises of mean zero obeying and ξ j (t) ξ k (t ) s = 0, where the notation · s denotes stochastic averaging over the available samples in the limit S → ∞. Moreover, ξ j dt andξ j dt are standard Wiener increments, which are implemented by Gaussian random variables of variance 1/∆t at each time step of length ∆t. The equations (5) are the ones to be solved numerically and our subsequent analysis in this paper is based upon them. They contain the full quantum mechanics of the system, provided that the noise amplification catastrophe alluded to above does not occur (The useful simulation time t sim beforehand is estimated in Appendix B).

IV. SINGLE-MODE PERFORMANCE
Let us start with the baseline single mode case, because it is very revealing regarding the capabilities of the method, and allows us to easily compare to the exact solution, as was given by Drummond and Walls [61]. It also turns out to be an excellent guide for assessing which many-site systems can be simulated, and lets us understand more involved multi-mode systems that will follow. We omit the site indices j in this Section. The observables of most interest have the following stochastic estimators in the positive-P calculations: For all results we present in this work, we begin simulations in vacuum (α = α = 0) and evolve until the steady state is reached (or until excessive noise amplification makes further simulation pointless). Appendix B gives a perspective on other initial states.

A. Regimes of usefulness
A basic starting question is whether the stationary state can be reached. For many-site systems, a rough minimum requirement is that single site simulations can do so -under all the local conditions found in the large system. Hence, the fundamental importance of determining the conditions under which a single site system can reach the stationary state. We have carried out positive-P calculations across the whole spectrum of parameters for the single site system, and assessed them according to whether a stationary state with useful signal-to-noise ratio is reached. That is, whether for practical numbers of realisations, the values of the observables we consider in the steady state are not masked due to the selfamplification of the noise. U was chosen as an arbitrary energy scale. Fig. 2 presents the results of this benchmarking, over many orders of magnitude of the parameters F , U , and γ, when ∆ = 0. This is one of the main results of the paper. Green square: numerics reaches the stationary state and remains stable; Yellow square: remains stable, but poor signal-to-noise ratio makes accurate determinations intractable (especially for g2); Blue square: numerics reaches the stationary state but does not remain stable later; Open circle: numerics becomes unstable before reaching the stationary state. The broad grey lines indicate crossovers between physical regimes listed in Sec. II; the red dashed line shows the empirical estimate of the usability region (8).
The stable region in which numerical integrations reach the steady state and remain well behaved is shown in green, and is attained for all parameters F, U, ∆ when the damping γ becomes sufficiently large. Examples of such calculations are shown in Fig. 3(a, b, c). Dynamics that do not reach the steady state before the noise instability occurs, such as Fig. 3(d), are shown as a small open circle. The blue squares are on the edge of stability, such that a stationary state is reached, but noise instability similar to that shown in Fig. 3(d) sets in some time after. The yellow square cases are stable, but mode occupation is too low compared to the vacuum noise, and useful information cannot be extracted. We find that for nonzero ∆, the regime of stability is qualitatively almost identical to that in Fig. 2, particularly on a log-log scale (see Appendix C for details). Dependence on ∆ is investigated further in Sec. IV C.
An empirical rule that largely captures the regime of usability, based on the data shown in Fig. 2, is: The uncertainty is about ±0.01 on the exponent, and 10% on the prefactor. For very low driving, a more appropriate rule is In the usable regime, numerical effort scales linearly with the number of sites, and quadratically with the precision (according to the central limit theorem, since all samples have independent noise input). Much lower damping may become accessible through the use of stochastic gauges, particularly in the high occupation regime where they were shown to be effective for this Hamiltonian [68].

B. Typical behavior
Here, we now look in more detail at specific examples presented in Fig. 3. The case in Fig. 3(a) with F = 1000U , γ = 31.6U and ∆ = 0 is representative of the strongly driven regime, with a few oscillations before settling down to a steady state with high occupation and almost perfect coherence (g 2 ≈ 1). The crossover regime that mixes all three regimes mentioned in Sec. II, and is often studied [32,34,62,63,65,69], is shown in Fig. 3(b). There F = U , γ = 3.16U , ∆ = 0, and occupation is O(1). This case is notable in that we can obtain large antibunching, indicating strong quantum effects, while remaining stable and despite rather strong dissipation.
Getting into lower occupations and stronger antibunching, Fig. 3(c) shows the case of F = 0.1U , γ = 2U , ∆ = 0. Notice that the statistical error in g 2 is becoming more pronounced, despite averaging over 10 6 trajectories. This is still a well behaved simulation, however, without significant noise amplification. The fairly low signal to noise ratio is a consequence of low occupation. When damping is insufficient to stabilize the long time behavior, a case like Fig. 3(d) occurs, here with F = 0.1U , γ = U , ∆ = 0. The exact stationary value is approached, but the evolution does not convincingly stabilize before noise amplification appears (first spiking near U t ≈ 5) and leads to an instability (U t ≈ 6.2). As is common for higher order moments, the g 2 estimation becomes too noisy to be useful some time before.

C. Nonzero detuning
While detuning ∆ does not appreciably influence the regime of stability, the physics is significantly affected. Fig. 4(a,c) shows the variation of occupation and bunching when F = U , γ = 3.16U , close to the mixing region of all three regimes. This is quite a strongly damped case compared to many theoretical studies, but still shows strong bunching and antibunching. The form of this variation of g 2 with ∆ is qualitatively consistent with the phase diagram calculated in the weakly damped γ = 0.05U regime [34], just with a reduced degree of bunching/antibunching owing to the stronger dissipation relative to U . Comparison with exact results shows that all detunings can be reliably and stably simulated, even close to the ∆ = 0 limit of stability (8).
The fact that those simulations remain stable for ∆ = 0 can be attributed to the decreasing occupation. To understand this, note first that for the undamped system, single particle energy shifts of the kind represented by ∆ were shown not to affect the stable simulation time given by (B1), for set values of U and mean particle number N [59]. In the damped system, a similar indirect-only dependence is expected, but N does depend on ∆. Since |∆| > 0 generally reduces the particle number (Fig. 4), the estimate (B1) indicates increased t sim , so one expects increased stability and smaller γ values than in Fig. 2 to become accessible. This is borne out in Fig. 5, which shows a study of this stability dependence at F = U , the typical case of interest. Near the edge of the stable region, only very rare trajectories exhibit instability, such that small ensembles are usually still well behaved (cyan color in Fig. 5). Such a trade-off between better precision in larger ensembles, but encountering instability if one generates too many trajectories, is typical for the positive-P method in borderline stable/unstable regimes.
The unstable region is more asymmetric around ∆ = 0 than the density in Fig. 4. Bunching correlates with increased number fluctuations, such that maximum excursions of occupation are larger for ∆ > 0 than for ∆ < 0, making instability persist at ∆ > 0 for larger damping.  Fig. 2, with the addition of cyan cases when the instability after the steady state is reached was seen only for very large ensembles (S = 10 6 ) but not seen in S = 10 5 ensembles. Fig. 4 compares the positive-P and exact results to those of a leading competitor for scalable quantum simulations -the truncated Wigner approach (TW). This method's equations are described in Appendix D. Unlike the positive-P method, truncated Wigner involves an approximation, namely that the exact evolution equation for the Wigner distribution (the equivalent of (A5) for that representation) contains third order derivative terms, which must be neglected in order to obtain a stochastic differential equation from the resulting Fokker-Planck equation. Physically, this means that some quantum correlations are not included in the description, which becomes an issue when looking at problems with a higher degree of entanglement and low mode occupations. In other words, the more semiclassical the problem is, the better it is described with the truncated Wigner approach, which fails for very quantum cases. It is therefore useful to show that the positive-P may be applicable in situations where the truncated Wigner approximation fails to give accurate results, as well as compare their properties under conditions where either method would be viable.

D. Comparison to truncated Wigner
One can see that while truncated Wigner gives a qualitatively good description of the occupation (though with some deviations), two-body correlations g 2 are on the whole completely inaccurate. Unphysical predictions of g 2 < 0 with the truncated Wigner method are also seen. These are typical known problems with the truncated Wigner approach when occupations are low. The method gives much more accurate results for high occupations, such as in the strongly driven regime. Table I gives examples of this behavior for some other values of the parameters.
For both methods, the statistical uncertainty on the steady state result is obtained by partitioning the S ≈ 10 6 trajectories into (roughly) s ∼ 100 subensembles, each containing S/s trajectories. For each subensemble i, we extract the steady state value O i of a given observableÔ. We then consider these s values as independent measurements, so that our best estimate is The positive-P and truncated Wigner also differ with regard to the signal-to-noise ratio (SNR). At low occupations, the SNR is far superior in positive-P, while at high occupations it is comparable. This can be seen very clearly in Table. I. On the other hand, despite systematic errors and SNR issues, the truncated Wigner never suffers from the noise catastrophe of Fig. 3(d).  (2) .01015 (15) .00010 (8) 99.3309(5) g2: exact [61] . 86243 .90930 .799984 .9966697 Positive-P .8628 (6) .9093 (5) .801 (5) .996675(9) tr. Wigner .779(2) -12(2)  A systematic comparison of the applicability of the two methods is made in Fig. 6, using the ∆ = 0 case. We assess the accuracy of the truncated Wigner by calculating both the systematic and statistical relative errors, defined as |O − O ex |/|O ex | and δ stat O/|O|, respectively, for each of the four observables shown in Fig. 3: N , g 2 , | a | 2 /N and phase arg a . Here O ex is the corresponding value of the observable obtained from the exact solution of Ref. [61]. We then define ∆ T W as the maximum relative error out of the entire set (see Appendix D for  Fig. 6 correspond to values of ∆ T W = 0.01, 0.03, 0.1, 0.3 from top to bottom. We take ∆ T W = 0.03 as the nominal limit of sensible applicability of the truncated Wigner. This curve determines the upper filled region in Fig. 6, corresponding to the model parameters that can be accurately simulated by the method.

details). The blue contours in
Explicit conditions for the TW accuracy region can be obtained by fitting the ∆ T W = 0.03 curve, in the asymptotic regimes of weak and of strong dissipation, to a straight line. The obtained results are shown as red dashed lines in Fig. 6, from which we find The lower filled region in Fig. 6 refers instead to the regime of sensible applicability of the positive-P, where the empirical limits (IV A) are used. The bottom line of this comparison is that the regimes of applicability of the positive-P and truncated Wigner methods are mostly complementary. TW is sufficient for small damping, high driving (alternatively -large N ), where the system behaves largely semi-classically and quantum correlations are small, while positive-P should be a method of choice for low driving, appreciable damping (low and moderate N ) where quantum correlations are significant. Both are good in the high damping, high occupation regime. Together, these two phase space approaches cover the vast majority of the parameter space. What is left is the low occupation (the limit in (10a) with N ≈ (F/U ) 2/3 gives N ≈ 2.5), low damping regime, which fortunately suits tensor network methods best.
1,1 (τ ) for the driven site in the steady state ( the same example as in Fig. 7). Positive-P (red) is compared with numerical exact solution of the master equation (dashed black).
For the first many-mode example, we consider a situation where nontrivial behavior can occur in a system of only two sites. Strong two-particle interference effects leading to g 2 → 0 pose no problem to simulate. A calculation of the so-called unconventional photon blockade [71] (using parameters from [70]) proceeds easily, as shown in Fig. 7. The steady state value obtained with 10 6 realizations is −0.001 ± 0.004. This system consists of two sites "1" and "2", in which only site 1 is driven. Destructive two-photon interference leads to the effect seen in Fig. 7, which demonstrates that two photons never occur together in this site in the steady state, giving an excellent single photon source.
Using this example, we can also show how to calculate multi-time correlations with positive-P. Any multitime correlation function that is normally-and timeordered can be calculated in the positive-P representation in a simple way, by averaging the corresponding product of phase space variables over the trajectories [72,73]. This follows by straightforward extension of the derivation found in Gardiner [73] for the Glauber-P representation. Such is not the case in the Truncated Wigner approach, which is based on symmetrically ordered operators, making computing useful time correlations challenging [74,75]. As an example, in Fig. 8, we show the two time density correlations g (2) 1,1 (τ ) of the driven site in the steady state: In the positive-P representation this can be calculated as , (12) where N 1 (t) = Re α 1 (t) α * 1 (t) s as defined in (6), and the factors inside the numerator average in (12) are constructed using different time values coming from the same realization. The form of g (2) 1,1 (τ ) shows the characteristic oscillations with the delay τ , as seen in previous literature on the unconventional photon blockade [70,71].
For this two site system, it is possible for us to compare to exact numerical solutions of the master equation. It can be seen in both Figs. 7 and 8, that there is a strong agreement between the positive-P and more direct numerical integration of the master equation.

B. Lieb Lattices
Going beyond the two mode case to more complicated systems, we begin by considering the much studied case of a Lieb lattice [21,23,25,26], which exhibits frustration and a flatband structure. The unit cells contain 3 sites (labeled A,B,C), and only some connections allow tunneling between cells, as per the schematic shown in Fig. 9. A 1d Lieb lattice has been implemented e.g. by polaritons in an array of micropillars [21,26]. A Lieb lattice pumped locally only on the C sites has dark B sites that have far more striking departures from coherence than the single sites of Sec. IV [69]. In directly pumped sites, the field is usually close to being pinned by the coherent pump, whereas the dark sites are free to evolve to a much less classical stationary state.  Table II. Comparison between positive-P and corner space normalization calculations for the stationary state of 1d (top) and 2d (bottom) Lieb lattices. Here, Jij = J for all connected sites, U = 0.3γ, ∆ = 0, and Fc = 0.1γ in the driven sites (C sites only) with periodic boundary conditions. nA and nB are the occupations of A and B sites, respectively, while g (2) B is the on-site two body correlation on the B sites. g (2) B,nn = a † B,j a † B,k aB,j a B,k /n 2 B is the normalized density correlation between B sites in nearest neighbor unit cells j and k.
The study of [69] used the corner space renormalization method [32] to obtain accurate predictions for small lattice sizes, which provide a convenient benchmark for the precision and accuracy of the positive-P approach. Table II shows a comparison. There is excellent agreement, and similar precision. Due to the much more favorable scaling of the positive-P method, huge lattices are easily accessible. Results for lattices with up to 100 × 100 unit cells are shown. They indicate that for 1d, the 12 unit cell lattice saturates the infinite size limit. However, for the 2d system, the 4 × 4 lattice that was achievable by corner space renormalization does not yet reach the macroscopic limit in terms of density correlations g (2) .

C. Uniform square lattices
A uniform square lattice with tunneling between all nearest neighbor sites is also of much current interest. Here, we will use the notation of [62], where J ij = J/z between nearest neighbor sites. The lattices have periodic boundary conditions and M = m × m sites in total. The coordination number is z = 4 when m > 2, and z = 2 for the special case of m = 2 in which left/right connections are to the same site.
The homogeneous case with uniform F , U , J, ∆ has been studied using a self-consistent mean field (SCMF) approach pioneered by LeBoite et al. [62]. They found a flat-band, collective excitations, and a tunneling induced transition to bistability. Later work has also shown bimodality in the photon number distribution and a hysteretic cycle around a 1st order phase transition at higher tunneling [34]. The idea behind the SCMF is that the tunneling terms in the Hamiltonian can be expressed in the mean-field picture as which is equivalent to an effective coherent driving of One then self-consistently solves for the exact quantum expressions from [61] for a in a single mode while using F eff ( a ) from (14) as the coherent driving. This can be done by iteration, starting with the bare F . It is a similar approach to the self consistent mean-field widely used for conservative Bose-Hubbard models. Eq. (14) also lets one see that it may be useful to approximate coherent transport into the region of interest with an effective driving F ≈ −J a in some systems. For cases with negligible quantum depletion, a symmetry broken "Gross-Pitaevskii" (GP) approach can also be used. This is equivalent to setting the quantum noises ξ and ξ in (5) to zero: (15) so that the observable predictions (6) and (7) reduce to N = |α| 2 and g 2 = 1. Both approaches have evident gaps in the description: the SCMF assumes a uniform system (or potentially, a local density approximation), and does not take into account any spatial correlations. The GP approach can treat inhomogeneities properly, but does not take into account quantum depletion at all. How does the full quantum approach of the positive-P compare? Figure 10. Square Lattices. Simulations with different lattice sizes for F = U , γ = 3.16 U , ∆ = 0 are shown: a 100×100 lattice with J = 2U (violet), a 2 × 2 lattice with J = 2U (green), and a single site with J = 0 (yellow). Panels show: g 2 the density-density correlation g2 averaged over all sites (top-left), the average occupation per site N (top-right), the average amplitude a (bottom-left). In both the latter, the J = 2U cases overlap. Bottom-right: the nearest neighbor 1st order coherence (normalized) as given by (16). Also shown are 1-site exact values [61] (dotted) and the SCMF predictions [62] for many modes (dashed).
First, we consider the uniformly driven case with periodic boundary conditions. F, U, ∆, and J are identical at all sites. In Fig. 10, positive-P simulations are shown for the crossover regime case of F = U , γ = 3.16U , ∆ = 0 studied in Fig. 3(b), now with a nonzero tunneling J = 2U on small and large lattices. The 1 site case is also shown for reference in yellow. The move from single site to modest lattice to huge lattice is basically effortless in terms of calculation difficulty. Numerical comparisons with standard estimates are shown in Table III. The quantity estimates full quantum single-mode [61] SCMF [62] positive-P gives the average 1st order coherence between nearest neighbor sites, where N is the average occupation.
The first thing to note is that there is a significant influence of J: basically none of the observables agree between the 1-site model shown in yellow and the lattice calculations. Furthermore, the 2 × 2 lattice is not sufficient to reach the asymptotic behavior, as seen in both 1st and 2nd order correlations. However, the mean amplitude and occupation can mislead one into thinking that the asymptotic limit has been reached. Since huge lattices of 100 × 100 are easily accessible, a positive-P calculation can be used to determine the size required to reach the asymptotic regime. In the case of the parameters of Fig. 10, a 5 × 5 lattice is needed for accurate convergence, as shown in Table.  An important question is whether the mean-field approach is faithful to the asymptotic limit of many sites. The SCMF is remarkably good for N and a , but poor for density fluctuations g 2 . This potentially sheds some doubt on past results obtained this way [34,62], at least in similar regimes. Going to even simpler approaches, a cut down version of the SCMF simply calculates the single-mode exact value with detuning modified as per ∆ → ∆ eff = ∆ + J. This is also shown in Table. III. We note that there is already a large improvement over the J = 0 estimate, except for g 2 , which is sensitive to quantum correlations. None of the estimates are able to give any information about g 1nn .

D. Non-uniform pumping
A situation where particularly large lattices are necessary is when the system is nonuniform, or excitations involve many sites collectively [8,34,60]. Systems of 10 4 − 10 6 sites pose no problem for the positive-P approach, potentially allowing for complicated geometries, extensive transport, or simulations of emergent phenomena. Fig. 11 shows results for a truly large 256 × 256 system with complicated geometry. Spreading of N (x, y) away from the pumped area is observed simultaneously with coherent spatial oscillations as a surface effect around the pumping zone. These behaviors are captured well by the Gross-Pitaevskii equation (15). The SCMF approach does not replicate the emergent local structure due to tunneling, though the bulk density is properly described.
On the other hand, the density fluctuations are not well described by either of the approximate methodsonly the positive-P gives an accurate description, even in the bulk. This last aspect is consistent with what we saw in Fig. 10 and Table. III. At the points furthest from the driven region, g 2 seems to tend towards the SCMF estimate, though it becomes very noisy, as one would also expect in experiment, due to the very low density (e.g. observe the regions around x = −70, 5, 25, 50, the furthest points from the driven region for which the occupation is still sufficient to have g 2 measurable beyond the noise). Notably, the positive-P calculation allows one to predict the spatial variation of g 2 in the vicinity of the surface, which is not possible either accurately or even qualitatively by the approximate approaches.
One realization of the simulation shown in Fig. 11 took 80s on a single PC processor (Intel Xeon E5645, 2.40GHz). Calculations on a 1000 × 1000 lattice took 1h per realization under the same fairly basic conditions. The calculation time grows approximately linearly with J for these parameters due to time step requirements.

VI. NONZERO TEMPERATURE
The master equation (3) assumes dissipation into empty modes. A more general form is which can be used to model systems coupled to baths with finite occupations N B j . The correction to the FPE of (A5) is then with complex white noises η j of mean zero that obey In Fig. 12, we give an example of a single mode with coherent drive and decay into occupied modes. Compared to the vacuum bath case, the steady state value of g 2 falls much closer to the value g 2 = 2 that would occur for thermal states; meanwhile, the coherence | â | 2 falls much lower than in previous examples with dissipation to empty modes, showing that the positive-P method still works well for less coherent states. We thus expect the method to apply to condensates with low condensate fractions, materialization and other problems with no or weak coherence.
In the absence of coherent driving F , interactions U , and tunneling J ij , Eq.(18) with the bath coupling leads to a stationary distribution of P ( v) = ⊗ j P j with This is a thermal ensemble with occupations n j = |α j | 2 and on average N B j quanta at site j. The thermal occupation of each mode with energy E j = −∆ j can be considered as Bose distributed N B j = {exp[(E j −µ)/k B T ]−1} −1 in which T and µ are resultant effective parameters of the reservoir, and in equilibrium -also of the system.
When both the tunneling and the temperature are appreciable, so that density fluctuations become important, a proper treatment of the coupling of the system to the reservoir should involve the extended single particle states, instead of the local (site) basis. For a Markovian reservoir, where particle and energy exchange takes place through interaction between system and reservoir quanta, a model that has often been used [45,[76][77][78][79][80] replaces local N B j in (17) with an effective Bose-Einstein A thermal bath of this kind in the relatively high temperature limit exp[(E − µ)/k B T ] → 1 + (E − µ)/k B T has been implemented using the positive-P method for the closely related continuum ultracold atom systems [81]. They differ from our Hamiltonian by having kinetic energy rather than site-to-site tunneling, and lacking coherent driving. Since these terms describe one-particle processes, their contribution to the stochastic differential equations is obtained by simply replacing a j → α j and a † j → α * j in the Heisenberg equations of motion. Hence for our driven-dissipative model the corresponding equa-tions become: with a reservoir at temperature T (k B = 1), and coupling constant Γ. The conditions of applicability of (19) and (22) are different, though there is an overlap regime.
In that regime one can identify the correspondences , with µ incorporated into the ∆ j . Space-dependent temperature profiles can be included through a site-dependence of T .
A rigorous derivation and consideration of applicability criteria for the equations (22) goes beyond the scope of the article, but we include them for completeness of the picture regarding thermal effects. We also mention that thermal baths that take into account the quantum particle-like nature at high energies have been implemented in [45,80].

VII. CONCLUSIONS
We have described the essential elements for applying the positive-P method to the driven dissipative Bose-Hubbard model, and benchmarked its accuracy -confirming lack of systematics down to the 4th significant digit in our test cases. The method appears to be a versatile and robust way to describe the full quantum mechanics of even very large systems, allowing for the study of the stationary state -and its time correlations -provided that the dissipation is sufficiently strong.
In particular, one can reach strong antibunching (Figs. 3(c)), high driving and occupation (Figs. 3(a)), and crossover regimes, and even a strong photon blockade with perfect antibunching and destructive interference (Fig. 7). Large non-uniform systems with 256 × 256 sites (Fig. 11) and more are easily accessible, opening the way to the study of much uncharted territory: e.g. dissipative phase transitions in a flat-band with large scale fluctuations [8,34,36], or the point at which bimodality predicted by semiclassical approaches [62,64] morphs to a 1st order phase transition [7,63]. It can also be readily used to determine the minimal sizes of systems required to reach the asymptotic regime. We have studied Lieb lattices (Table II) and simple orthogonal 2d lattices (Table. IV, Fig. 10) in this regard, showing that 4×4 systems e.g. tend not to be in the asymptotic limit.
The method exhibits clear superiority over various mean field approaches and the truncated Wigner approximation in the difficult regime when occupation is small, and provides the ability to simultaneously and accurately study correlations, interference, tunneling and nonlocal effects. Due to the numerical instabilities, the positive-P approach cannot describe strongly driven weakly damped conditions. However, this more semiclassical regime is extremely well approximated by the related truncated Wigner method with (from the technical point of view) very similar stochastic equations to be solved. Thus, between truncated Wigner, which gives very accurate results for large occupations, and positive-P, we have a viable method for all conditions where either drive or dissipation are significant effects. Notably, the positive-P approach gives full quantum results in the medium to large dissipation regime, whereas most other full quantum approaches such as DMRG, tensor networks etc. work more easily under the opposite, low dissipation, conditions.
In the usable regime, numerical effort scales merely linearly with the number of sites, and quadratically with the precision. Space and time-dependence of all parameters in the model is easily incorporated with no extra numerical effort. Nonlocal interactions can also be efficiently treated [82]. Thus we suggest positive-P as the method of choice to access large systems in the very regions that are currently experimentally relevant, especially in drivendissipative but correlated photonic platforms.
Due to the additional stability provided by dissipation, the positive-P is applicable to a much wider range of problems in open dissipative systems than in closed systems. Its great success in describing the archetypal driven dissipative Bose-Hubbard model shown here implies that positive-P may be an ideal method for simulating various kinds of open quantum systems that either consist of or can be mapped onto bosons. Such promising extensions, such as incoherent driving or systems with coupled spins and bosons, will be the subject of future work. the stronger restriction in practically all cases. In particular, there are very large regions in which this error is seen without being masked by attendant statistical error. We also notice that at low γ the highest systematic error comes from the coherent amplitude characterized by | a | 2 /N , whereas at high γ the limiting systematic error is from g 2 , though the corresponding errors in N and | a | are also comparable. We finally take the largest relative error as the overall assessment of the errors expected from the truncated Wigner. This quantity is displayed in Fig. 6. From the above discussion, its behavior closely follows the systematic relative errors in either the coherent amplitude | a | 2 /N or in the g 2 .