Persistent spin dynamics in the pressurized spin-liquid candidate YbMgGaO$_4$

TL;DR

This study investigates how pressure affects the crystal structure and spin dynamics of the spin-liquid candidate YbMgGaO$_4$, revealing persistent spin fluctuations up to high pressures and low temperatures.

Contribution

It provides new insights into the pressure stability of the crystal structure and the persistent spin dynamics in YbMgGaO$_4$, a key candidate for quantum spin liquids.

Findings

Crystal structure remains stable up to 10 GPa.

Persistent spin dynamics observed up to 2.6 GPa and 250 mK.

Spin-spin autocorrelation exhibits power-law behavior under pressure.

Abstract

Single-crystal x-ray diffraction, density-functional band-structure calculations, and muon spin relaxation ($\mu$SR) are used to probe pressure evolution of the triangular spin-liquid candidate YbMgGaO$_4$. The rhombohedral crystal structure is retained up to at least 10 GPa and shows a nearly uniform compression along both in-plane and out-of-plane directions, whereas local distortions caused by the random distribution of Mg$^{2+}$ and Ga$^{3+}$ remain mostly unchanged. The $\mu$SR data confirm persistent spin dynamics up to 2.6 GPa and down to 250 mK with no change in the muon relaxation rate. Longitudinal-field $\mu$SR reveals power-law behavior of the spin-spin autocorrelation function, both at ambient pressure and upon compression.

Figures

Tables

Table 1. Table 1: Details of the single-crystal XRD data collection and structure refinement at 10 GPa.

$T$ (K)	298
$a$ (Å)	3.34272(13)
$c$ (Å)	24.6033(7)
wavelength (Å)	0.41114
$h_{\min }$, $h_{\max }$	$-4\le h\le 3$
$k_{\min }$, $k_{\max }$	$-5\le k\le 5$
$l_{\min }$, $l_{\max }$	$-39\le l\le 38$
No. of reflections	126
No. of refinable parameters	9
$R_{\mathrm{int}}$	0.0249
$R_I$	0.0310

Table 2. Table 2: Comparison of the relaxed geometrical parameters for the ambient-pressure crystal structure in the primitive cell with the R ​ 3 ​ m 𝑅 3 𝑚 R3m symmetry. Different band-structure codes ( FPLO vs. VASP ), exchange-correlation potentials (LDA  Perdew and Wang ( 1992 ) vs. GGA  Perdew  et al. ( 1996 ) ), and approaches to the treatment of the Yb 4 ​ f 4 𝑓 4f shell (valence vs. core) have been used.

Code	FPLO	FPLO	FPLO	FPLO	VASP	Experiment
$V_{\mathrm{xc}}$	LDA	LDA	GGA	GGA	GGA
Yb $4f$	valence	core	valence	core	core
$d_1/d_2$ (Å)	2.229/2.248	2.214/2.233	2.231/2.257	2.215/2.243	2.258/2.270	2.234(4)
${\beta }_1/{\beta }_2$ (deg)	98.62/99.73	99.52/100.68	98.09/99.60	98.88/100.62	97.31/98.02	99.44(13)

Table 3. Table 3: Details of the μ 𝜇 \mu SR experiments.

$P$ (GPa)	Spectrometer	Pressure cell	$T$ (K)	Mode
0	Dolly	none	0.26–200	ZF/LF
0	GPD	MP35N+MP35N	0.25–10	ZF
1.9	GPD	MP35N+CuBe	0.25–10	ZF
2.6	GPD	MP35N+MP35N	0.25–10	ZF/LF

Equations

$$p(V)=\frac{3B_0}{2}\left[{\left(\frac{V_0}{V}\right)}^{\frac{7}{3}}-{\left(\frac{V_0}{V}\right)}^{\frac{5}{3}}\right],$$

$$A(t)=f_{\mathrm{PC}}{A}_{\mathrm{PC}}+(1-f_{\mathrm{PC}})e^{-\lambda t}$$

$$M=\frac{1}{N}\sum _i^N\frac{(A_i-A_{i+1}{)}^2}{({\delta }_i-{\delta }_{i+1}{)}^2}$$

