$1\leq i\leq N$ and $\mathbf{\Omega}_{R}$ is diagonal with entries $\left(\omega_{R_{1}},\ldots,\omega_{R_{M}}\right)$.

We now block-diagonalize $\mathbf{M}_{1}$ by applying a Schrieffer-Wolff transformation to get

$$\widetilde{\mathbf{M}}_{1}=\exp\left(-\mathbf{S}\right)\mathbf{M}_{1}\exp\left(\mathbf{S}\right)\label{eq:block-diagonal-M1}$$

where $\mathbf{S}$ is skew-symmetric and $\widetilde{\mathbf{M}}_{1}$ block-diagonal which can be computed up to desired order in the bare couplings using Eqs. (B.12) and (B.15) in [@Winkler]. We note that since this transformation is unitary it will keep the capacitance matrix identity such that we have the following block-diagonal Hamiltonian in the final frame

$$\mathcal{H}=\frac{1}{2}\boldsymbol{q}^{T}\boldsymbol{q}+\frac{1}{2}\boldsymbol{\phi}^{T}\widetilde{\mathbf{M}}_{1}\boldsymbol{\phi}+\mathcal{O}(\varphi_{J}^{4})\label{eq:block-diagonal-hamiltonian}$$

where the final coordinate fluxes $\boldsymbol{\phi}$ are related to the initial coordinates $\mathbf{\Phi}$ by

$$\mathbf{\Phi}=\left(\begin{array}{c}
\mathbf{\Phi}_{J}\\
\mathbf{\Phi}_{R}
\end{array}\right)=\left(\begin{array}{cc}
\mathbf{C}_{0}^{-1/2} & \mathbf{0}\\
\mathbf{0} & \mathbf{1}
\end{array}\right)\mathbf{T}\exp(\mathbf{S})\boldsymbol{\phi}$$

and $\mathcal{O}(\varphi_{J}^{4})$ term standing for higher order nonlinear corrections giving anharmonicities and dispersive shifts between modes calculated in Appendix .

Using Eq. (B.15c) in [@Winkler], to second order in the bare couplings

$$\left(\widetilde{\mathbf{M}}_{1}\right)_{ij}=\frac{1}{2}\underset{k}{\sum}\left(\mathbf{M}_{1}\right)_{ik}\left(\mathbf{M}_{1}\right)_{kj}\left[\frac{1}{\omega_{i}^{2}-\omega_{R_{k}}^{2}}+\frac{1}{\omega_{j}^{2}-\omega_{R_{k}}^{2}}\right]\label{eq:M1-tilde-ij}$$

where again $i$ and $j$ are qubit labels and $k$ labels internal modes. $\left(\mathbf{M}_{1}\right)_{ik}$ is the $\left(i,k\right)$-th entry of the matrix $\mathbf{M}_{1}$ and from Eq. we have

$$\left(\mathbf{M}_{1}\right)_{ik}=\omega_{i}^{2}C_{i}^{1/2}r_{ki}\label{eq:Mik}$$

Noting again $r_{k}^{T}r_{k}=\mathbf{A}_{k}$ we can write

$$\begin{aligned}
\left(\mathbf{M}_{1}\right)_{ik}\left(\mathbf{M}_{1}\right)_{kj} & = & \omega_{i}^{2}\omega_{j}^{2}C_{i}^{1/2}C_{j}^{1/2}r_{ki}r_{kj}\nonumber \\
 & = & \omega_{i}^{2}\omega_{j}^{2}C_{i}^{1/2}C_{j}^{1/2}\left(\mathbf{A}_{k}\right)_{ij}\end{aligned}$$

Hence we can re-write Eq.

$$\begin{aligned}
\left(\widetilde{\mathbf{M}}_{1}\right)_{ij} & = & \frac{1}{2}\underset{k}{\sum}\omega_{i}^{2}\omega_{j}^{2}C_{i}^{1/2}C_{j}^{1/2}\left[\frac{\left(\mathbf{A}_{k}\right)_{ij}}{\omega_{i}^{2}-\omega_{R_{k}}^{2}}+\frac{\left(\mathbf{A}_{k}\right)_{ij}}{\omega_{j}^{2}-\omega_{R_{k}}^{2}}\right]\nonumber \\
 & = & -\frac{1}{2}\omega_{i}^{2}\omega_{j}^{2}C_{i}^{1/2}C_{j}^{1/2}\mathrm{Im}\left[\frac{Z_{ij}(\omega_{i})}{\omega_{i}}+\frac{Z_{ij}(\omega_{j})}{\omega_{j}}\right]\end{aligned}$$

Quantizing the system by introducing 