Introduction

A quantum spin liquid (QSL) is an exotic quantum state in which spins are highly entangled and remain disordered down to zero-temperature limit1,2,3,4. It has attracted intensive interest because of its potential relevance to high-temperature superconductivity and quantum information applications. The most fascinating feature of certain gapless QSL is that its fermionic-like spin excitations, or spinons, behave like mobile charge carriers in a paramagnetic metal with a Fermi surface1,2,3,4. Since the first reported gapless QSL candidate EtMe3Sb[Pd(dmit)2]25, a spin-1/2 triangular lattice antiferromagnet (TAF) in 2010, the search for other candidates remains as a hot topic during the last decade. One commonly used experimental probe for QSL is a broad continuous magnetic excitation, or continuum mode, in the inelastic neutron scattering (INS) spectrum6,7,8.

Meanwhile, the community recently started to pay attention to the chemical disorder effects on quantum magnetism. For example, some recent studies proposed that the randomness in a quantum magnet can induce spin-singlet dimers of varying strengths in a spatially random manner, which can account for the continuum mode due to its widely distributed binding energy. Indeed, the disorder is unavoidable in most of the studied gapless QSL candidates. For instance, the kagome lattice herbertsmithite ZnCu3(OH)6Cl2 has Zn2+/Cu2+ site mixture9 (whether there is a small gap in this QSL candidate is still under debate); the LiZn2Mo3O8 with breathing kagome lattice has Li+/Zn2+site mixture10; the Ca10Cr7O28 with bilayer kagome lattice has disorder among the two different Cr3+positions11,12, and the H3LiIr2O6 with honeycomb lattice has mobile Hydrogen ions13. Therefore, this so-called random-singlet (RS) or valence bond glass state14,15,16,17,18,19,20,21 seriously prompts re-consideration of the intrinsic magnetic ground state of them: whether they are true gapless QSL or just spin liquid-like?

YbMgGaO4 (YMGO), a TAF with the effective spin-1/2 Yb3+ ions22,23, is at the center of this controversy. On one hand, the observed continuum mode7,8,24,25, the Cp ~ T0.7 behavior for the specific heat22, the temperature-independent plateau for the Muon spin relaxation (MuSR) rate26, and the saturated DC susceptibility below 0.1–0.2 K27 all suggest a gapless QSL state. Its origin has been interpreted as a U(1) QSL state with spinon Fermi surface8,25,27,28, or a resonant valence bond-like state24, or a J1–J2 exchange interactions resulted in QSL7,29. On the other hand, no residual κ0/T term on the thermal conductivity has been observed30. It needs to be emphasized that since the itinerant spin excitations, spinons, carry entropy5,31,32,33, a nonzero residual linear κ0/T term at ultralow temperature for thermal conductivity should exist. It is a well-recognized signature for spinon. The absence of the κ0/T term in YMGO makes it difficult to reconcile with spinon physics. The frequency-dependent AC susceptibility peak34 further suggests that the disordered occupancy of the inter-triangular layers by Mg2+and Ga3+ions23,35 lead to a glassy ground state. By including this disorder, other types of ground states have been proposed, such as the mimicry of a spin liquid36,37, RS led spin-liquid-like behavior16,38, and randomness induced spin-liquid-like phase in J1–J2 model39.

In this paper, with the motivation to settle down this dispute, which is certainly significant for QSL studies, we performed thermal conductivity, specific heat, magnetic torque, DC magnetization, and AC susceptibility measurements on high-quality YMGO single crystals at ultralow temperatures with two crystallographic axes and detailed field scans to study its intrinsic magnetic ground state. We observed a residual κ0/T term and series of quantum spin state transitions (QSSTs) in the zero-temperature limit, which strongly suggests that a QSL state with itinerant excitations and quantum spin fluctuations survives disorder in YMGO.

Results and discussions

Thermal conductivity

