I found the answer to my question and I'm going to share it here for future learners, despite going against some haughty moderators' criteria, I think it constitutes a good question.
First of all, indeed $Q_g = Q_1 - Q_2$ or any similar reasoning is incorrect. However, the rest of the equations were well posed. Let me reproduce them here in addition with the final missing step:
I want to find the charge $Q_1$ and $Q_2$ of $C_1$ and $C_2$ to quantise the problem. My guess for the Kirchhoff voltage loops and current equations would be
\begin{align}
\text{Voltage loops: }&\begin{cases}
V_{C_1} + V_{L_1} = \frac{1}{C_1} Q_{1} + L_1 \dot{I}_1 = 0 \\
V_{C_2} + V_{L_2} = \frac{1}{C_2} Q_2 - L_2\dot{I}_4 = 0
\end{cases}\notag\\
\text{Current junctions: }&\begin{cases}
I_1 = I_2 + I_3 \\
I_2 = I_4 + I_5.
\end{cases}\notag
\end{align}
And by using the current junction conditions we can write $I_1$ and $I_4$ in terms of the capacitors charge:
\begin{equation}
\begin{cases}
\frac{1}{C_1} Q_{1} + L_1 (\dot{I}_2 + \dot{I}_3) = \frac{1}{C_1} Q_{1} + L_1 (\ddot{Q}_g + \ddot{Q}_1) = 0\\
\frac{1}{C_2} Q_2 + L_2(\dot{I}_2 - \dot{I}_5) = \frac{1}{C_2} Q_2 - L_2(\ddot{Q}_g - \ddot{Q}_2) = 0.
\end{cases}\notag
\end{equation}
Now, I want to eliminate $Q_g$ since I'm only interested in $Q_1$ and $Q_2$. For this, we can read off one more loop equation for the loop concerning only the capacitors. This is
$$\frac{Q_2}{C_2} + \frac{Q_g}{C_g} - \frac{Q_1}{C_1} = 0, $$
which is the correct expression for the charge conservation. With this, we yield
\begin{equation}
\begin{cases}
\frac{1}{C_1} Q_{1} + L_1 (\ddot{Q}_g + \ddot{Q}_1)=\frac{1}{C_1} Q_{1} + L_1\ddot{Q}_1 + C_g L_1(\frac{\ddot{Q}_1}{C_1} -\frac{\ddot{Q}_2}{C_2} ) = 0\\
\frac{1}{C_2} Q_2 - L_2(\ddot{Q}_g - \ddot{Q}_2) =\frac{1}{C_2} Q_{2} + L_2\ddot{Q}_2 + C_g L_2(\frac{\ddot{Q}_2}{C_2} -\frac{\ddot{Q}_1}{C_1} )= 0,
\end{cases}\notag
\end{equation}
which can of course be symmetrised a bit by dividing
\begin{equation}
\begin{cases}
\frac{1}{L_1 C_1} Q_{1} + \ddot{Q}_1 + C_g (\frac{\ddot{Q}_1}{C_1} -\frac{\ddot{Q}_2}{C_2} ) = 0\\
\frac{1}{C_2 L_2} Q_{2} + \ddot{Q}_2 + C_g(\frac{\ddot{Q}_2}{C_2} -\frac{\ddot{Q}_1}{C_1} )= 0.
\end{cases}\notag
\end{equation}
It will be easier to work with the variable $\Phi_n = Q_n/C_n $ instead so that we can read off the lagrangian easily. With that we have
\begin{equation}
\begin{cases}
\frac{1}{L_1} \Phi_{1} + C_1 \ddot{\Phi}_1 + C_g(\ddot{\Phi}_1 -\ddot{\Phi}_2 ) = 0\\
\frac{1}{ L_2} \Phi_{2} + C_2 \ddot{\Phi}_2 + C_g(\ddot{\Phi}_2 -\ddot{\Phi}_1 )= 0,
\end{cases}\notag
\end{equation}
with which we reverse engineer the Lagrangian to be
$$ \mathcal{L}(\Phi_1,\Phi_2;\dot{\Phi}_1,\dot{\Phi}_2) = \frac{C_1}{2}\dot{\Phi}_1^2 - \frac{1}{2}\frac{\Phi_1^2}{L_1} + \frac{C_2}{2}\dot{\Phi}_2^2 - \frac{1}{2}\frac{\Phi_2^2}{L_2} + \frac{C_g}{2}\left(\dot{\Phi}_1 - \dot{\Phi}_2 \right)^2, $$
as desired.
Quantasing the problem is now trivial. First, we perform a Legendre transformation to obtain the classical Hamiltonian
$$ \mathcal{H} = \sum_n\left( \dot{\Phi_n} \frac{\partial \mathcal{L}}{\partial \dot{\Phi_n}}\right) - \mathcal{L} = \frac{C_1}{2}\dot{\Phi}_1^2 + \frac{1}{2}\frac{\Phi_1^2}{L_1} + \frac{C_2}{2}\dot{\Phi}_2^2 + \frac{1}{2}\frac{\Phi_2^2}{L_2} + \frac{C_g}{2}(\dot{\Phi}_1 - \dot{\Phi}_2 )^2$$
which written in terms of canonical conjugates, $ p_n \doteq \frac{\partial \mathcal{L}}{\partial \dot{\Phi_n}} = C_n \dot{\Phi}_n $, is simply
$$ \mathcal{H} = \frac{1}{2C_1}p_1^2 + \frac{1}{2}\frac{x_1^2}{L_1} + \frac{1}{2C_2}p_2^2 + \frac{1}{2}\frac{x_2^2}{L_2} + \frac{C_g}{2}\left(\frac{p_1}{C_1} - \frac{p_2}{C_2} \right)^2.$$
We could finish here and conclude these are two harmonic oscillators coupled via quadratic terms and call it a day. However, we can simplify the coupling term a bit further by grouping the terms proportional to momentum squared:
$$ \mathcal{H} = \left(\frac{1}{2C_1} + \frac{C_g}{2C_1^2}\right)p_1^2 + \frac{1}{2}\frac{x_1^2}{L_1} + \left(\frac{1}{2C_2} + \frac{C_g}{2C_2^2}\right)p_2^2 + \frac{1}{2}\frac{x_2^2}{L_2} - C_g\frac{p_1}{C_1} \frac{p_2}{C_2},$$
this was at the expense of making the oscillators dependent on the coupling $C_g$ which isnt necessary. We now define bosonic creation and annihilation operators,
$$ x_n = \sqrt{\frac{\lambda_n}{2}}(a^\dagger_n + a_n), \quad p_n = i\sqrt{\frac{1}{2 \lambda_n}}(a^\dagger_n - a_n), \quad \text{together with} \quad [x,p] = i \iff [a_n, a_n^\dagger]= 1.$$
with all other commutators vanishing, and leave us the freedom to choose $\lambda_n $ wisely. A wise choice for the diagonalisation is $\lambda_n = \sqrt{2L_n \left(\frac{1}{2C_n} + \frac{C_g}{2C_n^2}\right)} $. With this, we have
$$ \mathcal{H} = \sum_n \sqrt{\frac{1}{L_n}\left(\frac{1}{C_n} + \frac{C_g}{C_n^2}\right)} \left(a^\dagger_n a_n +\frac{1}{2}\right) + \frac{C_g}{2\sqrt{\lambda_1\lambda_2}}\left(a^\dagger_1 - a_1\right)\left(a^\dagger_2 - a_2\right).$$
It is clear in the limit $C_g=0$ this reduces to two independent LC oscillators.
