A cluster algorithm for Monte Carlo simulations of spin ice

We present an algorithm for Monte Carlo simulations of a nearest-neighbor spin ice model based on its cluster representation. To assess its performance, we estimate a relaxation time, and find that, in contrast to the Metropolis algorithm, our algorithm does not develop spin-freezing. Also, to demonstrate the efficiency, we calculate the spin and charge structure factors, and observe pinch points in a high-resolution color map. We then find that Debye screening works among defects and brings about short-range correlations, and that the deconfinement transition triggered by a fugacity of defects $z$ is dictated by a singular part of the free-energy density $f_{\rm s}\propto z$.