  The efficient simulation of quantum systems is a primary motivating factor
for developing controllable quantum machines. For addressing systems with
underlying bosonic structure, it is advantageous to utilize a naturally bosonic
platform. Optical photons passing through linear networks may be configured to
perform quantum simulation tasks, but the efficient preparation and detection
of multiphoton quantum states of light in linear optical systems are
challenging. Here, we experimentally implement a boson sampling protocol for
simulating molecular vibronic spectra [Nature Photonics $\textbf{9}$, 615
(2015)] in a two-mode superconducting device. In addition to enacting the
requisite set of Gaussian operations across both modes, we fulfill the
scalability requirement by demonstrating, for the first time in any platform, a
high-fidelity single-shot photon number resolving detection scheme capable of
resolving up to 15 photons per mode. Furthermore, we exercise the capability of
synthesizing non-Gaussian input states to simulate spectra of molecular
ensembles in vibrational excited states. We show the re-programmability of our
implementation by extracting the spectra of photoelectron processes in H$_2$O,
O$_3$, NO$_2$, and SO$_2$. The capabilities highlighted in this work establish
the superconducting architecture as a promising platform for bosonic
simulations, and by combining them with tools such as Kerr interactions and
engineered dissipation, enable the simulation of a wider class of bosonic
systems.