Figure 1a shows the ultralow-temperature thermal conductivity (κ) of YMGO measured along either the a or c-axis, plotted with κ/T vs. T. Apparently, the κa displays a T2 behavior in a rather broad temperature range (200–600 mK) with a slope change at T < 200 mK. As shown in Fig. 1a, the κc also exhibits a T2 behavior at T < 600 mK with κc/T = 0 for zero-temperature limit. Since YMGO is a two-dimensional spin system with negligibly weak spin interaction along the c-axis, the κc should represent a pure phonon heat transport of this material. This T2 behavior for phonon heat transport is further verified by the high-magnetic-field data of κa. As shown in Fig. 1a, b, with applying 14 T (or 10 T) field along the a-axis (or the c-axis), the κa displays a good T2 behavior at T < 700 mK accompanied with κa/T = 0 for zero-temperature limit. Since these fields are high enough to polarize all the spins, below the magnon gap, the high-field thermal conductivity should be a purely phononic term without any contribution (carrying heat or scattering phonon) from magnetic excitations. It is therefore reasonable that in high fields there is no residual term of κa/T at T → 0. Based on the above comparisons, it is obvious that in zero field the κa behaves as T2 at 200 mK < T < 600 mK with a residual thermal conductivity of κ0/T = 0.0058 W K−2 m−1 and thereafter, the slope change leads to a smaller κ0/T = 0.0016 W K−2 m−1 at T < 200 mK. Moreover, the larger κa/T in high fields indicates that, in zero magnetic field, the phonons are rather strongly scattered by magnetic excitations that are be gapped out in high magnetic fields and suppressed at low temperatures. It firmly suggests the presence of magnetic excitations in the thermal conductivity result at zero magnetic field, and the magnetic excitation certainly does not have a large gap and could be gapless.

Fig. 1: Ultralow-temperature thermal conductivity of YbMgGaO4 single crystals.
figure 1

a Data measured along the a-axis (κa) and the c-axis (κc), plotted as κ/T vs. T2. The solid lines are some linear fitting results. The κc display a T2 behavior at T < 600 mK with zero intercepts at T = 0. The κa behave as T2 at 200 mK < T < 600 mK with a residual thermal conductivity of κ0/T = 0.0058 W K−2 m−1 and a smaller κ0/T = 0.0016 W K−2 m−1 at T < 200 mK. b The a-axis thermal conductivity in zero field and in 10 T (14 T) magnetic field applied along the a- (the c-) axis. The high-field κc data display a T2 behavior at T < 700 mK with κ0/T = 0.

Full size image

It is a bit strange that the phonon thermal conductivity at ultralow temperature has a T2 behavior. Usually, the phonon thermal conductivity of a high-quality insulating crystal at the boundary scattering limit should have a T3 temperature dependence. One possible explanation is the surface reflection effect, which could weaken the temperature dependence and gives a T-power-law behavior with a power smaller than 3. However, the calculated phonon means free path (see Supplementary Information) is much smaller than the sample width, at least at several hundreds of millikelvins. Therefore, the surface reflection may not be the origin of the T2 behavior. Here, we would like to leave it as an open question.

The nonzero residual thermal conductivity at zero field immediately implies that YMGO is a QSL with itinerant gapless spin excitations having a long-range algebraic (power-law) temperature dependence. In reality, it is very rare to observe this nonzero κ0/T term in the studied QSL candidates. So far only the organic EtMe3Sb[Pd(dmit)2]25,33 and the inorganic 1T-TaS240 exhibit a nonzero κ0/T term, both of which are spin-1/2 TAFs. But it is also notable that some recent studies reported a zero κ0/T term in these two materials41,42. By the following ref. 5’s method, we estimate the mean free path (ls) of the spin excitations in YMGO by calculating

$$\frac{{\kappa }_{0}}{T}=\frac{\pi {k}_{B}^{2}}{9{{\hslash }}}\frac{{l}_{s}}{{ad}}=\frac{\pi }{9}{\left(\frac{{k}_{B}}{{{\hslash }}}\right)}^{2}\frac{J}{d}{\tau }_{s}$$
(1)

