\delta \sum_{j=1}^L \hat{n}_j + \frac{U}{2} \sum_{j=1}^L \hat{n}_j(\hat{n}_j-1)\\ &&+ \sum_{j=1}^L\Omega_j ( \hat{a}_j + \hat{a}^\dagger_j).
\label{eq:Hamiltonian}\end{aligned}$$ The operators $\hat{a}^\dagger_j(\hat{a}_j)$ create (annihilate) a particle at site $j$ of the chain and $\hat{n}_j = \hat{a}^\dagger_j\hat{a}_j$ is the local number operator. Photons in the chain can tunnel to their neighbouring sites with a hopping strength given by $J$. The notation $\langle j,l\rangle$ means that the summation runs over all neighbouring sites. The two-body interaction strength for photons confined inside the same cavity is given by $U$. The amplitude of the driving field at each site is $\Omega_j$, while its detuning from the onsite single-photon resonance is given by $\delta=\omega_L-\omega_c$. For simplicity we impose periodic boundary conditions, such that $\hat{a}_1=\hat{a}_{L+1} $.

The dissipative nature of the setup implies that injected photons have a finite lifetime inside the cavity array before they escape. In the Born-Markov approximation, the coupling of the system to its environment at zero temperature is described by the dissipator in the Lindblad form [@gardiner_quantum_noise; @open_quantum_systems] $$\label{eq:dissipator}
\mathcal{D}[\hat{\rho}]=\frac{\gamma}{2}\sum_j(2\hat{a}_j\hat{\rho}\hat{a}^\dagger_j - \hat{\rho}\hat{n}_j - \hat{n}_j\hat{\rho}).$$ The full dynamics of the density matrix $\hat{\rho}$ is then governed by a master equation, which includes both the unitary evolution under $\hat{H}$ and the photonic losses $$\label{eq:q_master}
\partial_t \hat{\rho} = -i[\hat{H},\hat{\rho}] 