Persistent spin dynamics and absence of spin freezing in the $H$-$T$ phase diagram of the 2D triangular antiferromagnet YbMgGaO₄

Plate-like high-quality single crystals of YbMgGaO₄ were grown by a floating-zone method and characterized as reported previously Ma et al. (2018). Figure 1

shows two-thirds of a hexagonal unit cell of YbMgGaO₄ (space group $R\bar{3}m$, No. 166) Li et al. (2015a, b). Mg and Ga ions are distributed randomly on the indicated sites, leading to unavoidable local disorder in the structure. There are two inequivalent oxygen locations (O1 and O2) in the unit cell. The triangular 2D Yb³⁺ layers are well separated by two Mg/Ga planes and four O1/O2 planes (Yb³⁺ layer spacing 8.377 Å), and inter-layer magnetic interaction is weak.

Implanted positive muons in oxides are likely to stop near O^2- ions. To date there are no reported calculations of $\mu^{+}$ site locations in YbMgGaO₄ [Forareview; see]Bonfa16. In lieu of such information, we consider candidate sites that are centered with respect to oxygen nearest neighbors. There are two such crystallographically inequivalent sites ($\mu$1 and $\mu$2), shown in Fig. 1. Dipolar magnetic fields from neighboring Yb³⁺ ions and nuclei are unlikely to be very different if the actual $\mu^{+}$ sites are closer to O^2- ions.

In a typical $\mu$SR experiment Schenck (1985); *Blundell99; *Brewer03; *Yaouanc11, $\sim$100% spin-polarized positive muons are implanted into the sample and stop at interstitial sites. Each $\mu^{+}$ precesses in the (in general fluctuating) local field at its site, and at the time of its decay (mean lifetime ${\sim}2.2~{}\mu$s) emits a positron preferentially in the direction of the $\mu^{+}$ spin. Positrons are detected by an array of scintillation counters, and a histogram $N(t)$ of count rate vs time after implantation is obtained from each counter.

For two identical counters on opposite sides of the sample, the $\mu$SR asymmetry time spectrum $A(t)=[N_{1}(t)-N_{2}(t)]/[N_{1}(t)+N_{2}(t)]$ is proportional to the component $P_{\mu}(t)$ of the ensemble $\mu^{+}$ spin polarization $\mathbf{P}_{\mu}(t)$ along the axis joining the counters and the sample. Corrections for differences in counter geometries and efficiencies are straightforward Schenck (1985); Blundell (1999); Brewer (2003); Yaouanc and Dalmas de Réotier (2011). Unlike its cousin NMR, which, apart from NQR and field-cycling techniques, requires a strong applied magnetic field to polarize the nuclear spins, the polarized $\mu^{+}$ beam allows $\mu$SR to be carried out in arbitrary fields.

Experiments may be considered in three classes:

• •

  In high TF-$\mu$SR, the precession frequency in a paramagnet is shifted from its value in vacuo. The fractional frequency shift $K$ (often called the Knight shift after its discoverer Knight (1949)) and its inhomogeneous distribution (often the main contribution to the linewidth) provide a microscopic characterization of the field-induced static magnetization.

• •

  In sufficiently strong LF-$\mu$SR, the resultant of applied and any local static fields is essentially parallel to the initial $\mu^{+}$ spin orientation. This results in “decoupling” of the static fields, and the remaining relaxation is purely dynamic in origin Hayano et al. (1979).

• •

  In ZF- and weak LF-$\mu$SR, the time evolution of $P_{\mu}(t)$ is due to either or both of two mechanisms: (1) decay due to dynamic fluctuations of the local fields, and (2) $\mu^{+}$ spin precession in the resultant of a distribution of (quasi)static local fields and the applied field, if any (static Kubo-Toyabe (KT) relaxation Hayano et al. (1979)). The latter are normally due to nuclear dipolar fields.

  This situation can be modeled by a dynamically-damped static KT function $G_{\mathrm{KT}}(t)$:

  $$P_{\mu}(t)=G_{d}(t)\,G_{\mathrm{KT}}(t)\,,$$

  (1)

  where the dynamic damping function $G_{d}(t)$ is often a simple exponential but can be more complex (as in the present work). Choices of $G_{\mathrm{KT}}(t)$ are discussed below in Sec. III.2.2. It is, however, often difficult to disentangle static and dynamic relaxation contributions from ZF-$\mu$SR data alone; this is one reason why strong LF-$\mu$SR is valuable.

A mosaic of YbMgGaO₄ single crystals (total mass $\approx 0.3$ g) was mounted on a pure silver sample holder using diluted GE varnish, with the crystal $c$ axes parallel to the $\mu^{+}$ beam direction. $\mu$SR experiments were carried out using the dilution refrigerator spectrometer at the M15 beam line of TRIUMF, Vancouver, Canada. Field zeroing to ${\lesssim}10$ mOe was achieved using a variation of the method of Morris and Heffner Morris and Heffner (2003).