$$A_{\mathrm{PC}}(t)=\left[\frac{1}{3}+\frac{2}{3}(1-{\sigma }_{\mathrm{PC}}^2{t}^2)e^{-{\sigma }_{\mathrm{PC}}^2{t}^2/2}\right]\exp (-{\lambda }_{\mathrm{PC}}t).$$

Full text

I Introduction

Spin-liquid states in frustrated magnets are nowadays actively studied as hosts for unconventional excitations representing magnetic monopoles Balents (2010); Gingras and McClarty (2014) and other exotic quasiparticles Savary and Balents (2017); Hermanns et al. (2018). One relatively unexplored aspect in this field is the evolution of spin-liquid materials under pressure and the change in spin dynamics caused by tuning magnetic interactions via lattice compression.

Here, we focus on the spin-liquid candidate YbMgGaO4 Li et al. (2015a, b, 2020) that recently evolved as a unique triangular antiferromagnet with the robust three-fold symmetry, persistent spin dynamics, and a broad continuum of (potentially fractionalized) magnetic excitations. The crystal structure of this compound features triangular layers of the pseudospin-$\frac{1}{2}$ Yb3+ ions that are well separated by slabs of non-magnetic Mg2+ and Ga3+ (Fig. 1a). Thermodynamic measurements Li et al. (2015a, b) and muon spin relaxation ($\mu$SR) Li et al. (2016) suggest the absence of magnetic order down to at least 50 mK. Weak spin freezing around 100 mK was indicated by ac-susceptibility data Ma et al. (2018), although it involves only a tiny amount of the magnetic entropy Li et al. (2015a) and leaves no signatures in either dc-susceptibility Li et al. (2019) or $\mu$SR Li et al. (2016).

Magnetic excitations of YbMgGaO4 form a broad continuum that can be interpreted in terms of gapless spinons Shen et al. (2016); Li et al. (2017a); Li and Chen (2017); Shen et al. (2018) or as arising from short-range valence bonds Li et al. (2017b, 2019), the latter suggestion being particularly interesting, as it makes a direct link to Anderson’s pioneering work at the outset of the spin-liquid research Anderson (1973). Continuum features were also reported for other Yb3+ triangular antiferromagnets Ding et al. (2019); Ma et al. ; Dai et al. , but they show a different distribution of the spectral weight and probably have a separate origin. Indeed, YbMgGaO4 is known to be strongly affected by structural randomness that arises from the random distribution of Mg2+ and Ga3+ between the Yb3+ layers and modulates magnetic interactions via random crystal electric fields acting on Yb3+ Li et al. (2017c); Paddison et al. (2017). No such randomness occurs in other Yb-based triangular antiferromagnets.

In the following, we explore the effect of hydrostatic pressure on the structure and magnetism of YbMgGaO4, and juxtapose the behavior of this material with the evolution of other spin-liquid candidates upon compression. We show that in YbMgGaO4 spins remain dynamic up to at least 2.6 GPa, and quantify associated structural changes for an eventual comparison across different classes of spin-liquid materials.

II Crystal structure

II.1 Average structure

Room-temperature x-ray diffraction (XRD) data were collected at the ID15B beamline of the European Synchrotron Radiation Facility (ESRF) between ambient pressure and 10 GPa. A diamond anvil cell loaded with He gas and a small single crystal of YbMgGaO4 from the batch reported in Ref. Li et al., 2015b were used for the experiment.

Three positional parameters ($z$-coordinates) and three independent atomic displacement parameters were refined for Mg/Ga, O1, and O2. The position of Yb was fixed at $(0,0,0)$, with its thermal ellipsoid refined anisotropically to account for local displacements caused by the random distribution of Mg and Ga. Altogether, 9 structural parameters were refined from about 120 symmetry-independent reflections collected at each pressure point sup . Details of the data collection and structure refinement are given in Appendix A.

