+ \mathcal{D}[\hat{\rho}]$$

A mean-field description of the problem can be derived in terms of coherent fields $\psi_j = \langle \hat{a}_j \rangle$ by assuming that all normal-ordered operator products factorize [@Iacopo_QFL]. This leads to the following equations of motion $$\label{eq:MF}
i \dot{\psi}_j = -\Big(\delta + i\frac{\gamma}{2}\Big)\psi_j -J(\psi_{j+1} + \psi_{j-1}) + U|\psi_j|^2\psi_j + \Omega_j$$

For this work we assume a uniform drive field $\Omega_j=\Omega$ at all sites in the chain. In the steady state we thus find one or two stable homogeneous density solutions in the mean-field description, depending on the amplitude and detuning of the pump. They are found as solutions of $$n_0((\delta - Un_0 + 2J)^2 + \gamma^2/4) = |\Omega|^2,$$ where $n_0=|\psi_0|^2$ in terms of the spatially uniform mean-field steady state $\psi_j = \psi_0$. The parameter $$\label{eq:delta}
 \Delta= \delta - Un_0 + 2J$$ is the renormalized laser detuning from the interaction-blueshifted optical resonance. We restrict our analysis to the case $\Delta < 0$, such that the system is in the optical limiter regime or in the high-density branch of the hysteresis loop of a bistable regime [@Iacopo_QFL]. As we will briefly review later, this restriction asserts a gapped spectrum of excitations [@Iacopo_probing]. In the numerical analysis that follows, the drive 