We demonstrate mutual synchronization of Josephson oscillations in multiple stacks of intrinsic Josephson junctions of the cuprate superconductor Bi$_2$Sr$_2$CaCu$_2$O$_{8+delta }$. Detailed analysis of the full polarization parameters allows the determination of a phase correlation between the stacks: a simultaneous emission state is described by a linear combination of individual emission states with a phase retardation. This proves that the stacks are coupled via a Josephson plasma in a superconducting substrate and the coupling matrices can be extracted from polarization analyses. Our findings suggest a route towards the realization of high-power terahertz sources based on the synchronization of a large number of intrinsic Josephson junctions.