Two- and three-point functions at criticality: Monte Carlo simulations of the three-dimensional $(q+1)$-state clock model

We simulate the improved $(q+1)$-state clock model on the simple cubic lattice at the critical point on lattices of a linear size up to $L=960$. We compute operator product expansion (OPE) coefficients for the three-dimensional XY universality class. These are compared with highly accurate estimates obtained by using the conformal bootstrap method. We find that the results are consistent.