Magnetic-charge ordering and phase transitions in monopole-conserved square spin ice

Magnetic-charge ordering and corresponding magnetic/monopole phase transitions in spin ices are the emergent topics of condensed matter physics. In this work, we investigate a series of magnetic-charge (monopole) phase transitions in artificial square spin ice model using the conserved monopole density algorithm. It is revealed that the dynamics of low monopole density lattices is controlled by the effective Coulomb interaction and the Dirac string tension, leading to the monopole dimerization which is quite different from the dynamics of three-dimensional pyrochlore spin ice. The condensation of the monopole dimers into monopole crystals with staggered magnetic-charge order can be predicted clearly. For the high monopole density cases, the lattice undergoes two consecutive phase transitions from high-temperature paramagnetic/charge-disordered phase into staggered charge-ordered phase before eventually toward the long-range magnetically-ordered phase as the ground state which is of staggered charge order too. A phase diagram over the whole temperature-monopole density space, which exhibits a series of emergent spin and monopole ordered states, is presented.


Simulation details
We employ the Monte Carlo method to track the magnetic states and associated phase transitions in such dipolar square spin ice model with different monopole density. In our simulations, the long-range dipolar interaction is treated using the Ewald summation scheme, avoiding the errors induced by the conventional finite truncation method. In order to track the monopole-ordered phases, we employ the conserved monopole algorithm (CMA) 1 , which applies to a statistical ensemble with a fixed monopole density over the whole T-range.
According to the CMA, any spin flip creating or destroying monopole is strictly forbidden.
Considering a neighboring Vi-Vj vertex-pair, which is denoted as the Vi,j-pair, any single spin flip will carry a monopole from one vertex to the other. Therefore, relationship E1,3=E1,4<E2,3=E2,4 is satisfied, where Ei,j is the energy of the vertex-pair Vi-Vj. It is clear that the spin flip from the V2,3-pair to the V1,3-pair or from the V2,4-pair to the V1,3-pair has higher probability than that for the flip from the V1,3-pair to the V2,3-pair or from the V1,3-pair to the V2,4-pair, because we have E2,3-E1,3=E2,4-E1,4=E2-E1=>0. However, the single spin flip scheme alone is very inefficient to reach the ground state and one may include other schemes which can accelerate the ground state searching in addition to the single spin flips. Here, the loop spin flip method 2 is used. This method will not move the monopoles but is very efficient in bypassing the energy barriers that separate degenerate states and allows efficient tracking the long-range ordered states without violating the ice rule. Besides, the loop spin flip is a prominent conserved monopole algorithm because it strictly sticks to the flip rules of CMA.
The simulation starts from an initial state with fixed monopole defects (i.e. fixed and equivalent numbers of V3 vertex and V4 vertex). Such initial state can be obtained by the following procedure. We first fill up the lattice with the V1 vertices, and then run the spin flips which can create new monopoles until a desired monopole density () in the lattice is reached.
Subsequently, we allow sufficient loop spin flip events to reach a random state as the initial state. In the following Monte Carlo steps, we employ the parallel temperature Monte Carlo scheme 3 which is very efficient for complicated spin systems and frequently used in parallel algorithm. Each Monte Carlo step in our simulation contains N/2 single spin flips and N/2 loop spin flips, where N is the lattice size (L 2 ). At each temperature point, we complete 2×10 5 MC steps for relaxation process and perform 1×10 5 sampling averages in the next 5×10 5 MC steps.

Physical properties of the system in the  =0.5 case
Here we present some of our simulation results at the case  =0.5. Figure S1 shows the specific heat and the order parameters of the system. There are three order parameters in our system, and they are the vacuum order parameter (MV), charge crystal order parameter (MC) and the staggered charge order parameter (Q). Three peaks in the specific heat curve indicate that the system goes through three phase transitions in an annealing process. At the temperature range T/D>0.4, the system is in the monopole fluid state. As the temperature decreases, MV quickly rises to 1, indicating that the system is in the vacuum ordered state. Furthermore, a phase transition occurs at about 0.23D leading the system to a staggered charge ordered state. When the temperature decreases blow 0.1D, the charge crystal order parameter rapidly rises to maximum, corresponding to a phase transition to the monopole crystal phase.