Data were taken over the temperature range 22 mK–8 K. Due to the large low-temperature specific heat of YbMgGaO₄ Xu et al. (2016), temperatures were stabilized for 15 minutes before taking data. ZF- and LF-$\mu$SR data were taken with the initial $\mu^{+}$ ensemble spin polarization parallel to the crystal $c$ axes. For TF-$\mu$SR the $\mu^{+}$ spins were rotated by ${\sim}90^{\circ}$ before implantation. Data were analyzed using the Paul Scherrer Institute musrfit fitting program Suter and Wojek (2012) and the TRIUMF physica programming environment ¹¹1http://computing.triumf.ca/legacy/physica/.

TF-$\mu$SR experiments were carried out in YbMgGaO₄ over the temperature range 26 mK–0.825 K in a transverse applied field $H_{T}=10$ kOe. Representative Fourier transforms of TF asymmetry time spectra are shown in Fig. 2.

They exhibit three lines: line 1 [Fig. 2(a)], which is very broad and shifted to higher frequency, and lines 2 and 3 [Fig. 2(b)], which are narrow and overlap considerably.

The original TF asymmetry time spectra $A^{\mathrm{TF}}(t)$ (not shown) were therefore fit by a sum of three Gaussian-damped oscillating components:

where a common phase $\varphi$ was assumed. Parameter values for $T=26$ mK are listed in Table 1.

Amplitudes $A_{1}^{\mathrm{TF}}$ and $A_{2}^{\mathrm{TF}}$ from the fits are nearly equal and considerably larger than $A_{3}^{\mathrm{TF}}$. The curves in Fig. 2 are Fourier transforms of the Gaussian asymmetry terms in Eq. (2), with parameters from the fits: frequencies $\omega_{i}/2\pi$, widths $\sigma_{i}/2\pi$, and areas proportional to $A_{i}$. Their sum is normalized to the experimental Fourier transform in Fig. 2(b).

One of the three lines arises from muons that miss the sample and stop in the silver sample holder, and the other two are from the sample. Line 3, which is the narrowest, is assumed to be the “silver” signal, and $\omega_{3}/2\pi$ is used as the reference frequency for the Knight shifts of the other two components. Additional evidence for this attribution is discussed below in Sec. III.2.1. The sample lines are numbered in order of their shifts (Fig. 2), consistent with the $\mu^{+}$ site numbering in order of their distance from the nearest Yb³⁺ layer (Fig. 1).

In order for Knight shifts to reflect local magnetism, a correction must be made for the contributions to the raw shift $K_{\mathrm{raw}}$ from macroscopic Lorentz and demagnetizing fields  Carter et al. (1977); *Akashin92:

where $A_{\mathrm{Ld}}=4\pi\left(\textstyle{\frac{1}{3}}-D\right)/v_{\mathrm{mol}}$, $\chi_{\mathrm{mol}}$ is the molar bulk susceptibility, $D$ is the sample demagnetization coefficient, and $v_{\mathrm{mol}}$ is the molar volume (50.6 cm³/mol for YbMgGaO₄). For field normal to the thin mosaic samples, $D$ is estimated to be 0.8(1), so that $A_{\mathrm{Ld}}=-0.12(2)$ mol emu^-1. This correction can be important in $\mu$SR shift measurements.

The rates $\sigma_{i}$ in Table 1 are faster than the dynamic rates discussed below in Secs. III.2 and III.3, particularly for line 3. $\sigma_{i}$ is therefore attributed to static distributions of TF frequencies, most likely related to the crystalline-electric-field randomness observed in neutron scattering Li et al. (2017a), which lead to rms widths $\delta K=\sigma/\omega_{3}$ of Knight shift distributions. We do not fully understand the large difference between $\sigma_{1}$ and $\sigma_{2}$.

The temperature dependencies of $K_{\mathrm{raw}}$, $K_{\mathrm{corr}}$, and $\delta K$ for lines 1 and 2 are given in Fig. 3.

The error in $K_{\mathrm{corr}}$ is largely due to uncertainty in $A_{\mathrm{Ld}}$ [Eq. (3)]. Over the temperature range of the data both the Knight shifts and $\chi_{\mathrm{mol}}(T)=M(T)/H$ in 10 kOe Li et al. (2019) are practically constant.

The relation between the Knight shift and $\chi_{\mathrm{mol}}$ is often analyzed using the Clogston-Jaccarino relation Clogston and Jaccarino (1961)

