  We describe an efficient numerical method for simulating the dynamics of
interacting spin ensembles in the presence of dephasing and decay. The method
builds on the discrete truncated Wigner approximation for isolated systems,
which combines the mean-field dynamics of a spin ensemble with a Monte Carlo
sampling of discrete initial spin values to account for quantum correlations.
Here we show how this approach can be generalized for dissipative spin systems
by replacing the deterministic mean-field evolution by a stochastic process,
which describes the decay of coherences and populations while preserving the
length of each spin. We demonstrate the application of this technique for
simulating nonclassical spin-squeezing effects or the dynamics and steady
states of cavity QED models with hundred thousand interacting two-level systems
and without relying on any symmetries. This opens up the possibility to perform
accurate real-scale simulations of a diverse range of experiments in quantum
optics or with solid-state spin ensembles under realistic laboratory
conditions.