The $R\bar{3}m$ symmetry of YbMgGaO4 keeps all Yb–O distances equal but allows a trigonal distortion of the YbO6 octahedra with the O–Yb–O angle $\beta$ deviating from $90^{\circ}$. Incidentally, this angle is equal by symmetry to the Yb–O–Yb bridging angle $\alpha$, which is responsible for superexchange interaction (Fig. 1b). Two structural parameters, one distance and one angle, are thus sufficient to characterize both the local geometry of Yb3+ and the geometry of the nearest-neighbor exchange pathway within the average structure.

No changes in the crystal symmetry were observed in our experiment. Pressure evolution of the $a$ and $c$ lattice parameters demonstrates similar compressibility along the different crystallographic directions (Fig. 2a). The fit with the second-order Birch-Murnaghan equation of state,

[TABLE]

using the EosFit routine Angel et al. (2014) returned the bulk modulus of $B_{0}=142(2)$ GPa and the unit-cell volume of $V_{0}=253.12(7)$ Å3. A similar fit for individual lattice parameters using $a^{3}$ or $c^{3}$ in place of $V$ yields $B_{0,a}=151(3)$ GPa and $B_{0,c}=130(2)$ GPa suggesting a nearly isotropic compression of the structure. YbMgGaO4 is more compressible than Yb2Ti2O7 ($B_{0}=219$ GPa Mishra et al. (2012)) and other rare-earth pyrochlores that typically feature the bulk moduli in excess of 200 GPa Turner et al. (2017).

The Yb–O distances shrink by about 0.6% at 2.6 GPa (the highest pressure of our $\mu$SR experiment, Sec. III) and by 1.7% at 10 GPa (Fig. 2b). The $\alpha=\beta$ angle shows a weak downward trend only, with a larger error bar caused by the lower accuracy in the determination of oxygen position due to the low scattering power of oxygen. In order to confirm this downward trend, we relaxed the experimental structures using density-functional (DFT) band-structure calculations. As DFT can not treat the Mg/Ga disorder explicitly, ordered structural models have to be used, as explained in Appendix B. This approximation leads to a constant offset between the DFT results and experiment. Nevertheless, not only the qualitative trends but also the slope are well reproduced (Fig. 2). We thus confirm that the Yb–O–Yb angle decreases under pressure. Compared to the ambient-pressure value, it changes by $0.07^{\circ}$ at 2.6 GPa and by $0.2^{\circ}$ at 10 GPa.

II.2 Local structure

Moderate changes of the average crystal structure are accompanied by a strong elongation of the Yb thermal ellipsoid. The in-plane displacements characterized by $U_{11}$ are not affected by pressure, whereas the out-of-plane displacement component increases by 70% (Fig. 1d). This out-of-plane displacement has been previously linked to the local distortions of the YbO6 octahedra caused by the random (and, generally, asymmetric) distribution of the differently charged Mg2+ and Ga3+ ions around Yb3+ Li et al. (2017c). At first glance, the increase in $U_{33}$ implies a strong enhancement of the structural randomness under pressure, but the actual situation is more complex.

We analyze pressure evolution of the local structure by DFT relaxations for several ordered models of YbMgGaO4. First, we consider the models of Ref. Li et al., 2017c that were previously used to interpret crystal-field excitations of Yb3+, and evaluate $\Delta z_{\rm Yb}$ as the displacements of the Yb atoms relative to each other. These displacements indicate the spread of the Yb electron density, which is gauged by the $U_{33}$ parameter of the average structure.

Three scenarios visualized in Fig. 3a,b are considered, with the Yb3+ ions sandwiched between: i) two Ga3+ ions on one side and two Mg2+ ions on the other side (octahedra A); ii) one Ga3+ and one Mg2+ ion on each side (octahedra B); iii) two Ga3+ ions on one side and a combination of Ga3+ and Mg2+ on the other side (octahedra C). Pressure has a strong effect on the Yb position in the octahedra A and nearly no effect on the Yb atoms in B and C (Fig. 1e). Surprisingly, this effect is a shift of the whole octahedron A along the $c$ direction without any change in the octahedron itself. Using the geometrical parameters shown in Fig. 1c, we determine that for the octahedra A the deformation expressed by $(d_{2}-d_{1})$ changes from 1.65 % at 0 GPa to 1.66 % at 10 GPa, whereas $(\beta_{2}-\beta_{1})$ evolves from $4.29^{\circ}$ to $4.32^{\circ}$ (Fig. 3d,e). The changes in the octahedra B and C are equally small. Since $\alpha=\beta$ in our material, the distribution of the Yb–O–Yb angles also does not change with pressure. Therefore, neither local randomness at each Yb site nor randomness of the exchange couplings should be affected by pressure.

