Crystalline Electric Field Randomness in the Triangular Lattice Spin-Liquid YbMgGaO$_4$

TL;DR

This study investigates the crystalline electric field levels in the spin-liquid candidate YbMgGaO$_4$ using inelastic neutron scattering, revealing broadening effects due to site disorder that influence magnetic excitations.

Contribution

It demonstrates that CEF level broadening in YbMgGaO$_4$ arises from Mg/Ga site mixing, affecting magnetic properties and excitation spectra.

Findings

CEF excitations at 39, 61, and 97 meV with broadening

Site disorder causes distribution of Yb--O distances and CEF parameters

Broadening explains low-energy excitation features in the polarized phase

Abstract

We apply moderate-high-energy inelastic neutron scattering (INS) measurements to investigate Yb$^{3+}$ crystalline electric field (CEF) levels in the triangular spin-liquid candidate YbMgGaO$_4$. Three CEF excitations from the ground-state Kramers doublet are centered at the energies $\hbar \omega$ = 39, 61, and 97\,meV in agreement with the effective \mbox{spin-1/2} $g$-factors and experimental heat capacity, but reveal sizable broadening. We argue that this broadening originates from the site mixing between Mg$^{2+}$ and Ga$^{3+}$ giving rise to a distribution of Yb--O distances and orientations and, thus, of CEF parameters that account for the peculiar energy profile of the CEF excitations. The CEF randomness gives rise to a distribution of the effective spin-1/2 $g$-factors and explains the unprecedented broadening of low-energy magnetic excitations in the fully polarized ferromagnetic phase of YbMgGaO$_4$, although a distribution of magnetic couplings due to the Mg/Ga disorder may be important as well.

Tables

Table 1. Table 1: Combined CEF fits for YbMgGaO 4 , as described in the main text. Fit #1 is performed against the whole dataset. For the fit #2, the region between 73 and 90 meV is excluded. For the fit #3, the energy of the 87 meV shoulder is used as the third CEF excitation. For the fit #4, the energy of the 139 meV mode is used as the third CEF excitation, see text for further details.

fits	#1		#2		#3		#4
	$X_i^{\mathrm{obs}}$	$X_i^{\mathrm{cal}}$	$X_i^{\mathrm{obs}}$	$X_i^{\mathrm{cal}}$	$X_i^{\mathrm{obs}}$	$X_i^{\mathrm{cal}}$	$X_i^{\mathrm{obs}}$	$X_i^{\mathrm{cal}}$
$\hslash {\omega }_1$(meV)	39.36(8)	39.36	39.31(8)	39.31	39.37(7)	39.24	39.36(8)	39.50
$\hslash {\omega }_2$(meV)	61.3(1)	61.3	61.3(1)	61.1	61.3(1)	60.8	61.3(1)	62.18
$\hslash {\omega }_3$(meV)	96.6(2)	96.4	97.0(2)	97.1	87(1)	89	138.7(3)	136.1
$I_2/I_1$	0.72(2)	0.70	0.72(2)	0.78	0.72(2)	0.84	0.72(2)	0.41
$I_3/I_1$	0.71(2)	0.91	0.71(2)	0.82	0.71(2)	0.69	0.71(2)	0.79
$g_{⟂}$ Li et al. (2015b)	3.06(2)	3.21	3.06(2)	3.19	3.06(2)	3.18	3.06(2)	3.02
$g_{\parallel }$ Li et al. (2015b)	3.72(3)	3.68	3.72(3)	3.67	3.72(3)	3.44	3.72(3)	4.06
$B_2^0$(meV)	$-0.130$		$-0.141$		$-0.108$		$0.571$
$B_4^0$(meV)	$0.00179$		$0.0020$		$-0.0013$		$0.0091$
$B_4^3$(meV)	$-0.93$		$-0.96$		$-0.92$		$-1.34$
$B_6^0$(meV)	$-0.00095$		$-0.00096$		$-0.00082$		$-0.0012$
$B_6^3$(meV)	$-0.0247$		$-0.0239$		$-0.0186$		$-0.0455$
$B_6^6$(meV)	$0.0190$		$0.0169$		$0.0151$		$0.00179$

Table 2. Table 2: Fitted CEF wavefunctions.

Fit #1 :
$|{\omega }_{0,\pm }⟩$ = $\pm 0.71|\pm 7/2⟩\mp 0.36|\mp 5/2⟩+0.60|\pm 1/2⟩$
$|{\omega }_{1,\pm }⟩$ = $0.97|\pm 3/2⟩\pm 0.24|\mp 3/2⟩$
$|{\omega }_{2,\pm }⟩$ = $\mp 0.21|\pm 7/2⟩+0.02|\mp 7/2⟩+0.07|\pm 5/2⟩\mp 0.92|\mp 5/2⟩-0.31|\pm 1/2⟩\mp 0.02|\mp 1/2⟩$
$|{\omega }_{3,\pm }⟩$ = $0.67|\pm 7/2⟩+0.09|\mp 5/2⟩\mp 0.74|\pm 1/2⟩$
Fit #2 :
$|{\omega }_{0,\pm }⟩$ = $\pm 0.71|\pm 7/2⟩\mp 0.35|\mp 5/2⟩+0.61|\pm 1/2⟩$
$|{\omega }_{1,\pm }⟩$ = $|\pm 3/2⟩$
$|{\omega }_{2,\pm }⟩$ = $\mp 0.25|\pm 7/2⟩-0.02|\mp 7/2⟩-0.09|\pm 5/2⟩\mp 0.93|\mp 5/2⟩-0.25|\pm 1/2⟩\pm 0.02|\mp 1/2⟩$
$|{\omega }_{3,\pm }⟩$ = $0.66|\pm 7/2⟩+0.03|\mp 5/2⟩\mp 0.75|\pm 1/2⟩$
Fit #3 :
$|{\omega }_{0,\pm }⟩$ = $\pm 0.70|\pm 7/2⟩\mp 0.35|\mp 5/2⟩+0.62|\pm 1/2⟩$
$|{\omega }_{1,\pm }⟩$ = $0.71|\pm 3/2⟩\mp 0.71|\mp 3/2⟩$
$|{\omega }_{2,\pm }⟩$ = $\mp 0.29|\pm 7/2⟩+0.03|\mp 7/2⟩+0.11|\pm 5/2⟩\mp 0.93|\mp 5/2⟩-0.21|\pm 1/2⟩\mp 0.02|\mp 1/2⟩$
$|{\omega }_{3,\pm }⟩$ = $0.66|\pm 7/2⟩-0.03|\mp 5/2⟩\mp 0.75|\pm 1/2⟩$
Fit #4 :
$|{\omega }_{0,\pm }⟩$ = $\pm 0.66|\pm 7/2⟩+0.08|\mp 7/2⟩-0.02|\pm 5/2⟩\mp 0.14|\mp 5/2⟩+0.73|\pm 1/2⟩\mp 0.09|\mp 1/2⟩$
$|{\omega }_{1,\pm }⟩$ = $|\pm 3/2⟩$
$|{\omega }_{2,\pm }⟩$ = $\mp 0.15|\pm 7/2⟩-0.03|\pm 5/2⟩\mp 0.99|\mp 5/2⟩-0.05|\pm 1/2⟩$
$|{\omega }_{3,\pm }⟩$ = $0.73|\pm 7/2⟩-0.08|\mp 5/2⟩\mp 0.68|\pm 1/2⟩$