Here, a (~3.40 Å) and d (~25.1 Å) are the in-plane and out-of-plane lattice constants, respectively. From the observed κ0/T = 0.0058 W K−2 m−1, the ls is obtained as 78.4 Å, indicating that the excitations are mobile to a distance 23 times as long as the inter-spin distance without being scattered. In comparison, the κ0/T value of YMGO is ~30 times smaller than that of EtMe3Sb[Pd(dmit)2]2 (~0.19 W K−2 m−1), and ~10 times smaller than that of 1T-TaS2 (~0.05 W K−2 m−1), which may be attributed to two possible reasons. First, the spinon velocity is much smaller in YMGO due to the small J value. Second, the ls is much shorter in YMGO than that in EtMe3Sb[Pd(dmit)2]2 (~1.20 μm)5, which is related to stronger scattering between phonon and spinon. Note that although 1T-TaS2 has a very large J value40, its ls (~50 Å) is also not so long, which may be due to the spin–phonon scattering.

Accordingly, the smaller κ0/T = 0.0016 W K−2 m−1 below 200 mK leads to a ls ~21.6 Å, which still covers 6 times of the inter-spin distance. While the exact origin for this slope change (or reduction of the κ0/T) is not clear, we suspect it should be related to the effect of disorder. Here, 200 mK is comparable to the AC susceptibility peak position observed at 80 mK. This result strongly suggests that despite the heavily chemical disorder on the Mg/Ga sites, the itinerant spin excitations of YMGO survive when approaching zero temperature.

It is necessary to point out that in an earlier study, a T2 behavior of κ measured in the ab plane at T < 300 mK with a negative intercept at T = 0 was observed and was interpreted as the non-existence of magnetic heat transport30. However, if one compares those data with our κa data, it can be found that those data also exhibit a slope change at ~300 mK; the data above 300 mK exhibit a T2 behavior with a nonzero κ0/T = 0.0018 W K−2 m−1. In the present work, first, we measured the thermal conductivity along both the a- and c-axis. As discussed above, the comparison between the κa and κc clearly show the existence of κ0/T term for κa. Second, the phonon mean free path of our sample is much larger than that of ref. 30 (see Fig. S2 in Supplementary Information). Therefore, it is very clear that our samples display better thermal conductivity, indicating higher sample quality, and should exhibit more intrinsic physical properties of YMGO. Nevertheless, the ultralow-temperature thermal conductivity data from different groups all indicate that at temperatures above 200 or 300 mK, there are itinerant spin excitations associated with the QSL state, while the disorder effect suppresses the transport of spin excitations at lower temperatures.

Figure 2 shows the magnetic-field dependence of κa at various temperatures and with B // a or B // c. First, the κa is significantly enhanced at high magnetic fields, indicating the existence of magnetic scattering of phonons that can be smeared out in high fields. Second, the κa exhibits several structures at the low-field region. For B // a, the κa isotherm measured at 92 mK shows a minimum at Ba1 = 0.5 T and another anomaly at Ba3 ~ 3 T. For B // c, two more obvious minima are observed for the 92 mK data at Bc1 = 0.5 T and Bc2 = 1.6 T, respectively. These structures under magnetic fields of both directions disappear gradually with increasing temperature.

Fig. 2: Magnetic field dependence of the a-axis thermal conductivity of YbMgGaO4.
figure 2

a, b Magnetic field dependence of κa at various temperatures with field-applied along either the a-axis or the c-axis. c, d A magnified view of κa(B) data at low fields. At 92 mK and with B // a, there is a minimum at 0.5 T and a slope change at 3 T, indicated by arrows. While at 92 mK and with B // c, there are two minima at 0.5 and 1.6 T. With increasing temperature, these anomalous minima become weaker and disappear above 850 mK.

