Crystal field splitting, local anisotropy, and low energy excitations in the quantum magnet YbCl$_3$

TL;DR

This study investigates YbCl3's magnetic properties, revealing a well-isolated ground state doublet with easy-plane anisotropy and low-energy excitations indicative of reduced-dimensional antiferromagnetic correlations.

Contribution

The paper determines the crystal field Hamiltonian and characterizes the low-energy magnetic excitations in YbCl3, providing insights into its anisotropic spin-1/2 behavior.

Findings

Ground state doublet is well isolated from other levels.

YbCl3 exhibits local easy-plane anisotropy.

Low energy excitations suggest reduced dimensionality of antiferromagnetic correlations.

Abstract

We study the correlated quantum magnet, YbCl$_3$, with neutron scattering, magnetic susceptibility, and heat capacity measurements. The crystal field Hamiltonian is determined through simultaneous refinements of the inelastic neutron scattering and magnetization data. The ground state doublet is well isolated from the other crystal field levels and results in an effective spin-1/2 system with local easy plane anisotropy at low temperature. Cold neutron spectroscopy shows low energy excitations that are consistent with nearest neighbor antiferromagnetic correlations of reduced dimensionality.

Tables

Table 1. Table 1: Refined crystal field parameters in units of meV determined as described in the text. Each coefficient is presented divided by the corresponding Steven’s parameter α J subscript 𝛼 𝐽 \alpha_{J} , β J subscript 𝛽 𝐽 \beta_{J} and γ J subscript 𝛾 𝐽 \gamma_{J}   Hutchings ( 1964 ) .

$B_2^0$	$B_2^2$	$B_4^0$	$B_4^2$	$B_4^4$
$-2.820$	$-25.956$	$6.170$	$42.336$	$-36.335$
$B_6^0$	$B_6^2$	$B_6^4$	$B_6^6$
$-3.004$	$10.764$	$7.482$	$49.327$
$B{\left(i\right)}_4^2$	$B{\left(i\right)}_4^4$	$B{\left(i\right)}_6^2$	$B{\left(i\right)}_6^4$	$B{\left(i\right)}_6^6$
$-6.42\times {10}^{-3}$	$-0.015$	$8.56\times {10}^{-3}$	$-0.067$	$-0.036$

Table 2. Table S2: Summary of Rietveld refinement parameters for YbCl 3 . Data were refined in the space group C ​ 2 / m 𝐶 2 𝑚 C2/m (12).

	T=10 K  $R_w$=2.42  a=6.7291(3) Å, b=11.6141(3) Å, c=6.3121(3) Å $\alpha =\gamma ={90}^{\circ }$, $\beta =110.5997{\left(7\right)}^{\circ }$				T=100 K  $R_w$=2.93  a=6.7326(5) Å, b=11.6201(1) Å, c=6.3280(5) Å $\alpha =\gamma ={90}^{\circ }$, $\beta =110.551{\left(1\right)}^{\circ }$
Atom	x	y	z	Occ	x	y	z	Occ
Yb	0	0.1663(6)	0	1.0	0	0.1665(1)	0	1.0
Cl(1)	0.2589(1)	0.3214(1)	0.2431(1)	0.992(3)	0.2587(1)	0.3209(1)	0.2422(1)	0.997(4)
Cl(2)	0.2188(1)	0	0.2493(1)	0.997(4)	0.2179(2)	0	0.2487(3)	0.991(6)

Table 3. Table S3: Tabulated wave functions of the crystal field states in YbCl 3 obtained within the LS-coupling approximation. The crystal-field energies (in meV) are tabulated horizontally, the m J subscript 𝑚 𝐽 m_{J} -values of the ground-state multiplet vertically. Only coefficients of the wave functions > 10 − 3 absent superscript 10 3 >10^{-3} are shown. The imaginary parts of the eigenvectors are omitted because they are ≈ 2 absent 2 \approx 2 order of magnitude smaller than the real ones.

$m_J$	$0.0$	$0.0$	$21.08$	$21.08$	$32.03$	$32.03$	$39.28$	$39.28$
$-7/2$		-0.695	-0.196			0.689	-0.062
$-5/2$	-0.318			-0.249	-0.169			0.899
$-3/2$		-0.546	-0.428			-0.691	-0.205
$-1/2$	-0.343			-0.846	-0.139			-0.383
$1/2$		0.343	-0.846			0.139	0.383
$3/2$	0.546			-0.428	0.691			0.205
$5/2$		0.318	-0.249			0.169	-0.899
$7/2$	0.695			-0.196	-0.689			0.062

