  There has been a problem of gauge ambiguities with the Rabi Hamiltonian due
to the fact that it can be derived from two formally different but physically
equivalent fundamental Hamiltonians. This problem has recently been resolved
for models with single quantized electromagnetic mode. In this work, we
mathematically and numerically verify this for multimode models. With this
established, we combine the numerical methods, matrix product states (MPS) and
numerical mode decomposition (NMD), for analyzing cavity QED systems. The MPS
method is used to efficiently represent and time evolve a quantum state.
However, since the coupling structure of the Rabi Hamiltonian is incompatible
with MPS, it is numerically transformed into an equivalent Hamiltonian that has
a chain coupling structure, which allows efficient application of MPS. The
technique of NMD is used to extract the numerical electromagnetic modes of an
arbitrary environment. As a proof of concept, this combined approach is
demonstrated by analyzing 1D cavity QED systems in various settings.