Table 3. Table 3: Crystal structure of YbMgGaO 4 at 100 K determined from single-crystal XRD data ( a = 3.4061 Å and c = 25.1297 Å). In contrast to Ref.  Li et al. ( 2015b ) , we considered the fully position ordered model with the Yb atoms at (0,0,0), resulting in the elongated thermal ellipsoid for Yb. Alternatively, this elongation can be described by the split position of Yb, thus producing the disordered structural model reported in Ref.  Li et al. ( 2015b ) . However, our current analysis invalidates such a model, because the Yb atoms are not simply split into two positions, as the model would suggest, and rather show a continuous spread along the c-axis, following the disordered arrangement of Mg and Ga in the structure.

Space group: $R\overline{3}m$
Yb	z	0
	U11 = U22 = 2U12	0.0127(3)
x = y = 0	U33	0.0240(4)
	U33/U11	1.89(7)
U13 = U23 = 0	Occupancy	1
Mg/Ga	z	0.21442(7)
	U11 = U22 = 2U12	0.0184(6)
x = y = 0	U33	0.0134(8)
	U33/U11	0.73(7)
U13 = U23 = 0	Occupancy	0.5
O1	z	0.2912(3)
	U11 = U22 = 2U12	0.020(2)
x = y = 0	U33	0.015(4)
	U33/U11	0.8(3)
U13 = U23 = 0	Occupancy	1
O2	z	0.1284(4)
	U11 = U22 = 2U12	0.027(3)
x = y = 0	U33	0.018(4)
	U33/U11	0.7(2)
U13 = U23 = 0	Occupancy	1
Number of variables		12
Residuals	R(F) (I $>$ 3${\sigma }_I$)	0.049
	Rw(F) (I $>$ 3${\sigma }_I$)	0.063

Table 4. Table 4: Optimized geometries of the YbO 6 octahedra for different stacking sequences of Mg/Ga. Here, R O ​ 1 subscript 𝑅 𝑂 1 R_{O1} is the distance between Yb and its neighboring O1, and θ O1 subscript 𝜃 O1 \theta_{\rm O1} is the angle between the Yb–O1 bond and the c -axis. For undistorted YbO 6 octahedra, there is only one set of R O1 subscript 𝑅 O1 R_{\rm O1} and θ O1 subscript 𝜃 O1 \theta_{\rm O1} , whereas two sets of R O1 subscript 𝑅 O1 R_{\rm O1} and θ O1 subscript 𝜃 O1 \theta_{\rm O1} are obtained for the distorted YbO 6 octahedra.

YbO6 octahedrons	Local stackings	Probability	$R_{O1}$	${\theta }_{O1}$
1	MgGaYbGaMg/	1/7	2.2615	60.3953
	GaMgYbMgGa		undistorted	undistorted
2	MgGaYbMgGa/	1/7	2.2685	60.1025
	GaMgYbGaMg		2.2570	60.6052
3	GaMgYbMgMg/	1/7	2.260	60.456
	MgMgYbMgGa		2.251	60.865
4	GaGaYbGaMg/	1/7	2.321	57.927
	MgGaYbGaGa		2.243	61.244
5	MgGaYbMgMg/	1/7	2.258	60.552
	MgMgYbGaMg		2.240	61.393
6	GaGaYbMgGa/	1/7	2.316	58.105
	GaMgYbGaGa		2.247	61.054
7	GaGaYbMgMg/	1/7	2.321	57.904
	MgMgYbGaGa		2.228	61.972

Table 5. Table 5: Relative CEF parameters, B n m superscript subscript 𝐵 𝑛 𝑚 B_{n}^{m} (YbO 6 )/ ⟨ B n m ⟩ delimited-⟨⟩ superscript subscript 𝐵 𝑛 𝑚 \langle B_{n}^{m}\rangle , obtained from the point-charge model for different Yb-O geometries from Table  4 .

YbO6 octahedra	Probability	$n=2$; $m=0$	$n=4$; $m=0$	$n=4$; $m=3$	$n=6$; $m=0$	$n=6$; $m=3$	$n=6$; $m=6$
1	1/7	1.0000	1.0000	1.0000	1.0000	1.0000	1.0000
2	1/7	0.9916	1.0005	0.9973	0.9951	0.9766	0.9942
3	1/7	1.0529	0.9897	1.0130	1.0240	1.1605	1.0352
4	1/7	0.8510	1.0071	0.9570	0.9153	0.6276	0.9138
5	1/7	1.1154	0.9749	1.0271	1.0496	1.3554	1.0768
6	1/7	0.8481	1.0116	0.9578	0.9171	0.6064	0.9106
7	1/7	0.9233	0.9880	0.9731	0.9429	0.8615	0.9646
Expected deviations		0.0792	0.0095	0.0221	0.0433	0.2063	0.0470

Table 6. Table 6: Relative CEF parameters, B n m superscript subscript 𝐵 𝑛 𝑚 B_{n}^{m} (Mg/Ga)/ ⟨ B n m ⟩ delimited-⟨⟩ superscript subscript 𝐵 𝑛 𝑚 \langle B_{n}^{m}\rangle , obtained from the point-charge model for different distributions of Mg 2+ and Ga 3+ around Yb 3+ (total 6 nearest-neighbor Mg/Ga sites around each Yb 3+ ). All ions are kept at their ideal positions obtained from XRD (see Table  3 ).

