A neutron scattering measurement of crystalline-electric fields in magnesium rare-earth selenide spinels

TL;DR

This study uses inelastic neutron scattering to determine the crystal electric field excitations in magnesium rare-earth selenide spinels, revealing their ground state properties and magnetic behaviors relevant to frustrated magnetism.

Contribution

The paper provides the first detailed CEF analysis for MgRE$_2$Se$_4$ spinels, experimentally confirming ground state characteristics and magnetic anisotropies for multiple rare-earth variants.

Findings

MgTm$_2$Se$_4$ has a non-magnetic ground state

MgYb$_2$Se$_4$ exhibits effective $S=1/2$ spins with specific g-factors

MgHo$_2$Se$_4$ has a ground state doublet with Ising spins and accessible low-lying CEF levels

Abstract

The symmetry of local moments plays a defining role in the nature of exotic grounds states stabilized in frustrated magnetic materials. We present inelastic neutron scattering (INS) measurements of the crystal electric field (CEF) excitations in the family of compounds MgRE$_2$Se$_4$ (RE $\in$ $\{$Ho, Tm, Er and Yb$\}$). These compounds form in the spinel structure, with the rare earth ions comprising a highly frustrated pyrochlore sublattice. Within the symmetry constraints of this lattice, we fit both the energies and intensities of observed modes in the INS spectra to determine the most likely CEF Hamiltonian for each material and comment on the ground state wavefunctions in the local electron picture. In this way, we experimentally confirm MgTm$_2$Se$_4$ has a non-magnetic ground state, and MgYb$_2$Se$_4$ has effective $S=\frac{1}{2}$ spins with $g_\parallel = 5.188(79)$ and $g_\perp = 0.923(85)~\mu_B$. The spectrum of MgHo$_2$Se$_4$ indicates a ground state doublet containing Ising spins with $g_\parallel = 2.72(46)$, though low-lying CEF levels are also seen at thermally accessible energies $\delta_E = 0.591(36)$, 0.945(30) and 2.88(7)~meV, which can complicate interpretation. These results are used to comment on measured magnetization data of all compounds, and are compared to published results on the material MgEr$_2$Se$_4$.

Figures

Tables

Table 1. Table 1: Table listing the fit cubic lattice parameter ( a 𝑎 a ), and the partial coordinate position of Se ( x Se subscript 𝑥 Se x_{\text{Se}} ) extracted from Rietveld refinements of XRD data.

Material	$a$ (Å)	$x_{\text{Se}}$
MgHo2Se4	11.5508(2)	0.2466(1)
MgEr2Se4Reig-i Plessis et al. (2019)	11.5207(14)	0.2456(9)
MgTm2Se4	11.48493(5)	0.24614(7)
MgYb2Se4	11.45591(3)	0.24595(8)

Table 2. Table 2: A list of all observed phases and their percent masses in the samples used in this paper. Next to each compound name is an example of the mark showing their peak positions in Fig.  2 .

MgHo2Se4	( )	71.15( 0.68)	MgTm2Se4	( )	90.31( 0.40)	MgYb2Se4	( )	68.11(0.86)	MgYb2Se4	( )	92.01( 0.52)
Ho2Se3	( )	8.84( 0.22)	Tm2O2Se	( )	5.57( 0.06)	Yb2O2Se	( )	3.92(0.28)	Yb2O2Se	( )	3.02( 0.06)
Ho2O2Se	( )	7.62( 0.14)	TmSe	( )	4.12( 0.08)	YbSe	( )	14.88(0.33)	YbSe	( )	4.97( 0.13)
HoSe	( )	12.38( 0.24)				Yb7.24Se8	( )	13.09(0.21)

Table 3. Table 3: A list of the expected energies and neutron scattering cross sections of transitions from the the i t ​ h superscript 𝑖 𝑡 ℎ i^{th} to j t ​ h superscript 𝑗 𝑡 ℎ j^{th} CEF levels in MgHo 2 Se 4 . Listed are both expectations from fitted B m ​ n subscript 𝐵 𝑚 𝑛 B_{mn} presented in Table  6 , and from rescaling the CEF potential of MgEr 2 Se 4 .

MgHo2Se4
$i=0$  ( )
best fit				rescaled
$j$	$E\left(meV\right)$	$I_{5K}$	$I_{100K}$	$E\left(meV\right)$	$I_{5K}$	$I_{100K}$
1	0.59	36.29	14.47	0.76	51.61	20.78
2	0.95	52.00	20.74	0.87	35.71	14.38
3	2.88	0.30	0.12	2.69	0.81	0.32
4	17.74	0.64	0.26	20.02	0.61	0.25
5	19.20	0.69	0.28	20.29	0.54	0.22
6	20.71	4.48	1.79	21.92	4.59	1.85
7	22.56	2.76	1.10	23.54	3.15	1.27
8	24.27	0.05	0.02	26.07	0.00	0.00
9	26.08	0.74	0.29	27.63	0.34	0.14
10	34.87	0.09	0.03	31.54	0.13	0.05
$i=1$   ( )
2	0.35	0.01	0.02	0.11	0.04	0.09
3	2.29	1.93	2.84	1.93	5.40	11.70
4	17.15	0.00	0.00	19.26	0.93	2.02
5	18.61	1.14	1.67	19.53	0.32	0.70
6	20.12	0.07	0.10	21.16	0.01	0.03
7	21.97	0.07	0.10	22.78	0.48	1.04
8	23.68	0.05	0.08	25.31	0.03	0.06
9	25.49	0.15	0.22	26.87	0.13	0.29
10	34.28	0.06	0.08	30.78	0.03	0.07
$i=2$    ( )
3	1.94	3.44	11.02	1.82	0.98	2.70
4	16.80	0.74	2.38	19.15	0.00	0.00
5	18.26	0.22	0.69	19.41	0.71	1.97
6	19.77	0.01	0.04	21.05	0.07	0.21
7	21.61	0.21	0.66	22.66	0.04	0.12
8	23.32	0.03	0.10	25.20	0.01	0.03
9	25.13	0.18	0.57	26.75	0.03	0.10
10	33.93	0.02	0.06	30.67	0.01	0.04

Table 4. Table 4: A list of the expected energies and neutron scattering cross sections of transitions from the the i t ​ h superscript 𝑖 𝑡 ℎ i^{th} to j t ​ h superscript 𝑗 𝑡 ℎ j^{th} CEF levels in MgTm 2 Se 4 . Listed are both expectations from fitted B m ​ n subscript 𝐵 𝑚 𝑛 B_{mn} presented in Table  6 , and from rescaling the CEF potential of MgEr 2 Se 4 .

