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.