We also consider a different scenario where two adjacent Yb atoms within the same layer have a different local environment. In layer D (Fig. 3c), Yb1 is surrounded by 3 Ga atoms on one side and 3 Mg atoms on the other side. In contrast, Yb2 and Yb3 are surrounded by 2 Ga and 4 Mg atoms each. These dissimilar local environments lead to different $\alpha_{1}$ and $\alpha_{2}$ for the Yb1–Yb2 and Yb2–Yb3 superexchange pathways and reflect the effect of structural randomness on nearest-neighbor magnetic interactions in YbMgGaO4. However, the difference between these two pathways appears to be only weakly pressure-dependent and decreases upon compression (Fig. 3e).

Our results suggest that the visible enhancement of $U_{33}$ under pressure is not related to the increased distortions around the Yb3+ ions. It rather indicates a change in the overall position of the octahedra A that may be explained by the accumulation of different charges above (Ga3+) and below (Mg2+) the Yb layer. But the key result at this juncture is that such a change reflects a re-arrangement within the Mg/Ga slabs and does not affect the structure of the magnetic layer itself. The YbO6 octahedra undergo a uniform compression and simply keep the deformation that they had at ambient pressure. The randomness effect on the exchange couplings should be largely unchanged or even decrease under pressure.

III Spin dynamics

Pressure evolution of spin dynamics is probed by muon spin relaxation ($\mu$SR). The experiments were performed at the GPD and Dolly spectrometers at the Paul Scherrer Institute (PSI) at ambient pressure and at the GPD spectrometer under hydrostatic pressure Khasanov et al. (2016) down to 250 mK on a polycrystalline sample of YbMgGaO4. The data were collected in zero field (ZF) and in the longitudinal-field (LF) mode, where the magnetic field was applied parallel to the spin of the implanted muons.

III.1 ZF $\mu$SR

ZF $\mu$SR time spectra measured at ambient pressure and at 2.6 GPa, the highest pressure of our experiment, are compared in Figure 4. The similar behavior of the $\mu$SR time spectra indicates no change of the magnetic ground state, and the absence of oscillations excludes pressure-induced magnetic ordering in YbMgGaO4 up to at least 2.6 GPa. To estimate the temperature dependence of the relaxation rate $\lambda$, we fitted the ZF-$\mu$SR time spectra by

[TABLE]

where $A_{\rm PC}$ is the pressure-cell contribution described in Appendix C. The fraction of signal coming from the pressure cell was $f_{\rm PC}=0.7$ at 0 GPa and 2.6 GPa, and $f_{\rm PC}=0.5$ at 1.9 GPa, where a different cell was used in an attempt to reduce the background.

The zero-field muon relaxation rate $\lambda$ is temperature-independent between 4 and 40 K (Fig. 4). Its increase below 4 K indicates the onset of spin-spin correlations that fully develop around 0.8 K, where $\lambda$ flattens out and remains temperature-independent upon further cooling. This temperature evolution is essentially similar to the ambient-pressure $\mu$SR data reported in Ref. Li et al., 2016 and remains unchanged at 2.6 GPa (Fig. 4). Even absolute values of $\lambda$ are the same as at ambient pressure within the error bar.

Above 40 K, $\lambda$ shows a steady decreasing trend described by an activated behavior $\lambda=A+B\,e^{-\Delta/T}$ with $\Delta=320\pm 20$ K, which is reminiscent of the lowest crystal-field excitation energy of about 450 K Li et al. (2017c). This observation suggests that at high temperatures the relaxation is governed by an Orbach process Orbach (1961) involving the excited crystal-field doublets of Yb3+.

