In circuit QED, a Rabi model can be derived by truncating the Hilbert space of the anharmonic qubit coupled to a linear, reactive environment. This truncation breaks the gauge invariance present in the full Hamiltonian. We analyze the determination of an optimal gauge such that the differences between the truncated and the full Hamiltonian are minimized. Here, we derive a simple criterion for the optimal gauge. We find that it is determined by the ratio of the anharmonicity of the qubit to an averaged environmental frequency. We demonstrate that the usual choices of flux and charge gauge are not necessarily the preferred options in the case of multiple resonator modes.Circuit QED [1, 2] is a central subject of quantum information science that has deepened our understanding of light-matter interaction [3][4][5]. Most implementations consist of a two-level system (qubit) that is coupled to a linear environment. The qubit is formed by the two lowest energy levels of an anharmonic multilevel-system. For the physics of interest only the qubit subspace is important. The Schrieffer-Wolff (SW) transformation [6,7] is the standard method to perturbatively derive an effective Hamiltonian description within this subspace. For most purposes, it is sufficient to consider the effective Hamiltonian only to first order, yielding the well known quantum Rabi model (QRM). However, since the Hamiltonian of the non-truncated system is unique only up to a unitary transformation, the effective description is gauge dependent to every finite order [8,9]. This gauge ambiguity becomes particularly important in the (ultra) strong coupling regime. It has been found that the QRM derived in a gauge where the qubit-resonator coupling is mediated by the flux variables leads to different predictions than the one where the coupling is mediated by the charge variables [10][11][12].In this work, we look at the issue from a different perspective. We use the gauge degree of freedom to find an optimal gauge such that the results of the effective model are as close as possible to full model. Importantly, we take account of the need for a multimode description [13][14][15][16] in the quest for achieving the ultra-strong coupling regime [17][18][19][20]. To increase the flexibility, we not only consider the extremal cases of purely flux or charge mediated coupling but perform a general gauge transformation that smoothly interpolates between the two. A similar transformation has been used in [21] to extend the Jaynes-Cummings model into the ultra-strong coupling regime.We find that already the second order term of the effective Hamiltonian within the SW method is a good indicator of the validity of the QRM. Based on this observation, we derive a simple analytical criterion for the optimal gauge and benchmark it against numerical sim- FIG. 1. (a) Circuit diagram of a qubit with potential U (φq), capacitance Cq, and inductance Lq that is coupled to a general reactive environment. In the Forster form, the latter is represented by N resonators with capacitance...