  NV centers in diamond are a promising platform for highly sensitive quantum
sensors for magnetic fields and other physical quantities. The quest for high
sensitivity combined with high spatial resolution leads naturally to dense
ensembles of NV centers, and hence to strong, long-range interactions between
them. Hence, simulating strongly interacting NVs becomes essential. However,
obtaining the exact dynamics for a many-spin system is a challenging task due
to the exponential scaling of the Hilbert space dimension, a problem that is
exacerbated when the system is modelled as an open quantum system. In this
work, we employ the Matrix Product Density Operator (MPDO) method to represent
the many-body mixed state and to simulate the dynamics of an ensemble of NVs in
the presence of strong long-range couplings due to dipole-dipole forces. We
benchmark different time-evolution algorithms in terms of numerical accuracy
and stability. Subsequently, we simulate the dynamics in the strong interaction
regime, with and without dissipation.