III.2 LF $\mu$SR

LF measurements complement the data obtained in zero field. Even a longitudinal field, which is more than 10 times higher than the local static field estimated from the low-temperature value of $\lambda$, does not decouple the muon relaxation, suggesting that the spins are dynamic in nature as has been seen at ambient pressure Li et al. (2016).

The LF data measured in different fields follow universal scaling initially proposed by Keren et al. Keren et al. (1996) for systems with glassy dynamics. This scaling manifests itself in the power-law behavior of the muon relaxation rate $\lambda(H)$ extracted by fitting individual LF curves with a stretched exponential $e^{-[\lambda(H)t]^{\beta}}$ supplied with the time-independent and field-independent background, which originates from the sample holder and pressure cell. The ensuing values of $\lambda$ follow $\lambda(H)\sim H^{-\gamma}$ with $\gamma\simeq 0.3$ at 0 GPa and 0.8 at 2.6 GPa (Fig. 5).

This behavior is confirmed by the scaling of muon asymmetry plotted against $t/H^{\gamma}$. To this end, the data points from all fields up to $t=6$ $\mu$s 111The data above $t=6$ $\mu$s were excluded due to their large error bars. are arranged with increasing $t/H^{\gamma}$ for every value of $\gamma$. An empirical data mismatch function is calculated by taking the difference between the neighboring points and weighing them by the corresponding error bars. This mismatch function is defined as

[TABLE]

where $N$ is the number of data points, and $A_{i}$ and $\delta_{i}$ correspond to the asymmetry and error bar of the $i$-th data point, respectively. The lowest value of the mismatch function is obtained at $\gamma=0.35$ (0 GPa) and $\gamma=0.85$ (2.6 GPa) that produce the universal scaling over at least three orders of magnitude in $t/H^{\gamma}$ and show excellent agreement with the $\gamma$ values from the analysis of $\lambda(H)$.

The scaling may not hold in the whole time and field range – e.g., at very short times comparable to the width of the muon pulse Keren et al. (1996) or in very high fields that affect spin dynamics – but the scaling over three orders of magnitude (Fig. 5) is congruent with all earlier observations Keren (2004); MacLaughlin et al. (2004) and serves as a robust evidence of the power-law behavior $q(t)\sim t^{\gamma-1}$ of the dynamic spin-spin autocorrelation function $q(t)=\langle S_{i}(t)\cdot S_{i}(0)\rangle$. It indicates collective dynamics, as for example in spin-glass systems where $\gamma$ increases toward 1.0 upon approaching the freezing temperature from above Keren et al. (1996).

IV Discussion and summary

YbMgGaO4 is remarkably insensitive to pressure and thus different from other spin-liquid candidates. For example, herbertsmithite Mendels and Bert (2010); Norman (2016) becomes magnetically ordered at 2.5 GPa Kozlenko et al. (2012), whereas stoichiometric Yb2Ti2O7 Rau and Gingras (2019) shows pressure-induced magnetic order already at 0.1 GPa Kermarrec et al. (2017). These materials have different compressibilities, with the bulk modulus of $B_{0}=70$ GPa in herbertsmithite Kozlenko et al. (2012) and $B_{0}=219$ GPa in Yb2Ti2O7 Mishra et al. (2012). YbMgGaO4 shows an intermediate value of 142 GPa and is certainly comparable to other spin-liquid candidates as far as elastic properties are concerned. Its structural parameters are affected by pressure. The persistence of spin dynamics will then indicate that structural changes influence exchange couplings only weakly, whereas the material lies sufficiently far away from the phase boundary between the spin-liquid and magnetically ordered states. Alternatively, spin dynamics may not be caused by frustration that controls this phase boundary, and originate from structural randomness, which is not affected by pressure.

