  We study the ground state phase diagram of a nonlinear two-photon
Rabi-Hubbard (RH) model in one dimension using quantum Monte Carlo (QMC)
simulations and density matrix renormalization group (DMRG) calculations. Our
model includes a nonlinear photon-photon interaction term. Absent this term,
the RH model has only one phase, the normal disordered phase, and suffers from
spectral collapse at larger values of the photon-qubit interaction or
inter-cavity photon hopping. The photon-photon interaction, no matter how
small, stabilizes the system which now exhibits {\it two} quantum phase
transitions: Normal phase to {\it photon pair} superfluid (PSF) transition and
PSF to single particle superfluid (SPSF). The discrete $Z_4$ symmetry of the
Hamiltonian spontaneously breaks in two stages: First it breaks partially as
the system enters the PSF and then completely breaks when the system finally
enters the SPSF phase. We show detailed numerical results supporting this, and
map out the ground state phase diagram.