Mg/Ga configurations	Probability	$n=2$; $m=0$	$n=4$; $m=0$	$n=4$; $m=3$	$n=6$; $m=0$	$n=6$; $m=3$	$n=6$; $m=6$
6Mg	1/64	0.8990	1.0113	1.0107	0.9891	0.9702	1.0005
5Mg+Ga	3/32	0.9327	1.0075	1.0071	0.9927	0.9801	1.0004
4Mg+2Ga	15/64	0.9663	1.0038	1.0036	0.9964	0.9901	1.0002
3Mg+3Ga	5/16	1.0000	1.0000	1.0000	1.0000	1.0000	1.0000
2Mg+4Ga	15/64	1.0337	0.9962	0.9964	1.0036	1.0099	0.9998
Mg+5Ga	3/32	1.0673	0.9925	0.9929	1.0073	1.0199	0.9996
6Ga	1/64	1.1010	0.9887	0.9893	1.0109	1.0298	0.9995
Expected deviations		0.0315	0.0035	0.0033	0.0034	0.0093	0.00017

Equations

$$R_p=\sqrt{\frac{1}{7}\sum _{i=1}^7{\left(\frac{X_i^{\mathrm{obs}}-X_i^{\mathrm{cal}}}{{\sigma }_i^{\mathrm{obs}}}\right)}^2},$$

$$C_m^{\mathrm{CEF}}=\frac{1}{k_B{T}^2}\frac{{\partial }^2\ln [{\sum }_{k=0}^{3}2\exp (-\frac{\hslash {\omega }_k}{k_{B}T})]}{\partial (\frac{1}{k_{B}T}{)}^2}.$$

$$\frac{d^2\sigma }{d\Omega d\omega }\propto |F(\mid \text{Q}\mid ){|}^2\times [\frac{Q_y^2{S}^{xx}(\text{Q},E)+Q_x^2{S}^{yy}(\text{Q},E)}{\mid \text{Q}{\mid }^2}+S^{zz}(\text{Q},E)],$$

$$H_0=B_2^0{O}_2^0+B_4^0{O}_4^0+B_4^3{O}_4^3+B_6^0{O}_6^0+B_6^3{O}_6^3+B_6^6{O}_6^6,$$

$$|{\omega }_{k,\pm }⟩=\sum _{m_J=-7/2}^{7/2}{C}_{m_J}^{k,\pm }|J=7/2,m_J⟩,$$

$$\frac{d^2\sigma }{d\Omega d\omega }\propto |F(Q){|}^2\exp \left(\frac{-\hslash {\omega }_{k1}}{k_{B}T}\right)\times \sum _{\alpha =x,y,z}|⟨{\omega }_{k2}|J_{\alpha }|{\omega }_{k1}⟩{|}^2\delta ({\omega }_{k1}-{\omega }_{k2}+\omega ),$$

$$g_{\perp }=g_J|⟨{\omega }_{0,\pm }|J_{\pm }|{\omega }_{0,\mp }⟩|$$

$$g_{\parallel }=2g_J|⟨{\omega }_{0,\pm }|J_z|{\omega }_{0,\pm }⟩|$$

$$B_n^m=C_{\mathrm{Yb}}(n,m)\sum _i\frac{Z_i{Y}_n^m({\theta }_i,{\phi }_i)}{R_i^{n+1}},$$

Full text

Introduction.—Quantum spin liquid (QSL) is a novel state of matter with zero entropy and without conventional symmetry breaking even at zero temperature. Such states were proposed to host ’spinons’, exotic spin excitations with fractional quantum numbers Balents (2010); Lee (2008); Wen (2004). Although many candidate QSL materials with two-dimensional or three-dimensional interaction topologies on the triangular, kagome, and pyrochlore lattices were reported Shimizu et al. (2003); Yamashita et al. (2008, 2009); Itou et al. (2008, 2010); Yamashita et al. (2011); Shores et al. (2005); Han et al. (2012); Fu et al. (2015); Fåk et al. (2012); Clark et al. (2013); Li et al. (2014); Hodges et al. (2002); Ross et al. (2011), they typically suffer from magnetic or non-magnetic defects Lee et al. (2007); Freedman et al. (2010); Li and Zhang (2013); Li et al. (2013); Gomilšek et al. (2016), spatial anisotropy Shimizu et al. (2003); Itou et al. (2008); Li et al. (2014), antisymmetric Dzyaloshinsky-Moriya anisotropy Moriya (1960); Zorko et al. (2008); Fortune et al. (2009), and (or) interlayer magnetic couplings Fortune et al. (2009); Doi et al. (2004); Shirata et al. (2012) that reduce or even completely release magnetic frustration Fortune et al. (2009); Shirata et al. (2012); Dommange et al. (2003); Rigol and Singh (2007); Rousochatzakis et al. (2009).

Many of the aforementioned shortcomings can be remedied in a new triangular antiferromagnet YbMgGaO4 that was recently reported by our group Li et al. (2015a, b, 2016). No spin freezing was detected down to at least 0.048 K, which is about 3% of the nearest-neighbor interaction $J_{0}\sim 1.5$ K Li et al. (2016). Residual spin entropy is nearly zero at 0.06 K, excluding any magnetic transitions at lower temperatures Li et al. (2015a). Below 0.4 K, thermodynamic properties evidence the putative QSL regime with temperature-independent magnetic susceptibility $\chi=\text{const}$ Li et al. (2016) and power-law behavior of the magnetic heat capacity, $C_{m}\sim T^{2/3}$ Li et al. (2015a), the observations that are consistent with theoretical predictions for the U(1) QSL ground state (GS) on the triangular lattice Motrunich (2005); Lee and Lee (2005); Motrunich (2006).

Very recently, two inelastic neutron scattering (INS) studies of YbMgGaO4 Shen et al. (2016); Paddison et al. (2017) reported continuous excitations at transfer energies of $0.1-2.5$ meV extending well above the energy scale of the magnetic coupling $J_{0}\sim 0.13$ meV. These spectral features were identified as fractionalized excitations (’spinons’) from the QSL GS Shen et al. (2016). Surprisingly, though, magnetic excitations remain very broad in both energy and wave-vector ($Q$) even in the almost fully polarized state at 7.8 T, where only narrow spin-wave excitations of an ordered ferromagnet are expected Paddison et al. (2017). This persistent broadening of magnetic excitations may be related to a very inconspicuous intrinsic structural disorder that we uncover and quantify by INS measurements at high energies, where crystalline electric field (CEF) excitations of Yb3+ ions can be probed.

In this Letter, we report a comprehensive investigation of the CEF excitations in YbMgGaO4. They are observed at the energies of $\hbar\omega$ = 39, 61, and 97 meV and show not only a pronounced broadening, but also a very peculiar energy profile with a shoulder around 87 meV on the side of the 97 meV excitation. These peculiarities are rationalized by considering the frozen Mg/Ga disorder that affects the local environment of Yb3+ and, thus, the CEF parameters. Their randomness gives rise to a distribution of electronic $g$-factors and explains the broadening of low-energy magnetic excitations, thus rendering structural randomness an important ingredient of the QSL physics in YbMgGaO4.