MgTm2Se4
$i=0$  ( )
best fit				rescaled
$j$	$E\left(meV\right)$	$I_{5K}$	$I_{100K}$	$E\left(meV\right)$	$I_{5K}$	$I_{100K}$
1	0.88	28.92	11.12	0.57	25.62	12.27
2	12.26	6.56	2.52	15.49	5.67	2.71
3	12.55	0.00	0.00	16.66	0.00	0.00
4	18.12	0.64	0.25	20.34	0.54	0.26
5	27.92	0.01	0.00	33.31	0.19	0.09
6	37.57	0.97	0.37	42.51	1.09	0.52
7	46.76	0.04	0.02	46.04	0.00	0.00
8	48.41	0.00	0.00	47.18	0.00	0.00
$i=1$    ( )
2	11.38	0.11	0.28	14.92	0.24	0.40
3	11.68	0.33	0.87	16.09	0.44	0.73
4	17.25	0.40	1.07	19.77	0.97	1.63
5	27.04	0.00	0.00	32.74	0.00	0.00
6	36.70	0.21	0.56	41.94	0.18	0.30
7	45.88	0.02	0.07	45.47	0.01	0.01
8	47.54	0.00	0.00	46.62	0.18	0.30
$i=2$  ( )
3	0.29	0.00	2.14	1.17	0.00	1.61
4	5.86	0.00	3.59	4.85	0.00	2.79
5	15.66	0.00	0.00	17.82	0.00	0.01
6	25.31	0.00	0.21	27.02	0.00	0.10
7	34.50	0.00	0.26	30.55	0.00	0.01
8	36.16	0.00	0.00	31.69	0.00	0.21

Table 5. Table 5: A list of the expected neutron scattering cross sections of transitions from the first 3 CEF levels in MgYb 2 Se 4 , given by the best fit of the data, as well as what is expected from rescaling the CEF potential of MgEr 2 Se 4 .

MgYb2Se4
$i=0$  ( )
best fit				rescaled
$j$	$E\left(meV\right)$	$I_{5K}$	$I_{100K}$	$E\left(meV\right)$	$I_{5K}$	$I_{100K}$
1	26.01	4.10	3.78	19.02	5.73	4.98
2	29.13	5.99	5.52	28.28	5.70	4.96
3	54.95	0.13	0.12	54.12	0.07	0.06
$i=1$
2	3.12	0.00	0.04	9.26	0.00	0.17
3	28.93	0.00	0.07	35.10	0.00	0.24
$i=2$
3	25.81	0.00	0.20	25.84	0.00	0.20

Table 6. Table 6: The CEF parameters (in meV) of the compounds MgRE 2 Se 4 for RE = Ho, Tm and Yb from the best fits of the INS data shown in this paper, and results for RE = Er taken from our previous paper Reig-i Plessis  et al. ( 2019 ) .

	MgHo2Se4	MgEr2Se4Reig-i Plessis et al. (2019)	MgTm2Se4	MgYb2Se4
B20	$7.10\left(78\right)\cdot {10}^{-2}$	$-4.214\left(63\right)\cdot {10}^{-2}$	$-0.1707\left(5\right)$	$-0.507\left(61\right)$
B40	$5.52\left(24\right)\cdot {10}^{-4}$	$-6.036\left(30\right)\cdot {10}^{-4}$	$-1.880\left(6\right)\cdot {10}^{-3}$	$2.89\left(4\right)\cdot {10}^{-2}$
B43	$9.07\left(13\right)\cdot {10}^{-3}$	$-1.3565\left(67\right)\cdot {10}^{-2}$	$-4.533\left(18\right)\cdot {10}^{-2}$	$0.310\left(50\right)$
B60	$-2.44\left(18\right)\cdot {10}^{-6}$	$3.264\left(16\right)\cdot {10}^{-6}$	$-4.59\left(11\right)\cdot {10}^{-6}$	$1.98\cdot {10}^{-4}$
B63	$2.94\left(83\right)\cdot {10}^{-5}$	$-3.791\left(75\right)\cdot {10}^{-5}$	$2.48\left(2\right)\cdot {10}^{-4}$	$-2.30\cdot {10}^{-3}$
B66	$-2.09\left(17\right)\cdot {10}^{-5}$	$2.194\left(65\right)\cdot {10}^{-5}$	$-1.441\left(59\right)\cdot {10}^{-4}$	$1.33\cdot {10}^{-3}$

Table 7. Table 7: Calculated ground state wavefunctions for MgRE 2 Se 4 with RE = Ho, Er, Tm, and Yb.

Ho	$|{\psi }_0^{\pm }⟩$	=	$0.323\left(66\right)|\pm 7⟩\mp 0.711\left(11\right)|\pm 4⟩+0.0797\left(45\right)|\pm 1⟩\mp 0.317\left(51\right)|\mp 2⟩-0.516\left(17\right)|\mp 5⟩\mp 0.131\left(42\right)|\mp 8⟩$
ErReig-i Plessis et al. (2019)	$|{\psi }_0^{\pm }⟩$	=	$\pm 0.9165\left(7\right)|\pm 15/2⟩+0.3600\left(11\right)|\pm 9/2⟩\pm 0.1581\left(16\right)|\pm 3/2⟩-0.0731\left(15\right)|\mp 3/2⟩\pm 0.0036\left(7\right)|\mp 9/2⟩+$
			$0.0034\left(14\right)|\mp 15/2⟩$
Tm	$|{\psi }_0⟩$	=	$0.6700\left(5\right)|6⟩+0.1468\left(7\right)|3⟩+0.2431\left(19\right)|0⟩-0.1468\left(7\right)|-3⟩+0.6700\left(5\right)|-6⟩$
Yb	$|{\psi }_0^{\pm }⟩$	=	$-0.9684\left(48\right)|\pm 5/2⟩\pm 0.218\left(16\right)|\mp 1/2⟩+0.1204\left(99\right)|\mp 7/2⟩$

Table 8. Table 8: Calculated values for the g-factors of the ground state of the rare earth ions in MgRE 2 Se 4 , relative to local ⟨ 111 ⟩ delimited-⟨⟩ 111 \langle 111\rangle directions.

	$g_{\parallel }$	$g_{⟂}$
Ho	$2.72\left(46\right)$	$0.000\left(63\right)$
ErReig-i Plessis et al. (2019)	$16.591\left(12\right)$	$0$
Tm	$0$	$0$
Yb	$5.188\left(79\right)$	$0.923\left(85\right)$

Equations

$$H_{CEF}=\sum _{nm}{B}_{nm}{\hat{O}}_{nm},$$

