We study the two-dimensional two-component Coulomb gas in the canonical ensemble and at inverse temperature $β>2$. In this regime, the partition function diverges and the interaction needs to be cut off at a length scale $λ\in (0,1)$. Particles of opposite charges tend to pair into dipoles of length scale comparable to $λ$, which themselves can aggregate into multipoles. Despite the slow decay of dipole--dipole interactions, we construct a convergent cluster expansion around a hierarchical reference model that retains only intra-multipole interactions. This yields a large deviations result for the number of $2p$-poles as well as a sharp free energy expansion as $N\to\infty$ and $λ\to0$ with three contributions: (i) the free energy of $N$ independent dipoles, (ii) a perturbative correction, and (iii) the contribution of a non-dilute subsystem. The perturbative term has two equivalent characterizations: (a) a convergent Mayer series obtained by expanding around an i.i.d.\ dipole model; and (b) a variational formula as the minimum of a large-deviation rate function for the empirical counts of $2p$-poles. The Mayer coefficients exhibit transitions at $β_p=4-\tfrac{2}{p}$, that accumulate at $β=4$, which corresponds to the Berezinskii-Kosterlitz-Thouless transition in the low-dipole-density limit. At $β=β_p$ the $p$-dipole cluster integrals switch from non-integrable to integrable tails. The non-dilute system corresponds to the contribution of large dipoles: we exhibit a new critical length scale $R_{β, λ}$ which transitions from $λ^{-(β-2)/(4-β)}$ to $+\infty$ as $β$ crosses the critical inverse temperature $β=4$, and which can be interpreted as the maximal scale such that the dipoles of that scale form a dilute set.