Triangular antiferromagnets develop $120^{\circ}$ order or stripe order depending on the extent of exchange anisotropy Luo et al. (2017) and on the ratio $J_{2}/J_{1}$ of the second-neighbor to nearest-neighbor exchange coupling Zhu et al. (2018); Maksimov et al. (2019). The formation of a spin liquid is mostly controlled by $J_{2}/J_{1}$, with the complete suppression of magnetic order between $J_{2}/J_{1}\simeq 0.07$ and 0.15 Zhu et al. (2018); Maksimov et al. (2019); Li et al. (2020). Experimental estimates of this parameter vary in a broad range between 0.0 and 0.26 Li et al. (2020) and would place YbMgGaO4 in any of the regions on the phase diagram: $120^{\circ}$ order Li et al. (2015b), stripe order Paddison et al. (2017); Zhu et al. (2017), or the boundary between the latter and the spin-liquid phase Zhang et al. (2018).

External pressure of 2.6 GPa leads to a shortening of the Yb–O distance $d$ by 0.013 Å and a reduction in the Yb–O–Yb angle $\alpha$ by $0.07^{\circ}$ (Fig. 2). Using superexchange theory of Ref. Rau and Gingras, 2018, we estimate that the change in the angle increases $J_{1}$ by 6 %, whereas the shortening of the Yb–O distance increases $J_{1}$ by another 5 % assuming exponential dependence of Slater-Koster parameters on the interatomic distance. Overall, we expect that $J_{1}$ increases by 11 %, which is similar in magnitude to pressure-induced changes in Cu-based quantum magnets Zayed et al. (2017); Zvyagin et al. (2019). No significant changes are expected for $J_{2}$, because long-range couplings are less sensitive to the structural geometry. In YbMgGaO4, the Yb–O–O angles decrease by as little as $0.11^{\circ}$ upon compression to 2.6 GPa. In contrast, in a Cu-based quantum magnet even the $2^{\circ}$ change of the respective angle reduces the coupling by 20 % only Prishchenko et al. (2017).

We conclude that $J_{1}$ increases, $J_{2}$ is roughly unchanged, and the $J_{2}/J_{1}$ ratio of YbMgGaO4 should be reduced under pressure. Assuming the parametrization of Ref. Zhang et al., 2018 with $J_{2}/J_{1}=0.18$, it would imply that pressure pushes the material into the spin-liquid region, in agreement with the persistence of spin dynamics observed experimentally. Alternatively, and perhaps even more likely, the lack of significant changes under pressure may indicate that spin dynamics is triggered by structural randomness and thus unaffected by pressure.

The only change we observe is the increase in the $\gamma$ parameter that describes scaling behavior in longitudinal fields. For a system with glassy dynamics the increase in $\gamma$ would reflect an evolution toward a frozen state Keren et al. (1996) that was indeed proposed for YbMgGaO4 below 100 mK (at ambient pressure) based on the ac-susceptibility data Ma et al. (2018). However, neither dc-susceptibility Li et al. (2019) nor $\mu$SR Li et al. (2016) support bulk spin freezing at ambient pressure. The scaling behavior may be in fact unrelated to glassy dynamics, and indeed it was also reported for a variety of systems with critical fluctuations influenced by structural disorder MacLaughlin et al. (2004); Keren (2004). On the other hand, the very presence of this scaling serves as an additional argument for collective spin dynamics caused by the structural randomness, and the importance of this randomness for the spin-liquid behavior of YbMgGaO4.

Altogether, hydrostatic pressure leads to a uniform compression of the YbMgGaO4 structure with the reduction in the Yb–O distances and Yb–O–Yb angles, whereas local distortions of the YbO6 octahedra and consequent randomness effects are nearly unchanged. Spin dynamics is not affected by pressure and appears to be collective yet influenced by the structural randomness. This puts YbMgGaO4 into the group of materials, where randomness effects can be integral to the spin-liquid formation, but collective spin dynamics is nevertheless observed.

Acknowledgements.