$$\begin{gathered} H_{CEF}=B_{20}{\hat{O}}_{20}+B_{40}{\hat{O}}_{40}+B_{43}{\hat{O}}_{43}+B_{60}{\hat{O}}_{60} \\ +B_{63}{\hat{O}}_{63}+B_{66}{\hat{O}}_{66}, \end{gathered}$$

$$I_{ij}={\Sigma }_{\alpha }{\left|⟨{\psi }_j|{\hat{J}}_{\alpha }|{\psi }_i⟩\right|}^2$$

$$B_{nm}=\frac{{\theta }^{(n)}⟨r^n⟩}{{\theta }_0^{(n)}⟨r^n{⟩}_0}{\left(\frac{a}{a_0}\right)}^{-n+1}{B}_{nm}^0,$$

$$\frac{1}{2}{g}_{ii}{\sigma }_i=\left[\begin{array}{cc}⟨+|j_i|+⟩& ⟨+|j_i|-⟩\\ ⟨-|j_i|+⟩& ⟨-|j_i|-⟩\end{array}\right]$$

Full text

I Introduction

The strategic combination of frustrated lattice geometries and strong local-ion anisotropy is a well-established route for stabilizing novel spin states in quantum materialsMoessner (2001); Greedan (2006); Castelnovo et al. (2012). This fact endows a special significance to the local crystal electric field (CEF) Hamiltonian of magnetic momentsBertin et al. (2012), which dictates size, dimensionality, and allowed interactions in effective spin systems at low temperatures. In $f$-electron systems in particular, spin and orbital degrees of freedom are strongly coupled, leading to dramatic changes in the CEF splittings depending on the number of valence electrons. Additionally the small radius and shielding of the $f$-electron orbitals leads to small CEF splittings, necessitating full consideration of excited states. As a result, even within a family of closely related structures, each change in the number of valence electrons creates an entirely new effective spin system and leads to a wide range of interesting behaviors.

A particularly important example from recent years is the so-called ‘227’ family of pyrochlore oxides, A2RE2O7 (A = cation, RE = rare earth), where rare-earth moments occupy a frustrated sublattice of corner-sharing tetrahedra Gardner et al. (2010). The phase diagram of these materials is exceedingly rich, and includes such diverse states as non-collinear orderChampion et al. (2001), spin glassGreedan et al. (1986), classical spin liquid Gardner et al. (1999, 2001), and both classicalHarris et al. (1997); Ramirez et al. (1999); Bramwell and Gingras (2001) and quantumApplegate et al. (2012); Ross et al. (2011); Pan et al. (2014); Sibille et al. (2016); Wen et al. (2017); Savary and Balents (2017) variants of the “spin ice” states. This variety mirrors the number of different local symmetries selected through the interaction between the valence shell and the CEFBertin et al. (2012), which spans possibilities from strongly Ising-like Harris et al. (1997); Bramwell and Gingras (2001) to XY Hallas et al. (2018a, b) to Heisenberg Champion et al. (2001); Stewart et al. (2008); Wills et al. (2006) moments. Virtual transitions associated with low-lying CEF states have further been credited with inducing quantum fluctuationsMolavian et al. (2007); Rau et al. (2016), while effects of multipolar local ion symmetries are suggested to lead to unexpected spin ordersLhotel et al. (2015), quantum spin ice Huang et al. (2014); Li and Chen (2017), or other enriched spin liquid statesSibille et al. (2015).

This wide variety of exotic states in the single family of isostructural 227 oxides has generated strong interest in other materials in which rare earth moments comprise pyrochlore sublattices, with potentially new CEF environments. Of these, perhaps most prominent have been the rare-earth (RE) spinel chalcogenides: ARE2X4, with X $\in$ {S, Se}. Both spinel and 227 pyrochlore families have global $Fd\bar{3}m$ symmetry and the magnetic cations comprise identical frustrated sublattices. The local environment of the moments are substantially different, however, as demonstrated in Fig 1. The A-site cations in 227 pyrochlores are surrounded by heavily-distorted cubes of O2- atoms, with a large trigonal distortion along the $\langle 111\rangle$ directions. In contrast, the moments in spinels lie at the center of nearly perfect octahedral cages of X2- anions, with trigonal fields arising from both the compression or expansion of the REX6 octahedra and the antiprism of neighboring cation sitesMun et al. (2014); Yaresko (2008); Reig-i Plessis et al. (2019). This substantial difference in local environment allows specific rare earth ions to adopt a drastically different symmetry in the two material families. For example, Er3+ has XY-like moments in the ‘227’ pyrochloresPoole et al. (2007) and Ising-like moments in the spinelsReig-i Plessis et al. (2019).

Among ternary rare-earth chalcogenides, the spinel phases have been confirmed for compounds with A $\in$ {Cd, Mg} and RE $\in$ {Ho, Er, Tm, Yb}Flahaut et al. (1965, 1966). Earliest measurements of material properties were performed in the 1960’s-1980’s, and employed mostly X-ray diffraction (XRD) Flahaut et al. (1965, 1966); Fujii et al. (1972); Ben-Dor and Shilo (1980), magnetizationFujii et al. (1972); Pokrzywnicki et al. ; Pokrzywnicki ; Pokrzywnicki et al. (1977); Ben-Dor and Shilo (1980); Pawlak et al. (1988) and Mössbauer spectroscopyBen-Dor et al. (1981). X-ray diffraction measurementsFlahaut et al. (1965, 1966); Fujii et al. (1972); Ben-Dor and Shilo (1980) confirmed early on the ideal cubic $Fd\bar{3}m$ structure for the entire series, and further indicated that this high symmetry persists to the lowest measured temperatures. This observation stands in contrast to the symmetry-lowering cooperative structural transitions that are typically observed in spinel oxidesLee et al. (2010); Suzuki et al. (2007); Garlea et al. (2008); Yamada et al. (1999).

Original magnetization studies reported no order in the compound CdHo2X4 above 2 K Fujii et al. (1972), and at least one argued for a spin singlet ground state on the basis of an observed temperature independent paramagnetic signal Ben-Dor and Shilo (1980). Early work on CdEr2X4 claimed the onset of magnetic order in the temperature region $T$ = 4 - 10 K based on magnetization Fujii et al. (1972) and Mössbauer spectroscopy Ben-Dor et al. (1981), though these reports stood in conflict with one another and failed to appreciate the consequences of local spin anisotropy on their data. Early reports on CdYb2Se4 provided a more comprehensive analysis, and determined the CEF excitation energies of 20.6 meV and 40.7 meV using a model which accounted for an octahedral environment of Se-anionsPokrzywnicki et al. ; Pokrzywnicki , though important contributions from neighboring cations were ignored. The same study estimated a nearest-neighbor exchange energy to be around $J\approx$ 2.2 K. Studies of CdTm2Se4 concluded having a spin singlet ground state, consistent with expectationsBen-Dor and Shilo (1980); Pokrzywnicki et al. (1977).