Full size image

Specific heat

Figure 3 shows the magnetic-field dependence of specific heat (Cp) at 300 mK with B // a or B // c. For B // a, the Cp isotherm shows an obvious slope change at Ba3 ~ 3.0 T. Its derivative additionally shows another change around Ba2 = 1.5 T, which corresponds to a weaker slope change of Cp. For B // c, one clear slope change occurs at Bc2 = 1.7 T, which is also highlighted as a peak on its derivative. These features, both shown by the κ and Cp, pointing to some field-induced crossovers or transitions, or some special field-dependent magnetic scatterings.

Fig. 3: Low-temperature specific heat of YbMgGaO4.
figure 3

a, b The magnetic-field dependence of specific heat, Cp, measured at 300 mK and with B // a or B //c. c, d The derivative (dCp/dB) for B // a or B //c. The dashed lines are the guide for eyes and indicate the anomalies associated with some transition fields Ba2, Ba3, and Bc2.

Full size image

Magnetic torque

To further investigate the possible transitions, magnetic torque (τ) measurements were performed. Figure 4a, b shows the calculated τ/B with applied field along the a- or c-axis. At 30 mK, the data shows a weak hump around 1–2 T for both directions. In principle, torque magnetometry measures magnetic anisotropy. Accordingly, the d(τ/B)/dB for B // a (Fig. 4c) clearly shows a valley at Ba2 = 1.2 T, a peak at Ba3 = 2.7 T, and another valley at Ba4 = 7.0 T. With increasing temperature, the lower field valley, and peak both become weak and disappear with T ≥ 1.5 K. For B // c, the d(τ/B)/dB (Fig. 4d) does not show a peak but a flat regime starting around Bc2 = 1.7 T, and a valley at Bc3 = 5.0 T. The Ba2, Ba3, and Bc2 observed here are closely corresponded to the anomaly observed from the κ(B) and Cp(B) data. This further confirms the existence of field-induced magnetic crossovers or transitions. It is noticed that the magnetic torque was also measured for YMGO at 350 mK by Steinhardt et al.43 Its derivative shows broad maxima near 3.5 T for B // c and 5.5 T for B // a, which is different from our results obtained at 30 mK while approaching zero temperature.

Fig. 4: Ultralow-temperature magnetization of YbMgGaO4.
figure 4

a, b The calculated magnetization (M) from the measured magnetic torque as torque/B with field-applied along either the a-axis or the c-axis (see Supplementary Fig. S1). c, d The derivative (dM/dB) for B // a- or c-axis. At 30 mK, the derivative shows two anomalies at Ba2 and Ba3 for B // a, which are related to 1/3Ms and \(\sqrt{3}/3\)Ms, respectively. Meanwhile, only one anomaly at Bc2 was observed from the derivative for B //c, which is related to 1/2Ms.

Full size image

It is noticed that a recent study reported the DC magnetization measured at 500 mK with B // ab plane44 or B // c. We also measured the 500 mK magnetization with B // a or B // c. Our results are the same as the reported ones, as shown in Fig. S4 in Supplementary Information. For B // a, the saturation field Bs is around 7.0 T with saturation moment Ms = 1.45 μB. For B // c, Bs = 5.0 T and Ms = 1.8 μB. Therefore, the Ba4 and Bc4 represent the saturation fields. It is also noted that the magnetization at Ba2 = 1.5 T, Ba3 = 3.0 T, and Bc2  = 1.7 T is around 1/3, \(\sqrt{3}/3\), and 1/2 of the saturation value, respectively, with the assumption that the critical field values are similar between 300 and 500 mK. In ref. 44, they further measured magnetization at 40 mK with B // c. Its derivative shows a flat regime around 1.8 T, which is related to 1/2Ms. This result is consistent with our observation here. For their data measured at 500 mK with B // ab plane, its derivative shows a flat regime around 3.5 T, which again was ascribed to 1/2Ms. This explanation is different from our observation for B // a case, in which both phase boundaries at 1/3Ms and \(\sqrt{3}/3\)Ms were observed by the combination of κ(B), Cp(B), and τ(B) data. One reason for this discrepancy could be the fact that the DC magnetization was measured at 500 mK, a relative high temperature. As shown in our data, Fig. 4c, the anomalies on the d(τ/B)/dB curve at 30 mK, an ultralow temperature, become weak or flattened pretty quickly with increasing temperature. The Cp(B) data at 200 mK was also reported in ref. 44. However, its Cp(B) data with B // ab plane shows no distinct feature around 3.0 T, which our B // a data clearly shows. The slope change around 1.5 T, as we observed, also has not been discussed in ref. 44. Meanwhile, we also noticed that an early study on the DC magnetization measured on powder sample at 500 mK reported two critical fields at 1.6 and 2.8 T22, which were suggested as the phase boundaries of an up–up–down (UUD) phase. This study supports our observation with B // a.