Experimental technique.—Moderate-high-energy INS data and low-energy INS data were collected, respectively, on the MERLIN Bewley et al. (2006) and LET Bewley et al. (2011) spectrometers at the ISIS pulsed neutron facility, Rutherford Appleton Laboratory, U.K. sup . Several incident energies $E_{i}$ were used at MERLIN. In the following, we focus on the data obtained with $E_{i}=153.5$ meV that provides the best trade-off between the resolution and energy coverage 111Note that the data obtained with $E_{i}=307$ meV extend to much higher energies and reveal a weak excitation, presumably of CEF origin, at 139 meV. This feature has only 20% of the intensity of the 97 meV excitation and is likely due to a multiple scattering excitation, see Supplemental material for details. Our YbMgGaO4 ($m_{\rm Yb}$ = 14.03 g) and LuMgGaO4 ($m_{\rm Lu}$ = 6.36 g) powder samples for the MERLIN experiment were prepared using solid-state reactions Li et al. (2015a). Single crystals of YbMgGaO4 for the LET experiment were grown by the floating zone technique Li et al. (2015b).

CEF excitations.—According to the Hund’s rules, free Yb3+ (4$f^{13}$) ion has the spin angular momentum $s=1/2$ and orbital angular momentum $L=3$ resulting in the eight-fold-degenerate ground state with the total angular momentum $J=7/2$ and Landé $g$-factor $g_{J}=8/7$ for the GS multiplet. In the idealized YbMgGaO4 structure, Yb3+ has trigonal local symmetry with the point group D3d (see Fig. 1(a)) that splits this multiplet into four Kramers doublets sup (see Fig. 1(b)).

Raw INS spectra for both YbMgGaO4 and its nonmagnetic analog LuMgGaO4 are shown in panels (a) and (b) of Fig. 2, respectively. Their comparison reveals three features that are identified as CEF excitations of Yb3+ based on the following observations: (1) These excitations are absent in the non-magnetic reference compound LuMgGaO4 (Fig. 2 (b)). (2) The lowest-lying excitation at around $\sim$ 39.4 meV (see Fig. 2 (c)) is consistent with the energy separation $\Delta$ $\sim$ 36.5(1) meV between the ground-state Kramers doublet and the first excited state, as found in our previous heat capacity measurements Li et al. (2015b). (3) No systematic anharmonic effect is observed, thus excluding phonon origin of the excitations sup ; Balkanski et al. (1983). (4) $Q$-independent excitation energies (see Fig. 2 (a)) suggest their single-ion nature. (5) The intensities decrease with $Q$ following the square of the magnetic form factor of the Yb3+ ion (see Fig. 3 (b)). (6) The lowest excitation at around 39 meV is far above $J_{0}\simeq 0.13$ meV and, therefore, unrelated to the spin-spin correlations in YbMgGaO4. All these facts indicate that three spectral features are single-ion CEF excitations. For an Yb3+ ion with $J=7/2$ in the D3d symmetry, we indeed expect four CEF doublets and thus three CEF excitations from the ground-state doublet.

A closer inspection of these CEF excitations reveals two unusual features, though. First, all excitations are much broader than the instrumental resolution. For example, at $E_{i}$ = 153.5 meV the total FWHMs (full width at half maximum obtained from the convoluted Lorentzian and Gaussian peak profiles) are 10.1(4) meV ($\hbar\omega_{1}\sim 39$ meV), 10.9(4) meV ($\hbar\omega_{2}\sim 61$ meV), and 12.2(7) meV ($\hbar\omega_{3}\sim 97$ meV), much larger than the instrumental resolutions (Gaussian component) of 6.7, 5.6, and 4.3 meV, respectively. Through convolution calculations, we determine the additional broadening (Lorentzian component) of 5.5, 7.9, and 10.7 meV, respectively. Given the high quality of our sample sup and the low temperature of the measurement ($T=5$ K), we conclude that this broadening is intrinsic.

Another peculiar feature is the shape of the highest CEF excitation that shows the main peak around 97 meV and a shoulder at $\sim$ 87 meV (see Fig. 2 (d)). The $Q$-dependence of the intensity at $\sim$ 87 meV follows the square of the magnetic form factor of the Yb3+ ion sup , thus proving the CEF origin of this spectral feature. It contributes about 40% of the overall intensity of the highest-energy excitation and is clearly intrinsic. Further, there are no phonon modes observed between $\sim$ 70 and 120 meV in LuMgGaO4 (see Fig.  2(b)) and hence the $\sim$ 87 meV shoulder could not be due to CEF-phonon coupling Adroja et al. (2012).

Combined CEF fit.—To determine the CEF parameters, we fit energy dependence of the experimental intensity in four fashions (see Fig. 3 (a) and sup ). Fit #1 is performed against the whole dataset with a three-peaks fit (see Fig. 3 (a)), for fit #2 we excluded the region between 73 and 90 meV, where the additional shoulder is observed, while fit #3 is performed against the whole dataset with a four-peaks fit and uses the additional shoulder energy ($\sim$ 87 meV) as $\hbar\omega_{3}$, and fit #4 uses the additional mode energy (138.7(3) meV) as $\hbar\omega_{3}$. All fits share the same measured relative INS intensities, $I_{2}/I_{1}$ and $I_{3}/I_{1}$ (see Fig. 3 (b)), and the same measured effective spin- 1/2 g-factors Li et al. (2015b). Through combined fits to these seven observables – $\hbar\omega_{1}$, $\hbar\omega_{2}$, $\hbar\omega_{3}$, $I_{2}/I_{1}$, $I_{3}/I_{1}$, $g_{\perp}$, and $g_{\parallel}$ – we obtain all six CEF parameters, $B_{n}^{m}$ sup , by minimizing the following deviation function,

[TABLE]

where $X_{i}^{\rm obs}$ and $\sigma_{i}^{\rm obs}$ are the experimental value and its standard deviation, respectively, whereas $X_{i}^{\rm cal}$ is the calculated value. Qualitatively similar CEF parameters and wavefunctions are found from fits #1, #2, #3 and #4, as shown in sup , respectively. Magnetic part of the experimental heat capacity ($C_{m}$) is very well reproduced with the first three parameter sets (see Fig.  3 (c)):

[TABLE]

