  Analog quantum simulation - the technique of using one experimentally
well-controlled physical system to mimic the behavior of another - has quickly
emerged as one of the most promising near term strategies for studying strongly
correlated quantum many-body systems. In particular, systems of interacting
photons, realizable in solid-state cavity and circuit QED frameworks, for
example, hold tremendous promise for the study of nonequilibrium many-body
phenomena due to the capability to locally create and destroy photons. These
systems are typically modeled using a Jaynes-Cummings-Hubbard (JCH)
Hamiltonian, named due to similarities with the Bose-Hubbard (BH) model. Here,
we present a non-perturbative procedure for transforming the JC Hamiltonian
into a dressed operator representation that, in its most general form, admits
an infinite sum of bosonic $k$-body terms where $k$ is bound only by the number
of excitations in the system. We closely examine this result in both the
dispersive and resonant coupling regimes, finding rapid convergence in the
former and contributions from $k\gg1$ in the latter. Through extension to a
two-site JCH system, we demonstrate that this approach facilitates close
inspection of the analogy between the JCH and BH models and its breakdown for
resonant light-matter coupling. Finally, we use this framework to survey the
many-body character of a two-site JCH for general system parameters,
identifying four unique quantum phases and the parameter regimes in which they
are realized, thus highlighting phenomena realizable with finite JCH-based
quantum simulators beyond the BH model. More broadly, this work is intended to
serve as a clear mathematical exposition of bosonic many-body interactions
underlying JC-type systems, often postulated through analogy to Kerr-like
nonlinear susceptibilities or by matching coefficients to obtain the
appropriate eigenvalue spectrum.