In recent years, interest in RE spinel chalcogenides has seen a revival, with a sharper focus on the frustrated nature of interactionsLau et al. (2005) and the resultant potential for novel forms of magnetismLago et al. (2010); Legros et al. (2015); Yaouanc et al. (2015); Higo et al. (2017); Dalmas de Réotier et al. (2017); Gao et al. (2018); Reig-i Plessis et al. (2019). Indeed, both MgEr2Se4 Reig-i Plessis et al. (2019) and CdEr2Se4 Lago et al. (2010); Gao et al. (2018) have independently been identified as strong candidates for a classical spin ice state. Ordered states have been observed in both MgYb2X4Higo et al. (2017) and CdYb2X4Dalmas de Réotier et al. (2017) X $\in$ {S, Se}, but are notable for their highly renormalized moments and the existence of persistent spin dynamics at low temperaturesHigo et al. (2017); Dalmas de Réotier et al. (2017). Similar anomalous fluctuations are reported in CdHo2S4 Yaouanc et al. (2015) below a reported ordering transition, along with several features which draw parallels to the “partially ordered” pyrochlore system Tb2Sn2O7Rule et al. (2007, 2009). The presence of an ordered state and of a local moment size of $8~{}\mu_{B}$ in CdHo2Se4Yaouanc et al. (2015) are in direct contradiction to the singlet ground state predicted from magnetization measurements Ben-Dor and Shilo (1980).

In each of the above cases, the exact nature of the magnetic ground state is closely entwined with the local CEF environment of the constituent RE moments, as has been acknowledged on several occasionsLago et al. (2010); Reig-i Plessis et al. (2019); Gao et al. (2018); Higo et al. (2017); Yaouanc et al. (2015); Rau and Gingras (2018). A necessary condition for realizing spin ice physics in MgEr2Se4Reig-i Plessis et al. (2019) and CdEr2Se4Lago et al. (2010); Gao et al. (2018) is the presence of a ground state Kramers doublet with Ising symmetry. The ordered states in MgYb2X4 and CdYb2X4 have been discussed in the context of frustrated anisotropic exchange models, in which the particular choice of CEF parameters can select from a variety of distinct ordered or spin liquid phasesRau and Gingras (2018). Material specific calculations predict the existence of several low energy CEF states in CdHo2S4Gao et al. (2018), which draw even stronger parallels between this material and Tb2Sn2O7 and lend special significance to the low-temperature fluctuationsMolavian et al. (2007); Mirebeau et al. (2005, 2007).

There thus exists a strong motivation for systematic and high precision measurements of the CEF Hamiltonian across this family of compounds. Some early studies of cadmium spinels acquired this information through careful fits of magnetization data, however the stated results are broadly inconsistent with conclusions from modern studiesPokrzywnicki et al. ; Pokrzywnicki ; Pawlak et al. (1988); Ben-Dor and Shilo (1980); Pokrzywnicki et al. (1977). In the current paper, we instead measure crystal field excitations directly using inelastic neutron scattering (INS), and use fits of both the energy and intensity of observed transitions to determine the most likely CEF Hamiltonian consistent with the symmetry of the $Fd\bar{3}m$ space group. This spectroscopic analysis is the de-facto choice for CEF measurements in rare earth systemsBertin et al. (2012); Siddharthan et al. (1999); Champion et al. (2003); Rosenkranz et al. (2000) due to the method’s precision and symmetry-driven modeling, which is largely insensitive to the presence of impurity phases, defects and other mechanisms which adversely affect bulk thermodynamic data. In a recent publication, we showed how a similar analysis of INS data can be used to confirm the Ising-like effective spins in the material MgEr2Se4, and additionally identified several low-lying CEF excitations capable of inducing quantum fluctuationsReig-i Plessis et al. (2019). Below, we extend this analysis to three other closely related spinel selenide systems. In MgTm2Se4, we confirm the ground state is well characterized as a spin singlet, with the first excited state at $E=0.876(16)$ meV – low enough to thermally excite non-zero local moments at temperatures of only a few Kelvin. In MgHo2Se4, we identify 10 separate CEF excitations, and determine an Ising-like ground state Kramers doublet with multiple low-lying excited states, drawing intriguing parallels to the 227 pyrochlores Tb2Ti2O7 and Tb2Sn2O7. Fits of MgYb2Se4 were underconstrained, but we determine a best fit Hamiltonian which is capable of reproducing both INS and magnetization data at a variety of fields. Compared to previous estimatesHigo et al. (2017); Rau and Gingras (2018), our analysis is notable for the much stronger inferred easy-axis anisotropy. Results for each material are compared to measured magnetization data, and the implications for spin-spin interactions and magnetic ground states are discussed.

II Experimental methods

Polycrystalline samples of RE spinel chalcogenides were prepared via a two-step solid state reaction at Illinois using the same method described in detail in Ref. Reig-i Plessis et al., 2019, and sample quality was confirmed via powder x-ray diffraction (XRD) using a PANalytical X’Pert3 powder diffractometer at the Center of Nanophase Materials Science at Oak Ridge National Laboratory (ORNL). INS measurements were performed using the SEQUOIA Fine-Resolution Fermi Chopper Spectrometer at the ORNL’s Spallation Neutron Source. Spectra were measured with a variety of initial neutron energies, $E_{i}$ and temperatures, $T$, as dictated by the relevant energy scales of the CEF transitions predicted by point charge calculations. Specifically, measurements were taken with $E_{i}=6$, 11, 30, and 50 meV and $T=5$ and 100 K for MgHo2Se4, with $E_{i}=30$, 50, and 100 meV and $T=5$ and 100 K for MgTm2Se4, and with $E_{i}=30$, 50, and 100 meV and $T=5$ and 250 K for MgYb2Se4. Magnetization measurements were taken on a Quantum Design MPMS3 instrument in the Illinois Materials Research Laboratory, utilizing the DC measurement mode. Measurements were performed at temperatures of $T$ = 2, 5, 10, 20, and 40 K for all samples, with additional measurements at 80 and 120 K for MgHo2Se4. Supplementary measurements were performed on MgYb2Se4 as a function of temperature at a constant field of 100 Oe, as laid out below.

III X-ray diffraction