Inherent structural disorder and CEF randomness.—Peculiar shape of the CEF excitations is rooted in subtle details of the YbMgGaO4 crystal structure. Our single-crystal x-ray diffraction study excludes any global symmetry reduction or a site mixing between Yb and Mg/Ga sup . On the other hand, Mg and Ga share one crystallographic site, thus forming different local configurations around each Yb3+ ion. The most obvious effect of this Mg/Ga disorder is the variation of the electrostatic potential imposed on Yb3+. We assess it by calculating CEF parameters using the point-charge model sup and find that, as long as Yb occupies its ideal position at (0,0,0), the random distribution of Mg and Ga gives rise to only a weak broadening of the CEF excitations, $\Delta(\hbar\omega_{1})$ = 0.27 meV, $\Delta(\hbar\omega_{2})$ = 0.26 meV, and $\Delta(\hbar\omega_{3})$ = 0.39 meV, and does not account for our experimental observations. Moreover, all three CEF excitations remain symmetric.

A further effect of the Mg/Ga disorder is local charge misbalance that may push Yb out of its ideal position, as reflected by the enhanced values of the Yb atomic displacement parameter, with the thermal ellipsoid elongated along the $c$-direction sup ; Cava et al. (1998). We probed this effect quantitatively by constructing several representative Mg/Ga configurations and optimizing their geometry using density-functional calculations sup . We indeed observed that exact positions of both Yb and its neighboring oxygens are affected by the local distribution of Mg2+ and Ga3+. The resulting distortions of the YbO6 octahedra give rise to a pronounced distribution of the CEF parameters and render the highest-lying CEF excitation asymmetric sup . Both the $\sim$ 87 meV shoulder and the overall broadening of the CEF excitations can be well reproduced (see Fig. 3 (d) and sup ). Consequently, the effective spin-1/2 $g$-factors show a maximum distribution as follows: $\Delta g_{\perp}$ $\sim$ 0.3 and $\Delta g_{\parallel}$ $\sim$ 1.2 sup . This distribution of the $g$-values has immediate ramifications for low-energy excitations, as we show below. Eventually, intersite magnetic couplings should be random too Ross et al. (2012), although in a much more complicated manner.

Prior to discussing the low-energy excitations, let us note that the broadening of CEF excitations is not uncommon for rare-earth compounds Princep et al. (2015); Babkevich et al. (2015); Wen et al. (2016). This effect is usually ascribed to antisite defects, as in the “stuffed” quantum spin ice, Yb2Ti2-xYbxO7-x/2 with $x\sim 0.01-0.02$ Ross et al. (2012); Gaudet et al. (2015, 2016). On the other hand, site mixing beyond the rare-earth site, as in Tb2Sn2-xTixO7, is believed to merge the CEF excitations into a broad continuum Gaulin et al. (2015). Interestingly, YbMgGaO4 with its complete Mg/Ga disorder and without any detectable Yb antisite defects retains well-defined CEF excitations, albeit with a peculiar energy profile that can be reproduced, perhaps for the first time, by considering local atomic relaxation depending on the distribution of Mg and Ga around Yb3+.

Spin-wave continuum.—Under the field of 8.5 T applied along the c-axis, the spin system is fully polarized according to the static magnetization measurement at 1.9 K sup . Therefore, at 8.5 T and 0.1 K the fully polarized (ferromagnetic) state should give rise to narrow spin-wave excitations having the width of about 0.16 meV according to the instrumental resolution of LET. In contrast, the measured INS signals are still broadly distributed in the energy space with a width of more than 0.5 meV (see Fig. 4 (a)). This width is obviously larger than the instrumental resolution Bewley et al. (2011) and than the width of spin-wave excitations in single crystals of a similar Yb3+ material, Yb2Ti2O7 Ross et al. (2011). Similarly broad excitations were observed in the recent INS measurement performed in the applied field of 7.8 T at 0.06 K Paddison et al. (2017). Our magnetization data sup confirm that at 8.5 T the Yb3+ spins are fully polarized along the c-axis. Nevertheless, the excitations remain very broad. Thus, the natural explanation to this observed spin-wave continuum is the aforementioned randomness of the effective spin-1/2 $g$-factors and (or) couplings, since linear spin-wave theory should be applicable to YbMgGaO4 at 8.5 T and 0.1 K.

We model the spin-wave excitations in the ab-plane under the high applied field along the c-axis using the expression Xu et al. (2013)

[TABLE]

where $F(\mid\textbf{Q}\mid)$ is the magnetic form factor of Yb3+ and $S^{\alpha\alpha}(\textbf{Q},E)$ ($\alpha$ = $x$, $y$ or $z$) is the dynamic spin structure factor calculated by the Spinw-Matlab code based on the linear spin-wave theory Toth and Lake (2015) and coupling parameters reported earlier Li et al. (2015b). By considering the maximum broadening of $g_{\parallel}$, $\Delta g_{\parallel}$ = 1.2, and the instrumental Gaussian broadening, $\sigma$ = 0.16 meV, we are able to reproduce the broadening of about 0.7 meV at the zone center. However, the signal is even broader at the zone boundary. While the on-site randomness is clearly very important for the low-energy physics, a distribution of magnetic couplings is relevant too and requires further investigation.

Conclusions.—The CEF excitations of Yb3+ in the triangular QSL YbMgGaO4 have been studied by moderate-high-energy INS measurements. Large broadening and peculiar energy profile of the CEF excitations is observed and ascribed to the structural randomness, namely, the random distribution of Mg and Ga that affects local coordination of the Yb3+ ions. We propose that this inherent structural disorder results in the distribution of the effective spin-1/2 $g$-factors that is responsible for the persistent broadening of low-energy magnetic excitations in the fully polarized ferromagnetic state, although the distribution of magnetic couplings seems to be relevant too. Our results put forward structural randomness as an important ingredient in the spin-liquid physics, an observation that goes hand in hand with the recent report on the suppression of thermal conductivity at low temperatures Xu et al. (2016). In 4$f$-based materials, the randomness of CEF levels can be easily introduced without generating strong structural disorder, thus opening interesting prospects for the design of new spin-liquid materials Sanders et al. (2016).

Note added.—After the submission of our paper, the results of Ref. Paddison et al., 2017 were updated as it progressed from a preprint to a published article incorporating high-energy INS data that are largely consistent with our results.