Phase diagram

Two magnetic phase diagrams are constructed by using these critical fields, as shown in Fig. 5. Since no anomaly was observed on τ/B around 0.2–0.5 T, the Ba1 and Bc1 are unlikely related to spin-state transitions. Accordingly, besides the paramagnetic phase at high temperatures, QSL at low temperature and zero field, and fully polarized phase at high fields, with increasing field, there are three phases for B // a (Fig. 5a) and two phases for B // c (Fig. 5b).

Fig. 5: Magnetic phase diagrams of YbMgGaO4.
figure 5

a For B // a. b For B // c. The data points are obtained from the magnetic torque (τ) and field dependence of thermal conductivity (κ) measurements. The dashed lines are phase boundaries. For B // a, there are three phases (I: thee canted 120° spin structure, II: the UUD phase, and III: the oblique phase) in the low-temperature and low-field region. Whereas, there are two phases (I: unknown phase and II: UUUD phase) at low temperature and low field for B // c. The dashed lines are the guide to the eyes.

Full size image

For a spin-1/2 TAF with 120° spin ordering and easy-plane anisotropy, the strong quantum spin fluctuations can stabilize a so-called UUD phase within a certain magnetic field regime applied in the triangular plane while approaching zero temperature45,46,47,48,49, which has been observed in Ba3CoSb2O950,51,52,53,54, a TAF with effective spin-1/2 Co2+ ions. This UUD phase exhibits a magnetization plateau with 1/3 saturated moment (Ms). Recently, even in QSL candidates AYbCh2 (A = Na and Cs, Ch = O, S, Se), one TAF family with effective spin-1/2 Yb3+ ions and easy-plane anisotropy, the UUD phase has been proposed under applied fields55,56,57,58. Experimentally, more complex QSSTs have been observed for Ba3CoSb2O959,60,61,62,63,64,65. For example, with increasing field along the ab plane, its 120° spin structure at zero field is followed by a canted 120° spin structure; the UUD phase; a coplanar 2:1 canted oblique phase with one spin in the 120° spin structure rotated to be parallel with another spin, which gives \(\sqrt{3}/3\)Ms; and another coplanar phase before entering the fully polarized state.

Since YMGO also has the easy-plane anisotropy7,8,29, its field-induced anomaly at ultralow temperatures could also be related to these QSSTs. By comparing YMGO’s phase diagram with B // a to that of Ba3CoSb2O9, we propose a phase I as the canted 120° spin structure, phase II as the UUD phase, and phase III as the oblique phase. The observation of the QSSTs strongly supports the existence of quantum spin fluctuations approaching zero temperature in YMGO, which clearly differentiate it from a conventional spin glass state. It needs to be pointed out is that while the phase boundaries of the UUD phase were observed, the 1/3Ms plateau does not exist for YMGO, which could be due to the disorder effect. In another TAF Rb1−xKxFe(MoO4)266, the disorder introduced by the K-doping also weakens the magnetization plateau feature related to the UUD phase.

