Partial up-up-down order with the continuously distributed order parameter in the triangular antiferromagnet TmMgGaO4

Frustrated quasidoublets without time-reversal symmetry can host highly unconventional magnetic structures with continuously distributed order parameters even in a single-phase crystal. Here, we report the comprehensive thermodynamic and neutron diffraction investigation on the single crystal of TmMgGaO$_4$, which entails non-Kramers Tm$^{3+}$ ions arranged on a geometrically perfect triangular lattice. The crystal electric field (CEF) randomness caused by the site-mixing disorder of the nonmagnetic Mg$^{2+}$ and Ga$^{3+}$ ions, merges two lowest-lying CEF singlets of Tm$^{3+}$ into a ground-state (GS) quasidoublet. Well below $T_c$ $\sim$ 0.7 K, a small fraction of the antiferromagnetically coupled Tm$^{3+}$ Ising quasidoublets with small inner gaps condense into two-dimensional (2D) up-up-down magnetic structures with continuously distributed order parameters, and give rise to the \emph{columnar} magnetic neutron reflections below $\mu_0H_c$ $\sim$ 2.6 T, with highly anisotropic correlation lengths, $\xi_{ab}$ $\geq$ 250$a$ in the triangular plane and $\xi_c$ $<$ $c$/12 between the planes. The remaining fraction of the Tm$^{3+}$ ions remain nonmagnetic at 0 T and become uniformly polarized by the applied longitudinal field at low temperatures. We argue that the similar model can be generally applied to other compounds of non-Kramers rare-earth ions with correlated GS quasidoublets.


