  We provide the necessary framework for carrying out stochastic positive-P and
gauge-P simulations of bosonic systems with long range interactions. In these
approaches, the quantum evolution is sampled by trajectories in phase space,
allowing calculation of correlations without truncation of the Hilbert space or
other approximations to the quantum state. The main drawback is that the
simulation time is limited by noise arising from interactions.
  We show that the long-range character of these interactions does not further
increase the limitations of these methods, in contrast to the situation for
alternatives such as the density matrix renormalisation group. Furthermore,
stochastic gauge techniques can also successfully extend simulation times in
the long-range-interaction case, by making using of parameters that affect the
noise properties of trajectories, without affecting physical observables.
  We derive essential results that significantly aid the use of these methods:
estimates of the available simulation time, optimized stochastic gauges, a
general form of the characteristic stochastic variance and adaptations for very
large systems. Testing the performance of particular drift and diffusion gauges
for nonlocal interactions, we find that, for small to medium systems, drift
gauges are beneficial, whereas for sufficiently large systems, it is optimal to
use only a diffusion gauge.
  The methods are illustrated with direct numerical simulations of interaction
quenches in extended Bose-Hubbard lattice systems and the excitation of Rydberg
states in a Bose-Einstein condensate, also without the need for the typical
frozen gas approximation. We demonstrate that gauges can indeed lengthen the
useful simulation time.