It is surprising to see that the phase boundary between phase I and II for B // c case is related to 1/2 Ms. As learned from Ba3CoSb2O9, while for B // c, the 120° spin structure will be followed by an umbrella spin structure and an oblique phase, between which the phase boundary is related to \(\sqrt{3}/3\)Ms. Meanwhile, Ye and Chubukov calculated the phase diagram of a 2D isotropic triangular Heisenberg antiferromagnet in a magnetic field and predicted a novel up–up–up–down (UUUD) spin sate with a 1/2Ms magnetization plateau for J2/J1 > 0.125 (J1: nearest-neighbor exchange interaction, J2: next nearest-neighbor exchange interaction)67. The reported J2/J1 values by simulating the spin excitations obtained from INS data7 and terahertz spectroscopy data29 are 0.22 and 0.18, respectively. Therefore, we tend to assign phase II as the UUUD phase, while the nature of Phase I is not clear at this stage. Again, the chemical disorder could be the reason for the disappearance of the 1/2Ms plateau.

In the case that the minimum observed at 0.5 T for κ(B) belong to some special magnetic scatterings, the scenario of a spinon Fermi surface QSL, that supports gapless magnetic excitations discussed above, may give a possible understanding from the conventional wisdom of Kohn anomaly. Since the charge degrees of freedom in YMGO are frozen out, only Zeeman coupling can be included as the magnetic field is applied. In the weak-field regime, the field does not destroy the QSL ground state and the spinon remains to be a valid description of the magnetic excitation28. The magnetic fields would modify the spinon Fermi surface, and the Fermi surfaces of spin-up and spin-down spinons expand and shrink with increasing field, respectively. Just like the electron-phonon coupling68, the spin–lattice coupling may enhance the phonon scatterings with modified spinon Fermi surfaces, resulting in the thermal conductivity modulation. In the high-temperature regime, the QSL breaks down and the structures disappear as in the experiment.

AC susceptibility

Finally, we revisited the AC susceptibility (χ′) to check the possibility for a spin glass state in YMGO, although the observation of the residual κ0/T and QSSTs clearly disputes this scenario. Compared to the reported χ′ data performed at zero DC magnetic field and without anisotropic information, our χ′ data was measured with AC field both along the a- and c-axis, and with applied DC field. Figure 6a, b shows the χ′(T) measured with applied DC field along either the a- or c-axis. With B = 0 T, the χ′ shows a peak around 80 mK, which is lower than the reported data exhibiting a peak around 100 mK34. With increasing B, the intensity of χ′ below 80 mK increases and eventually becomes flat or saturated with B ≥ 0.05 T. Figure 6c, d shows the frequency dependence of the peak position with the AC excitation field either along the a- or c-axis. For both of them, the peak’s position (T0) shifts to higher temperatures with increasing frequency. As shown in Fig. 6e, f, its frequency dependence can be fit to an Arrhenius law f = f0 exp[−E/(kBT0)], which yields an activation energy Ea = 3.8(6) K and Ec = 2.5(8) K for the excitation field along the a- and c-axis, respectively. While the frequency-dependent peak of χ′ normally represents a spin glass transition as discussed for YMGO and YbZnGaO434, the saturation of the χ′ below this peak position under a small DC field and the anisotropic activation energy both indicate that it should not be treated as a conventional spin glass69.

Fig. 6: AC susceptibility of YbMgGaO4.
figure 6

a, b Temperature dependence of the real part of the AC susceptibility, χ′, measured with different DC magnetic fields along either the a-axis or the c-axis. Here, the maximum value of each data was scaled to 1 to clearly show the DC magnetic field effects on AC susceptibility. c, d The Frequency dependence of the χ′ peak. Here, the absolute value of the AC susceptibility was obtained by rescaling it to the DC susceptibility (see Supplementary Information). e, f Arrhenius law fit of the χ′ peak position, T0. The applied AC excitation field (~1 Oe) is along either the a-axis for (c, e) or the c-axis for (d, f).