Figure 2 shows the results of powder XRD measurements, along with solid curves representing the best refinements using the FULLPROF software suiteRodríguez-Carvajal (1993). Refinements assumed a majority phase with the $Fd\bar{3}m$ symmetry expected for a normal spinel structure, while accounting for the possibility of several common impurity phases. The positions of associated Bragg peaks are indicated by sets of lines in Figure 2, with the majority phase indicated at the top in black and the impurity phases below in different colors. The brightest reflections were reliably fit to the spinel MgRE2Se4, with the resulting cubic lattice parameters and fractional coordinate of the Se anions listed in Table 1. The values for lattice parameters are $\sim 10\%$ larger than those typically observed in 227 pyrochlore oxidesGardner et al. (2010) and, combined with a larger RE-anion distance, result in a significantly lower energy scale for CEF excitations in the spinel selenides. The fractional coordinate of the Se anions in the $Fd\bar{3}m$ space group quantify the trigonal distortion of RESe6 octahedra, with measured values showing minimal deviation from the undistorted case at $x=0.25$. Accordingly, subsequent point charge calculations demonstrate that the dominant non-cubic contribution to the CEFs at the RE site comes from neighboring cations, and not distortions of the local chalcogen environment. This observation is in conflict with previously used models of CEFs for these compoundsPokrzywnicki et al. ; Pokrzywnicki ; Higo et al. (2017), and provides a further point of contrast between spinels and 227 pyrochlores.

The composition of the prepared samples varied with each synthesis, but all contained the same limited number of impurity phases. The exact distribution of phases in the large volume samples studied with INS are listed in Table 2, along with one higher purity sample of MgYb2Se4 which was prepared for follow-up studies of magnetization. In addition to the majority spinel phase, all samples investigated in this study had sizable fractions of binary rare earth selenide compounds. This is consistent with the high vapor pressure of Mg, resulting in loss during reaction steps. All samples were further seen to contain between 3-8% of the rare earth oxiselenide, which is consistent with the strong tendency for precursor metals to oxidize before forming selenium binaries. The tendency towards metallicity and, in the case of YbSe the lack of local moments, minimize the contribution of the RE monoselenide (RESe) impurities to the INS spectrum Reid et al. (1964). The RE-sesquiselenides (RE2Se3)Prokofiev et al. (1995) and oxiselenides (RE2O2Se)Quezel et al. (1972) are insulators with known structures and were accounted for in subsequent analysis. It is worth noting at this point, however, that the RE-oxiselenides have diffraction patterns which overlap substantially with peaks of the predicted spinel patterns, and have further been reported to have antiferromagnetic ordering transitions at temperatures below 5 KQuezel et al. (1972). The existence of previously unappreciated volume fractions of oxiselenide impurities is thus a leading contender to explain reports of unindexable magnetic Bragg peaks in a number of published neutron powder diffraction studies of RE spinel chalcogenidesReig-i Plessis et al. (2019); Gao et al. (2018); Dalmas de Réotier et al. (2017).

IV Inelastic neutron scattering

Typical INS spectra for MgHo2Se4, MgTm2Se4 and MgYb2Se4 are shown in Fig. 3, which for each material reveal the existence of multiple dispersionless modes at finite energy transfer. The scattering intensity has contributions from both CEF transitions and phonons, in addition to various sources of background. The variation of the scattering intensity as a function of momentum transfer $(Q)$ was used to determine whether observed scattering modes originate from phonons or are magnetic. Each spectrum was measured with multiple incident neutron energies, as a means of separating intrinsic and spurious sources of scattering and to balance energy range and resolution.

Fig. 4 and Fig. 5 show the variation of neutron scattering intensity versus energy, extracted from data in Fig. 3 by integrating over finite regions in $Q$ at positions chosen to maximize the available fit range at each incident neutron energy, $E_{i}$. Data is represented by blue dots, whereas solid curves represent the results of fits described in following sections. The solid red curves are estimates of the slowly varying contributions to background scattering, obtained by performing a cubic spline interpolation between points chosen away from obvious peaks. For all materials, the scattering above the slowly varying background takes the form of well defined peaks which are largely described by the CEF transitions laid out below. Tick marks below the data show the energy positions of thermally accessible transitions between CEF states, described in more detail in the figure caption.

As shown in Fig. 4(a) and (b), the magnetic and phonon excitations are well-separated in the material MgTm2Se4 and therefore the identification of the CEF transitions is most straightforward. The best fit curve shows excellent agreement with the data at both base temperature (5K) and at 100 K, where transitions from thermally populated levels contribute significantly to the scattering pattern. The INS data for MgHo2Se4 presented in Fig. 4(c)-(f) shows multiple overlapping peaks below 30 meV, but they are still clearly above background and mostly captured by the CEF fits. The only exceptions are observed excesses of scattering at energies $E\approx 1.3$ meV and 16 meV, which we respectively associate with an impurity phase discussed below and with a phonon mode also seen in MgYb2Se4.

In MgYb2Se4, the CEF excitations overlap appreciably with dispersionless optical phonon modes, which mildly complicates analysis. The constant-E cuts for MgYb2Se4 shown in Fig. 5(a) and (b) reveal four different modes – a distinct peak near 17 meV, and three closely grouped peaks between 22 and 38 meV. Integrating over a finite energy range in the relevant spectra allows us to compare the $Q$-dependence of these excitations with the expectations for magnetic and phonon modes. The cuts presented in Fig. 5(c) and (f) clearly reveal the modes at $E=16.75$ meV and 43.75 meV to be phonons, as the $Q$-dependence of the intensity varies as $I\propto Q^{2}$. On the other hand, the cuts shown in Fig. 5(d) and (e) have $I\propto f(Q)^{2}$, where $f(Q)$ is the magnetic form factor for Yb3+, and therefore these excitations are identified as CEF levels. The phonon mode near 17 meV is consistent with the excess scattering in both MgTm2Se4 and MgHo2Se4 spectra at the same energy. The two identified CEF excitations constitute two of three predicted transitions for MgYb2Se4 at $T=5$ K, which is a $J=7/2$ systemKramers (1930), though the absence of the third transition in the measured spectra is significant in that it places a upper bound on its scattering intensity.

IV.1 CEF model fitting

The INS constant-Q cut data shown in Fig. 4(a)–(f) and Fig. 5(a) and (b) was fit was fit using the Stevens operator approach, which considers only the ground state J-multiplet determined by Hund’s rules. This assumption is what is known as the LS coupling model Ruminy et al. (2016), which has been demonstrated to be valid for rare earth atoms heavier than DyGaudet et al. (2018a); Princep et al. (2013); Rosenkranz et al. (2000) due to the large energy of the next J-multiplet Elliott et al. (1953); Blume et al. (1964). The associated CEF Hamiltonian is given by

[TABLE]

