validity of the two-level approximation, which is an essential step in this derivation, depends explicitly on the choice of gauge once the system enters the ultrastrong coupling regime. In particular, while in the electric dipole gauge the two-level approximation can be performed as long as the Rabi frequency remains much smaller than the energies of all higher-lying levels, it can dramatically fail in the Coulomb gauge, even for systems with an extremely anharmonic spectrum. We extensively investigate this phenomenon both in the single-dipole (Rabi) and multi-dipole (Dicke) case, and considering the specific examples of dipoles confined by double-well and by square-well potentials, and of circuit QED systems with flux qubits coupled to an LC resonator.