where $A_{\mathrm{hf}}$ is the hyperfine coupling constant and temperature is an implicit variable. In the present case there is not enough temperature dependence to either quantity for this procedure to be reliable. Instead, we assume $K_{0}$ in Eq. (4) is negligible, as would be expected in a magnetic insulator, so that $A_{\mathrm{hf}}=K_{\mathrm{corr}}/\chi_{\mathrm{mol}}$. The hyperfine field $H_{\mathrm{hf}}$ (the local field at a $\mu^{+}$ site per $\mu_{B}$ of paramagnetic moment on the local-moment sites) is given by $H_{\mathrm{hf}}=N_{A}\mu_{B}A_{\mathrm{hf}}$, where $N_{A}$ is Avogadro’s number.

Values of $A_{\mathrm{hf}}$ and $H_{\mathrm{hf}}$ for the two lines are given in Table 1. Also included in Table 1 are calculated Yb³⁺ dipolar hyperfine fields at each of the two candidate $\mu^{+}$ sites discussed in Sec. II.1 (Appendix A); local-moment/$\mu^{+}$ coupling fields are expected to be predominantly dipolar in insulators Schenck (1985). The calculated and measured values are of the same order of magnitude, but with a discrepancy in sign for line 1. This is not understood, but it may mean that the location of site $\mu$1 is incorrect.

$\mu$SR experiments are extremely sensitive to static magnetism, ordered or disordered Schenck (1985), in the present case with an upper limit ${\sim}10^{-3}\mu_{B}$ per Yb ion on any sample-average static Yb³⁺ moment down to 22 mK (Sec. III.3.2). This is a strong constraint on proposals involving magnetic ground states Luo et al. (2017); Zhu et al. (2017); Parker and Balents (2018); Kimchi et al. (2018), although in the disordered local-singlet valence solid Kimchi et al. (2018) the “ultimate fate of the …defect spins is unknown.”

The lack of spin freezing is in apparent contradiction with the observation of a cusp in the low-field ac susceptibility of YbMgGaO₄ at $\sim$0.1 K Ma et al. (2018), and the conclusion that a spin-glass-like transition occurs there. We note that suppression of static order by the muon itself, a possible explanation of this discrepancy, has been shown to be unlikely Keren (2004). We also note that large static Yb³⁺ moments in a macroscopic spurious phase would result in a fraction of the ZF asymmetry spectrum with oscillations or rapid relaxation. This is not seen at the level of a few percent. The resolution of these issues is not clear, and more work will be required to elucidate the impact of the present results on the QSL debate.

Our $\mu$SR results suggest reinterpretation of some of the previous work. The observation of two equal-amplitude lines in the TF-$\mu$SR data is evidence for the two $\mu^{+}$ stopping sites in the two inequivalent O^2- lattice planes. This, rather than, for example, Mg²⁺/Ga³⁺ site mixing Li et al. (2017a), is the main source of inhomogeneity in the dynamic relaxation, and strongly suggests a two-exponential analysis of LF- and ZF-$\mu$SR data [Eqs. (5) and (6)] rather than a stretched-exponential analysis.

The fact that good fits can be obtained with either two-exponential or stretched-exponential functions is not surprising. In general, the relaxation function $P(t)$ associated with a spatially inhomogeneous distribution $P(\lambda)$ of rates $\lambda$ is given by the Laplace transform

Two-exponential and stretched-exponential functions are effectively two choices of $P(\lambda)$.

Solving Eq. (9) for $P(\lambda)$ is a mathematically ill-posed problem Press et al. (1992); $P(t)$ is not sensitive to the form of $P(\lambda)$. The difficulty is greater when the range of the independent variable is limited, as it is for $\mu$SR asymmetry spectra by the muon lifetime ⁴⁴4This is the problem in the present case. At late enough times, the sum of two exponentials with very different rates differs significantly from a stretched exponential. We also note that for $P(t)=\exp[-(\lambda^{\ast}t)^{\beta}]$, the relation of $\lambda^{\ast}$ to $P(\lambda)$ depends on $\beta$; $\lambda^{\ast}$ is not simply an average rate, and one should not compare values of $\lambda^{\ast}$ for different values of $\beta$ Johnston (2006).

The temperature dependence of $\beta$ in fits to ZF data is also consistent with the two-exponential scenario. (Quasi)static relaxation in ZF and weak LF by nuclear dipolar fields should not be ignored, although the lack of structure in the $\mu^{+}$ asymmetry spectra and the limitation imposed by the muon lifetime make the analysis difficult. We argue that the reported increase of $\beta$ with increasing temperature in ZF stretched-exponential fits is not a temperature-dependent effect of the spin dynamics Li et al. (2016a). Instead, the decrease of ZF dynamic relaxation with increasing temperature (Fig. 8) increases the importance of the temperature-independent “super-exponential” (Gaussian) nuclear dipolar field distribution. This renders the overall ZF relaxation function less sub-exponential at high temperatures, leading to $\beta\approx 1$ at $\sim$10 K Li et al. (2016a).