  We study the dynamics of lattice models of quantum spins one-half, driven by
a coherent drive and subject to dissipation. Generically the meanfield limit of
these models manifests multistable parameter regions of coexisting steady
states with different magnetizations. We introduce an efficient scheme
accounting for the corrections to meanfield by correlations at leading order,
and benchmark this scheme using high-precision numerics based on
matrix-product-operators in one- and two-dimensional lattices. Correlations are
shown to wash the meanfield bistability in dimension one, leading to a unique
steady state. In dimension two and higher, we find that multistability is again
possible, provided the thermodynamic limit of an infinitely large lattice is
taken first with respect to the long time limit. Variation of the system
parameters results in jumps between the different steady states, each showing a
critical slowing down in the convergence of perturbations towards the steady
state. Experiments with trapped ions can realize the model and possibly answer
open questions in the nonequilibrium many-body dynamics of these quantum
systems, beyond the system sizes accessible to present numerics.