I. INTRODUCTION
Geometrical frustration can render the ground state(s) of the correlated spin system macroscopically degenerate and completely disordered in the classical Ising case [1][2][3][4][5][6], or trigger strong quantum fluctuations that prevent the conventional symmetry breaking even down to ∼ 0 K in the quantum case [7][8][9][10][11][12][13][14][15]. In most of the previously proposed frustrated magnets, the (effective) S = 1/2 dipole moments of the ground-state (GS) doublets are protected either by time-reversal symmetry in the case of Kramers ions with an odd number of electrons per site [4,[12][13][14][15][16] or by the local crystal electric field (CEF) symmetry in the case of non-Kramers ions with an even number of electrons per site [17][18][19][20]. In non-Kramers ions without symmetry-protected doublets, two close-lying singlets will typically occur [21], as in the three-dimensional dipolar Ising ferromagnet LiTbF 4 [22][23][24] and in the kagome magnet Pr 3 Ga 5 SiO 14 [25]. However, geometrically frustrated non-Kramers magnets with correlated GS quasidoublets are still rare to date. Here, the complete site-mixing disorder between two nonmagnetic ions with different valences significantly distributes the energies of the two lowest-lying and nearly degenerate * yuesheng li@hust.edu.cn † hao.deng@frm2.tum.de ‡ philipp.gegenwart@physik.uni-augsburg.de CEF singlets. And the intersite spin interactions combined with the single-ion terms and the randomness can lead to exotic quantum phases at low temperatures [26].
In a search for such a kind of material with the twodimensional (2D) triangular arrangement of the 4f ions, we explored structural siblings of YbMgGaO 4 , which we recently characterized as a quantum spin liquid (QSL) candidate with the experimental evidence for breaking and re-arrangement of uncorrelated/resonating valence bonds [27][28][29]. Despite the disorder effect on the CEF and the putative QSL state [29][30][31] caused by the sitemixing disorder between nonmagnetic Mg 2+ and Ga 3+ ions in YbMgGaO 4 , the GS CEF doublet of the Kramers Yb 3+ ion always gives rise to the effective S = 1/2 magnetic moment at T ∆ CEF /k B ∼ 460 K, where ∆ CEF is the energy gap to the first excited level [16]. While, the many-body correlated physics may be significantly changed when Yb 3+ is replaced by the non-Kramers Tm 3+ (4f 12 ) ion on the triangular lattice, as the time reversal symmetry is no longer preserved, and the previously symmetry-protected degeneracy of the GS CEF doublet may get "lifted". And thus exotic quantum phases may emerge in the new frustrated magnet, TmMgGaO 4 . Recently, single crystals of TmMgGaO 4 were successfully synthesized by Cevallos et al. [32], which provides an opportunity to study its exotic correlated magnetism experimentally.
In this paper, we report a thorough single-crystal investigation of the low-temperature magnetism of TmMgGaO 4 , including heat capacity, Faraday force magnetization (susceptibility), magnetocaloric effect, and neutron diffraction measurements, down to 30 mK. The Mg 2+ /Ga 3+ disorder significantly distributes the energies of the two lowest-lying CEF singlets, thus mixing them into a GS quasidoublet. At low temperatures and in small longitudinal fields, a fraction of the Tm 3+ ions -those characterized by a small gap between the two CEF singlets -give rise to the novel 2D Ising up-up-down (uud) phase with the continuously distributed order parameter, under the frustrated intersite couplings on the triangular lattice. The remaining large fraction with the large inner gaps remains nonmagnetic at 0 T and becomes uniformly polarized by the applied longitudinal field. Using the random many-body-correlated model of the GS quasidoublets, we can naturally interpret most of the low-T magnetic properties. A similar model can be generally applied to other non-Kramers rare-earth magnets with correlated GS quasidoublets.

II. TECHNICAL DETAILS
High-quality single crystals (∼ 1 cm) of TmMgGaO 4 , Tm 0.04 Lu 0.96 MgGaO 4 , and Yb 0.04 Lu 0.96 MgGaO 4 were grown by the floating zone technique (Appendix A) [11,32]. The Faraday force magnetization [33], heat capacity, magnetocaloric effect (magnetic Grüneisen ratio) [34,35] down to 30 mK were measured in a 3 He-4 He dilution refrigerator (Appendix B). The neutron diffraction experiments were carried out in the ab plane (L = 0) and along the c axis (L = 0), on the CEA-CRG single crystal diffractometer D23 [36] of Institut Laue-Langevin (ILL) in France and on the single crystal diffractometer POLI [37] of Heinz Maier-Leibnitz Zentrum (MLZ) in Germany, respectively, down to 60 mK and up to 5 T. Using the Matlab codes, we performed CEF, exact diagonalization (ED), and spin-wave calculations for a model spin Hamiltonian. Then we simultaneously fit this model to the temperature dependence of direct-current (dc) susceptibility and heat capacity measured at ∼ 0 T, as well as the field dependence of magnetization at 40 mK, by minimizing the following function, where N 0 , X obs i and σ obs i are the number of the data points, the observed value and its standard deviation, respectively, whereas X cal i is the calculated value.

III. SINGLE-ION PHYSICS
Generally, the CEF of Tm 3+ with the D 3d point-group symmetry of TmMgGaO 4 splits the 13-degenerate GS of the free Tm 3+ ion with the total angular momentum J = 6, |m J (m J = 0, ±1..., ±J), into five sin-glets (3A 1g +2A 2g ) and four doublets (4E g ), according to the symmetry analysis. In the following, we will further determine the low-lying CEF states of Tm 3+ in TmMgGaO 4 by thermodynamical measurements.
At ∼ 1.9 K, the effective spin-1/2 moments of Tm 3+ can be fully polarized above ∼ 8 T applied along the c axis [see Fig. 1 (a)]. Through a linear fit to this highfield magnetization data, we obtain the fitted intercept that measures the saturated effective spin-1/2 magnetic moment g /2, where g = 13.18 (1), and the small slope that corresponds to the van Vleck susceptibility, χ vv = 0.003(1) cm 3 /mol [10,11]. While, along the ab plane the magnetization shows a linear field dependence between 0 and 12 T with the nearly zero intercept [see Fig. 1 (a)], suggesting the strict Ising anisotropy of the Tm 3+ magnetic moments (Appendix C). Between ∼ 30 and 60 K, the magnetic entropy of TmMgGaO 4 is measured to be constant, S m ∼ Rln2 per mole [see Fig. 2 (e)], confirming the GS CEF (quasi)doublets, and thus the formation of the effective Ising spin-1/2 moments of Tm 3+ below 60 K (Appendix C).
This Ising nature is rooted in the GS CEF quasidoublets [38]. To better understand its nature, we first prepared the highly diluted samples of Tm x Lu 1−x MgGaO 4 (x = 0.04), where no intersite interactions occur, and single-ion physics of Tm 3+ can be probed. The diluted Yb x Lu 1−x MgGaO 4 sample with the Yb 3+ Kramers ion was also studied as reference. In both cases, the dilution eliminates any intersite magnetic couplings, as confirmed by the diminutively small Curie-Weiss temperatures, θ w (x = 0.04) ∼ 0.16θ w (x = 1) [see Fig. 1 (a) for Tm x Lu 1−x MgGaO 4 and Ref. [10] for Yb x Lu 1−x MgGaO 4 ] . The small Curie-Weiss temperature is obtained by the fit to χ measured in the temperature range where S m ∼ Rln2, and thus the CEF effect due to excitations to higher CEF levels is negligible, and we obtain θ w = -3(J zz 1 + J zz 2 )/2 (see below) [11,39]. Therefore, the magnetic ions should be almost homogeneously distributed in both diluted samples, otherwise significant Curie-Weiss temperatures should be expected. The difference between the Kramers Yb 3+ and non-Kramers Tm 3+ cases is clearly seen in C m /T , where the signal of the diluted Yb 3+ sample diverges at low temperatures, whereas the diluted Tm 3+ sample reveals a finite zerotemperature value of C m /T . This finite value indicates a distribution of the energy splitting between the two lowest-lying CEF singlets, |E 1 and |E 2 [39]. This transforms two singlets into a quasidoublet and gives rise to the Ising anisotropy [38,39]. A similar single-ion scenario was recently reported for the Ising spin chain compound PrTiNbO 6 [39]. Here, we use the same approach and model the distribution with a Lorentzian function centered at ∆ = E 2 -E 1 and having the full width at half maximum (FWHM) ω. The non-zero ω arises from the site mixing of Mg 2+ and Ga 3+ that, with their different charges, generate random CEF on the rare-earth site. By fitting C m /T of the diluted samples [ Fig. 1 (b)], we find ∆ = 5.9 K and ω = 5.3 K for the Tm 3+ compound to be Tm0.04Lu0.96MgGaO4, YbMgGaO4, and Yb0.04Lu0.96MgGaO4 at 0 T. The red and blue lines show, respectively, the fits to the data for Tm0.04Lu0.96MgGaO4 and Yb0.04Lu0.96MgGaO4 with the Lorentzian distributions of E2-E1. (c) Temperature dependence of susceptibilities measured in the field of 0.05 T applied both parallel and perpendicular to the c axis. The inset shows the CEF levels from the combined CEF fit, with the black and violet lines for the CEF doublets and singlets, respectively. The GS quasidoublet is shown by the black-violet line. (d) Temperature dependence of the magnetic heat capacity measured at 0 T. The inset presents sketch of the 2D three-sublattice magnetic dipole structure with the green and red lines showing the triangular lattice and the magnetic unit cell, respectively. The lines show the combined CEF fit to the magnetic data of TmMgGaO4 above 90 K in both (c) and (d). (e) Magnetic neutron diffraction of TmMgGaO4 measured on D23 at 60 mK and 0 T in the ab plane (L = 0). (f) L dependence of selected static structure factors measured on POLI at 60 mK and 0 T. The magnetic structure factors are normalized by the magnetic form factor of Tm 3+ . compared with ∆ = 0 K and ω = 0.19 K for Yb 3+ , where the GS doublet is protected by time-reversal symmetry. Whereas this protection does not occur in the case of Tm 3+ , a robust GS quasidoublet can still form, because ω is comparable to ∆.
Above 90 K, both the CEF randomness and intersite couplings, with the energy scales of ∼ 10 K [see Fig. 1 (a) and (b), see also below], can be neglected, and the combined CEF fit can be carried out for both magnetic susceptibilities and heat capacity measured on the single crystal of TmMgGaO 4 [see Fig. 1 (c) and (d)] (Appendix C). Fitting thermodynamic data for Tm 0.04 Lu 0.96 MgGaO 4 leads to less accurate results, owing to the large error bar for the specific heat at high temperatures. Moreover, above ∼ 90 K the normalized susceptibilities for TmMgGaO 4 and Tm 0.04 Lu 0.96 MgGaO 4 nearly match, suggesting that intersite couplings play no significant role in this temperature range.
The average inner gap between the two lowest-lying CEF singlets is fitted to be (E 2 -E 1 ) CEF = 6.3 K (Appendix C), very similar to that in Tm 0.04 Lu 0.96 MgGaO 4 , ∆ = 5.9 K (see above). ∆ CEF ∼ 450 K is also obtained for TmMgGaO 4 , which is close to that of YbMgGaO 4 [11,16]. At low temperatures (T ∆ CEF ), the components of the pseudospin-1/2 magnetic moment tensor are calculated as m α ij = µ B g J E i |J α |E j (i, j = 1, 2, and α = x, y, z), where g J = 7/6 is the Landé g factor and J α is the component of the total angular momentum operator. And we obtain the general form of the tensor as under the subspace of |E 1 and |E 2 , where |G| = g CEF ∼ g [see Fig. 1 (a)]. The eigenstates of Eq. (2) are |σ = ± = 1 √ 2 ( G |G| |E 1 ± |E 2 ) with the Ising eigen-moments, m x = m y = 0 and m z ∼ ±µ B g /2. Therefore, in the dipole approximation both |E 1 and |E 2 are nonmagnetic, but their linear superposition |σ = ± become magnetic. Un-der the subspace of |σ = ± , we get Therefore, by resetting (E 2 +E 1 )/2 = 0 K the low-T single-ion CEF term should be taken into account in the effective spin-1/2 Hamiltonian (see below), with the random inner gap ∆ = E 2 -E 1 . Different local environments at the Tm 3+ sites, with different distributions of Mg 2+ /Ga 3+ , give rise to the different CEF parameters, thus leading to different values of ∆ as well as g [16]. Eq. (4) introduces transverse magnetic field term into the spin Hamiltonian [26]. Neutrons have a magnetic moment that is sensitive to the Tm 3+ dipole moments only. Therefore, our data directly evidence that an ideal 2D three-sublattice magnetic component of the dipole moments [see Fig. 1 (d)] forms in the 3D crystal structure of TmMgGaO 4 below 0.7 K and 2.6 T (see below), but the interlayer spin-spin correlations are negligible. A similar conclusion has been derived in Ref. [42] based on the time-of-flight inelastic neutron scattering (INS) data. To the best of our knowledge, experimental examples of ideal 2D magnetic structures in a real 3D material are highly rare to date, as the interlayer magnetic interactions are always present. The negligible interlayer couplings/correlations are likely caused by the extremely large interlayer distance of c/3 ∼ 8.4Å in TmMgGaO 4 [32].
We can't assign unique values or detailed distributions (in the case of randomness) of S z 1 , S z 2 , S z 3 , only based on the conventional magnetic structure refinement, because any three-sublattice structure illustrated in the inset of Fig. 1 (d) with arbitrary S z 1 = S z 2 = S z 3 gives the additional magnetic reflections sharing the equal structure factor, except for small differences [see Fig. 1 (f)] caused by the sample shape, crystal extinction, and similar effects [43]. Here, S z i = GS|S z i |GS and |GS is the GS of the effective spin-1/2 system at ∼ 0 K, which is applicable to the following ED calculations. The neutron diffraction data do not contain enough information for the refinement of the 2D magnetic structure.
To explore the 2D correlated magnetism of TmMgGaO 4 , as well as the detailed GS magnetic structures, one has to turn to the following disordered effective spin-1/2 Hamiltonian on the triangular lattice, with the longitudinal field of H applied along the c axis, Different from the 2D Ising Hamiltonian reported in Ref. [26], we should further consider the second-neighbor interaction in Eq. (5), owing to the large observed g factor, g = 13.18. In the limit of the magnetic dipole-dipole interaction, the average second-neighbor coupling is estimated to be non-negligible, J zz 9Å. The dipolar interaction is relatively long-ranged, and interactions beyond second neighbors can also be envisaged. However, the presence of these further-neighbor interaction terms significantly complicates the calculation and the modelling. In all of the existing references on TmMgGaO 4 , other groups also restrict themselves to the second-neighbor interaction [3,42,44]. Even if the couplings beyond second neighbors are non-negligible, we do not see strong reasons to include them into the effective spin Hamiltonian, in contrast to other perturbations, such as the randomness of the CEF gap (energy scale ∼ 8 K). And the present model of TmMgGaO 4 with only the first-and secondneighbor interactions should be just effective. The real situation may be more complicated, but we seek to explain the bulk of experimental observations within the minimum model that captures the essential physics.
Between 30 and 60 K, S m is a constant of Rln2 [see Fig. 2 (e)], both the CEF excitations to higher levels and the intersite spin-spin correlations have marginal effect [39], and the mean-field approximation of the effective spin-1/2 system, χ = C /(T -θ w ), is applicable. Here, C = N A µ 0 µ 2 ef f /k B is the Curie constant, and the Curie-Weiss temperature of θ w = -3(J zz 1 + J zz 2 )/2 reflects the intersite magnetic couplings on the triangular lattice along the c axis (in the spin space) [11]. Through the Curie-Weiss fit to the susceptibility measured along the c axis between 30 and 60 K, we obtain an effective mo-  ment of µ ef f = g µ B /2 = 6.5(1)µ B and θ w = -16.44 (3) K. And we further get J zz 1 + J zz 2 ∼ -2θ w /3 ∼ 10 K. In the above Curie-Weiss fit, we neglected the small van Vleck susceptibility (the CEF effect to the susceptibility), χ vv < 0.2%χ at T ≤ 60 K [see Fig. 1 (c)]. As g ∼ 2Jg J , |m J =±J dominate the GS CEF quasidoublet, and thus the non-Ising intersite coupling terms should be neglected in TmMgGaO 4 [19,20].

V. FITS TO THERMODYNAMIC DATA
Since only the three-sublattice magnetic structure is observed by neutron diffraction below T c ∼ 0.7 K and below µ 0 H c ∼ 2.6 T, we carry out ED calculations using the 9-site and 12-site clusters with different periodic boundary conditions (PBC). No significant finite-size effects have been observed (Appendix E). We use five different models: In #1, we fix the parameters, ∆ = 9.01 K, J zz 1 = 6.61 K, J zz 2 = 0.30 K, and g = 12.11, reported in Ref. [42] without any distribution, and get R p = 147. In #2, we refine the above four parameters, and obtain ∆ = 5.71(6) K, J zz 1 = 10.9(1) K, J zz 2 = 1.11(2) K, g = 13.6(1), and get the least-R p = 61.7 [45]. In #3 and #4, we further induce the Gaussian and Lorentzian dis-tributions to ∆, g , J zz 1 , J zz 2 , respectively, due to the Mg/Ga site-mixing disorder. Each local chemical environment of Tm 3+ has a definite ∆, and thus a definite g [16]. Moreover, the intersite couplings should also be distributed around their average values due to the CEF randomness, via both f -p virtual electron hopping processes [19,20,46] and magnetic dipole-dipole interactions (∝ g 2 ). Indeed, #3 and #4 fit the thermodynamic data much better, with much smaller least-R p = 20.5 and 17.8, respectively. #3 gives ∆ = 5.66 (6) (1)]. Finally, we also try to fit the thermodynamic data without any intersite couplings (J zz 1 = J zz 2 = 0, in the #0 model), but the quality of the fit is very low with a large least R p = 197 (see Fig. 2). Moreover, the four fitted parameters, ∆ = 10 , are inconsistent with the aforementioned values. Since only the single-ion terms are considered, the finite-size effects are completely excluded in this case. Therefore, we conclude that the low-T physics of TmMgGaO 4 goes well beyond single-ion CEF effects, and the antiferromagnetic intersite couplings (J zz 1 and J zz 2 ) are critically important. Besides the too large observed R p (see Fig. 2), model #1 [42] can't well explain the measured magnetic properties of TmMgGaO 4 for the following reasons: First, the mean-field approximation, J zz 1 +J zz 2 ∼ -2θ w /3, must be fulfilled at high temperatures (30 ≤ T ≤ 60 K). The least-R p fitted results of #2, #3, and #4 obey the above relationship very well. In contrast, model #1 gives J zz 1 + J zz 2 = 6.9 K, much smaller than the reported value of -2θ w /3 = 12.7 K measured at 1 T in Ref. [42], where this model was used. Second, the refined values of ∆ obtained from #2, #3, and #4 models are much closer to 5.9 K measured on Tm 0.04 Lu 0.96 MgGaO 4 [ Fig. 1 (b)], than to ∆ = 9.01 K reported in Ref. [42]. Third, model #1 essentially fails to reproduce the anomalies of the magnetization around µ 0 H c measured at low temperatures [see in 0 T (see below). Therefore, we don't exclude from the fit any T -dependent data around T c . And the deviation between the experimental data and least-R p #3 (#4) calculation is relatively large only in C m around T c [see Fig. 2  Here, we define the order parameter for this phase as where N is the number of the triangular sites. The order parameter is strongly dependent on ∆ and, therefore, continuously distributed. It takes the maximum of 4/9 at ∆ = 0 K, and gradually vanishes at the boundaries, |∆| ∼ δ (see Fig. 3). Here, δ(H ) is strongly dependent on the applied longitudinal magnetic field [see with the intensity of ∼ 11000 (see Appendix E) at |Q| = 3.5276π/a, which is obviously larger than the average value of ∼ 1900 measured at 1.5 T (see Appendix D). Therefore, the fraction of the uud compound can be estimated to be ∼ 1900/11000 ∼ 18%, which is slightly larger than the above value of ∼ 11% obtained from the thermodynamic data. In this case, the fraction should be in a range of ∼ 11−20% at 1.5 T.
The above calculations are based on the hypothesis that local symmetries are preserved and same values of ∆ occur within each cluster. However, the real situation in TmMgGaO 4 is much more complicated. Therefore, we also perform calculations in the presence of the spatially randomly distributed inner gap, ∆. The average value of ∆, ∆, is comparable to its FWHM, according to the magnetic heat capacity measured on the highly diluted sample of Tm 0.04 Lu 0.96 MgGaO 4 (see above). For simplicity, we assume that 1, 2, 3, 2, 1 pseudospins feature ∆ = 0, 0.5∆, ∆, 1.5∆, 2∆, respectively, in order to mimic the Lorentzian distribution. These pseudospins are randomly arranged on the 9-site cluster by Matlab [see Fig. 4 (a)]. Considering the possible size effect, we also performed the ED calculation on the 12-site cluster with PBC. Similarly, we assume that 1, 3, 4, 3, 1 pseudospins feature ∆ = 0, 0.5∆, ∆, 1.5∆, 2∆, respectively. The 12 pseudospins are randomly arranged in the cluster [see Fig. 4 (b)].
The calculations with the randomness inside the cluster largely confirm the main conclusions drawn from the previous calculations where local symmetries were kept. The local order parameter shows a pronounced variation (see Fig. 4). Around the pseudospin with the smaller inner gap, the local order parameter is much larger than that around the pseudospin with the larger inner gap. Moreover, the uud and uniformly polarized components appear in different parts of the cluster depending on the local value of ∆ [see Fig. 4 (b)]. These effects become even more obvious on larger clusters, but do not differ qualitatively from the results for the "homogeneous" clusters with same value of ∆. Therefore, even the model without internal randomness within the cluster should capture the essential physics of TmMgGaO 4 .

VII. PHASE DIAGRAM AND DISCUSSION
Around the critical points, such as T = T c at 0 T and H = H c at the low temperatures, the integral intensities of the magnetic reflections just completely disappear (see Appendix D). The correlation length in the ab-plane (ξ ab ) can be estimated from the intrinsic broadening of the magnetic reflections along H and K. Unlike the conventional long-range magnetic transition where the magnetic Bragg peaks keep coherent at all temperature below T c [40], TmMgGaO 4 shows diffuse magnetic Bragg peaks with short correlation lengths, ξ ab ∼ 200Å around T c and H c [see Fig. 5 (a) and (b)], well consistent with the absence of the sharp λ peaks in the temperature dependence of the magnetic heat capacity at T c . While, at the phase space well below the above critical points the magnetic Bragg peaks become sharp with (quasi-)long correlation lengths of ≥ 1000Å, which is more than two orders of magnitude larger than the lattice constant of the triangular lattice, a = 3.4097Å. Moreover, the longitudinal magnetic field applied up to ∼ 1.5 T along the c axis gradually shifts the critical point to a higher temperature, and the peak of the magnetic heat capacity becomes sharper and sharper [please see C m /T data in Fig. 8 (b) of Appendix B]. These observations are in contradiction to the formation of the conventional short-range spin-glass GS [29,47,48].
To obtain the detailed low-T phase diagram for TmMgGaO 4 , we further measured temperature dependence of the magnetic heat capacity (C m ) at different magnetic fields, as well as the field dependence of C m , the first derivative of magnetization (dM /dH ), and magnetic Grüneisen ratio (Γ m ) (see Fig. 8 in Appendix B). The phase diagram is shown in Fig. 5 (c). Above T c =  Fig. 12 (e) in Appendix D], no magnetic neutron reflections with fractional H and K are observed, suggesting the paramagnetic phase with a large frustration factor, |θ w |/T c ∼ 23. At low temperatures, the spin system of TmMgGaO 4 is fully polarized by high longitudinal fields, along with a very weak van Vleck susceptibility caused by excitations to higher-lying CEF levels. As the applied field decreases, both dM /dH and Γ m show a broad hump at ∼ 3.5 T, indicating a crossover from the fully polarized phase to the uniformly partially polarized phase. No magnetic neutron reflections are observed above µ 0 H c = 2.61(2) T [see Fig. 12 (f) in Appendix D], and both the optimized models #3 and #4 indeed produce δ = 0 K above ∼ 3 T, such that the uud order vanishes.
At µ 0 H c , the emergence of the additional magnetic reflections clearly indicates the field-induced magnetic transition along with δ > 0 K, confirmed by the relatively narrowed peak observed in the field dependence of dM /dH , Γ m , and C m . Narrowed peaks are observed at ∼ 1 K in the temperature dependence of C m /T at applied fields below µ 0 H c , consistent with the phase transition toward the (quasi-)long-range uud order. At µ 0 H ∼ 1.5 T, the measured C m /T peak becomes sharpest and λ-shape at 1.61(7) K, and the magnetic neutron reflections take the maximum intensities with long-range correlations [ξ ab ∼ 4000Å, see Fig. 5 (b)] at 60 mK, which can be well interpreted by #3 and #4 models with the maximum δ(H ) [see Fig. 2 (d) and Fig. 3 (f)]. As µ 0 H further decreases, a broad peak is observed in the field dependence of dM /dH , Γ m , and C m , at ∼ 0.3 T, possibly suggesting some very delicate transition or crossover from the long-range (ξ ab ∼ 4000Å) to quasi-long-range (ξ ab ∼ 1000Å) orders along with the decrease of δ(H ). The above information of the correlation length is also roughly marked in the phase diagram [ Fig. 5 (c)].
Despite the success of the models #3 and #4 in simultaneously reproducing the low-energy thermodynamic properties and magnetic order probed by neutron diffraction, we emphasize that a more sophisticated model would be required to describe all experimental data, including the spin-wave excitations reported in Ref. [42]. First, the asymmetry of the peak is clearly observed in the field dependence of the structure factor on the magnetic reflection measured at 60 mK [see Fig. 2 There is still about 30% of the maximum intensity (at 1.5 T) observed at 0 T, and the structure factor quickly disappears at µ 0 H c = 2.61(2) T. In contrast, both models #3 and #4 give the symmetric peak profile centered at ∼ 1.5 T [see Fig. 2 (d)]. Second, the calculated INS excitations using #3 and #4 indeed get broader, but may still deviate from the reported spin-wave result [42] (see Appendix E). The asymmetric distribution functions of ∆, g , J zz 1 , and J zz 2 by considering the detailed Mg/Ga arrangements, as well as the inherent correlations among these Hamiltonian parameters, would be required to reach a consistent interpretation for all observations.
A similar effective spin-1/2 Ising Hamilitonian with a continuous distribution of the microscopic parameters can be applied to other non-Kramers rare-earth magnets with correlated GS quasidoublets, such as the Pr 3+ effective spin-1/2 chain compound, PrTiNbO 6 , with the similar site-mixing disorder between nonmagnetic Ti 4+ and Nb 5+ ions [39]. Therefore, the correlated magnetism of the GS quasidoublets should be general as well in condensed matter physics, as the structural disorder is usually inevitable in a real material. Our present work paves the road to understanding this kind of novel many-body physics.

