  Strong coupling of quantum emitters with confined electromagnetic modes of
nanophotonic structures may be used to change optical, chemical and transport
properties of materials, with significant theoretical effort invested towards a
better understanding of this phenomenon. However, a full theoretical
description of both matter and light is an extremely challenging task. Typical
theoretical approaches simplify the description of the photonic environment by
describing it as a single or few modes. While this approximation is accurate in
some cases, it breaks down strongly in complex environments, such as within
plasmonic nanocavities, and the electromagnetic environment must be fully taken
into account. This requires the quantum description of a continuum of bosonic
modes, a problem that is computationally hard. We here investigate a compromise
where the quantum character of light is taken into account at modest
computational cost. To do so, we focus on a quantum emitter that interacts with
an arbitrary photonic spectral density and employ the cumulant or cluster
expansion method to the Heisenberg equations of motion up to first, second and
third order. We benchmark the method by comparing with exact solutions for
specific situations and show that it can accurately represent dynamics for many
parameter ranges.