Acknowledgment—We thank Jinchen Wang and Yu Li for helpful discussion. YL thanks Duc Le for introducing him to the Mantid (3.7.1) software for CEF fits. AT acknowledges Tom Fennell for pointing our recent examples of CEF excitations broadening. This work was supported by the NSF of China and the Ministry of Science and Technology of China (973 projects: 2016YFA0300504). Y.S.L. was supported by the start-up funds of Renmin University of China. The work in Augsburg was supported by the German Federal Ministry for Education and Research through the Sofja Kovalevskaya Award of Alexander von Humboldt Foundation. Q.M.Z. was supported by the Fundamental Research Funds for the Central Universities, and by the Research Funds of Renmin University of China.

*Supplementary material

*Crystalline Electric Field Randomness in the Triangular Lattice Spin-Liquid YbMgGaO4

We present here:

1. General techniques.

1.1 Moderate-high-energy MERLIN inelastic neutron scattering (INS) experiments.

1.2 Low-energy LET INS experiments.

1.3 Crystalline electric field (CEF) Hamiltonian and its relationship to observations.

1.4 CEF fits.

1.5 Calculation of the CEF parameters using the point-charge model.

2. X-ray diffraction (XRD) patterns for YbMgGaO4 and LuMgGaO4 powders used in the MERLIN INS experiments.

3. First-derivative absorption electron spin resonance (ESR) spectra for (Yb0.04Lu0.96)MgGaO4 and YbMgGaO4 Li et al. (2015a).

4. Shape of the 97 meV CEF excitation.

5. Energy dependence of the CEF excitations measured for YbMgGaO4 from 5 to 300 K.

6. Combined CEF fit results for YbMgGaO4.

7. Fitted CEF wavefunctions.

8. Crystal structure of YbMgGaO4 at 100 K determined from single-crystal XRD data.

9. Local structural configurations of Yb3+ obtained from density-functional band-structure calculations.

10. Calculated CEF randomness in the framework of the point-charge model.

10.1 Calculated CEF randomness caused by the distribution of nearest-neighbor oxygen environments.

10.2 Calculated CEF randomness caused solely by the site-mixing between Mg2+ and Ga3+.

11. INS spectra measured at $E_{i}$ = 307 meV, and the 138.7(3) meV mode.

12. Elastic neutron diffraction pattern for the single crystal sample of YbMgGaO4 in the ab-plane used in LET INS experiments.

13. Magnetic field dependence of the susceptibility ($dM/dH$) measured at 1.9 K along the c-axis.

I 1. general techniques

I.1 1.1 Moderate-high-energy Merlin INS experiments

The MERLIN measurements between 5 and 300 K were performed on both YbMgGaO4 and LuMgGaO4 samples with the incident beam energies $E_{i}$ = 92.5, 153.5, and 307 meV. At $E_{i}$ = 92.5 meV, the instrumental resolutions (full width at half maximum, FWHM) are 2.9 and 2.3 meV at the transfer-energies of 39 and 61 meV, respectively, whereas at $E_{i}$ = 153.5 meV, the resolutions are 6.7, 5.6, and 4.3 meV at 39, 61, and 97 meV, respectively Bewley et al. (2006). And at $E_{i}$ = 307 meV the instrumental resolution is 22.5 meV at $\hbar\omega$ $\sim$ 0 meV. Phase purity of both samples was confirmed by X-ray diffraction prior to the MERLIN INS measurements (see Fig. 5). All MERLIN data were processed and cut using the Mslice-Matlab code.

I.2 1.2 Low-energy LET INS experiments

Incident energies of 26.8 and 5.5 meV were chosen for both elastic and inelastic scattering with the energy resolution of 1.4 and 0.16 meV, respectively Bewley et al. (2011). The as-grown rod ($\sim$ 50 g) was cut into slices along the ab-plane. Ten best-quality ab-slices of the single-crystal (total mass $\sim$ 10 g) were selected for the neutron scattering experiment on LET by Laue X-ray diffraction on all surfaces. The slices were fixed to the copper base by CYTOP glue to avoid any shift in applied magnetic fields up to 8.5 T. The sample temperature of 0.1 K was achieved using a dilution refrigerator. The elastic neutron diffraction (see Fig. 14) showed that the alignment of the single crystals was sufficient for the INS study of the continuous excitations.

All INS data shown in the main text have been integrated over the momentum space, $-0.9\leq\eta\leq 0.9$ in $[0,0,-\eta]$. In Fig. 4 (a) of the main text, the INS data have been also integrated over $-0.03\leq\xi\leq 0.03$ in $[\xi/2,-\xi,0]$. In Fig. 4 (c), (e) and (g) of the main text, the INS data have been integrated over small energy ranges, (c) 1.15 $\leq$ $E$ $\leq$ 1.25 meV, (e) 1.65 $\leq$ $E$ $\leq$ 1.75 meV and (g) 2.05 $\leq$ $E$ $\leq$ 2.15 meV, respectively. All LET data were processed and cut using the Horace-Matlab on the ISIS computers.

I.3 1.3 Crystalline electric field (CEF) Hamiltonian and its relationship to observations

Under zero applied magnetic field, the generic CEF Hamiltonian that is invariant under the D3d point group symmetry is given by Bertin et al. (2012); Mirebeau et al. (2007); Ruminy et al. (2016)

[TABLE]

where $B_{n}^{m}$ ($n,m$ are integers and $n\geq m$) are CEF parameters that will be determined experimentally, and the Stevens operators $O_{n}^{m}$ are polynomial functions of the components of the total angular momentum operator $J_{z}$, $J_{+}$, and $J_{-}$ ($J_{\pm}=J_{x}\pm iJ_{y}$). For the moderate-high-energy INS measurements ($>10$ meV), it is sufficient to consider only the single-ion CEF Hamiltonian [Eq. (4)] and neglect the couplings between the Yb3+ ions ($J_{0}\sim 0.13$ meV) Li et al. (2015b). Under the eight-dimensional representation space $|J=7/2,m_{J}\rangle$ ($m_{J}$ = $\pm$1/2, $\pm$3/2, $\pm$5/2 and $\pm$7/2), we diagonalize the single-ion CEF Hamiltonian and obtain four ascending energy levels, $\hbar\omega_{k}$ (k = 0, 1, 2 and 3). Each energy level is two-fold-degenerate with the eigenfunctions (Kramers doublet)

[TABLE]

where $|C_{m_{J}}^{k,\pm}|=|C_{-m_{J}}^{k,\mp}|$.

In a paramagnetic system (single-ion approximation), the differential neutron cross section of the powder sample for the CEF Kramers doublet transition $|\omega_{k1}\rangle\rightarrow|\omega_{k2}\rangle$ is given in the dipole approximation by Bertin et al. (2012); Mirebeau et al. (2007); Ruminy et al. (2016); Furrer et al. (1988); Mazumdar et al. (2008); Chi et al. (2008); Babkevich et al. (2015):