VIII. CONCLUSIONS
We performed an extensive single-crystal study on the 2D frustrated magnetism of TmMgGaO 4 , with the perfect triangular lattice of non-Kramers rare-earth Tm 3+ ions. The distribution of two nearly degenerate GS CEF singlets (quasidoublet) caused by the Mg/Ga disorder is clearly evidenced by the magnetic heat capacity of highly diluted Tm 0.04 Lu 0.94 MgGaO 4 with the negligible intersite couplings, as well as the combined CEF fits to the high-T thermodynamic data of TmMgGaO 4 . At low temperatures, the effective spin-1/2 Hamiltonian of the correlated quasidoublets is experimentally determined. It gives rise to the small fraction of the 2D uud phase of the Ising magnetic dipoles with small inner gaps (|∆| ≤ δ), as well as the main nonmagnetic phase at 0 T with large inner gaps (|∆| > δ), which become uniformly polarized at a finite longitudinal applied field of H . Our correlated quasidoublet model naturally explains the strongest magnetic reflections observed at µ 0 H ∼ 1.5 T, as well as the vanishing intensity with increasing or decreasing H . The similar effective spin-1/2 model with a distribution of the microscopic parameters should be applied to other non-Kramers rare-earth magnets with the disorder-induced GS CEF quasidoublets.   Fig. 6) were grown in a high-temperature optical floating zone furnace (FZ-T-10000-H-VI-VPM-PC, Crystal Systems Corp.), using 53.0%, 60.7%, and 60.9% of the full power of the four lamps (the full power is 1.5 kW for each lamp), respectively [11,32,39]. The single crystals were oriented by the Laue x-ray diffraction, and were cut consequently by a line cutter along the crystallographic ab plane. The cut planes were cross-checked by both Laue (see Fig. 6) and conventional x-ray diffraction (see Fig. 6). No significant impurity phase of the TmMgGaO 4 sample was observed by the single-crystal x-ray and neutron diffraction, consistent with the previously reported work [32]. These single-crystal samples are well transparent, and thus we have full confidence in the absence of any impurity phases, also from a visual inspection of the crystals with a microscope. We also show the x-ray diffraction data measured on the TmMgGaO 4 powder in  The dc magnetization (1.8 ≤ T ≤ 400 K and 0 ≤ µ 0 H ≤ 7 T) was measured by a magnetic property measurement system (MPMS, Quantum Design) using single crystals of ∼ 100 mg. The dc magnetization up to 14 T was measured by a vibrating sample magnetometer (VSM) in a physical property measurement system (PPMS, Quantum Design). The heat capacity (1.8 ≤ T ≤ 400 K and 0 ≤ µ 0 H ≤ 12 T) was measured using single crystals of ∼ 10 mg in a PPMS. N-grease was used to facilitate thermal contact between the sample and the puck below 210 K, while H-grease was used above 200 K. The sample coupling was better than 99%. The contributions of the grease and puck under different external fields were measured independently and subtracted from the data. It is very difficult to precisely measure the magnetic heat capacity of Tm 0.04 Lu 0.96 MgGaO 4 and Yb 0.04 Lu 0.96 MgGaO 4 above 10 K, due to the high dilution of the magnetic ions and the inevitable thermal disturbance.
The ac susceptibility (1.8 ≤ T ≤ 30 K and 0 T)  [39]. The two-τ model [39] is applicable in most cases. No signatures of the poor thermal contact between the sample and holder were observed.
In Fig. 7 (c), we show a typical zero-field relaxation curve and its fit with the two-τ model. In applied magnetic fields, the data may deviate from the two-τ model at short times even at higher temperatures [as in Fig. 7  (d)], similar to our previous report on the spin-chain compound PrTiNbO 6 [39]. This slight deviation may be caused by the thermal decoupling between the phonon (lattice) and electronic/nuclear subsystems [49]. Similar to this previous work [39], we chose to exclude the heat capacity data with the adj. R 2 smaller than 0.9995. The magnetic heat capacity (C m ) of TmMgGaO 4 was obtained by subtracting C p of LuMgGaO 4 from C p of TmMgGaO 4 [see Fig. 8 (a)]. We fitted the 0, 0.2, and 0.5 T heat capacities using the function, C n ( 169 ∆/T )+Aexp(-∆/T ), from the lowest temperature up to the temperature of the minimum in C m /T [see Fig. 8 (b)]. Here C n ( 169 ∆/T ) is the nuclear heat capacity expressed by a two-level model, 169 ∆ and ∆ are the nuclear and electronic spin gaps, respectively, and A is a pre-factor [39].
The dc magnetization (M ) of TmMgGaO 4 between 0.024 and 2.0 K at magnetic fields up to 8 T applied along the c axis, was measured by a high-resolution capacitive Faraday force magnetometer in a 3 He-4 He dilution refrigerator [33] [see Fig. 8 (c) for dM /dH ]. The magnetic Grüneisen ratio or magnetocaloric effect, Γ m = (dT /dH)/(µ 0 T ) = -(dM /dT )/C p , was measured by the alternating field technique (ν = 0.02 and 0.04 Hz) in   [34,35] [see Fig. 8 (d)]. Field dependence of the heat capacity was also measured at 0.3 K [see Fig. 8 (e)], where both the nuclear and lattice contributions are negligible. The peak positions correspond to the magnetic field induced transitions or crossovers (see Fig. 8).
rameters that will be determined experimentally, and the Stevens operators O m n are polynomial functions of the components of the total angular momentum operators J z , J + , and J − (J ± = J x ± iJ y ). The eigenvalues and eigenvectors of Eq. (C1) are given by E j and | E j (j = 1−13), respectively. Under an external magnetic field of H along the x-, y-or z-direction (z is along the c axis, and x is along the a axis), the CEF Hamiltonian can be expressed as, with α = x, y, and z respectively. The eigenvalues and eigenvectors of Eq. (C2) are given by E α j and | j, α , respectively. The single-ion dc magnetic susceptibility can be calculated by, , (C3) and the single-ion magnetic heat capacity under 0 T can be calculated by, ular and parallel to the c axis, respectively. Through the combined fit to the high-T magnetic susceptibilities and heat capacity measured above 90 K [see Fig. 1 (c) and (d) in the main text], all of the six CEF parameters (median values), B m n , can be determined experimentally (see Table I). All of the thirteen eigenvalues (the relative values) and eigenvectors of Eq. (C1) are then obtained (see Table II). The resulting GS g tensor naturally features the strict Ising anisotropy, g CEF ⊥ = 0 and g CEF = 12.5 (see main text).
Finally, we checked the single-ion physics by measuring specific heat of the strongly diluted sample, Tm 0.04 Lu 0.96 MgGaO 4 , and observed the finite zerotemperature value, C m /T ∼ 0.65 JK −2 per mol Tm. This indicates the mixing of the two lowest-lying CEF singlets and the formation of a quasidoublet, which renders the Ising anisotropy [38]. . The high-temperature background was measured on the sample at 2 K and 0 T, where the spin system is paramagnetic [see Fig. 11 (a)], and subtracted from the low-T data.
In order to evaluate the correlation length of the threesublattice magnetic order, we choose a broad magnetic  is the intrinsic scattering signal from the three-sublattice magnetic order, and I bckgr , I 0 , ω L , Ω 0 are fitting parameters for the background, integral intensity, intrinsic reflection width, peak center, respectively. We obtained ω L ∼ 0.34 o and 0.08 o at 0 and 1.5 T, respectively. If we fit the magnetic reflections using a single Lorentzian function, FWHM = 0.87(2) o and 0.69 (1) o are obtained at 0 and 1.5 T, respectively, with FWHM < ω L +σ Ω . Therefore, at 60 mK the correlation length of the threesublattice magnetic order can be estimated as, [40,41]. With λ i = 2.364Å and θ = 13.4 o , we obtain ξ ab ∼ 850 and 3800 A at 0 and 1.5 T, respectively. The measured ξ ab is more than two orders of magnitude larger than the lattice parameter, a = 3.4097Å, but it is still much smaller than the crystal size. Along the c axis, the interlayer correlation length, ξ c , can be estimated as, ξ c ∼ 2π/FWHM L < c/12, where FWHM L > 12 2π c is the broadening of the magnetic reflections along L [see Fig. 11 (c) and (d)].
The magnetic field dependence of the intensity on the nuclear reflections measured at ∼ 60 mK up to 5 T is shown in Fig. 11 (e). A rapid increase of the intensity is observed, especially on the reflections of (1,1,0) and (2,-1,0), at ∼ 3 T with increasing the applied field, which seems consistent with the #2, #3, and #4 models [ Fig. 11 (e)]. However, the measured increase of the intensity (with large error bars) up to 5 T seems smaller than the expected values. The increase of the intensity mainly reflects the uniform/bulk magnetization process [ Fig. 2 (c)]. The spin system of TmMgGaO 4 is almost fully polarized at 5 T, and thus the overall increase of the magnetic intensity with integer indexes can't be tiny [see Fig. 11 (e)]. One possible explanation is that the dominant nuclear part may weakly depend on the applied magnetic field through the magnetostriction effect owing to the spin-lattice coupling.
At ∼ 60 mK, the magnetic reflections show the maximum intensities at µ 0 H ∼ 1.5 T, while completely disappear at µ 0 H ∼ 3 T. Well below the critical points, T c = 0.7 K and µ 0 H c = 2.6 T, the magnetic reflections are coherent with a correlation length of ≥ 1000Å [see Fig. 12 (c) and (d)].
We performed the combined fit to the intensities of the magnetic reflections measured at 0.1−3.5 K in the field of 0 T, by sharing the same fitting parameters, T c and β. We fixed T 0 ≡ 0.001 K to ensure the conditional function, whereas A HK0 were the fitted pre-factors for the reflections (H,K,0). Through the combined fit [see Fig. 12 (e)], the critical temperature and exponent, T c = 0.70(5) K and β = 0.103(3), were obtained. Similarly, we also fitted the intensities measured in the fields of 2−5 T applied along the c axis at 60 mK, by sharing the same fitting parameters, H c and β . We  Around (just below) the critical points, the (quasi-)long-range spin order is replaced by the short-range one, and thus the FWHM of the magnetic reflections increases quickly [see Fig. 12 (g) and (h)] [53].
Appendix E: Exact diagonalization calculations and simulations.
Including the precise single-ion and bond disorder effects into the many-body correlated model of TmMgGaO 4 is a challenging problem. For simplicity, we kept all symmetries of the system (space group: R3m), and assumed the distributions Y -Y = K Y (∆-∆) around the average value in both #3 and #4 models (see main text), where Y = J zz 1 , J zz 2 , g is the Hamiltonian parameter, and K Y is the fitting linear parameter proportional to the FWHM of Y . Each set of the Hamiltonian parameters corresponds to one local CEF environment (Mg 2+ /Ga 3+ arrangement), and we perform the ED calculations for each of these sets. To facilitate the calculations, we truncated the distribution function at |∆-∆|/FWHM(∆) ≥ 1.3 and 2.0 with P (∆)/P (∆) ≤ 0.9% and ≤ 5.9%, for #3 (Gaussian) and #4 (Lorentzian) models, respectively, and then normalized the numerical distribution function by γ P (∆ γ ) = 1.
For each set of the Hamiltonian parameters, we calcu-late the magnetization [see Fig. 2 (c) in main text] using where E j (H ) and |j, H are the eigenvalue and eigenstate of Eq. (5) (see main text), after the ED calculation. The dc magnetic susceptibility is obtained as χ = N A M /H [see Fig. 2 (a) in main text]. The zerofield heat capacity [see Fig. 2 (b) in main text] can be calculated as And the static structure factor of the Ising dipole moment [see Fig. 2
We also perform the ED calculation on the 12-site cluster with different PBC [see Fig. 4 (b) for the geometry]. The calculated thermodynamic data are shown in Fig. 13 (f), with the previous 9-site ED result for comparison. Although certain differences are observed at low temperatures, the overall trend is similar. The measured signal to noise ratio (the standard deviation) of the magnetic heat capacity is much larger than that of the magnetization (susceptibility) data due to the technical difference (Fig. 2), and thus the slight difference [ Fig. 13 (f)] in the magnetic heat capacity obtained on different clusters won't significantly affect the final (fitted) result [please see Eq. (1)]. Indeed, our best parameterization (∆ = 5.7 K, J zz 1 = 10.9 K, g = 13.6, J zz 2 = 1.1 K) shows excellent agreement with the theoretical result reported in the recent preprint (after the initial submission of our present work), where the authors used quantum Monte Carlo (QMC) method and arrived at ∆ = 0.54J zz 1 = 6.2 K, J zz 1 = 0.99 meV = 11 K, g = 1.101×12 = 13.2, J zz 2 = 0.05J zz 1 = 0.6 K [44], using the same model (#2). At relatively high temperatures and/or in high longitudinal magnetic fields, the size effect is relatively small. On the other hand, at low temperatures (≤ 1 K) and at ∼ 0 T our ED calculation becomes semi-quantitative. Therefore, we get a relatively large deviation from the QMC result on J zz 2 [44], because this coupling mostly affects the lowenergy part of the spectrum. Therefore, we mainly focus on the low-T (∼ 60 mK) physics of TmMgGaO 4 in magnetic fields ∼ 1.5 T (µ 0 H g µ B /k B ∼ 13 K), where the up-up-down order is most stable and the ED calculations should be accurate enough.
The spin-wave excitations can be calculated by the Spinw-Matlab code based on the linear spin-wave theory [54] [see Fig. 14 (a) -(d)]. With the above code, the #1 model largely reproduces the spin-wave excitation measured at 0 T, which is sensitive to the main nonmagnetic phase with large inner gaps, while it completely fails to explain the thermodynamic properties and magnetic neutron diffraction measured under the longitudinal field (see Fig. 2 in main text). For example, this calculated spin-wave excitations of the #1 model show a full gap of ≥ 0.4 meV, and obviously can't account for the highly enhanced reflections at K points anymore, at ∼ 1.5 T [see Fig. 14 (a)]. Moreover, the measured width of the spin-wave excitation seems much wider than the reported instrumental resolution (σ E = 0.114 meV) at E i = 4.8 meV [42].
Very recently, Ref. [44] pointed out that the linear spin-wave approximation may get invalid in TmMgGaO 4 . Similarly, we also calculate the INS spectra using the ED results as Here, k i and k f are the incident and final neutron wavevectors, and Z(H , T ) is the partition function. By setting E = 0 meV (k i = k f ) and σ E = 0 meV, Eq. (E5) is equivalent to Eq. (E3) (the integral structure factor of the Ω-scan of the neutron diffraction), after the normalization by the magnetic form factor. Due to the size effect of the ED calculation, the resulted resolution of the transfer momentum is very low [see Fig. 14 (e) -(h)]. While, at K points our calculated diffraction intensities using #1 and #2 models and at 0 T by the ED are well consistent with the recently reported QMC results equipped with stochastic analytical continuation [44] [see Fig. 14 (e) and (f), respectively]. Similarly, our ED calculation using the #1 model at K points clearly contradicts with the measured spin-wave excitations [42], as well as the Spinw calculation based on the linear spin-wave approximation [ Fig. 14 (a)]. As the #1 model gives an energy gap of ∼ 0.42 meV at K points similar to that reported in Ref. [44], and much larger than σ E at 0 T. On the other hand, our ED calculations using the random #3 and #4 models well reproduce the observed (quasi-)gapless (gap < σ E ) feature at K points in 0 T [42], as well as the enhanced magnetic reflection intensity at ∼ 1.5 T [see Fig. 14 (g) and (h)]. Therefore, we emphasize the impor-tant ingredient, the distribution of the effective spin-1/2 Hamilitonian parameters, on the correlated magnetism of non-Kramers GS quasidoublets. It is caused by the nonmagnetic Mg/Ga site-mixing disorder, which is expected to be uniformly distributed at the Mg/Ga sites in RMgGaO 4 (R = rare-earth), according to the diffraction and structure refinements [10,11,32]. Interestingly, the coherent magnetic reflections and (quasi-)long-range magnetic order, instead of the short-range spin-glass GS, are observed in TmMgGaO 4 , despite the above randomness. Ref. [3] also reported that the long-range order can survive in a triangular Ising antiferromagnet, in the presence of the uniform bond randomness.