FIG. 5: (Color online) Spin collinearity P [Eq. (19)] at h = 0 for b = 0 . 1 , showing the onset of nematic order at T Q ≈ 0 . 13 . P is measured for the farthest spin pair along the 〈 110 〉 chains in the py- rochlore lattice in each system size ranging from L = 4 to L = 16 , and averaged over the 〈 110 〉 chains running in different directions. 7 FIG. 4: (Color online) (a) Illustration of collinear E –symmetry ne- matic state at h = 0 , showing loop–like coordination of parallel spins associated with the “ice” manifold. (b) Four–sublattice lo ng range order with E symmetry induced by FM J$_{3}$ . (c) Canted nematic state with partial magnetization under applied field. of the tetrahedral symmetry group T $_{d}$, defined in Sec. II B, and that local fluctuations of order λ local E take on their maximal value λ local E = 16 / 3 . This concludes our survey of symmetry breaking for T = 0 , but it leaves open the question, what happens at finite tem- perature? By analogy with ordered systems where “order from disorder” is effective, thermal fluctuations might be e x- pected to select a single configuration from the “ice” mani- fold, and so restore N´eel order. To address this question, we have performed extensive Monte Carlo simulations using a local–update Metropolis algorithm to sample spin configura - tions. We typically perform 10 6 MC samplings for measure- ments after 10 5 steps for thermalization. We have checked the convergence of the results by comparing those for different initial spin configurations. In particular, to minimize thehys- teresis associated with first order transitions, we used mix ed initial conditions in which different parts of the system are assigned different ordered or disordered states$^{32}$. Where the acceptance rate in MC updates becomes extremely slow, we used the exchange MC method$^{33}$, to avoid local spin-freezing at low temperatures. Results are divided into five bins to es- timate statistical errors by variance of average values in the bins. The system sizes in the present work are up to L = 16 , where L is the linear dimension of the system measured in the cubic units shown in Fig. 1, i.e., the total number of spins N is given by 16 L $^{3}$. We show the results for b = 0 . 1 and b = 0 . 6 , case by case, both of which exhibit qualitatively the same behavior. In the absence of applied magnetic field, the spin collinear- ity P ( r ) grows sharply below a transition temperature T$_{Q}$ ≈ b , as illustrated in Fig. 5. As expected, for T → 0 , P ( r ) → 1 , implying that all spins have a single common axis. The nematic order parameter Q [Eq. (18)] is plotted in Fig. 6(a), together with the heat capacity in Fig. 6(b) for a range of system sizes from L = 4 to L = 16 . Here the heat capacity is calculated by the fluctuation of internal energyas dictated by this effective electrodynamics. This power–la w decay of spin correlations is a signal property of the “ice” manifold. However it should not be taken to imply that “up” and “down” spins are entirely uncorrelated. Within each uudd tetrahedron, each spin has twice as many AF aligned neigh- bors as FM aligned ones, and the net correlation on nearest neighbor bonds is Locally, order is well formed. More formally, we can state that these uudd tetrahedra belong to the two–dimensional E irrep C$_{v}$ = 〈H $^{2}$〉 − 〈H〉2 T $^{2}$N . (31) 〈 S$_{i}$ . · S$_{j}$ 〉 n . n . = − 1 3 . (30)