[TABLE]

where $F(Q)$ = $\langle j_{0}(\frac{Q}{4\pi})\rangle$ + $\frac{2-g_{J}}{g_{J}}\langle j_{2}(\frac{Q}{4\pi})\rangle$ is the magnetic form factor in the dipole approximation, and $|\langle\omega_{k2}|J_{\alpha}|\omega_{k1}\rangle|^{2}=|\langle\omega_{k2,+}|J_{\alpha}|\omega_{k1,+}\rangle|^{2}+|\langle\omega_{k2,-}|J_{\alpha}|\omega_{k1,+}\rangle|^{2}+|\langle\omega_{k2,+}|J_{\alpha}|\omega_{k1,-}\rangle|^{2}+|\langle\omega_{k2,-}|J_{\alpha}|\omega_{k1,-}\rangle|^{2}$ is proportional to the total transition probability. At 5 K, only three CEF excitations from the ground Kramers doublet are taken into account, and the CEF transitions from the excited doublets are ignored.

At low temperatures ($T\ll\hbar(\omega_{1}-\omega_{0})/k_{B}$), the effective spin-1/2 $g$-factors for the ground-state Kramers doublet ($|\omega_{0,\pm}\rangle$) are given by:

[TABLE]

for the field perpendicular and parallel to the $c$-axis, respectively

I.4 1.4 CEF fits

All CEF simulations and combined fits to the experimental data (excitation energies, relative INS intensities, and effective spin-1/2 g-factors) were performed using Matlab programs. Our Matlab code for the CEF calculations (excitation energies and relative INS intensities) is fully compatible with the widely used code Mantid (3.7.1).

I.5 1.5 Calculation of the CEF parameters using the point-charge model

CEF parameters can be calculated for a given crystal structure within the point-charge approximation using the expression Baldoví et al. (2013)

[TABLE]

