Phase boundary near a magnetic percolation transition
EEuropean Physical Journal B manuscript No. (will be inserted by the editor)
Phase boundary near a magnetic percolation transition
Gaurav Khairnar a,1 , Cameron Lerch b,1,2 , Thomas Vojta c,1 Department of Physics, Missouri University of Science and Technology, Rolla, Missouri 65409, USA Department of Mechanical Engineering and Materials Science, Yale University, New Haven, Connecticut 06520, USAReceived: date / Accepted: date
Abstract
Motivated by recent experimental observa-tions [Phys. Rev. , 020407 (2017)] on hexagonal fer-rites, we revisit the phase diagrams of diluted magnetsclose to the lattice percolation threshold. We performlarge-scale Monte Carlo simulations of XY and Heisen-berg models on both simple cubic lattices and latticesrepresenting the crystal structure of the hexagonal fer-rites. Close to the percolation threshold p c , we find thatthe magnetic ordering temperature T c depends on thedilution p via the power law T c ∼ | p − p c | φ with expo-nent φ = 1 .
09, in agreement with classical percolationtheory. However, this asymptotic critical region is verynarrow, | p − p c | (cid:46) .
04. Outside of it, the shape of thephase boundary is well described, over a wide rangeof dilutions, by a nonuniversal power law with an ex-ponent somewhat below unity. Nonetheless, the perco-lation scenario does not reproduce the experimentallyobserved relation T c ∼ ( x c − x ) / in PbFe − x Ga x O .We discuss the generality of our findings as well as im-plications for the physics of diluted hexagonal ferrites. Disordered many-body systems feature three differenttypes of fluctuations, viz., static random fluctuationsdue to the quenched disorder, thermal fluctuations, andquantum fluctuations. Their interplay can greatly affectthe properties of phase transitions, with possible con-sequences ranging from a simple change of universalityclass [1] to exotic infinite-randomness criticality [2, 3],classical [4] and quantum [5, 6] Griffiths singularities,as well as the destruction of the transition by smearing a e-mail: [email protected] b e-mail: [email protected] c e-mail: [email protected] [7–10]. Recent reviews of some of these phenomena canbe found in Refs. [11–13]. Randomly diluted magneticmaterials are a particularly interesting class of systemsin which the above interplay is realized. Here, the dis-order fluctuations correspond to the geometric fluctua-tions of the underlying lattices which can undergo a ge-ometric percolation transition between a disconnectedphase and a connected (percolating) phase [14].Recently, the behavior of diluted magnets close tothe percolation transition has reattracted attention be-cause of the unexpected shape of the phase bound-ary observed in the diluted hexagonal ferrite (hexafer-rite) PbFe − x Ga x O [15]. Pure PbFe O orders fer-rimagnetically at temperatures below about 720 K [16].The ordering temperature T c can be suppressed by ran-domly substituting nonmagnetic Ga ions for Fe ions inPbFe − x Ga x O . It vanishes when x reaches the crit-ical value x c ≈ .
6. This value is very close the perco-lation threshold x p = 8 .
846 of the underlying lattice ,suggesting that the transition at x c is of percolationtype [15]. Remarkably, the phase boundary follows thepower law T c ( x ) = T c (0)(1 − x/x c ) φ with φ = 2 / x -range from 0 to x c . This disagrees withthe prediction from classical percolation theory [14, 17]which yields a crossover exponent of φ > x c . In this paper, we therefore reinvestigate the phaseboundary close to the percolation transition of dilutedclassical planar and Heisenberg magnets by means oflarge-scale Monte Carlo simulations. The purpose ofthe paper is twofold. First, we wish to test and verifythe percolation theory predictions, focusing not only on The lattice in question is the lattice of exchange interactionsbetween the Fe ions. a r X i v : . [ c ond - m a t . d i s - nn ] N ov the asymptotic critical behavior but also on the widthof the critical region and the preasymptotic properties.Second, we wish to explore whether the classical per-colation scenario can explain the experimental observa-tions in PbFe − x Ga x O [15].Our paper is organized as follows. In Sec. 2, we in-troduce the diluted XY and Heisenberg models and dis-cuss their qualitative behavior. Section 3 summarizesthe predictions of percolation theory. Our Monte Carlosimulation method is described in Sec. 4. Sections 5.1and 5.2 report our results for model systems on cubiclattices and for systems defined on the hexagonal ferritelattice, respectively. We conclude in Sec. 6. Consistent with the dual purpose of studying the crit-ical behavior of the phase boundary close to a mag-netic percolation transition and of addressing the ex-perimental observations in diluted hexaferrites [15], weconsider two models, viz., (i) site-diluted classical XYand Heisenberg models on simple cubic lattices and (ii)a classical Heisenberg Hamiltonian based on the hexa-ferrite crystal structure using realistic exchange inter-actions. Comparing the results of these different modelswill also allow us to explore the universality of the crit-ical behavior.2.1 Site-diluted XY and Heisenberg models on cubiclatticesWe consider a simple cubic lattice of N = L sites. Eachsite is either occupied by a vacancy or by a classicalspin, i.e., an n -component unit vector S i ( n = 2 forthe XY model and n = 3 for the Heisenberg case). TheHamiltonian reads H = − J (cid:88) (cid:15) i (cid:15) j S i · S j . (1)Here, the sum is over pairs of nearest-neighbor sites,and J > J to unity for thecubic lattice simulations.) The quenched independentrandom variables (cid:15) i implement the site dilution. Theytake the values 0 (vacancy) with probability p and 1(occupied site) with probability 1 − p . We employ pe-riodic boundary conditions. Magnetic long-range ordercan be characterized by the order parameter, the totalmagnetization m = 1 N (cid:88) i S i . (2) Fig. 1
Double unit cell of PbFe O . 24 Fe ions are lo-cated on five distinct sublattices. The qualitative behavior of this model as a functionof temperature T and dilution p is well understood (see,e.g., Ref. [18] for an overview). For sufficiently small di-lution, the system orders magnetically below a criticaltemperature T c ( p ). The critical temperature decreasescontinuously with p until it reaches zero at the perco-lation threshold p c of the lattice. For dilutions beyondthe percolation threshold, magnetic long-range order isimpossible because the system breaks down into finitenoninteracting clusters. The point p = p c , T = 0 is amulticritical point at which both the geometric fluctua-tions of the lattice and the thermal fluctuations becomelong-ranged.2.2 Hexaferrite Heisenberg modelPbFe O crystallizes in the magnetoplumbite struc-ture, as illustrated in Fig. 1. A double unit cell con-tains 24 Fe ions in five distinct sublattices; they arein the spin state S = 5 /
2. Below a temperature of about720 K , the material orders ferrimagnetically, with 16 ofthe Fe spins pointing up and the remaining 8 Fe ionspointing down [16]. Note that the high critical temper-ature and the high spin value suggest that a classicaldescription should provide a good approximation.In PbFe − x Ga x O , the randomly substituted Gaions, which replace the Fe ions, act as quenched spinlessimpurities. To model this system, we start from thehexaferrite crystal structure and randomly place eithera vacancy (with probability p ) or a classical Heisenbergspin S i (with probability 1 − p ) at each Fe site. Thedilution p is related to the number x of Ga ions in the unit cell by p = x/
12. The Hamiltonian reads H = − (cid:88) i,j J ij (cid:15) i (cid:15) j S i S j . (3)The quenched random variables (cid:15) i distinguish vacan-cies and spins, as before. The values of the exchangeinteractions J ij stem from the density functional cal-culation in Ref. [19]; they are scaled by a common fac-tor to approximately reproduce the critical temperature T c = 720 K of the undiluted material. In most of ourMonte Carlo simulations, we include only the leading(strongest) interactions which are between the followingsublattice pairs: 2a-4f IV , 2b-4f V I , 12k-4f IV , 12k-4f V I .These interactions are non-frustrated and establish theferrimagnetic order. We also perform a few test calcula-tions to explore the effects of additional couplings whichare significantly weaker but frustrate the ferrimagneticorder.The qualitative features of the phase diagram of themodel (3) are expected to be similar to those discussedin the previous section. With increasing dilution p , thecritical temperature T c ( p ) is continuously suppressedand reaches zero at the site percolation threshold. Thevalue of the percolation threshold of the lattice spannedby the leading non-frustrated interactions between theFe ions was determined in Ref. [15] by means of MonteCarlo simulations. They yielded p c = 0 . x c = 8 . In this section, we briefly summarize the predictions ofclassical percolation theory for the shape of the phaseboundary T c ( p ) close to multicritical point p = p c , T =0 [14, 17, 20]. Close to this point, two length scalesare at play, the percolation correlation length, ξ p whichcharacterizes the size of finite isolated clusters of lat-tice sites and the magnetic thermal correlation lengthon the critical infinite percolating cluster at p c denotedby ξ T . The percolation correlation length ξ p divergesas ξ p ∼ | p − p c | − ν p as the percolation threshold is ap-proached. The magnetic thermal correlation length be-haves as ξ T ∼ T − ν T for continuous-symmetry magnetsdescribed by the n -vector model with n > m ( p − p c , T ) = | p − p c | β X ( ξ T /ξ p ) . (4)For p < p c , the magnetic phase transition occurs ata particular value x c of the argument of the scaling function X . At the magnetic transition, we thereforehave ξ T = x c ξ p . This yields the power law relation T c ( p ) ∼ | p − p c | φ . (5)The crossover exponent φ takes the value φ = ν p /ν T .(In contrast, ξ T diverges exponentially, ξ T ∼ ( e − J/T ) − ν T ,for Ising magnets, leading to a logarithmic dependence T c ( p ) ∼ ln − (1 / | p − p c | ).)Using a renormalization group calculation, Coniglio[17] established the relation ν T = 1 / ˜ ζ R . Here, ˜ ζ R char-acterizes the resistance R of a random resistor networkon a critical percolation cluster of linear size L via R ∼ L ˜ ζ R .The exponent ˜ ζ R can be related to the well-knownconductivity critical exponent t which describes howthe conductivity σ of the resistor network depends onthe distance from the percolation threshold, σ ∼ | p − p c | t . To do so, consider a resistor network on a perco-lating lattice close to p c but on the percolating side. Itsbehavior is critical for clusters of size less than ξ p andOhmic for sizes beyond ξ p . For a d -dimensional systemof linear size L (cid:29) ξ p , we can employ Ohm’s law tocombine blocks of size ξ p , yielding R ( L ) = R ( ξ p ) (cid:18) Lξ p (cid:19) (cid:18) Lξ p (cid:19) − ( d − ∼ ξ ˜ ζ R p ξ d − p L − d . (6)The conductivity on the percolating side thus behavesas σ ∼ ξ − ( d − ζ R ) p ∼ | p − p c | ν p ( d −
2+ ˜ ζ R ) . Thus, we obtainthe hyperscaling relation, t = ( d − ζ R ) ν p . The re-sulting crossover exponent reads φ = ν p /ν T = ν p ˜ ζ R = t − ( d − ν p . Using the numerical estimates t/ν p = 2 . ν p = 0 . ζ R = 1 . φ = 1 . Shender and Shklovskii [20] had established an upper bound T >c ( p ) = A | p − p c | t − ν p for the critical temperature in threedimensions, in agreement with this relation. E Hot startCold start0 20 40 60 80 100 number of MC sweeps m Fig. 2
Equilibration of the energy per site E and the mag-netization m for a cubic lattice XY model of size L = 56,dilution p = 0 .
66 , and temperature T = 0 .
156 averaged over20 disorder configurations. The comparison of hot and coldstarts shows that the system equilibrates after roughly 50Monte Carlo sweeps despite being close to the multicriticalpoint. equilibration of small isolated clusters of lattice siteswhich might form as a result of dilution.For the cubic lattice calculations, we consider sys-tem sizes ranging from L = 10 to L = 112 . Wehave simulated 4000 − to 40 double unitcells (each double unit cell contains 24 Fe sites) us-ing 100 −
300 independent disorder configurations foreach size. All physical quantities of interest, such asenergy, magnetization, correlation length, etc. are aver-aged over the disorder configurations. Statistical errorsare obtained from the variations of the results betweenthe configurations.Measurements of observables must be performed af-ter the system reaches thermal equilibrium. We deter-mine the number of Monte Carlo sweeps required forthe system to equilibrate by comparing the results ofruns with hot starts (for which the spins initially pointin random directions) and with cold starts (for whichall spins are initially aligned). An example of such atest for a cubic lattice XY system close to multicrit-ical point is shown in Fig. 2. The energy and orderparameter attain their respective equilibrium values af-ter roughly 50 Monte Carlo sweeps. Similar numericalchecks were performed for other parameter values aswell as for the cases of Heisenberg spins on cubic andhexaferrite lattices. Based on these tests, we have cho-sen 150 equilibration sweeps (using a hot start) and 500measurement sweeps per disorder configuration for thecubic lattice simulations. For the hexaferrite lattice, we perform 1000 equilibration sweeps and 2000 measure-ment sweeps (using a hot start). Note that the combi-nation of relatively short Monte Carlo runs and a largenumber of disorder configurations leads to an overallreduction of statistical error [25–27].4.2 Data analysisWe employ the Binder cumulant [28] to precisely esti-mate the critical temperature T c . It is defined as g = (cid:20) − (cid:104)| m | (cid:105) (cid:104)| m | (cid:105) (cid:21) dis (7)where (cid:104) ... (cid:105) denotes the thermodynamic (Monte Carlo)average and [ ... ] dis denotes the disorder average. TheBinder cumulant g is a dimensionless quantity, it there-fore fulfills the finite-size scaling form g ( t, L ) = F ( tL /ν ) . (8)Here, F is a dimensionless scaling function, t = ( T − T c ) /T c denotes the reduced temperature, and ν is thecorrelation length exponent of the (magnetic) finite-temperature phase transition.At the critical temperature, t = 0, the Binder cu-mulants corresponding to different system sizes havethe universal value f (0), i.e., the critical temperatureis marked by a crossing of all Binder cumulant curves.To determine the crossing temperature, we fit the g vs T data sets corresponding to different system sizeswith separate quartic polynomials.(Quartic polynomi-als provide reasonable fits within the temperature rangeof interest while avoiding spurious oscillations.) The in-tersection point of these polynomials yields the crossingtemperature T ∗ . To estimate the errors of the crossingtemperature we use an ensemble method. For each g ( T )curve, we create an ensemble of artificial data sets g a ( T )by adding noise to the data g a ( T ) = g ( T ) + ∆g ( T ) r . (9)Here, r is a random number chosen from a normal dis-tribution of zero mean and unit variance, and ∆g ( T ) isthe statistical error of the Monte Carlo data for g ( T ).Note that we use the same random number r for theentire g ( T ) curve, leading to an upward or downwardshift of the curve. This stems from the fact that the sta-tistical error ∆g ( T ) is dominated by the disorder noisewhile the Monte Carlo noise is much weaker. This im-plies that the deviations at different temperatures of theBinder cumulant from the true average are correlated.Repeating the crossing analysis with these ensemblesof curves, we get ensembles of crossing temperatures.Their mean and standard deviation yield T ∗ and theassociated error ∆T ∗ , respectively. T g L=20L=28L=40L=56L=80L=112
Fig. 3
Binder cumulant g vs temperature T for the cubiclattice XY model with dilution p = 0 .
10. The statistical errorsarising from the Monte Carlo simulation are smaller than thesymbol size. The inset show the intersection region of thecurves more closely. All curves cross at the same point withintheir statistical errors.
In this section we report the results of our simulationsfor cubic and hexaferrite lattices occupied by XY orHeisenberg spins.5.1 Cubic LatticesWe investigate the behavior of both XY and Heisenbergmodels on cubic lattices. To check the validity of oursimulations, we first consider clean (undiluted) lattices.We find critical temperatures of T c = 2 . T c = 1 . p c = 0 . p < . T c from the cross-ing of the g ( T ) curves of the two largest system sizes, L = 80 and L = 112 . The ensemble method is ap-plied to find the error of T c . Fig. 3 shows an exampleof this situation for dilution p = 0 . p ≥ .
64) in the vicinity ofthe percolation threshold p c , the crossing of the Bindercumulant vs. temperature curves is less sharp. Specifi-cally, the crossing temperature T ∗ ( L ) of the curves for T g L=20L=28L=40L=56L=80L=112
Fig. 4
Binder cumulant g vs temperature T for the XY modelon a cubic lattice for dilution p = 0 .
65, i.e. close to p c . Thecurves do not all cross at the same temperature. Instead, thecrossing progressively shifts as L increases. The statisticalerrors arising from the Monte Carlo simulation are smallerthan the symbol size. linear system sizes L and √ L shifts visibly towardshigher temperatures as the system sizes are increased.An example (for p = 0 .
65) is shown in Fig. 4. This shiftis caused by corrections to the leading finite-size scalingbehavior (8). It can be modeled as T ∗ ( L ) = T c + bL − ω (10)where ω is the leading irrelevant exponent and b is anon-universal constant. To find the asymptotic (infinitesystem size) value of T c , we fit the crossing temperature T ∗ ( L ) to Eq. 10. As ω is expected to be universal, i.e., totake the same value for all dilutions near p c , we performa combined fit for all dilutions p ≥ .
64 and treat ω asa fitting parameter. This combined fit produces ω =1 . ± .
4. An example of the resulting extrapolation ispresented in Fig. 5 for p = 0 .
65. The figure shows thatthe finite-size shifts of the crossing temperature are notvery strong. This is further confirmed in Fig. 6 whichpresents an overview of the fits for all dilutions from p = 0 .
64 to p = 0 . T c ( p ) of the site-dilutedXY model on a cubic lattice is shown in Fig. 7. Theoverview given in the inset demonstrates that T c ( p ) isindeed continuously suppressed with increasing p andapproaches zero as p → p c . To analyze the functionalform of T c ( p ) close to p c , the main panel of Fig. 7 showsa log-log plot of T c vs. | p − p c | . We observe that thephase boundary follows two different power laws, closeto the percolation threshold p c and further away from p c . The asymptotic value of φ is determined from a fit ofthe data closest to p c (viz. p between 0 .
678 to 0 . φ = 1 . L T * p=0.65 Fig. 5
Extrapolation to infinite system size of the crossingtemperature T ∗ of the Binder cumulant curves for systemsizes L and √ L using ω = 1 .
5. The dilution is p = 0 .
65. A fitto Eq. (10) gives T c = 0 . T ∗ havebeen determined using the ensemble method described in Sec.4.2. L T * p=0.64p=0.65p=0.66p=0.67p=0.675p=0.678p=0.68p=0.6815p=0.6825 Fig. 6
Overview of the extrapolations of the crossing tem-peratures T ∗ for several dilutions near p c using ω = 1 .
5. Theerror bars ∆T ∗ are smaller than the symbols. estimate is a combination of the statistical error fromthe fit and a systematic error estimated from the ro-bustness of the value against changes of the fit interval.The asymptotic value of φ agrees reasonably well withthe prediction of percolation theory. The asymptoticpower law describes the data for dilutions above about p = 0 .
65. The asymptotic critical region thus rangesfrom about p = 0 .
65 to p c = 0 . T c ( p ) for p between p = 0 to p = 0 .
64 also follows a power law in goodapproximation. However, the exponent is significantlybelow unity, φ = 0 . p ln | p p c | l n ( T c ) = 0.80(1)= 1.09(2) p T c Fig. 7
Phase boundary of the site-diluted XY model on acubic lattice. Main panel: Log-log plot of T c vs. | p − p c | . Thestraight lines are power-law fits, T c ∼ | p − p c | φ . They are shownas solid lines within the fit range. The dotted and dash-dottedlines are extrapolations. For details see text. Inset: Overviewpresented as linear plot of T c vs. p . All error bars of the datapoints are smaller than the symbol size. T g L=20L=28L=40L=56L=80L=112
Fig. 8
Binder cumulant g vs temperature T for dilution p =0 .
65 on cubic lattice and Heisenberg spins. All curves cross atthe same temperature. Error bars are smaller than the symbolsize. we gradually increase dilution and find T c ( p ). In thecase of Heisenberg spins, we find that the correctionsto finite-size scaling are weaker than in the XY case.Even in the vicinity of p c , all Binder cumulant curvesintersect in a single point within their statistical errors.As an example, the g vs T data for p = 0 .
65 are shownin Fig. 8. The critical temperatures T c ( p ) and its er-ror are therefore determined from the Binder cumulantcrossing for system sizes L = 80 and L = 112 , thelargest systems simulated.The phase boundary of the site-diluted Heisenbergmodel on a cubic lattice is constructed from these data p ln | p p c | l n ( T c ) = 0.86(1)= 1.08(2) p T c Fig. 9
Phase boundary of the site-diluted Heisenberg modelon a cubic lattice. Main panel: Log-log plot of T c vs. | p − p c | .The straight lines are power-law fits, T c ∼ | p − p c | φ . They areshown as solid lines within the fit range. The dotted and dash-dotted lines are extrapolations. For details see text. Inset:Overview presented as linear plot of T c vs. p . All error barsof the data points are smaller than the symbol sizes. and shown in Fig. 9. Similar to the XY case, we ob-serve two separate power law exponents governing thephase boundary. The dilutions p (cid:38) .
65 constitute theasymptotic critical region with crossover exponent φ =1 . p (cid:46) .
62 is again smallerthan unity, φ = 0 . p c = 0 . T c for a givendilution is determined from the Binder cumulant cross-ings. Corrections to the finite-size scaling were foundto be negligible within the statistical errors. Thus, weused the Binder cumulant crossing of the two largestsystem sizes (28 and 40 double unit cells) to find T c . The resulting phase boundary is shown in Fig. 10. p ln | p p c | l n ( T c / K ) = 0.88(2)= 1.12(3) p T c ( K ) Fig. 10
Phase boundary for the Heisenberg model on ahexagonal ferrite lattice. The main panel shows the log-logplot of T c vs. | p − p c | . The statistical errors of the data (de-termined by the ensemble method) are smaller than the sym-bol size. The straight lines are fits to T c ∼ | p − p c | φ . Theyare shown as solid lines within the fit range. The dotted anddash-dotted lines are extrapolations. For details see text. Theinset shows a linear plot the complete phase boundary T c ( p ). The behavior of this phase boundary is very similarto the cubic lattice results. High dilutions, p (cid:38) . φ = 1 . φ = 0 . T c ( x ) = T c (0)(1 − x/x c ) / .In the simulations, the transition temperature T c is sup-pressed more rapidly with x than in the experimentaldata (see Fig. 11). To explore possible reasons for thisdiscrepancy, we also perform test simulations that in-clude additional weaker exchange interactions [19] thatfrustrate the ferrimagnetic order. The results of thesesimulations, which are included in Fig. 11, show thatthese weaker frustrating interactions have little effectat low dilutions. At higher dilutions, when the ferri-magnetic order is already weakened, the frustrating in-teractions further suppress the transition temperature.They thus further increase the discrepancy between theexperimental data and the Monte Carlo results. To summarize, motivated by recent experimental ob-servations on hexagonal ferrites, we have studied clas- x T / c / ( K / ) Experimental DataMonte Carlo, leading + frustrated interactionsMonte Carlo, leading interactions only
Fig. 11
Comparison between the numerically determinedphase boundary T c ( x ) and the experimental data forPbFe − x Ga x O [15]. The tuning parameter x is related tothe dilution by x/
12 = p . The Monte Carlo simulations showa more rapid suppression of T c with x . Including additionalweak frustrated interactions increases the discrepancy. sical site-diluted XY and Heisenberg models by meansof large-scale Monte Carlo simulations, focusing on theshape of the magnetic phase boundary. We have ob-tained two main results.First, for high dilutions close to the lattice percola-tion threshold, the critical temperature depends on thedilution via the power law T c ∼ | p − p c | φ in all studiedsystems. In this asymptotic region, we have found thevalues φ = 1 . φ = 1 . φ =1 . φ thus appears to be super-universal, i.e., ittakes the same value not just for different lattices butalso for XY and Heisenberg symmetry.Interestingly, the asymptotic critical region of thepercolation transition is very narrow, as the asymptoticpower-laws only hold in the range | p − p c | (cid:46) .
04. Atlower dilutions, the phase boundary still follows a powerlaw in | p − p c | , but with an exponent that appears tobe non-universal and below unity (in the range between0.8 and 0.9).Our second main result concerns the origin of the2 / T c ( x ) = T c (0)(1 − x/x c ) / , that wasexperimentally observed in PbFe − x Ga x O over theentire concentration range between 0 and close to thepercolation threshold [15]. Neither the asymptotic northe preasymptotic power laws identified in the simu-lations match the experimental result. In fact, in allsimulations, the critical temperature is suppressed morerapidly with increasing dilution than in the experiment. The observed shape of the magnetic phase boundary inPbFe − x Ga x O thus remains unexplained.Potential reasons for the unusual behavior may in-clude the interplay between magnetism and ferroelec-tricity in these materials [32] or the presence of quan-tum fluctuations (arising from the frustrated magneticinteractions mentioned above), even though it is hardto imagine that these stay relevant at temperatures ashigh as 720 K. Another possible explanation could bea statistically unequal occupation of the different ironsites in the unit cell by Ga ions. Exploring these pos-sibilities remains a task for the future. Disentanglingthese effects may also require additional experimentsintroducing further tuning parameters such as pressureor magnetic field in addition to chemical composition. Acknowledgements
We acknowledge support from the NSFunder Grant Nos. DMR-1506152, DMR-1828489, and OAC-1919789. The simulations were performed on the Pegasus andFoundry clusters at Missouri S&T. We also thank MartinPuschmann for helpful discussions.
References
1. G. Grinstein and A. Luther, Phys. Rev. B , 1329(1976).2. D. S. Fisher, Phys. Rev. Lett. , 534 (1992).3. D. S. Fisher, Phys. Rev. B , 6411 (1995).4. R. B. Griffiths, Phys. Rev. Lett. , 17 (1969).5. M. Thill and D. A. Huse, Physica A , 321 (1995).6. A. P. Young and H. Rieger, Phys. Rev. B , 8486 (1996).7. T. Vojta, Phys. Rev. Lett. , 107202 (2003).8. R. Sknepnek and T. Vojta, Phys. Rev. B , 174410(2004).9. G. Schehr and H. Rieger, Phys. Rev. Lett. , 227201(2006).10. J. A. Hoyos and T. Vojta, Phys. Rev. Lett. , 240601(2008).11. T. Vojta, J. Phys. A , R143 (2006).12. T. Vojta, J. Low Temp. Phys. , 299 (2010).13. T. Vojta, Ann. Rev. Condens. Mat. Phys. , 233 (2019).14. D. Stauffer and A. Aharony, Introduction to PercolationTheory (CRC Press, Boca Raton, 1991).15. S. E. Rowley, T. Vojta, A. T. Jones, W. Guo, J. Oliveira,F. D. Morrison, N. Lindfield, E. Baggio Saitovitch, B. E.Watts, and J. F. Scott, Phys. Rev. B , 020407 (2017).16. G. Albanese, F. Leccabue, B. E. Watts, and S. D´ıaz-Casta˜n´on, J. Mat. Sci , 3759 (2002).17. A. Coniglio, Phys. Rev. Lett. , 250 (1981).18. T. Vojta and J. A. Hoyos, in Recent Progress in Many-Body Theories , edited by J. Boronat, G. Astrakharchik,and F. Mazzanti (World Scientific, Singapore, 2008) p.235.19. C. Wu, Z. Yu, K. Sun, J. Nie, R. Guo, H. Liu, X. Jiang,and Z. Lan, Scientific Reports , 36200 (2016).20. E. Shender and B. Shklovskii, Physics Letters A , 77(1975).21. B. Kozlov and M. Lagu¨es, Physica A: Statistical Mechan-ics and its Applications , 5339 (2010).22. J. Wang, Z. Zhou, W. Zhang, T. M. Garoni, and Y. Deng,Phys. Rev. E , 052107 (2013).23. U. Wolff, Phys. Rev. Lett. , 361 (1989).24. N. Metropolis and S. Ulam, Journal of the American sta-tistical association , 335 (1949).25. H. G. Ballesteros, L. A. Fern´andez, V. Mart´ın-Mayor,A. Mu˜noz Sudupe, G. Parisi, and J. J. Ruiz-Lorenzo,Phys. Rev. B , 2740 (1998).26. T. Vojta and R. Sknepnek, Phys. Rev. B , 094415(2006).27. Q. Zhu, X. Wan, R. Narayanan, J. A. Hoyos, and T. Vo-jta, Phys. Rev. B , 224201 (2015).28. K. Binder, Zeitschrift f¨ur Physik B , 119 (1981).29. A. P. Gottlob and M. Hasenbusch, Physica A: StatisticalMechanics and its Applications , 593 (1993).30. R. G. Brown and M. Ciftan, Phys. Rev. B , 224413(2006).31. J. Wang, Z. Zhou, W. Zhang, T. M. Garoni, and Y. Deng,Phys. Rev. E , 052107 (2013).32. S. E. Rowley, Y.-S. Chai, S.-P. Shen, Y. Sun, A. T. Jones,B. E. Watts, and J. F. Scott, Scientific Reports6