where $B_{nm}$ are the CEF parameters and $\hat{O}_{nm}$ are Stevens operator equivalentsStevens (1952), with the appropriate matrix elements for $\hat{O}_{nm}$ given by the software EasySpin Stoll and Schweiger (2006).

For the D3d point group symmetry of the RE-site, this reduces to:

[TABLE]

where we have chosen a quantization axis along the local $\langle 111\rangle$ directions.

For a given set of Bmn, CEF levels are found by direct diagonalization of Eq. 2, resulting in level schemes visualized in Fig. 6. The neutron scattering cross section of an excitation from the $i^{\text{th}}$ to the $j^{\text{th}}$ level is proportional to the matrix element given by

[TABLE]

where $\hat{J}_{\alpha}$ is the angular momentum operator in the $\alpha$ direction and $\ket{\psi_{i}}$ is the eigenket of the $i^{\text{th}}$ level. Total scattering intensity is modeled as the convolution of these matrix elements with a pseudo-Voigt instrument resolution function with a fitted width that was uniquely determined for each incident energy data set. Both excitations from the ground state and between excited states were considered, and each transition is weighted with the appropriate Boltzmann factors at a given temperature.

The initial fit parameters were found by rescaling the published CEF parameters from our previous work on the material or MgEr2Se4Reig-i Plessis et al. (2019) using

[TABLE]

as demonstrated in Refs. Bertin et al., 2012 and Rau and Gingras, 2018. For these rescalings, the lattice parameters $a$ and $a_{0}$ are taken from our XRD fits, we used $\left<r^{n}\right>$ found in Ref. Forstreuter et al., 1997, and we used the Stevens parameters, $\theta^{(n)}$ defined in Ref. Hutchings, 1964. Subsequent fits then represent improvements over the rescaling predictions for CEF parameters, as they properly account for differences in local structure and covalency between individual materials.

In Tables 3 , 4 and  5, we list the energies and predicted neutron scattering intensities of relevant excitations calculated using CEF parameters from both the scaling analysis and best fits of neutron data, discussed below. The positions of these transitions are indicated in Figs. 4 and 5 with vertical tick marks.

For both MgHo2Se4 and MgTm2Se4, several transitions contribute to each of the peaks observed in constant-Q cuts of scattering data, though the intensity was overwhelmingly dominated by excitations out of the ground state. In order to access more transitions by thermal population of excited levels, we also include measurements at $T=100$ K. Refinements of CEF parameters were performed via a global least squares minimization routine using the constant-Q cuts presented in Fig. 4(a)-(f) and the predicted scattering intensity from all CEF transitions expected in the measured temperature range. The best fits are shown as solid red curves in these figures, and with few exceptions reproduce both the magnitude and position of all major peaks while predicting no scattering intensity that was not observed by experiment.

For MgYb2Se4, only the ground state CEF level has appreciable occupation at temperatures below 200 K, simplifying the magnetic spectrum. However, the strong overlap between CEF and phonon peaks makes the above procedure untenable, as it fits raw neutron intensity and associates all non-background scattering with the Hamiltonian in Eq. 4. Instead, the constant-Q cut data in Fig. 5(a) and (b) were fit to multiple pseudo-Voigt peaks shown as solid lines, with resulting peak intensities listed in Table 5. The refinement of CEF parameters was subsequently performed by consideration of these fitted energies and intensities. To deal with the underconstrained nature of fitting 6 CEF parameters to only 4 pieces of information, we fixed the $B_{nm}$ values for $n=6$ to the initial rescaled values and only refined parameters $B_{20}$, $B_{40}$, and $B_{43}$. We subsequently verified that varying the parameters $B_{60}$, $B_{63}$ and $B_{66}$ had minimal impact on the predicted neutron peak intensity and associated analysis.

The $B_{nm}$ parameters resulting from fits are given in Table 6, along with uncertainties. For MgHo2Se4 and MgTm2Se4, uncertainties are determined by stepping out in one direction in parameter space while continually optimizing other parameters, until the reduced $\chi^{2}$ is increased by one. For the MgYb2Se4 fit, we again kept the $B_{nm}|n=6$ values fixed when finding uncertainties. The full implications of these fitted parameters for the CEF levels and low-temperature effective spin systems of the three materials are laid out more fully in the following sections.

IV.2 Potential effect of impurities

To consider the potential contribution to the CEF signal from impurity phases, we modeled the expected CEF scheme and the associated inelastic neutron scattering of relevant sesquiselenide and oxiselenide phases using a simple point charge model. For these, we assumed perfect ionic bonding, included all ions out to 7.5 Å, and used structures taken from the above XRD refinements. The potential was calculated in a tesseral harmonic expansion $\gamma_{nm}$. For the cosine ($m\geq 0$) and sine ($m<0$) components of the tesseral harmonics, we got the coefficients of the tesseral harmonics in Cartesian coordinates from Ref. Prather, 1961. The CEF parameters are calculated as $B_{nm}=\langle r^{n}\rangle\theta_{n}\left(1-\sigma_{n}\right)\gamma_{nm}$, where $\langle r^{n}\rangle$ is the radius term, and $\sigma_{n}$ is the shielding parameter; both values are taken from calculations in Ref. Edvardsson and Klintenberg, 1998. The CEF Hamiltonian is then constructed using Eq. 1, and predicted neutron intensities are calculated as laid out above. For the calculations, we used the software EasySpin Stoll and Schweiger (2006) to generate the matrix elements of the Stevens operators.

The calculated scattering from the CEF levels is scaled according to molar fraction of the ion in the sample and plotted in all of the constant-$Q$ cuts shown in Fig. 5(a)–(b) and Fig. 4(a)–(d) as solid green curves. Similarly calculated CEF parameters Edvardsson and Klintenberg (1998); Gupta and Sen (1973) are known to reproduce experimental values within about 20%Blok and Shirley (1966), and are sufficient to reproduce the general shape and integrated intensity of peaks in neutron scattering spectra. Within these bounds on uncertainty, inspection of the calculated spectra can potentially explain excess scattering in the MgHo2Se4 spectra at 1.3 meV and 18 meV, and may overlap with peaks in MgYb2Se4 at 17 meV and 35 meV. Overall however, the energies of CEF levels from the impurity phases seem to be well removed from levels associated with the majority phases, and are significantly less intense. We thus conclude that excitations from impurities have a minimal effect on the fits of CEF levels laid out above.

IV.3 Results and interpretation

In addition to producing the energy level schemes presented in Fig. 6, the fitted CEF parameters in Table 6 were used to calculate the associated eigenkets and, in particular, the ground state wavefunctions, which determine the size and anisotropy of moments in the low-temperature effective spin states. In Table 7, we list the resulting ground state wavefunctions for the three magnesium spinel compounds investigated in this paper, along with our previously determined results on the material MgEr2Se4Reig-i Plessis et al. (2019), included for comparison. The wavefunctions of degenerate doublets were determined with a small guide field artificially applied along the $\langle 111\rangle$ direction.