We acknowledge ESRF for provision of beamtime at ID15B. IEC thanks J. Jacobs for the He gas load. Part of this work is based on experiments performed at the Swiss Muon Source S$\mu$S, Paul Scherrer Institute, Villigen, Switzerland. Swiss National Science Foundation has supported the work of G.S. (Grants No. 200021-175935 and Mobility grant P2EZP2-178604). The work in Augsburg was supported by the German Research Foundation (DFG) via the Project No. 107745057 (TRR80) and by the Federal Ministry for Education and Research through the Sofja Kovalevkaya Award of Alexander von Humboldt Foundation (AAT).

Appendix A X-ray diffraction

Membrane-driven LeToullec-type diamond anvil cells (DACs) were used, equipped with Boehler-Almax anvils. Stainless steel was used as the gasket material, and helium was loaded as the pressure-transmitting medium. Diffraction patterns were collected with a Mar555 flat-panel detector using steps of $0.5^{\circ}$ oscillations over a total $\omega$ scan range of $76^{\circ}$ about the vertical axis. The pressures were measured using the ruby fluorescence method. Lattice parameter determination and integration of the reflection intensities were performed using the CrysAlisPro software cry . Structures were refined using ShelxL Sheldrick (2008) within the ShelXle Hübschle et al. (2011) graphical interface.

Details of the data collection and structure refinement at 10 GPa are listed in Table 1. The refinements at other pressures were similar and can be found in the cif-file provided as Supplemental Material. We note that the ambient-pressure experiment was performed before loading the cell with the He gas. This may be the reason for the abrupt change in the Yb displacement parameters $U_{11}$ and $U_{33}$ between 0 and 0.78 GPa (Fig. 1d).

Appendix B Computational results

B.1 General methodology

Non-magnetic DFT calculations were performed in the FPLO Koepernik and Eschrig (1999) and VASP Kresse and Furthmüller (1996a); *vasp2 codes. Atomic positions were optimized until the energy minimum was reached, and residual forces dropped below 0.005 eV/Å. The typical $k$-mesh included 64 points within the first Brillouin zone. Experimental lattice parameters were used at each pressure, and only atomic positions were relaxed.

The results of structure optimization may be affected by the computational methodology. In Table 2, we analyze the role of three factors: i) band-structure code; ii) exchange-correlation potential; and iii) treatment of the Yb $4f$ shell. The calculations are performed for the simplest $R3m$ ordered model of YbMgGaO4 described in App. B.2 below.

The choice of the band-structure code affects the results of the relaxation. VASP calculations predict the too long Yb–O distances and, consequently, the too low angles. On the other hand, the relaxations within FPLO produce, irrespective of the methodology, the average distance $(d_{1}+d_{2})/2$ and the average angle $(\beta_{1}+\beta_{2})/2$ within, respectively, 0.01 Å and $0.5^{\circ}$ from the experimental values. The larger deviations of the VASP results are probably related to the lower accuracy of the default pseudopotential for the Yb atoms, whereas FPLO does not rely on pseudopotentials and introduces no approximations to the crystal potential.

All choices mentioned here, including the VASP calculations, produce qualitatively similar results as a function of pressure. The trends obtained from DFT are thus robust. The results shown in the manuscript are obtained with FPLO, GGA functional, and Yb $4f$ states placed into the core, as this choice improved the convergence for larger unit cells.

B.2 Ordered structural models

As DFT can not treat mixed sites, ordered structural models are used for calculations. We adopted three types of such ordered models:

1. $R3m$ model obtained from the parent $R\bar{3}m$ structure by splitting the mixed Mg/Ga site into two, one fully occupied by Ga and the other one fully occupied by Mg. In this case, each Yb layer is sandwiched between the Mg and Ga layers leading to the off-center Yb displacement, two non-equivalent Yb–O distances $d_{1}\neq d_{2}$, and two O–Yb–O angles $\beta_{1}\neq\beta_{2}$. The $R3m$ structure was used in Table 2 for testing the effect of band-structure code, basis set, and exchange-correlation potential.

2. $P3m$ models that preserve three-fold symmetry and feature layer-by-layer ordering of Mg and Ga. We chose two of such ordered structures (Fig. 3) that correspond to the third and fourth structures from Fig. S5 of Ref. Li et al., 2017c.

