Existence of a Thermodynamic Spin-Glass Phase in the Zero-Concentration Limit of Anisotropic Dipolar Systems

The nature of ordering in dilute dipolar interacting systems dates back to the work of Debye and is one of the most basic, oldest and as-of-yet unsettled problems in magnetism. While spin-glass order is readily observed in several RKKY-interacting systems, dipolar spin-glasses are subject of controversy and ongoing scrutiny, e.g., in ${{\rm LiHo_xY_{1-x}F_4}}$, a rare-earth randomly diluted uniaxial (Ising) dipolar system. In particular, it is unclear if the spin-glass phase in these paradigmatic materials persists in the limit of zero concentration or not. We study an effective model of ${{\rm LiHo_xY_{1-x}F_4}}$ using large-scale Monte Carlo simulations that combine parallel tempering with a special cluster algorithm tailored to overcome the numerical difficulties that occur at extreme dilutions. We find a paramagnetic to spin-glass phase transition for all Ho ion concentrations down to the smallest concentration numerically accessible of 0.1%, and including Ho ion concentrations which coincide with those studied experimentally up to 16.7%. Our results suggest that randomly-diluted dipolar Ising systems have a spin-glass phase in the limit of vanishing dipole concentration, with a critical temperature vanishing linearly with concentration, in agreement with mean field theory.