Full size image

Meanwhile, the recent DC susceptibility measurements with B = 0.01–0.05 T for YMGO70 also revealed a saturation below 100–200 mK, which has been suggested as a signature for the presence of gapless spin excitations. Since the temperature dependence of the AC and DC susceptibility should behave similarly, the saturation we observed for χ′ could be due to the same origin. In fact, the reported INS measurement suggests that at most 16(3)% of the total spectral weight is elastic7. Correspondingly, the inelastic contribution (~84(3)%) is large compared with the 66% expected for a spin-1/2 glass71. Moreover, the DC susceptibility shows no divergence between the zero-field cooling and field cooling data, the μSR relaxation rate shows a plateau below 400 mK26, and the magnetic entropy has been already released by more than 99% at 80 mK from the specific heat measurement22. All these facts again rule out a frozen or glassy state for YMGO. This AC susceptibility peak could originate from the free impurity spins inside or attached to the system.

Conclusion

In summary, the observation of a residual κ0/T term and QSSTs strongly supports a QSL state with itinerant spin excitations and quantum spin fluctuations approaching zero temperature in YMGO, although its chemical disorder reduces the mean free path of the excitations and smears out the 1/3Ms plateau related to the UUD phase. This survival of QSL state in YMGO is surprising since the Mg/Ga site disorder is supposed to introduce random distribution of exchange interactions and therefore lead to an RS state or even a glassy state and any field-induced transitions should be expected to smear out completely. Therefore, the chemical disorder in YMGO must play a more complex role in the exchange interactions. Future studies on the local structure of YMGO to learn how exactly the Mg/Ga disorder affects the Yb–O local environment and the correlated exchange interactions are highly desirable to understand this survival of QSL in YMGO through chemical disorder.

Methods

Sample preparation and characterization

High-quality YMGO single crystals were grown by using the optical floating-zone technique8. By using the X-ray Laue photograph, the crystals were cut precisely along the crystallographic axes for the magnetic and thermal conductivity measurements.

Thermal conductivity measurements

Thermal conductivity was measured by using a “one heater, two thermometers” technique in a 3He/4He dilution refrigerator at 70 mK–1 K, equipped with a 14 T superconducting magnet72,73. The sample was cut precisely along the crystallographic axes with the longest and the shortest dimensions are along the a- or the c-axis. The magnetic fields were applied along either the a- or c-axis. Gold paint was used to make four contacts on each sample. The RuO2 chip resistors were used as heaters and thermometers and are connected to the gold contacts by using gold wires and silver epoxy. The RuO2 thermometers were pre-calibrated by using a RuO2 reference sensor (Lakeshore Cryotronics) mounted at the mixture chamber (the superconducting magnet was equipped with a cancellation coil at the height of the mixture chamber).

AC susceptibility measurements

The ac susceptibility was measured using the conventional mutual inductance technique (with a combination of ac current source and a lock-in amplifier) at SCM1 dilution fridge magnet of the National High Magnetic Field Laboratory, Tallahassee. The typical AC field strength is 1.1–1.6 Oe.

Torque measurements

Torque magnetometry was performed at the University of Michigan using a self-built capacitive cantilever setup mounted inside a dilution refrigerator74,75,76. The sample was mounted on a thin BeCu cantilever. To infer the magnetic torque \(\tau\), the cantilever deflection is measured by tracking the change of the capacitance between the cantilever and a gold thin film underneath. The report magnetization is the effective torque magnetization M = τ/B.

Specific heat measurements

The specific heat was measured with a Quantum Design physical property measurements system, equipped with a dilution refrigerator insert.

DC magnetization measurements

The DC magnetization curves were measured with a Quantum Design SQUID-VSM, equipped with a 3He refrigerator insert.