Equations

$$H=\sum _{n=2}^6\sum _{m=0}^{\le n}{B}_n^m{\hat{O}}_n^m+\sum _{n=4}^6\sum _{m=2}^{\le n}B(i{)}_n^m\hat{O}(i{)}_n^m$$

$$S(|Q|,\hslash \omega )=\sum _{i,i^{\prime }}\frac{({\sum }_{\alpha }|⟨i{|J_{\alpha }|i^{\prime }⟩|}^2){\mathrm{e}}^{-\beta E_i}}{{\sum }_j{\mathrm{e}}^{-\beta E_j}}L(\Delta E+\hslash \omega ,{\Gamma }_{i,i^{\prime }})$$

$$0.695|\pm \frac{7}{2}⟩-0.318|\mp \frac{5}{2}⟩+0.546|\pm \frac{3}{2}⟩-0.343|\mp \frac{1}{2}⟩.$$

$$\displaystyle\left\{\begin{array}[]{ll}\vec{a}_{mon}=\vec{a}_{hex}\\ \vec{b}_{mon}=\vec{a}_{hex}+2\vec{b}_{hex}\\ \end{array}\right.$$

$$\displaystyle\left\{\begin{array}[]{ll}\vec{a}_{hex}=\vec{a}_{mon}\\ \vec{b}_{hex}=(\vec{b}_{mon}-\vec{a}_{hex})/2\\ \end{array}\right.$$

$${\left(\begin{array}{c}h\\ k\end{array}\right)}_{mon}=\left(\begin{array}{cc}1& 0\\ 1& 2\end{array}\right){\left(\begin{array}{c}h\\ k\end{array}\right)}_{hex}$$

$${\left(\begin{array}{c}h\\ k\end{array}\right)}_{hex}=\left(\begin{array}{cc}1& 0\\ -1/2& 1/2\end{array}\right){\left(\begin{array}{c}h\\ k\end{array}\right)}_{mon}$$

$$|{\psi }_i⟩=\sum _m{a}_{im}|J,\mu +mp⟩$$

$${\rho }_e(\vec{r})=\sum _{k=0}^{min(2l,2J)}\sum _{q=-k}^{+k}{\rho }_{kq}(\vec{r})Y_k^q(\Omega )$$

$${\rho }_e(\vec{r})=R_{nl}^2(\vec{r})\sum _{k=0}^{min(2l,2J)}\sum _{q=-k}^{+k}{c}_{kq}{Y}_k^q(\Omega )$$

$${\rho }_e(\vec{r})=R_{nl}^2(\vec{r})\sum _{k=0}^{min(2l,2J)}\sum _{q=0}^{+k}{\zeta }_{kq}{Z}_k^q(\Omega ).$$

$${\zeta }_{kq}=\sqrt{\frac{2k+1}{4\pi }}{\theta }_k{b}_{kq}⟨{\psi }^{†}|O_k^q|\psi ⟩$$

$$⟨{\psi }^{†}|O_k^q|\psi ⟩=0$$

$$⟨{\psi }^{†}|O_k^q|\psi ⟩=\sum _m{a}_m^2⟨J,\mu +mp|O_k|J,\mu +mp⟩$$

$${\rho }_e(\vec{r})=R_{nl}^2(\vec{r})\sum _{k=0}^{min(2l,2J)}\sqrt{\frac{2K+1}{4\pi }}{\theta }_k$$

$$\sum _{q=mp}^k{b}_{kq}⟨{\psi }^{†}|O_k^q|\psi ⟩Z_{kq}(\Omega )$$

Full text

The Quantum Spin Liquid (QSL) is a state of matter hosting exotic fractionalized excitations and long range entanglement between spins with potential applications for quantum computing Savary and Balents (2016); Zhou et al. (2017); Takagi et al. (2019); Knolle and Moessner (2019). Since QSL physics relies on quantum fluctuations that are enhanced by low spin and low dimensionality, spin-1/2 systems on two-dimensional lattices provide a natural experimental platform for realizing a QSL phase. It has also been shown that an effective spin-1/2 system can be generated even in compounds with high-angular-momentum ions like Yb3+ and Er3+, where the combination of crystal-field effects and strong spin-orbit coupling lead to highly anisotropic interactions between effective spin-1/2 degrees of freedom Ross et al. (2014).

Magnetic frustration plays a central role in stabilizing QSL phases Balents (2010). While QSLs were traditionally associated with geometrically frustrated systems (e.g., triangular and kagome lattices), it has recently become well appreciated that exchange frustration due to highly anisotropic spin interactions can also stabilize QSL phases, even on bipartite lattices Witczak-Krempa et al. (2014); Rau et al. (2016a). Most famously, bond-dependent spin interactions on the honeycomb lattice give rise to the Kitaev model, an exactly solvable model with a gapless QSL ground state Kitaev (2006). A number of honeycomb materials, primarily containing 4d or 5d transition metals such as Ir or Ru have been put forth as realizations of the Kitaev model Jackeli and Khaliullin (2009); Chaloupka et al. (2010). Prominent examples include (Na,Li)2IrO3 Singh and Gegenwart (2010); Liu et al. (2011); Singh et al. (2012); Choi et al. (2012); Ye et al. (2012); Comin et al. (2012); Hwan Chun et al. (2015); Williams et al. (2016) and H3LiIr2O6 Kitagawa et al. (2018), as well as $\alpha$-RuCl3 Plumb et al. (2014); Sandilands et al. (2015); Sears et al. (2015); Majumder et al. (2015); Johnson et al. (2015); Sandilands et al. (2016); Banerjee et al. (2016); Sears et al. (2017); Banerjee et al. (2017); Baek et al. (2017); Do et al. (2017); Banerjee et al. (2018); Hentrich et al. (2018); Kasahara et al. (2018).

Recently, YbCl3 has been proposed to be a candidate material for Kitaev physics on a honeycomb lattice Xing et al. (2019); Luo and Chen (2019). YbCl3 crystallizes in the monoclinic space group $C12/m1$ ($\#$12). The crystal structure is composed of layers of Yb3+ ions coordinated by slightly distorted Cl octahedra as illustrated in Fig. 1. Despite being formally monoclinic at 10 K, the Yb-Yb distances of 3.864 Å and 3.886 Å and the Cl-Yb-Cl bond angles of 96.12*∘* and 96.73*∘* are nearly identicalsup (2019). The result of this atomic arrangement are well-separated, nearly-perfect honeycomb layers of Yb3+ ions in the $ab$-plane as shown in Fig. 1(a,c). The environment surrounding the Yb3+ cations depicted in Fig. 1(b) consists of $6$ Cl- anions arranged in distorted octahedra where the $b$-axis is the unique $C_{2}$ axis. Xing, et al. Xing et al. (2019) have reported that YbCl3 undergoes short range magnetic ordering at 1.2 K. A small peak in the heat capacity at $0.6$ K may indicate a transition to long range magnetic order. Moreover, the field dependence of the inferred ordering temperature suggests that the interactions in YbCl3 are 2-dimensional. On the other hand, Yb-based quantum magnets have been the subject of several recent investigations and, surprisingly, in many cases these materials have been found to possess strong effective Heisenberg exchange interactions Rau et al. (2016b); Haku et al. (2016); Park et al. (2016); Rau et al. (2018); Wu et al. (2016, 2019); Rau and Gingras (2018). Thus, key open questions for YbCl3 are the nature of the spin Hamiltonian and the role of potential Kitaev terms. It is likewise important to determine the single-ion ground state out of which the collective physics grows and additionally if the ground state doublet is well isolated and can be considered to be in the effective quantum spin-1/2 limit. In this paper we address these issues using inelastic neutron scattering and thermodynamic measurements to study the crystal field and low energy excitation spectrum in polycrystalline samples of YbCl3.

Anhydrous beads of YbCl3 and LuCl3 were purchased from Alfa Aesar and utilized in the experimental work presented here. Additional information and results of sample characterization are provided in the SIsup (2019). Refinements of neutron powder diffraction data did not reveal any significant chlorine deficiency or secondary phasessup (2019).

The crystal field excitations were measured with inelastic neutron scattering performed with the SEQUOIA spectrometer at the Spallation Neutron Source at Oak Ridge National Laboratory (ORNL) Granroth et al. (2010). Approximately 4.2 g of polycrystalline YbCl3 and 2.5 g of its non-magnetic equivalent LuCl3 were loaded into cylindrical Al cans and sealed under helium exchange gas. The use of the LuCl3 measurement as a background subtraction is described in the SIsup (2019). The samples and an equivalent empty can for Al background subtraction Kresch et al. (2008) were measured at $T=5$ K, $95$ K and $185$ K, with incident energies, $E_{i}=6$ meV, $45$ meV and $60$ meV with the high resolution chopper. Unless otherwise noted, all inelastic data presented here have had the measured backgrounds subtracted with the data reduced using the software packages Dave Azuah RT (2009) and MANTID Arnold et al. (2014).

Figures 2(a) and (b) show the inelastic neutron scattering spectra as a function of wave-vector transfer, $Q$, and energy transfer, $\hbar\omega$, measured at $T=5$ K and $95$ K respectively. Figure 2(c) is the wave-vector integrated scattering intensity from the $E_{i}=60$ meV measurements for $Q<3$ Å*-1*. The prominent higher energy modes are identified as crystal field excitations both from their $Q$-dependence and from comparison with the nonmagnetic analog LuCl3sup (2019). At $T=5$ K they are centered at energy transfers of $\hbar\omega=20.9$, $31.7$, and $39.5$ meV. The 20.9 meV mode is noticeably broadened toward higher energy transfers. Increasing temperature reduces intensity but does not appreciably shift or broaden these transitions, consistent with the behavior expected for crystal field excitations. Note, that there are some lower energy acoustic phonon modes in the data that are not well subtracted, particularly near 4 meV.

To understand the nature of the crystal field spectrum we analyze the energy levels following a formalism described by Wybourne Wybourne (1965); Judd (1956, 1959) and Stevens Stevens (1952). Given the $C_{2}$ site symmetry of the local Yb environment, the crystal field Hamiltonian consists of $14$ parameters Walter (1984). Prather’s convention Prather (1961) for the minimal number of crystal field parameters was achieved by rotating the environment by $\pi/2$ around the a-axis ($\hat{x}$-axis), i.e. the axis of quantization becomes the b-axis (the $\hat{z}$-axis in the rotated coordinate system). To constrain the parameters, we simultaneously fit the neutron scattering data at 5 and 95 K (Fig. 2(c)), the static magnetic susceptibility between 4 and 700 K (Fig. 3) and the field-dependent magnetization at 10 K (inset of Fig. 3).

Hund’s rules state that, for a $4f^{13}$ ion, $L=3$ and $S=1/2$, thus $J=\|L+S\|=7/2~{}$ Kittel (2005). Therefore the crystal field Hamiltonian can be written in terms of Steven’s operators as

[TABLE]

for $n$ even, where $B_{n}^{m}$ are the crystal field parameters, and $\hat{O}_{n}^{m}$ are the Steven’s operators Hutchings (1964) both in real and imaginary ($i$) form. Once Eq. 1 is diagonalized, the scattering function, $S(|Q|,\hbar\omega)$, can be written as

[TABLE]

where $\beta=1/k_{B}T$, $\alpha=x,y,z$, $\Delta E=E_{i}-E_{i^{\prime}}$, and $L(\Delta E+\hbar\omega,\Gamma_{i,i^{\prime}})$ is a Lorentzian function222As explained in greater detail below, a constrained two component Lorentzian has been used for the crystal field level at 21 meV with halfwidth $\Gamma_{i,i^{\prime}}$ that parameterizes the lineshape of the transitions between CEF levels (eigenfunctions of Eq. (1)) $i\rightarrow i^{\prime}$. We calculate the scattering function using this formalism, compare these values with the experimental data, and then vary the crystal field parameters to minimize the $\chi^{2}$ difference between the model and the data shown in Figs. 2(c) and Fig. 3.

The refinement of the Hamiltonian (Eq. 1) in the scattering function described in Eq. 2 yields the crystal field parameters presented in Tab. 1 and the set of eigenfunctions written in Tab. II of the SIsup (2019). The ground state eigenfunction is found to be

[TABLE]

The imaginary part of the eigenfunction is not shown because it was determined to be $\approx 2$ orders of magnitude smaller than the real part. The calculated $S(|Q|,\hbar\omega)$ is plotted at both temperatures and is shown in Fig. 2(c) as solid lines. The integrated intensity of the three crystal field excitations is reproduced as is the magnetic susceptibility (Fig. 3) and the field dependent magnetization at 10 K (inset of Fig. 3).

The CEF model demonstrates that the Yb3+ ions have a planar anisotropy and a calculated magnetic moment of $2.24(5)\mu_{B}$/ion is obtained for this ground state. The calculated components of the $g$-tensor for the ground state doublet, using the convention described above, for YbCl3 are $g_{z}=4.2(2)$, $g_{x}=3.8(2)$, and $g_{y}=2.2(2)$, which shows somewhat more anisotropy than Ref. [Xing et al., 2019]. Additionally, using the crystal field model derived here as a starting point, we calculated a magnetic torque diagram at 2.1 K for an applied field of 5 T (Fig. 3 inset). The result reproduces the data in Ref. [Xing et al., 2019] (note the difference in coordinate conventions), demonstrating that the crystal field ground state is anisotropic independent of any additional exchange anisotropy.

Despite the overall quality of the fits, one aspect of the crystal field excitation spectrum remains puzzling. The lineshape of the crystal field excitation centered at 21 meV extends toward higher energies. A similar broadening is not observed for the other crystal field excitations at 32 and 40 meV. Thus the broadening is a characteristic of the level at 21 meV and not of the ground state. To fully account for the spectral weight, we have modeled the lineshape for this excitation as two constrained Lorentzians with the widths fixed to be the same and the positions offset by a fixed amount. The lack of observable impurity peaks in the neutron diffraction datasup (2019) suggests that this effect is not due to an impurity phase. Deviations from ideal Cl stoichiometry are similarly hard to detect. Another possibility is that stacking faults result in a variation of the crystal field potential along the c-axis, though this was not evident in the diffraction data. The level at 21 meV would be more strongly affected by such stacking faults given the strong charge density out of the plane for this eigenfunction (see SIsup (2019) Fig. S4 for plots of the charge density for each eigenfunction). Additionally, first principles calculations of the phonon density of states suggests that this feature is not the result of hybridization of the crystal field level with a phononsup (2019). Careful studies of single crystals are required to further understand the origin of this broadening. Finally, we note that using a single Lorentzian in the CEF modeling does not significantly change the refined CEF parameters as the additional spectral weight is relatively small.

To probe for low-energy magnetic correlations, we performed inelastic neutron scattering measurements using the HYSPEC instrument Winn et al. (2015) at ORNL. For these measurements, the same sample used in the SEQUOIA measurements was cooled to $T=1.6$, $5$, and $10$ K and measured with $E_{i}=3.8$ meV at two different positions of the detector bank to cover a large range of $Q$. A $T=100$ K measurement of the YbCl3 sample was used as the background to isolate the magnetic scattering. Figure 4(a) shows the energy and wave-vector dependent magnetic spectrum in YbCl3. A broad dispersive mode with additional scattering is evident. The additional scattering may be due to a quantum continuum however, other explanations such as broadened excitations from a short ranged ordered state, magnon decay, etc, cannot be excluded with the data at hand. The $Q$ integrated scattering intensity in Fig. 4(c) shows a single peak at approximately 0.5 meV with no indication of a spin gap or additional scattering intensity above approximately 1.3 meV energy transfer. Given that long range magnetic order occurs at a maximum temperature of 0.6 KXing et al. (2019), the energy scale of the spin excitations suggests low dimensional and/or frustrated spin interactions in YbCl3. The $\hbar\omega$ integrated intensity in Fig. 4(b) is a broad function which peaks at approximately $Q=1.1$ Å*-1* likely corresponding to the reciprocal lattice points (1 1 0) and (0 2 0), which is consistent with spin correlations within the basal plane. The data in Fig. 4(c) were collected at $T$=1.6 K, which is at lower $T$ than the maximum in the specific heat capacityXing et al. (2019); sup (2019). Thus, the low temperature spin excitations may be responsible for a portion of the loss of entropy despite the lack of apparent long range order. Additional measurements using single crystals are required to fully understand the nature of the magnetic ground state and the spin excitation spectrum.

To investigate the spin-spin correlations in YbCl3, we performed Reverse Monte Carlo (RMC) calculations as implemented in Spinvert Paddison et al. (2013) (see SIsup (2019) for more details). Within this approximation, we fit the integrated intensity of the low energy excitation spectrum as a function of $Q$. Ising, XY, and Heisenberg types of spin correlations were all tried as initial starting points for the simulations (See the SIsup (2019) for further detail). The result of the RMC modeling is shown as a solid line in Fig. 4(b). The radial spin-spin correlation function was calculated for each final spin configuration as a means to investigate the orientation of the spins with respect to each other. Assuming a purely hexagonal geometry, nearest neighbor spins are antiferromagnetically correlated, second neighbor spins have weak ferromagnetic correlations, while there is a rapid decay of spin correlations at larger distances. This result is independent of the type of starting correlation used for the modeling.

We analyzed the spectroscopic properties of the interesting quantum magnet YbCl3. Our studies show that YbCl3 has crystal field excitations at $\hbar\omega=20.9$, $31.7$, and $39.5$ meV. The ground state is a well separated effective spin-1/2 doublet with easy plane ansisotropy and an average magnetic moment of $2.24(5)\mu_{B}$/ion. At $T$=1.6 K, where long range order is not believed to exist, the low energy dynamics of the YbCl3 are consistent with a low dimensional interacting spin system with antiferromagnetic nearest neighbor correlations.

Supplemental Information

I Sample Preparation

Anhydrous beads of YbCl3 (Alfa Aesar stock 40653) were utilized and care was taken to avoid exposure to air. The material was received in an argon-filled glass ampoule, which was opened inside a helium-filled glovebox. The material was stored in the glovebox prior to loading aluminum canisters for neutron scattering measurements. An indium seal was tightened inside the helium glovebox and checked with a helium leak detector. A powder sample of anhydrous LuCl3 (Alfa Aesar stock 35802) was handled in the same way so that a non-magnetic analogue could be examined.

II Heat Capacity

Specific heat measurements down to 0.4 K were performed in a Quantum Design Physical Property Measurement System equipped with a 3He insert. The low-temperature specific heat $C_{P}(T)$ (left axis) and $C_{P}(T)/T$ (right axis) of YbCl3 are shown in Fig. S1. Our $C_{P}(T)$ measurement shows a broad feature centered around 1.9 K, which corresponds to a maximum near 1.2 K in $C_{P}(T)/T$. Upon close inspection, a very weak feature can also be inferred near 0.6 K. These results are consistent with those reported by Xing et. al.[Xing et al., 2019]. As demonstrated by Xing et. al., the majority of the magnetic entropy is lost above 0.6 K.

III POWGEN refinement results

Neutron powder diffraction measurements were performed using the POWGEN time-of-flight diffractometer at the Spallation Neutron Source. Huq et al. (2011) Powder samples were measured in vanadium sample cans at $T=100$ K and $T=10$ K. The data were corrected for neutron absorption. The structural model was refined using the GSAS-IIToby and Von Dreele (2013) package. Data and refinement at 10 K are shown in Fig. S2. The refinement confirmed that our sample crystallizes into the monoclinic space group with $a=6.7291(3)$ Å, $b=11.6141(1)$ Å, $c=6.3120(3)$ Å and $\beta=110.5997(7)^{\circ}$. Alternately, as shown in the next section, neglecting the small monoclinic distortion, the in plance lattice parameters can be written in hexagonal notation as $a=b=6.720$ Å with $\gamma=120^{\circ}$. Table S2 shows a summary of our refinement parameters at 100 K and 10 K.

IV Transformation equations from Monoclinic to Hexagonal Notation

When the structure is weakly monoclinic it can be useful to transform the crystallographic parameters from monoclinic symmetry to hexagonal symmetry. Here we consider the transformation of a single layer from monoclinic to hexagonal notation. If we keep the $\vec{a}$ constant for the two systems, then the relations connecting the two symmetries are:

[TABLE]

or equivalently:

[TABLE]

where the subscripts “$mon$” and “$hex$” denote monoclinic and hexagonal notation respectively. From these vectors the reciprocal lattice vectors can be calculated using the standard formulasKittel (2005).

Therefore, in analogy, a monoclinic reflection $q=(h,k)$ can be transformed into its hexagonal equivalent as follows:

[TABLE]

[TABLE]

V Sequoia crystal field background analysis

In order to study the crystalline electric field (CEF) of YbCl3 we first investigated the background of our compound in the proximity of the three crystal field excitations. In principle, a good background subtraction to eliminate the phonon contamination can be done by measuring a non-magnetic equivalent of the main sample under the same experimental conditions. LuCl3 was used for this purpose.

We present in Fig. S3(a-d) a comparison of the data set collected with SEQUOIA for YbCl3 and LuCl3 with Ei = 60 meV at 5 K and 95 K. An empty can data set at the same temperature has been subtracted to eliminate the contribution due to the sample holder. The last column of Fig. S3(e-f) shows cuts in the momentum transfer range Q=[0,3] Å*-1* at both temperatures in analogy with our analysis in the main text, to compare the phonon density of states and the crystal field spectrum.

The energy transfer range below 15 meV is dominated by the phonons of Yb(Lu), in particular the acoustic phonon at Q=4 Å*-1* and the optical phonon at 9.8 meV have been used to normalize the intensities of the two data sets. Cuts were taken through the acoustic phonon away from the nuclear Bragg peak for both samples, and the integrated intensities were compared to find the proper normalizing factor. Then this quantity was cross-checked by taking cuts along the optic phonon at $Q\geq 4$ Å*-1*, where the magnetic form factor of Yb3+ is small. With the data sets normalized in this manner, the LuCl3 data set was then subtracted from the YbCl3 data set. The resulting data sets were used in the refinements of the crystal field model.

The wave functions of the refined crystal field model of YbCl3 are presented in Tab. S3. The complex parts are omitted because they are approximately two orders of magnitude smaller than the real ones. One of the possible ways to visualize the eigenfunctions is to calculate and plot the electronic charge density $\rho_{e}(\vec{r})=\langle\psi^{\dagger}|\psi\rangle$ of each of the crystal field levels. In general, under the influence of an external electric field, the crystal field states $|\psi_{i}\rangle$ can be written as:

[TABLE]

where $\mu$ denotes the CEF quantum number, $m$ is an integer and $J$ is the total angular momentum value. Here we choose, the axis quantization, $\hat{z}$, to be parallel to the $b$-axis, then following the procedure described in Ref. Walter (1984), $\rho_{e}$ can be expressed in terms of spherical harmonics as:

[TABLE]

where $l=3$ for rare earth ions, $Y_{k}^{q}(\omega)$ is a spherical harmonic, and $k$ is limited to even numbers due to the time reversal invariance of the charge density.

Now, an eigenfunction can always be partitioned into a radial part $R_{nl}(\vec{r})$ which depends only on the quantum number $n$ and $l$, and an angular part that is a linear combination of spherical harmonics:

[TABLE]

or alternatively using the tesseral harmonics $Z_{kq}$ as Hutchings (1964):

[TABLE]

Note that if a 2-fold axis perpendicular to the $p$-fold axis exists, then all the terms in the expansion are real.

The coefficients $\zeta_{kq}$ in Eq. S15 can be calculated using Steven’s operators as:

[TABLE]

where $\theta_{k}$ and $b_{kq}$ values can be found in Refs. Hutchings (1964); Walter (1986) respectively, and $O_{k}^{q}$ are the Steven’s operators. It can also be verified that:

[TABLE]

for $q$ not multiple of $p$ and:

[TABLE]

for $q$ multiple of $p$ and with m integer.

Finally, we can express the charge density using Eqs. S15-S18 as:

[TABLE]

for even $k$. The calculated angular part of the charge density for the 4 crystal levels is presented in Fig. S4. Note that we assumed $R_{nl}(\vec{r})=1$ for ease of visualization.

V.1 DFT Phonon Calculation on LuCl3

To confirm our understanding of the phonon density of states (DOS) and its relationship to the crystal field excitation spectrum, the phonon DOS of LuCl3 was calculated using Density Functional Theory(DFT). The calculations were performed using the Vienna Ab initio Simulation Package (VASP) Kresse and Furthmüller (1996). The calculation used Projector Augmented Wave (PAW) method Blöchl (1994); Kresse and Joubert (1999) to describe the effects of core electrons, and Perdew-Burke-Ernzerhof (PBE) Perdew et al. (1996) implementation of the Generalized Gradient Approximation (GGA) for the exchange-correlation functional. The energy cutoff was 600 eV for the plane-wave basis of the valence electrons. The lattice parameters and atomic coordinates determined by neutron diffraction at 10 K (see Tab. S2) were used as the initial structure. The electronic structure was calculated on a $6\times 3\times 6$ $\Gamma$-centered mesh for the unit cell, and a $3\times 3\times 3$ $\Gamma$-centered mesh for the $2\times 1\times 2$ supercell. The total energy tolerance for electronic energy minimization was $10^{-8}$ eV, and for structure optimization it was $10^{-7}$ eV. The maximum inter-atomic force after relaxation was below $0.001$ eV/Å, and the pressure after relaxation was below 2 MPa. The relaxed lattice parameters are $a=6.620$ Å, $b=11.437$ Å, $c=6.319$ Å, $\alpha=90^{\circ}$, $\beta=110.395^{\circ}$, $\gamma=90^{\circ}$. A Hubbard U term of 6.0 eV was applied on Lu for the localized $4f$ electrons Jiang et al. (2012). The optB86b-vdW functional Klimeš et al. (2009) for dispersion corrections was used to describe the van der Waals interactions between layers. Force constants were calculated using Density Functional Perturbation Theory (DFPT) as implemented in VASP, and the vibrational eigen-frequencies and modes were then calculated by solving the dynamical matrix using $Phonopy$ A Togo (2015). Non-analytical correction was applied to account for the LO-TO splitting Gonze and Lee (1997). The $OClimax$ software Cheng et al. (2019) was used to convert the DFT-calculated phonon results to the simulated INS spectra.

The calculated $S(\vec{Q},\omega)$, the relative DOS and the phonon dispersion across the highly symmetric directions of the Brillouin zone for LuCl3, are shown in Fig. S5. The DFT calculations are in very good agreement with the data displayed in Fig. S3(c) for LuCl3; in particular we can see that the first and third excited states (represented with dashed lines) are located in a range of energies where phonons are totally absent effectively ruling out a possible hybridization.

VI Reverse Monte Carlo Analysis

In order to extract information about the spin correlations in YbCl3 from the powder data set we collected at HYSPEC, we performed a Reverse Monte Carlo (RMC) calculation following the approach described in Ref. [Paddison et al., 2013]. The program $SpinVert$ was used to fit the temperature subtracted data set (as explained in the main text), integrated in the energy range $\hbar\omega$=[0.1,1.2] meV to avoid contamination from the elastic line.

A simulation super-cell of $6\times 6\times 6$ containing a total of $6^{4}=1296$ spins with periodic boundary conditions (PBCs) has been used to fit the spin excitation spectrum. For simplicity, the first attempt assumed all spins in the super-cell have Ising-like nature and are completely uncorrelated. The RMC has been done with a statistical average of 100000 configurations. The only fitting parameter in the code was a scaling factor to match the intensity of our calculation with the data set. The final spin configurations generated by the RMC were then averaged and the best fit with the HYSPEC data gave $R_{wp}=7.1$, and a $\chi^{2}=11.1$.

In order to verify our model we performed a campaign of RMC simulations changing the anisotropy of the spins to XY-like and/or Heisenberg, increasing the size of the simulation super cell to study possible size effects, and repeating the same type of calculations in monoclinic coordinates to check differences in the diffuse scattering pattern. We found that increasing the size of the super-cell above $6\times 6\times 6$ does not affect the quality of the fit, nor the calculated diffuse scattering; the calculation simply takes longer due to the higher number of spins in the super-cell. In the same way, performing the calculation in hexagonal or monoclinic notation does not affect the final result.

The lack of a spin Hamiltonian that describes the exchange interactions prevents us from fully understanding the dynamics in YbCl3 at this time. However, the analysis of the angular distribution of the spins in the Heisenberg model (see Fig. S6) indicates that the spins are preferentially laying in the $ab$ plane as suggested by the local anisotropy of the crystal field ground state. Finally, from the generated spin configurations, we can calculate and compare the radial spin-spin correlation function Paddison et al. (2013) for all three models (i.e., the Ising-like, XY-like, and Heisenberg models), to analyze the spin-spin correlations. Figure S6(b) shows that there is an antiferromagnetic nearest neighbor spin correlation. The second and third neighbors are ferromagnetically correlated. Furthermore, this behavior is qualitatively the same for all three models explored here.

The RMC calculations also provide a prediction of the diffuse scattering pattern which can be used to further constrain the type of spin-spin correlation when compared to single crystal neutron diffraction experiments. The calculated diffuse scattering patterns are shown in Fig. S7(b-c).