where $Z_{i}$ is the net charge (unit: $e$) of the lattice ion $i$ with the spherical coordinates ($R_{i}$,$\theta_{i}$,$\varphi_{i}$), $Y_{n}^{m}$ are spherical harmonics, and $C_{\rm Yb}(n,m)$ is a pre-factor that depends on $n$ and $m$ only. Since $B_{n}^{m}\propto\frac{1}{R_{i}^{n+1}}$ ($n$ $\geq$ 2), ions beyond nearest neighbors can be neglected. On the other hand, the point-charge calculation does not take into account effects of exchange fields and may not deliver exact values of $B_{n}^{m}$. Here, we used the point-charge calculations to obtain the distribution of individual CEF parameters for different scenarios of structural disorder, and eventually re-scaled average values of these distributions to the experimental (fit #1) CEF parameters listed in Table 1.

II 2. XRD patterns for YbMgGaO4 and LuMgGaO4 powders used in the MERLIN INS experiments

Fig. 5 shows the XRD patterns for the YbMgGaO4 and LuMgGaO4 powder samples used in the MERLIN INS measurements. No reflections of foreign phases are observed, suggesting that the concentration of possible Yb3+-related impurities, such as Yb3Ga5O12 and Yb2O3, should be well below $\sim$ 5%.

III 3. First-derivative absorption ESR spectra for (Yb0.04Lu0.96)MgGaO4 and YbMgGaO4 Li et al. (2015a)

If a fraction of the magnetic Yb3+ ions occupies the Mg2+/Ga3+ site, quasi-free spins of the Yb3+ defects should be detected by ESR. The narrow and strong hyperfine lines Misra and Isber (1998) from quasi-free defect Yb3+ spins observed in (Yb0.04Lu0.96)MgGaO4 completely disappear in YbMgGaO4 (see Fig. 6). From the comparison of the normalized signal intensities, we can estimate that the concentration of the Yb3+-related defects in YbMgGaO4 should be less than $\sim$ 0.04%.

IV 4. Shape of the 97 meV CEF excitation

The zoom-in view of the 97 meV CEF excitation is shown in Fig. 2 (d) of the main text. It has a local maximum of the intensity around 97 meV, and the one-peak fit to this CEF excitation gives an averaged peak center at 96.6(2) meV (see fit #1). However, the energy profile of this excitation is strongly asymmetric, with a clear shoulder at $\sim$ 87 meV.

The momentum transfer ($\mid$Q$\mid$) dependence of the integrated INS intensities at $\sim$ 87 meV and $\sim$ 97 meV are shown in Fig. 7. Both intensities follow the magnetic form factor of the free Yb3+ very well (see Fig. 7), suggesting the same CEF origin of these spectral features.

The INS intensity of the shoulder at $\sim$ 87 meV is about 40% of that around 97 meV (see Fig. 7). Neither Yb3+-related impurities (such as Yb3Ga5O12 and Yb2O3, $<$ 5%) nor defects (such as the site mixing between Yb3+ and Mg2+/Ga3+, $<$ 0.04%) can account for this shoulder at $\sim$ 87 meV. Thus, the $\sim$ 87 meV shoulder is intrinsic and reflects a peculiar distribution of the CEF parameters in YbMgGaO4.

V 5. Energy dependence of the CEF excitations measured for YbMgGaO4 from 5 to 300 K

No systematic anharmonic effect is observed excluding phonon origin of the excitations at $\hbar\omega$ = 39, 61, and 97 meV (see Fig. 8)  Balkanski et al. (1983).

VI 6. Combined CEF fit results for YbMgGaO4

VII 7. Fitted CEF wavefunctions

VIII 8. Crystal structure of YbMgGaO4 at 100 K determined from single-crystal XRD data

In Table 3, we provide average crystal structure of YbMgGaO4 determined from single-crystal XRD data at 100 K (the data are from the Ref. Li et al. (2015b)). When Yb3+ is placed into its ideal position at $(0,0,0)$, the displacement parameter U33 is nearly twice larger than U11 = U22. This reflects static disorder of Yb3+ along the c-axis Cava et al. (1998), because the large U33 of 0.0240 Å2 persists even at 100 K, where dynamic effects, such as phonons, should be minor. Mg2+ and Ga3+ are intermixed. On the other hand, no evidence for Yb antisite defects (for example, Yb atoms occupying the Mg/Ga position) is found.

IX 9. Local structural configurations of Yb3+ obtained from density-functional band-structure calculations

In order to go beyond the average structure from single-crystal XRD, we constructed several ordered Mg/Ga configurations and optimized their geometry using density-functional (DFT) band-structure calculations. All calculations were performed in the VASP code Kresse and Furthmüller (1996a, b) with the Perdew-Burke-Ernzerhof (PBE) flavor of the exchange-correlation potential Perdew et al. (1996) and the default Yb pseudopotential, where $4f$ states were placed into the core. The $4\times 4\times 4$ k-mesh was used, and residual forces were below 0.005 eV/Å in fully optimized structures. Lattice parameters were fixed to their experimental values.

The local configurations considered in our analysis are periodic and charge-neutral. Additionally, we restricted ourselves to those configurations that preserve the 3-fold symmetry of the Yb3+ site. While many other configurations violating this symmetry can be constructed, the 3-fold-symmetric configurations keep the problem tractable and are in fact sufficient to reproduce effects of CEF randomness in YbMgGaO4. We considered four different Mg/Ga distributions within the trigonal unit cell of YbMgGaO4 (Fig. 9). These four structures contain all 7 possible environments of Yb, where we define the environment with respect to the Mg/Ga distribution in the two neighboring slabs. This way, 7 A-B-Yb-C-D configurations, where each of the ABCD is either Mg or Ga, are obtained (Table 4). The restriction to only two neighboring slabs is well justified. For example, the MgMgYbGaGa configuration is independently probed in the structures shown in panels (b) and (d) of Fig. 9 resulting in less than 0.5 % deviation for the ensuing geometry of the YbO6 octahedron.

We find that the exact position of Yb as well as its oxygen coordination are strongly affected by the Mg/Ga distribution. The local configurations MgGaYbGaMg and GaMgYbMgGa retain inversion symmetry at the Yb site and keep the YbO6 octahedron undistorted. Any other configuration entails electrostatic field that pushes Yb atoms out of their ideal $(0,0,0)$ position. This effect is also detected experimentally by the largely elongated thermal ellipsoid in the averaged crystal structure of YbMgGaO4 (Table 3). Consequently, the YbO6 octahedra are distorted, as shown in Table 4.

X 10. Calculated CEF randomness in the framework of the point-charge model

X.1 10.1 Calculated CEF randomness caused by the distribution of nearest-neighbor oxygen environments

In Table 5, we list CEF parameters calculated using the point-charge model for 7 flavors of the YbO6 octahedra listed in Table 4. Relative distributions listed in Table 5 are then applied to obtain the CEF spectra by setting $\langle B_{n}^{m}\rangle$ = $B_{n}^{m}$(fit #1) (see Table 1 for each $B_{n}^{m}$(fit #1)). Each type of the YbO6 octahedron produces its three CEF excitations at slightly different energies (Fig. 10). The cumulative effect of this distribution is the symmetric broadening of the two lower CEF excitations and the complex profile of the third (highest-energy) CEF excitation, where the configurations 4 and 6 are responsible for the shoulder on the low-energy side (around 87 meV), whereas the configuration 5 would produce the shoulder on the high-energy side (around 105 meV).

The resulting energy profile is remarkably similar to our experimental observations, as further demonstrated by the 2D plot in Fig. 3 (d) in the main text. We note, however, that the much weaker (about a half of the intensity of the $\sim$ 87 meV shoulder) high-energy shoulder around 105 meV was not observed at its expected position. On the other hand, we do observe a weak excitation at a higher energy of 139 meV (see below), which may be related to the high-energy shoulder.

X.2 10.2 Calculated CEF randomness caused solely by the site-mixing between Mg2+ and Ga3+

We have also used the point-charge model to calculate the CEF parameters for different Mg/Ga distributions, while keeping all atoms in their ideal positions (see Table 3), i.e., neglecting effects of local structural relaxation. The distribution of the CEF parameters (see Table 6) appears to be minor compared to Table 5 discussed above. This minor distribution broadens the CEF excitations by $\Delta(\hbar\omega_{1})$ = 0.27 meV, $\Delta(\hbar\omega_{2})$ = 0.26 meV and $\Delta(\hbar\omega_{3})$ = 0.39 meV, which is much smaller than the broadenings of several meV observed experimentally. Therefore, the distribution of Mg/Ga does not explain the CEF randomness, and effects of local structural relaxation are crucial.

XI 11. INS spectra measured at $E_{i}$ = 307 meV: The 138.7(3) meV mode

The incident energy of $E_{i}$ = 307 meV gives access to a much broader energy range and reveals a weak excitation around 139 meV (see Fig. 11). After the subtraction by the corresponding INS signal of LuMgGaO4, this $\sim$ 139 meV mode becomes more clear, and its position is fitted to be 138.7(3) meV (see Fig. 12). Three main CEF excitations are confirmed, although they show significant broadening due to the poor instrumental resolution at $E_{i}$ = 307 meV ($\sigma$ $\sim$ 22.5 meV at $\hbar\omega$ $\sim$ 0 meV) (see Fig. 12 (b)). The intensity of the 138.7(3) meV mode is fitted to be $I_{139meV}/I_{3}$ = 0.21(1) (see Fig. 13), which is much smaller than that of the $\sim$ 87 meV shoulder, $I_{87meV}/I_{3}$ = 0.40(2) (see Fig. 7). Therefore, $I_{139meV}/I_{87meV}$ $\sim$ 0.5.

$|Q|$-dependence of the intensity of the 139 meV mode confirms its CEF nature, although a sizable background of unknown origin can be seen from the non-zero value of $b$ in Fig. 13. We assumed that this mode may be a trace of the third CEF excitations and used it for the CEF fit #4, where the energy of the third excitation was set to 139 meV (see Tables 1 and 2). Owing to its low intensity, this mode may be the third CEF excitation for a small fraction of the Yb3+ ions only (up to 20 %). Its intensity of one half of the 87 meV shoulder corresponds in fact to the 105 meV high-energy shoulder in our calculated spectra on Fig. 10. An alternative and perhaps even more plausible explanation for the 139 meV mode would be a multiple scattering CEF excitation, such as $\hbar\omega_{1}+\hbar\omega_{3}=136$ meV and $2\hbar\omega_{1}+\hbar\omega_{2}=139$ meV, because neutron can experience several scattering events before leaving the sample Babkevich et al. (2015).

We also note that our sequence of CEF excitations, including the weak 139 meV mode, is well in line with the recent report by Paddison et al. Paddison et al. (2017).

XII 12. Elastic neutron diffraction pattern for the single crystal sample of YbMgGaO4 in the ab-plane used in LET INS experiments

XIII 13. Magnetic field dependence of the susceptibility ($dM/dH$) measured at 1.9 K along the c-axis