I. INTRODUCTION
Dipolar interactions are ubiquitous in nature and often dominate other types of interactions in many simple systems, e.g., in insulating magnets. Additionally, dipolarlike couplings may also arise between defect states in crystals, where these are mediated by the phonon vacuum. Such interactions are known to induce magnetic and ferroic order in densely packed solids and liquids [1,2]. With dilution, the competing nature of the interaction and spatial disorder lead to a spin-glass (SG) order [3] at low temperatures. Mean-field theory suggests that the SG order is maintained at x ¼ 0 þ , with the critical temperature being linear in the concentration x [4,5]. However, at low concentrations, spatial inhomogeneities are large and could dominate the characteristics of the system. Because the dipolar nature of the interaction renders standard spatial renormalization group methods ineffective, rigorous analytic conclusions are currently beyond reach. Thus, the nature of anisotropic dipolar systems in general, and in the limit of low concentrations in particular, has been a longstanding controversy.
Experimentally, LiHo x Y 1−x F 4 is perhaps the best-studied dilute dipolar (strongly anisotropic) Ising system. This rareearth compound has attracted vast experimental, numerical, and theoretical interest in the past two decades. Its scrutinization has enhanced the understanding of many different magnetic phenomena, such as quantum phase transitions [6][7][8], large spin tunneling [9][10][11], quantum annealing [12], quantum entanglement [13], quantum domain-wall tunneling [14], random-field physics [15][16][17][18], and generic disordering mechanisms [19]. Thus, establishing the low-temperature phase of LiHo x Y 1−x F 4 at small Ho þ concentrations is not only of fundamental interest but is crucial for the further study of its characteristics. However, extremely long equilibration times have produced conflicting experimental results and a strong dependence on the used experimental protocol, with no clear evidence for the equilibrium phase of the system at low concentrations [13,[20][21][22][23][24][25][26]. Furthermore, where experimental data suggest the existence of a spin-glass phase, reported values for the critical temperature deviate markedly at low concentrations from the expected linear dependence on the Ho concentration [25].
Numerically, understanding the nature of dipolar Ising systems at small concentrations is notoriously difficult because spatial inhomogeneities are large and so are the required system sizes. Previous Monte Carlo simulations of the dilute dipolar Ising spin-glass model [27][28][29][30] showed no sign of a SG transition. More recent simulations where a better observable, namely, the finite-size two-point correlation function [31], was used suggest the existence of a SG phase down to a concentration of x ¼ 0.0625 [32]. However, the regime of much theoretical interest, where the typical distance between spins is much larger than interatomic distance, and thus fluctuations are large, could not be reached. Long equilibration times due to the slow dynamics of the system [33] limited the studied system sizes and concentrations, i.e., strong finite-size corrections in the data. As such, the nature of anisotropic dipolar systems at very low concentrations remains unclear.
Here, we present conclusive evidence for the existence of a SG phase in the dilute dipolar Ising model in the limit of x ¼ 0 þ , and for LiHo x Y 1−x F 4 for all experimentally relevant low concentrations. We use large-scale Monte Carlo simulations that combine parallel tempering [34] and a cluster algorithm [35,36] that allows us to efficiently handle the atypically large interactions stemming from rare nearby groups of spins and, at the same time, leaves the prevalent typical interactions for standard numerical treatment. We find clear evidence that the anisotropic dipolar glass has a SG phase at low temperatures for all studied concentrations down to x ¼ 10 −3 (almost 2 orders of magnitude smaller than the concentration reached in previous studies), with a critical temperature T c that is linear in the concentration x. Furthermore, our data show that for all x, the divergence of the correlation length at the transition is likely described by the same critical exponent ν. This strongly suggests that our results can be carried through to vanishing spin concentrations.
The paper is structured as follows. In Sec. II, we introduce a model Hamiltonian for LiHo x Y 1−x F 4 and outline our numerical approach to study the system. We present results in Sec. III, followed by a discussion and concluding remarks in Sec. IV.

II. MODEL AND NUMERICAL DETAILS
The tunnel splitting induced between the two polarized electronuclear single Ho ground states in the dilute LiHo x Y 1−x F 4 system by the crystal field and by the offdiagonal terms of the dipolar interaction is much smaller than the typical interaction down to extremely low Ho concentrations [37]. The same is true for the magnetic interaction between the nuclear spins of the F atoms. Thus, down to very low x, LiHo x Y 1−x F 4 is well described by a classical Ising spin model [19,25,32,33], i.e., Here, ϵ i ¼ f0; 1g is the occupation of the magnetic Ho 3þ ions on a tetragonal lattice (lattice constants a ¼ b ¼ 5.175 Å and c ¼ 10.75 Å) with four ions per unit cell [32,38], i.e., N ¼ 4L 3 spin sites. S i ∈ fAE1g are Ising spins. The magnetostatic dipolar coupling J ij between two Ho 3þ ions is given by Example finite-size scaling data collapse of the spin-glass two-point finite-size correlation function ξ L =L using an extended scaling form [39] for the extreme dilution x ¼ 0.001. The optimal scaling is accomplished via a Levenberg-Marquardt minimization procedure described in Ref. [40]. (b) Unscaled two-point finite-size correlation function divided by the linear system size ξ=L for the extreme dilution x ¼ 0.001 and system sizes L ¼ 36; …; 60. Data for different system sizes cross at the critical temperature T c . The vertical shaded area represents the estimated critical temperature T c ¼ 0.00052ð5Þ K; its width is the statistical uncertainty determined via an extended finite-size scaling [panel (a)]. The inset zooms into the critical temperature region showing that the data do indeed cross. r ij ¼ jr i − r j j, r i is the position of the i-th Ho 3þ ion and z ij ¼ ðr i − r j Þ ·ẑ is the component parallel to the easy axis. The dipolar constant is D=a 3 ¼ 0.214 K [29], and the antiferromagnetic nearest-neighbor exchange is set to J ex ¼ 0.12 K [38]. For the low concentrations x of interest to us here, this model is equivalent to the pure Ising dipolar model because the exchange interactions only slightly change the interaction strength of the rare nearby pairs, which, as we show below, do not affect the thermodynamics at and near the phase transition.
To determine the finite-temperature transition for a given value of x, we measure the two-point finite-size correlation function [31] Here, hÁ Á Ái T represents a thermal average, R i is the spatial location of the spin S i , and k min represents the smallest nonzero wave vector in the a-or c-axis direction, k min ¼ ð2π=L; 0; 0Þ or k min ¼ ð0; 0; 2π=LÞ, respectively. ξ L =L is dimensionless, and near the transition, it is expected to scale as ξ L =L ∼X½L 1=ν ðT − T c Þ. Because corrections to scaling are typically large for highly dilute systems, we use an extended scaling approach [39] that has proven to reduce scaling corrections and where ξ L =L∼ X½ðLTÞ 1=ν j1 − ðT=T c Þ 2 j; see Fig. 1(a). When T ¼ T c , the argument of the scaling function is zero (up to scaling corrections) and hence independent of L. As such, lines for different system sizes L cross [see Fig. 1(b)]. If, however, the lines do not meet, we know that no transition occurs in the studied temperature range. The best estimate of the critical temperature T c ðxÞ is determined by applying a Levenberg-Marquardt minimization combined with a bootstrap analysis to determine statistical error bars [40] to the aforementioned extended finite-size scaling analysis. An example of the resulting data collapse using the minimization is shown in Fig. 1(a) for x ¼ 0.001. Finally, we note that in Ref. [32] it was observed that the estimates of T c along the a axis tend to be systematically lower than the ones computed along the c axis of the material. A comparison with experimental results [25] showed an agreement between the experimental estimates and the numerical estimates along the c axis only. Similarly, in this work, our estimates of the transition temperatures from a paramagnetic (PM) phase to a SG phase computed along the c axis tend to be systematically higher for all studied dilutions and agree better with the experimental results of Quilliam et al. [25]. As such, all quoted results stem from simulation results with measurements along the c axis.
In the simulations, we use the Ewald summation method without a demagnetization factor to compute the periodic boundary conditions [32,41] for systems of up to approximately 7 × 10 6 lattice sites. To equilibrate the system at extreme dilutions, we use a combination of single spin-flip Monte Carlo dynamics and a cluster renormalization algorithm [35,36] combined with parallel tempering Monte Carlo [34]. The cluster renormalization algorithm is tailored to treat strongly coupled spins efficiently. It does not fully obey detailed balance; however, it has been successfully applied to different model systems [42]. The cluster renormalization technique works as follows: At the beginning of the simulation, set the random positions of the spins and search for clusters C i J i of spins coupled by an interaction of at least jJ i j. Once all clusters have been TABLE I. Simulation parameters for the different concentrations x and linear system sizes L studied. T min (T max ) is the lowest (largest) simulated temperature, and N T is the number of temperatures used in the parallel tempering Monte Carlo method. We thermalize and measure for 2 X Monte Carlo sweeps, and N sa is the number of disorder realizations. Note that to compute the spin-glass order parameter, we actually simulate 2N T replicas (N T for the tempering method and two per temperature to compute the spin overlap). Note that smaller system sizes have also been studied; however, for brevity, we do not list them because they share the same parameters as the smallest L listed for a given concentration x. labeled and recorded, search for sets of spins C iþ1 J iþ1 coupled by at least jJ iþ1 j ¼ jJ i j þ δ J and record these. Iteratively perform this renormalization i ¼ 0; 1; …; n times until all the sets of C n J n consist of only spin pairs. It is very important to carefully tune the renormalization procedure to the studied model. In this case, we use jJ 0 j ¼ x, where x is the concentration of the magnetic ions in the system. A suitable step value δ J is 2x. Spins in clusters are then flipped, regardless of their sign. Note that one Monte Carlo sweep in the simulation thus consists of the following procedure: For each spin in the system, we either perform a single-spin simple Monte Carlo flip with probability 0.75, or we flip a randomly selected cluster.
Finally, to verify that the data are properly thermalized, a logarithmic binning analysis is used: Observables are measured and averaged over an exponentially growing number of Monte Carlos sweeps 2 X , and their Monte Carlo time evolution is monitored. When at least three bins agree within error bars and are independent of Monte Carlo time, we deem the system to be in thermal equilibrium. Simulation parameters are listed in Table I.

III. RESULTS
Our main result is shown in the phase diagram depicted in Fig. 2, as well as Table II. The estimated critical temperatures T c ðxÞ show a clear linear behavior for 2 orders of magnitude, strongly suggesting that the SG phase extends to the zero concentration limit. Comparison to experiment shows that for x ¼ 0.167, we find T c ¼ 0.131ð3Þ, in excellent agreement with the experimental results of T c ¼ 0.133ð5Þ [45]. Similarly, for x ¼ 0.08, we find T c ¼ 0.048ð2Þ, close to the most recent experimental result, T c ¼ 0.065ð3Þ [25]. For lower concentrations, however, we find values of T c that agree well with the linear extrapolation of the experimental data from higher x but are lower than the experimentally obtained values for x ¼ 0.045 and x ¼ 0.018. We attribute this discrepancy to the microscopic time scale in LiHo x Y 1−x F 4 being very long at low temperatures and enhanced with the decrease of the Ho concentration. This results in the difficulty to equilibrate the system close to the critical temperature at the lowest experimentally studied concentrations [25,26]. Furthermore, it was argued that long equilibration times of small clusters lead to a quantum nonequilibrium state, whose nature depends on the degree of coupling to the environment [26].
In Fig. 3, we show a log-log plot of the high-dilution limit of the transition temperature T c versus concentration x phase diagram shown in the inset to Fig. 2. The solid line in the figure that separates the SG from the PM phase is a fit to T c ðxÞ ¼ ax with a ¼ 0.59ð1Þ. Allowing for a finite intercept, i.e., T c ðxÞ ¼ T 0 c þ ax, yields a ¼ 0.60ð2Þ and T 0 c ¼ −0.00097ð94Þ, which is statistically compatible with FIG. 2. Transition temperature T c versus concentration x phase diagram obtained from our simulations (circles), previous simulations from Refs. [19,32] (triangles), and experimental data from Refs. [6,18,25,43,44] (squares). The straight line between the spin-glass (SG) and paramagnetic (PM) phases is a linear fit to the data for small concentrations. Our results suggest that the spin-glass phase extends to the x ¼ 0 þ limit, as can be seen in more detail in the inset. For large concentrations and low temperatures, the system is a ferromagnet (FM).  (4) FIG. 3. Log-log plot of the high-dilution limit of the transition temperature T c versus the concentration x phase diagram obtained from our simulations. The solid line separating the spin glass (SG) from the paramagnetic (PM) phase is a fit to the numerical data of the form T c ðxÞ ¼ ax with a ¼ 0.59ð1Þ. The spin-glass phase therefore extends to the x ¼ 0 þ limit. a zero intercept. Therefore, we see strong evidence that the spin-glass phase extends to the x ¼ 0 þ limit. We note that despite the much-enhanced fluctuations in the distribution of interactions as x is reduced, equilibration times are similar for all system sizes. The simulations for the lowest concentration were limited by the time it takes to accurately compute the Ewald summation used to account for the periodic boundary conditions, and not the Monte Carlo simulation time, therefore showing the effectiveness of the implemented algorithm even for very high dilutions. This bottleneck is in fact much easier to overcome because the Ewald summation has to be performed only once for each system size at the beginning of the simulation.
In Fig. 4, we show the critical exponent ν as a function of x computed from a finite-size scaling of the two-point correlation function. For all concentrations depicted, data seem to agree approximately within error bars-an indicator for potential universal behavior. All critical parameters from the scaling analysis of the two-point correlation function are listed in Table II. Note, however, that to fully characterize a universality class, it is necessary to determine at least two critical exponents. We tried different analysis methods to extract a second critical exponent η from the spin-glass susceptibility. However, no robust estimate was possible, a common problem in spin-glass simulations (see also Ref. [40]).

IV. DISCUSSION
Dilute power-law interacting systems are natural candidates for emergent geometric similarity [46] whereby statistical mechanics of systems at different concentrations may be mapped onto each other [47]. The characteristic 1=r d falloff of the dipolar kernel implies linear scaling of typical interactions with concentration in any dimension d and suggests similar scaling of relevant temperature scales. However, dipolar systems at different concentrations are not quite geometrically similar. Rescaling of the interactions leaves the distribution of interactions practically unchanged at low and typical values but generates a progressively stronger tail at high values because the largest coupling is fixed by the lattice spacing, independent of the concentration. These large couplings produce physical correlation effects that impeded simulation progress in the past. Remarkably, focusing on and solving this relatively local high-energy bottleneck allows for essentially unimpeded progress on the rest of the problem. Our results provide strong support to the notion of emergent geometric similarity by locating and characterizing the spin-glass ordering transition over nearly 2 orders of magnitude in concentration, with the transition temperature scaling linearly with concentration.
We note here that this geometric similarity is destroyed by the application of a transverse field [48], which, in combination with the off-diagonal elements of the dipolar interaction, results in effective random fields in the longitudinal direction [15][16][17]. The emergent longitudinal fields have a large variance, are correlated with the interactions, and lead to a much more effective disordering of the spin-glass phase [48] than that predicted by the naive application of the Imry-Ma argument [49,50]. We emphasize that in the absence of an applied field, the effect of the off-diagonal dipolar terms on the thermodynamic phase of the system is negligible. The quantum fluctuations induced by these is much smaller than the interaction, and they do not change the Z 2 time-reversal symmetry (spin inversion) of the Hamiltonian.
Summarizing, using large-scale Monte Carlo simulations that combine parallel tempering with an innovative cluster renormalization algorithm [35,36], we have shown that the dilute dipolar Ising model has a spin-glass transition at low temperatures for concentrations down to x ¼ 10 −3 . Furthermore, a clear linear behavior of T c ðxÞ ∼ x is found in the highly dilute regime, strongly suggesting that the SG phase transition extends to the x ¼ 0 þ limit. The shaded area corresponds to their average value over all concentrations. The individual estimates come from an extended finite-size scaling analysis of the two-point correlation length; an example is shown in Fig. 1(a). All estimates agree within error bars, meaning that the critical exponent ν might be independent of the concentration x. This hints towards the possibility of a common universality class.