We study the dynamics of lattice models of quantum spins one-half, driven by a coherent drive and subject to dissipation. Generically the mean-field 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 mean field 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 mean-field 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 longtime 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.