Table Of Content

Microscopic Aspects of Magnetic Lattice Demagnetizing Factors M. Twengström,1 L. Bovo,2 M. J. P. Gingras,3,4,5 S. T. Bramwell,2 and P. Henelius1 1Department of Physics, Royal Institute of Technology, SE-106 91 Stockholm, Sweden 2London Centre for Nanotechnology and Department of Physics and Astronomy, University College London, 17-19 Gordon Street, London, WC1H OAH, U.K. 3Department of Physics and Astronomy, University of Waterloo, Waterloo, Ontario, N2L 3G1, Canada 4Canadian Institute for Advanced Research, 180 Dundas St. W., Toronto, Ontario, M5G 1Z8, Canada 5Perimeter Institute for Theoretical Physics, 31 Caroline St. N., Waterloo, Ontario, N2L 2Y5, Canada Consideringlocalizedmagneticmomentsonalattice,weshowthatthedemagnetizingfactorfora non-ellipsoidalsampledependsonthespindimensionality(Ising,XYorHeisenberg)andorientation, 7 as well as the sample shape and susceptibility. Our result indicates that the common practice of 1 approximating the true sample shape with ellipsoids for systems with long-range interactions may, 0 in certain cases, overlook important effects stemming from the microscopic aspects of the system. 2 n a Long-range interactions are fundamentally important ThedeterminationofN isafundamentalproblemand J in many areas of science, from cosmology, through the dates back to the work of Poisson and Maxwell [12]. In 6 gravitational interaction, to biology, through Coulomb’s the1940s,Osborn[2]andStoner[3]tabulatedN forgen- 2 law. A long-range interaction may be defined in d spa- eral ellipsoids, while more recently, Aharoni [13] treated tial dimensions by its two-body potential V(r) scaling cuboids in the χ 0 limit. These highly cited papers ] int → i with distance r as r−α where α d [1]. The paramount bearwitnesstotheimportanceofaccuratelycomputable c ≤ s probleminsuchsystemsishowtointegrateV(r)overan and easily accessible demagnetizing factors. Given that - extendedsystem. FollowingNewtonandEuler,theanal- i) it was realized already in the 1920s that N for a non- l r ysis of general systems has been largely based on the ex- ellipsoidal sample is a function not only of the sample t m act solutions for spheres and ellipsoids [2–6]. This raises shape,butalsoofχ itself[14,15]andthatii)manyex- int . the question of whether approximating other shapes to periments are routinely performed not on ellipsoids, but t a ellipsoids just neglects uninteresting details or whether on cuboids [16, 17], it is perhaps remarkable that it was m there are crucial properties that are lost in the approxi- only very recently that the χ-dependence of N was cal- - mation. Thedemagnetizingprobleminmagneticsystems culated for cuboids away from the χ 0 limit [18, 19]. int d → is a natural setting for exploring this question since it Theexistenceofdemagnetizingfactorsforcuboidssug- n o is accessible and of intrinsic importance in experiments, gests that their thermodynamics may be formulated in c andconstitutesaparagonfortheexplorationofthether- terms of an internal field, with corrections that become [ modynamics of systems with long-range interactions [1]. dependent on both shape and temperature [18] (through 1 Demagnetizing factors are also crucial in the experimen- χint). In this work we have found that, for magnetic lat- v tal investigation of superconducting systems. Analogues tices, the demagnetizing factor of cuboids depends also 8 occurinelectricalsystems[7](thedepolarizationfactor), on the local symmetry and allowed orientations of the 4 in the problem of strain fields around inclusions [8], and magneticmoments. Withreferencetothequestionposed 6 in the treatment of avalanching systems in confined ge- at the very beginning, our result illustrates a case where 7 0 ometries [9–11], just to name a few. a long-range interaction integrates in a qualitatively dif- . In an applied magnetic field, the thermodynamic en- ferent way for a cuboid and an ellipsoid, such that the 1 0 ergy of an ellipsoidal magnetic sample acquires a mag- discrete microscopic nature of the system matters in the 7 netostatic contribution E =(µ /2)Nm2/V, where N former case but not in the latter. We are aware of only mag 0 1 is the demagnetizing factor, V is the volume and m the a few previous studies where effects of such discreteness v: magneticmoment. AftersubtractingEmagfromthetotal have been discussed [7, 20–22]. In particular, Millev i energy, differentiation of the corrected energy amounts et al. [22] found a multiplicative separation of discrete X to an effective reduction of the applied field, B , by andcontinuumcontributionstothedemagnetizingeffects ext r µ NM, where M m/V is the magnetization. Viewed for systems with restricted geometries. In our work, we a 0 ≡ asacorrectiontoH B /µ ,thisisthedemagnetiz- demonstrate discreteness effects even in the 3d thermo- ext ext 0 ≡ ing field H = NM within the sample. With an inter- dynamic limit. Our interest in this problem was partly d − nalfieldH H +H ,onefindsanintrinsicmagnetic spurred by the perplexing recent experimental observa- int ext d ≡ susceptibilityχ =∂M/∂H ,ashape-independentma- tion of anomalous demagnetizing effects in the spin ice int int terialproperty,whichrelatestotheexperimentallydeter- material Dy Ti O [23]. 2 2 7 mined susceptibility, χ ∂M/∂H , through Onemayaskwhethersmalldifferencesinestimatedde- exp ext ≡ magnetizingfactorsreallymatterforexposingimportant 1 1 and interesting physics. The answer is found in Eq. (1). = N. (1) χint χexp − If χexp, and N are both small, then this transforma- 2 tion does not renormalize the sought χ appreciably, direction at site i, and m is dimensionless. We first int i being rather insensitive to the precise value of N. How- determinethecomponentofthelocalfieldalongtheIsing (cid:113) ever, in many physical systems displaying unusual mag- moment at site i, B = B ˆı, which is the sum of three i i· netic phenomena, χ can be large, and χ is, through contributions: exp int Eq. (1), a highly sensitive function of N. Examples include the rare-earth based spin ice materials Dy Ti O B(cid:113) =B(cid:113),dip+B(cid:113),ext+B(cid:113),self, (2) 2 2 7 i i i i and Ho Ti O , which support magnetic monopole ex- 2 2 7 which we discuss one by one. citations [24], and LiHo Y F which displays ultra- 1−x x 4 First, the dipolar field at site i produced by all the slow relaxational behavior [25]. Important demagnetiz- otherpointmagneticdipoleswithinthesample, B(cid:113),dip ingeffectsaremanifestinexperimentsrequiringanaccu- i ≡ ratelydirectedexternalfield, asinthecaseoftheelusive Bidip·ˆı, is given by the familiar form [32] Kasteleyntransition[26], sub-latticepinning[17,27–29], (cid:32) (cid:33) or multiple field-driven transitions [30], or disentangling B(cid:113),dip = µ0µB (cid:88) 3(ˆ·rîj)(ˆı·rîj)−ˆ·ˆı m . (3) the in-phase and out-of-phase response for non-zero fre- i 4π r3 j j(cid:54)=i ij quency stimuli [31]. In order to reach quantitative con- clusionsfromthededucedχ ,ortoaccuratelytesttheo- Second,weconsideranexternalfieldinthezˆdirection, int retical predictions, it is necessary to determine N highly Bext = Bextzˆ, with B(cid:113),ext = Bextcosθ , where cosθ i i i i ≡ accurately. Our work illustrates that one may indeed zˆ ˆı, the angle between the direction of the Ising axis at · achieve this (see Fig 1). site i and the direction of Bext, zˆ. i Third,isthecontributionfromtheself-field,B(cid:113),self. In i the classic case of a single point dipole [32, 33], a term 0.33 23µ0µBδ(r) must be added to ensure that the average magneticfieldinaspherecontainingthedipolegivesthe correct macroscopic field. Similarly, we add a self-field to ensure that the internal magnetic field in a uniformly 0.32 magnetized sample results in the expected value, for ex- N ample B=23µ0M for a uniformly magnetized sphere or cube[19]. Weareprimarilyconcernedwithparamagnetic 0.31 samplesinthelinearresponseregime,whereaweakmag- netic field induces a magnetization proportional to it, as in a typical χ measurement. For a non-ellipsoidal sam- 0.3 ple, the induced magnetization will in general be non- uniform, except in the χ 0 limit. In this limit, H d → vanishesand,asaconsequence,H andM areuniform. 0 10 20 30 int T [K] Our goal is to determine the self-field so that the magnetic field has the expected value in the χ 0 limit. We → demonstrate the basic idea with two examples. FIG. 1. N as a function of T for a cubic sample of the spin ice Dy Ti O . Shown are experimental results (blue open We first take a cubic sample with all moments aligned 2 2 7 circles) and our fitting parameter-free theory (red line). The intheglobalzˆdirection. InthiscaseB,M andH areall dashed reference line shows the exact N =1/3 value for a aligned with the zˆdirection for which the field equation sphere [2] B = µ (M +H) reduces to Bz = µ (Mz N Mz) = 0 0 0 2µ Mz, where N = 1 is the χ 0 limit−of N for a 3 0 0 3 → Determination of N via an iterative method In the cube [19]. If we consider a simple cubic lattice, it is well − spirit of problems (6.18) and (4.23) in Griffiths [32], we known that the lattice sum vanishes [34]. This implies use an iterative method to calculate the on-site field dis- that Bz,self=2µ Mz must be incorporated to ensure the 3 0 tribution inside a linear magnetic material placed in a expected net field value. uniform magnetic field. Initially, we assume that H Next, we consider the case of a lattice where the Ising int equals H and calculate the induced local magnetiza- axes are tilted by an angle θ with respect to the z- ext tion for an assumed χ . This magnetization generates axis, and half the spins tilt to the right and half to int a demagnetizing field that, in turn, modifies H . The the left so that there is no net magnetization in the int resulting field-magnetization equations are iterated until xˆ or yˆ directions. The total B, M and H fields are convergence. Withtheconvergedfieldandmagnetization in the zˆ-direction, but the question is what should the distributions in hand, one then computes N. B(cid:113) field parallel to the magnetic moments be? From Toproceed, weconsiderasampleofvolumeV with B=µ (M +H), it follows that B is generated by two 0 N magnetic moments. Initially, we focus on Ising moments terms, which we discuss separately. We begin with the m =m µ ˆı, whereˆı is the unit vector in the local Ising term generated directly by M, namely B1 =µ M, or i i B 0 3 B1,z=µ Mz=µ M(cid:113)cosθ, where M(cid:113) is the magnetiza- (MC) simulations for a few select physical cases to ver- 0 0 tion in the local Ising directions, M(cid:113)=V−1(cid:80)N m ˆı. ify the results of the iterative method. For a single data This equation is satisfied by B1,(cid:113)=µ M(cid:113). Thi=e1secio·nd point, the MC approach requires (105) core hours [36] 0 term, B2,z=µ Hz= µ N Mz= µ N M(cid:113)cosθ is gen- to reach the necessary precisionOfor (104) moments. 0 0 0 0 0 − − O erated by H . The field along the magnetic moment is With the MC method, we cannot set χ as done in the d loc B2,(cid:113)= µ N M(cid:113)cos2θ, and the net self-field becomes iterativemethod. Instead,weadjustthetemperature,T, 0 0 − in order to tune χ to the desired value. int B(cid:113),self =B1,(cid:113)+B2,(cid:113) =µ µ N (cid:2)1 N cos2θ(cid:3)m , (4) For definitiveness, we use the magnetostatic dipolar i 0 BV − 0 i Hamiltonian H = µ0µB (cid:80) Λ σ σ , where σ = 1 4π i>j ij i j i ± and Λ = [(ˆı ˆ) 3(ˆı rˆ )(ˆ rˆ )]/r3, and χzz is de- whichisvalidwhenthedipolarlatticesum,Eq.(3),van- ij · − · ij · ij ij ishes and when the average M is along Bext. In case termined according to the lattice sum does not vanish, it must be subtracted from the self-field in order to ensure the expected net χzz = ∂Mz = µ0µB (cid:42)(cid:32)(cid:88)N σ cosθ(cid:33)2(cid:43). (8) field value. ∂Hz k TV i B Eqs. (2, 3, 4) give the local field in terms of the set of i=1 local magnetizations, m . Next we “invert” this to get i Using Ewald boundary conditions we obtain χMC, while the m induced by B{(cid:113). }Using M=χH (linear media) int we g{etiB} =µ (M +Mi/χ)=µ χ+1M, or M= χ B, open boundary conditions yield χMexCp, and N is obtained 0 0 χ χ+1µ0 from Eq. (1). leading to (cid:18) (cid:19) (cid:113) 0.34 V χ B m = loc i , (5) i χ +1 µ µ loc 0 B 0.33 N where χ is the local susceptibility in theˆı direction. loc 0.32 (cid:113) We can now iterate the expressions for B in Eq. (2) i and m in Eq (5) until convergence, and then calculate i 0.31 N from Eq. (1), where χ is given by exp N SphereN=1/3 0.3 sccosθ=1iterative χexp =χzexzp =(cid:18)∂∂HMzz (cid:19) = VµB0µeBxt (cid:88)N micosθ. (6) 0.29 sbbccccccoccsooθssθθ===111M/M√CC2iterative ext T i=1 bcccosθ=1/√2MC pyrochlorecosθ=1/√3iterative The intrinsic susceptibility, χ , expresses the relation 0.28 pyrochlorecosθ=1/√3MC int Chen,Ref.[19] between Bext and induced M under so-called “Ewald”, Heisenberg,iterative 0.27 or “tin foil”, boundary condition [35] which eliminates 10 2 10 1 100 101 102 − − demagnetizing fields and corresponds to the N=0 limit. χ int Asaresult,bothχ andχ areresponsestoaninternal int loc field. Whileχ measurestheresponseinthedirectionof FIG. 2. N as a function of χ for cubic samples of various int int H ,χ measurestheresponsealongthelocalIsingaxis lattices(seemaintext). Thelinesaretheresultsoftheitera- ext loc ˆı. With zˆ ˆı = cosθ, Hz cosθ induces a magnetization tivecalculation,whilesymbolsareMCcheckpoints. Theblue M(cid:113) = χ ·Hz cosθ. Texhtis magnetization, in turn, has squares are from Chen et al. [19]. Maroon circles indicate an loc ext iterative calculation for spherically symmetric (Heisenberg) a component Mz = M(cid:113)cosθ = χ Hz cos2θ along zˆ, loc ext spins on an sc lattice. and therefore χ =χ cos2θ. int loc (cid:113) Tosumup,oncetheconvergedB andm distributions i i Experimental determination of N We determined N have been determined, N is calculated using Eq. (1), − for a cube of Dy Ti O spin ice using N = 1/χcs 2 2 7 exp − 1/χss +1/3, where χcs and χss are the experimen- (cid:34) N (cid:35)−1 exp exp exp N = µ0µB (cid:88)m cosθ 1 . (7) tally measured susceptibilities of a cubic and a spherical VBext i − χ cos2θ sample respectively. In order to match the level of the loc i=1 susceptibility of the cube and sphere in the high-T limit, Determination of N via Monte Carlo simulations χcs was shifted by about 1% (χcs χcs /1.0074) be- − exp exp → exp With the iterative method we reach system sizes of fore calculating N. (106) spins. However, it is not obvious that this The two samples were obtained from different larger O method, which is mean-field like and does not include crystals provided by D. Prabhakaran [37]. The spheri- fluctuations in the m ’s, would give the same result as a cal sample of diameter 4 mm was made by commercial i full statistical calculation for a given spin Hamiltonian. hand-cutting techniques similar to the ones employed in Therefore, we have also calculated N using Monte Carlo Ref. [23]. The cubic sample, of dimensions 2 2 2 × × 4 mm3,andwithalledgespreciselyorientedwiththecrys- tion[42,43]. TheresultsinFig.2indicatethatN clearly tallographic axes a, b, c, was cut and epi-polished on all depends on the symmetry and direction of the moments, sides [38]. We carefully controlled the crystal shape and and that the MC and iterative calculations agree. The orientationaswellasexperimentalsetuptominimizeall sc,bccandLiHoF latticeswithIsingspinsalignedalong 4 potential sources of uncertainty and systematic errors in B yield the same N, indicating that it is not directly ext the measurement. sensitive to the specific symmetry of the lattice itself. Results Our main experimental result is shown in TurningthelocalIsingaxesawayfromBextcausesamore − Fig. 1, where we display the experimentally determined rapiddecreaseofN withincreasingχint. Thepyrochlore N asafunctionofT forthecubicsinglecrystalsampleof lattice with tilt angle cosθ = 1/√3 yields a smaller N Dy2Ti2O7 spin ice. The main theoretical result is shown than the bcc lattice with cosθ = 1/√2 for χint > 0. It in Fig. 2, where N is displayed as a function of χ for is interesting to note that the spin ice pyrochlore lat- int cubic samples. We have considered the simple cubic (sc) tice yields the same result as the continuum method of and body centered cubic (bcc) lattices with the Ising di- Chen et al. [19], and the dipolar model with Heisenberg rection parallel to Bext. We also display results for two spins. It is therefore tempting to conjecture that mod- lattices with Ising spins tilted away from Bext. First, a els for which χ is independent of the direction of Bext bcc lattice, with spins pointing in the (101) and ( 101) will follow the behavior calculated by Chen et al. The L − directions. Second, a pyrochlore lattice built from the dependence in Fig. 3 shows that the iterative method is conventional cubic unit cell [39]. Finally, we include re- asymptotically exact in the limit of large L. sults for the dipolar model with spherically symmetric What are the lessons of experimental relevance per- Heisenberg spins on an sc lattice. Results (not shown) taining to the determination of N? If accurate measure- foratetragonalIsinglattice,relevanttoLiHoF [40],are ments of χ is the objective, then the corrections to 4 int found to be identical to the sc case. N(χ 0) identified here can become very important for χ (cid:38)→N. We emphasize that this is not simply a minor 0.33 quantitativecorrection, buttheeffectsmaybedramatic. For example, in the case of Dy Ti O , Tχ (T) features 2 2 7 int a peak, which is easily shifted outside the experimental temperature window by application of the ordinary 0.32 χ = 0) demagnetizing correction [23]. More generally, while the demagnetizing correction can in principal be controlled for highly elongated samples or ellipsoids, it N0.31 is not always easy to prepare real samples that approx- imate these idealizations. This is particularly true of non-metallic (and perhaps brittle) samples, which have 0.3 becomeofmuchgreaterinterestinrecentyears–spinice [44] and LiHoF [40] being well-known examples. There- 4 fore, insofarascuboidsareoftenthemostpracticalsam- ple shapes to accurately prepare, the best approach may 0.29 0 0.05 0.1 0.15 0.2 0.25 0.3 betousecuboidalsamplesandapplythetheoreticalcor- 1/L rections identified here. In conclusion, considering the demagnetizing problem FIG. 3. N as a function of inverse linear system size 1/L for as a paradigm for the study of effects caused by long- acubicsampleofthesclattice. Shownareχ =1.00(black), int range interactions, our results confirm that N may be 1.82 (red) and 7.53 (blue) for the iterative method (circles), precisely defined for cuboidal samples. This implies that and Monte Carlo method (squares). the free energy gains a non-additive part of the form F = (µ /2)VN(T)M2 in the approach to the ther- mag 0 All results are extrapolated to infinite system size and modynamic limit [45] and that the spatially averaged in- in Fig. 3 we compare the system-size dependence of the ternal field and magnetisation are conjugate thermody- twomethods. Extrapolationshavebeenperformedusing namic variables. By going beyond Maxwell’s continuum thefunctionalforma+b/L+c/L2fortheiterativedemag- theory, we show that N depends not only on the sample netizingfactorandχMexCp,whileweuseda+b/L3+c/L4for shape and χ,but also on microscopic factors: specifically χMC, where L is the linear system size of the cube [41]. the spin dimensionality and the effective magnetic sym- int Discussion The experimental result in Fig. 1 veri- metryofthecrystallographicbasis. Giventhatdetailsof − fies the explicit T dependence of N for a cuboid. Fur- the lattice and local environment may affect even such thermore, the theoretical prediction matches the experi- a fundamental and well-studied macroscopic property as mental data well, even though the short-range exchange N,itisinterestingtoaskwhethermicroscopicproperties interaction in Dy Ti O was not included in the simula- could be at play in other systems with long-range inter- 2 2 7 5 actions. Afullunderstandingofhowthethermodynamic [19] D.-X. Chen, E. Pardo, and A. Sanchez, IEEE Trans. limitisreachedforlong-rangedandconditionallyconver- Magn. 41, 2077 (2005). gent interactions continues to offer challenging problems [20] P.H.ChristensenandS.Mørup,J.Magn.Magn.Mater. 35, 130 (1983). to theorists and experimentalists alike. [21] E. Y. Vedmedenko, H. P. Oepen, and J. Kirschner, J. We thank D. Prabhakaran for providing crystals from Magn. Magn. Mater. 256, 237 (2003). which the samples were cut, and J. Rau for useful dis- [22] Y.T.Millev,E.Vedmedenko, andH.P.Oepen,J.Phys. cussions. The simulations were performed on resources D 36, 2945 (2003). provided by the Swedish National Infrastructure for [23] L. Bovo, L. D. C. Jaubert, P. C. W. Holdsworth, and Computing (SNIC) at the Center for High Performance S. T. Bramwell, J. Phys.: Condens. Matter 25, 386002 Computing (PDC) at the Royal Institute of Technology (2013). (KTH). MT and PH are supported by the Swedish Re- [24] C. Castelnovo, R. Moessner, and S. L. Sondhi, Nature 451, 42 (2008). search Council (2013-03968), MT is grateful for funding [25] A. Biltmo and P. Henelius, Nat. Comm. 3, 857 (2012). from Stiftelsen Olle Engkvist Byggmästare (2014/807), [26] T.Fennell,S.T.Bramwell,D.F.McMorrow,P.Manuel, LB is supported by The Leverhulme Trust through the and A. R. Wildes, Nat. Phys. 3, 566 (2007). Early Career Fellowship program (ECF2014-284). This [27] J. P.C. Ruff, R. G.Melko, and M. J. P. Gingras, Phys. researchwassupportedinpartbythePerimeterInstitute Rev. Lett. 95, 097202 (2005). forTheoreticalPhysics. ResearchatthePerimeterInsti- [28] H. Sato, K. Matsuhira, T. Tayama, Z. Hiroi, S. Takagi, tute issupported bythe Governmentof Canadathrough and T. Sakakibara, Journal of Physics: Condensed Mat- ter 18, L297 (2006). Innovation, Science and Economic Development Canada [29] H.Sato,K.Matsuhira,T.Sakakibara,T.Tayama,Z.Hi- and by the Province of Ontario through the Ministry of roi, and S. Takagi, Journal of Physics: Condensed Mat- Research, InnovationandScience. TheworkattheU.of ter 19, 145272 (2007). Waterloo was supported by the Canada Research Chair [30] P. Schiffer, A. P. Ramirez, D. A. Huse, and A. J. program (M.G., Tier 1). Valentino, Phys. Rev. Lett. 73, 2500 (1994). [31] J. A. Quilliam, L. R. Yaraskavitch, H. A. Dabkowska, B.D.Gaulin, andJ.B.Kycia,Phys.Rev.B83,094424 (2011). [32] D. J. Griffiths, Introduction to Electrodynamics, 3rd ed. [1] A. Campa, T. Dauxois, and S. Ruffo, Phys. Rep. 480, (Prentice-Hall, 1999). 57 (2009). [33] D. J. Griffiths, American Journal of Physics 50, 698 [2] J. A. Osborn, Phys. Rev. 67, 351 (1945). (1982). [3] E. C. Stoner, Philos. Mag. 36, 803 (1945). [34] J. D. Jackson, Classical Electrodynamics, 2nd ed. (John [4] M. Beleggia, M. De Graef, and Y. T. Millev, Philos. Wiley & Sons, 1975). Mag. 86, 2451 (2006). [35] P. P. Ewald, Ann. Phys. 369, 253 (1921). [5] M.Beleggia,M.DeGraef, andY.T.Millev,J.Phys.D [36] Intel Xeon E5-2698v3 processors with Cray Aries inter- 39, 891 (2006). connect between nodes. [6] G. Di Fratta, Proc. R. Soc. A 472, 0197 (2016). [37] D. Prabhakaran and A. Boothroyd, Crys. Growth 318, [7] I. Ponomareva, I. I. Naumov, I. Kornev, H. Fu, and 1053 (2011). L. Bellaiche, Phys. Rev. B 72, 140102 (2005). [38] SurfaceNet GmbH, http://www.surfacenet.de. [8] J. D. Eshelby, Proc. R. Soc. A 241, 376 (1957). [39] R. G. Melko and M. J. P. Gingras, J. Phys.: Condens. [9] G.DurinandS.Zapperi,J.Appl.Phys.87,7031(2000). Matter 16, R1277 (2004). [10] F. F. Csikor, C. Motz, D. Weygand, M. Zaiser, and [40] M. J. P. Gingras and P. Henelius, J. Phys.: Conf. Ser. S. Zapperi, Science 318, 251 (2007). 320, 012001 (2011). [11] L. Laurson, X. Illa, S. Santucci, K.-T. Tallakstad, K. J. [41] M. E. J. Newman and G. T. Barkema, Monte Carlo Måløy, and M. J. Alava, Nat. Comm. 4, 2927 (2013). Methods in Statistical Physics, 1st ed. (Oxford Univer- [12] D.-X. Chen, J. A. Brug, and R. B. Goldfarb, IEEE sity Press, 1999). Trans. Magn. 27, 3601 (1991). [42] T. Yavors’kii, T. Fennell, M. J. P. Gingras, and S. T. [13] A. Aharoni, J. Appl. Phys. 83, 3432 (1998). Bramwell, Phys. Rev. Lett. 101, 037204 (2008). [14] J. Würschmidt, Z. Phys. 12, 128 (1923). [43] P. Henelius, T. Lin, M. Enjalran, Z. Hao, J. G. Rau, [15] F.StäbleinandH.Schlechtweg,Z.Phys.95,630(1935). J. Altosaar, F. Flicker, T. Yavors’kii, and M. J. P. Gin- [16] J. A. Quilliam, S. Meng, C. G. A. Mugford, and J. B. gras, Phys. Rev. B 93, 024402 (2016). Kycia, Phys. Rev. Lett. 101, 187204 (2008). [44] M. J. P. Gingras, “Introduction to frustrated mag- [17] R. Higashinaka and Y. Maeno, Phys. Rev. Lett. 95, netism,” (Springer(Ed.Lacroix,Mendels&Mila),2011) 237208 (2005). Chap. Spin Ice. [18] D.-X. Chen, E. Pardo, and A. Sanchez, IEEE Trans. [45] P. Rhodes and G. Rowlands, Proc. Leeds Phil. Lit. Soc. Magn. 38, 1742 (2002). 6, 191 (1954).