In the structure from Fig. 3a, all Yb atoms develop a nearly symmetric local environment. This structure was used to evaluate pressure dependence of the Yb–O distances and Yb–O–Yb angles shown in Fig. 2. In this case, the octahedra are still deformed, but $d_{2}-d_{1}<0.005$ Å and $\beta_{2}-\beta_{1}<0.2^{\circ}$ indicate a negligibly small deformation. The “symmetric” structure (left panel of Fig. 3) is also used to obtain the displacements ($\Delta z_{\rm Yb}$) of Yb within the octahedra B (GaMgYbMgGa).

The structure from Fig. 3b covers the opposite scenario of the Yb atoms with a highly asymmetric local environment. It contains the octahedra A (GaGaYbMgMg) and C (GaGaYbMgGa) showing the largest values of $d_{2}-d_{1}$ and $\beta_{2}-\beta_{1}$. As explained in the main text, thermal displacement parameter of the average structure, $U_{33}$, does not reflect deformations of individual octahedra and shows instead relative displacements of the Yb atoms with respect to each other. In the ideal structure, the Yb atoms should be at $z=0,\frac{1}{3},\frac{2}{3}$. To assess the effect of pressure on $U_{33}$, we compare the positions of the Yb atoms within the octahedra A and C to the remaining Yb atom at $z\simeq\frac{2}{3}$ that shows a rather symmetric local environment and can be used as reference. The values of $\Delta z_{\rm Yb}$ are obtained as $\left|z_{A,C}-z_{\rm ref}-\frac{1}{3}\right|$, where $z_{\rm ref}$ is the $z$-coordinate of the reference Yb atom.

3. $P3m$ model with the four-fold $2\mathbf{a}\times 2\mathbf{b}\times\mathbf{c}$ supercell, where $\mathbf{a}$, $\mathbf{b}$, and $\mathbf{c}$ are lattice vectors of the parent $R\bar{3}m$ structure. This supercell allows to construct configurations, where adjacent Yb atoms have different local environment. We focus on the fragment shown in Fig. 3c, where Yb1–Yb2 and Yb2–Yb3 represent two non-equivalent superexchange pathways caused by the uneven distribution of Mg and Ga around the Yb layer.

Appendix C $\mu$SR experiment

$\mu$SR measurements were performed in double-walled MP35 pressure cells with Daphne oil 7373 as pressure-transmitting medium. The pressure value was determined by measuring the superconducting transition of a small piece of indium positioned next to the sample inside the pressure cell. Further experimental details are summarized in Table 3

Two different pressure cells were used. The 0 GPa and 2.6 GPa data were collected in the double-walled MP35N+MP35N type cell Khasanov et al. (2016), whereas the 1.9 GPa data were collected separately in the low-background double-walled MP35N+CuBe type cell Shermadini et al. (2017). For both cells, we employed Eq. (2) and expressed the pressure-cell contribution as the Kubo-Toyabe depolarization function multiplied by an exponential damping,

[TABLE]

In the case of the MP35N+MP35N type cell, $\sigma_{\rm PC}=0.29\,\mu\text{s}^{-1}$ was temperature-independent, whereas $\lambda_{\rm PC}$ remained constant down to 1 K and increased at lower temperatures, similar to Ref. Khasanov et al., 2016. The fraction of the signal coming from the pressure cell was $f_{\rm PC}=0.7$. In the case of the MP35N+CuBe type cell, both $\sigma_{\rm PC}=0.32\,\mu\text{s}^{-1}$ and $\lambda_{\rm PC}=0.025\,\mu\text{s}^{-1}$ are temperature-independent in agreement with those of Ref. Shermadini et al., 2017, and $f_{\rm PC}=0.5$. The usage of different cells may be the reason for a slight offset between the 1.9 GPa data and 2.6 GPa data in Fig. 4.

In longitudinal fields above 100 G, the contribution of the pressure cell should be mostly decoupled, and the data include only a constant, field-independent background. Fig. 6 shows the LF data fitted by the stretched exponentials to obtain $\lambda(H)$ presented in Fig. 5.