With no further analysis, one can immediately see that MgTm2Se4 has a ground state singlet with no net moment, in line with the previous conclusions of Ref. Ben-Dor and Shilo, 1980. For the other systems, the ground state is a doublet, which we can use to define a pseudo-spin-$\frac{1}{2}$ with effective up and down states (denoted $\ket{+}\text{and}\ket{-}$). With these states, we use

[TABLE]

to find the component of the moment parallel and perpendicular to the local $\langle 111\rangle$ directions. The results of these calculations are displayed in Table 8, where $\sigma_{i}$ are the Pauli matrices and $g_{zz}$ and $g_{xx}$ define $g_{\parallel}$ and $g_{\perp}$, respectively.

These values can be used to comment on the anisotropy of the effective spins. As an example, our previously determined results for MgEr2Se4 show $g_{\perp}=0$, indicating that material has fully Ising momentsReig-i Plessis et al. (2019). The current results imply that the moments in MgHo2Se4 also have perfect Ising symmetry, though one might expect deviations from this conclusion should one include spin-spin interactions which have the capacity to mix the CEF transitions.

For MgYb2Se4, we find $g_{\parallel}=5.188(79)$ and $g_{\perp}=0.923(85)$, implying an effective spin with strong anisotropy along the $\langle 111\rangle$ direction while falling far short of the Ising limit. These values imply significantly more anisotropy than the values of $g_{\parallel}=3.564$ and $g_{\perp}=2.204$ obtained from rescaling the CEF parameters from our previous MgEr2Se4 resultsReig-i Plessis et al. (2019); Rau and Gingras (2018), and are even farther removed from reports of nearly isotropic spins determined from fits of inverse magnetic susceptibility curvesHigo et al. (2017); Pawlak et al. (1988). Comparing the CEF parameters from these and the current study, the starkest contrast lay in the signs of $B_{43}$ and $B_{63}$ parameters and the magnitude of $B_{20}$. For MgYb2Se4, the parameters $B_{6m}|m=0,~{}3,~{}6$ have little consequence on the relative sizes of $g_{\parallel}$ and $g_{\perp}$, and on the goodness of the fit to our data, however both $B_{20}$ and the ratio between $B_{40}$ and $B_{43}$ are strongly associated with trigonal fields and thus the tendency of moments to point in the $\langle 111\rangle$ directions. The magnetic susceptibility studiesPawlak et al. (1988); Higo et al. (2017) underestimate the $B_{20}$ parameter and, more consequentially, fix the ratio between $B_{40}$ and $B_{43}$ to that expected in a perfect octahedral environment. This hugely underestimates the contribution to the potential from the next nearest neighbor atoms, and is likely responsible for the resultant underestimation of the anisotropy of the Yb3+ ions.

Finally, it is worth noting that the presence of optical phonons in several locations at the same energies as CEF levels raises the possibility that vibronic bound states might exist in these compounds. Vibronic bound states form as the result of strong magnetoelastic coupling, which produces significant hybridization/entanglement of CEF excitations and phononsThalmeier (1984), and have recently been observed in the related 227 pyrochlore compounds Ho2Ti2O7Gaudet et al. (2018b) and Tb2Ti2O7Ruminy et al. (2016). These bound states appear in INS spectra as a splitting of either CEF or phonon excitations, resulting in scattering from new modes with unusual momentum and temperature dependencesLoewenhaupt and Witte (2003). Of note in the current materials is the relatively intense optical phonon at E = 17 meV, which is close to the 17.7 meV CEF mode in MgHo2Se4 and 18 meV CEF mode in MgTm2Se4, as well as a possible 34 meV phonon mode which is near a 29.1 meV CEF mode in MgYb2Se4. Though initial inspection of our data reveals no clear evidence for vibronic modes, we suggest that the above materials and energies may be promising areas to search for bound states with follow-up higher resolution or polarized neutron scattering measurements.

V Magnetization

To check the CEF potential found by refinement of the INS data, we performed a series of magnetization measurements as a function of both temperature and applied field, with main results shown in Fig. 7. Symbols in this figure represent data, which is corrected for the demagnetizing field by assuming the powder sample is an oblate spheroid with a filling fraction of 60%. Solid lines in panels (a) – (c) represent predictions of a non-interacting model using the Hamiltonian in Eq. 2 with parameters Bnm from Table 6 and an additional Zeeman term to account for the role of applied field. Solutions of this modified Hamiltonian were found by direct diagonalization with the field pointing along $x$, $y$ and $z$ directions. For each direction, the expected moment is calculated using Boltzmann statistics before averaging to simulate a powder. This comprehensive approach is deemed more reliable than any that restricts attention to the ground state doublets only or treats the Zeeman term in the Hamiltonian perturbatively, as applied fields are known to both mix and shift the energy of low-lying excited CEF levels.

For all compounds, the measured and calculated magnetization show excellent agreement at high temperatures, as expected for strongly paramagnetic moments. This agreement extends to all temperatures for MgTm2Se4, which has a ground state composed of momentless singlets. For MgHo2Se4 and MgYb2Se4 however, the calculated curve begins to overestimate the measured values at the lowest temperatures. We attribute this discrepancy in MgYb2Se4 to the existence of net antiferromagnetic interactions, which are not accounted for in our independent moment CEF Hamiltonian. This conjecture is generally consistent with reports of negative Weiss constants in the literature on MgYb2Se4, which range from $\Theta_{CW}=-9.2$ KHigo et al. (2017) to $\Theta_{CW}=-44$ KPawlak et al. (1988), and reports of $\Theta_{CW}=-3.6(5)$ KYaouanc et al. (2015) and $\Theta_{CW}=-7.6(2)$ KLau et al. (2005) for CdHo2Se4, which is isostructural to MgHo2Se4. Though we caution against overinterpreting the results of Curie-Weiss fits in materials containing low-lying CEF transitions, these reports are sufficient to conclude antiferromagnetic interactions with an energy scale of a few Kelvin. For MgHo2Se4, we further note that the first two excited CEF levels (0.59 meV and 0.95 meV) are low enough in energy that interactions may mix these transitions with the ground state doublet and more fundamentally modify the effective spin state.

