Bosonization of primary fields for the critical Ising model on multiply connected planar domains

We prove bosonization identities for the scaling limits of the critical Ising correlations in finitely-connected planar domains, expressing those in terms of correlations of the compactified Gaussian free field. This, in particular, yields explicit expressions for the Ising correlations in terms of domain's period matrix, Green's function, harmonic measures of boundary components and arcs, or alternatively, Abelian differentials on the Schottky double. Our proof is based on a limiting version of a classical identity due to D.~Hejhal and J.~Fay relating Szeg\H{o} kernels and Abelian differentials on Riemann surfaces, and a systematic use of operator product expansions both for the Ising and the bosonic correlations.