As a first step in exploring the role of interactions in these compounds, we also performed a series of self-consistent calculations using a Weiss molecular field, $\lambda$, and compared results to magnetization data for MgYb2Se4. Specifically, for a series of temperatures and applied fields, magnetization was defined as the solution to the transcendental equation $M=M_{0}(H+\lambda M)$, where $M_{0}$ is the calculated moment in the non-interacting model and $\lambda$ was determined by fitting to the data. In our analysis, we found $\lambda$=-3.4 mol-Yb cm*-3* for MgYb2Se4. The curves associated with this analysis are shown as dashed curves in Fig. 7(c) and (d), and reveal a much improved match to both field and temperature data over the independent spin model. The current model is also greatly improved over calculations using CEF parameters of Ref. Higo et al., 2017, which concluded Heisenberg-like moments from fits of susceptibility vs temperature data. In particular, one can see that the more isotropic model, shown as a dotted red line for $T=2$ K, seems to be heading towards a much higher saturation moment than either the data or the predictions of the current paper.

The impact of interactions is further observed in the inverse susceptibility vs temperature curve, shown in Fig. 7(d) for MgYb2Se4. Here, we plot the calculated inverse susceptibility both without and with the interactions as solid and dashed curves respectively, along with data shown as blue circles. Whereas the non-interacting model prediction is systematically low, the curve including interactions matches the data quite well. In the same figure we also show the data from Ref. Higo et al., 2017 as red diamonds, demonstrating consistency between the two data sets aside from a small offset which can attributed to a small amount of disorder. This punctuates the fact that both isotropic and highly anisotropic models are capable of describing magnetization versus temperature data at low fields, and higher field and spectroscopic measurements are absolutely essential if one wishes to determine information about the local CEF environment of local moments.

VI Discussion and conclusions

The current manuscript outlines the determination with INS of the symmetry-allowed CEF parameters for three members of the RE-spinel selenide family MgRE2Se4 (RE = Ho, Tm, and Yb). The parameters obtained are substantially different and demonstrably more accurate than previous efforts to determine CEF schemes by fitting magnetic susceptibility curves at low applied fields. This can be seen in the inability of parameters determined by the latter methods to either reproduce higher-field magnetization data or to successfully predict the energies of excited CEF transitions, which can be measured directly with INSPawlak et al. (1988); Higo et al. (2017). In contrast, the parameters listed in Table 6 have been shown to largely reproduce the neutron scattering intensity as a function of both energy and temperature and magnetization data over a wide range of applied fields and temperatures. We can use these parameters to not only determine the ground state wavefunction of each material, as presented above, but also to revisit the role that ground state and exited levels have on low temperature magnetic properties.

For example, our measurements of MgHo2Se4 reveal a ground state Ising doublet with $m=1.36(23)\mu_{B}$ moments and antiferromagnetic interactions, which may make this material amenable to a long-ranged ordered state similar to the one determines for CdHo2Se4Yaouanc et al. (2015). Significantly however, we also observe several low lying excitations, including a singlet at 0.591(36) meV, a doublet at 0.945(30) meV, and a second singlet at 2.88(7) meV. This situation is reminiscent of the materials Tb2Ti2O7 and Tb2Sn2O7, where virtual transitions associated with low-lying CEF levels are strongly suspected to renormalize the effective HamiltonianMolavian et al. (2007) and induce quantum fluctuationsTakatsu et al. (2017, 2011). In the spinel selenides, the larger lattice parameters and rare earth to anion distances result in excited CEF energies even closer to the elastic channel, which implies an even faster timescale for quantum fluctuations.

In MgTm2Se4, the first excited CEF level contributes to the physics in a different way. Whereas our INS analysis concludes a simple singlet ground state, consistent with expectationsBen-Dor and Shilo (1980); Pokrzywnicki et al. (1977), our INS spectra also reveals the existence of a low lying singlet at $E=0.876(16)$ meV. The Brillouin-function-like field dependence of magnetization in Fig. 7(b) demonstrates that the appreciable occupation of this excited level by increasing either field or temperature endows the Tm3+ atoms with a considerable finite moment, raising the intriguing possibility of stabilizing ordering phenomena at finite temperature with applied field, even as the system strives toward a singlet ground state at $T=0$ K. Further, recent theoretical workLiu and Huang (2019) has pointed out that the lowest two CEF levels form a quasi-doublet, allowing the magnetism in this material to be described by a transverse-field Ising model with a potential for an exotic quantum spin ice phase.

Only in MgYb2Se4 is the effective spin system well-isolated from the lowest excited level, at $E=26.01(56)$ meV. A major insight of this work however is how highly anisotropic the Yb3+ moments are in this system, which we infer not just from the analysis of our INS data but also from the saturation magnetization, which falls far short of expectations for isotropic spins. The idea of strongly anisotropic Yb3+ effective spins stands in contrast to earlier predictions of isotropic moments from magnetizationHigo et al. (2017) or weaker anisotropy from scaling argumentsRau and Gingras (2018). The Yb pyrochlore-lattice materials stand out as rare examples where anisotropic interactions have been calculated and semiclassical phase diagrams have been produced as a function of material propertiesRau and Gingras (2018). Thus, follow-up neutron diffraction measurements of low temperature ordered phases in this material might provide an opportunity to immediately test the validity of our results, and perhaps contribute to the understanding of the larger family of Yb2M2O7 pyrochloresRau and Gingras (2019).

Discussion of these three materials may naturally be grouped with consideration of sulfur (MgRE2S4) and cadmium (CdRE2X4) analoguesFlahaut et al. (1965, 1966) and recent reports of classical spin-ice behavior in MgEr2Se4Reig-i Plessis et al. (2019) and CdEr2Se4Lago et al. (2010); Gao et al. (2018). Together these publications show growing interest in the chalcogenide spinels, as a relatively unexplored family of highly frustrated magnets with a diversity of exotic states that rivals that of the 227 pyrochlore oxides. Proper consideration of local CEF environments is the first necessary step in modeling and fully understanding the associated physics.

VII Acknowledgments

Acknowledgements.

This work was sponsored by the National Science Foundation, under grant number DMR-1455264-CAR. D.R. further acknowledges the partial support of by the U.S. D.O.E., Office of Science, Office of Workforce Development for Teachers and Scientists, Office of Science Graduate Student Research (SCGSR) program. The SCGSR program is administered by the Oak Ridge Institute for Science and Education for the DOE under contract number DE-AC05-06OR23100. R. D. M. acknowledges the support of the Research Experience for Undergraduates program funded by the NSF under Grant No. 1659598. Synthesis and magnetization measurements were carried out in the Materials Research Laboratory Central Research Facilities at the University of Illinois. X-ray scattering measurements were performed at the Center for Nanophase Materials Sciences at Oak Ridge National Lab. A portion of this work used resources at the Spallation Neutron Source, which is a DOE Office of Science User Facility operated by Oak Ridge National Laboratory.