Deviation from the dipole-ice model in the new spinel spin-ice candidate, MgEr$_2$Se$_4$

TL;DR

This study introduces MgEr$_2$Se$_4$ as a new spin-ice candidate with a normal spinel structure, revealing deviations from classical dipolar-ice models and highlighting the role of quantum fluctuations in its magnetic behavior.

Contribution

The paper reports the synthesis and characterization of MgEr$_2$Se$_4$, demonstrating its unique magnetic properties and deviations from traditional spin-ice models, emphasizing the importance of quantum effects.

Findings

MgEr$_2$Se$_4$ has a normal spinel structure with Er$^{3+}$ moments on a pyrochlore sublattice.

The ground state is dipolar-octupolar with significant multipolar character.

Heat capacity and neutron scattering show ice-like features but deviate from classical dipolar spin-ice models.

Abstract

In spin ice research, small variations in structure or interactions drive a multitude of different behaviors, yet the collection of known materials relies heavily on the `227' pyrochlore structure. Here, we present thermodynamic, structural and inelastic neutron scattering data on a new spin-ice material, MgEr$_2$Se$_4$, which contributes to the relatively under-explored family of rare-earth spinel chalcogenides. X-ray and neutron diffraction confirm a normal spinel structure, and places Er$^{3+}$ moments on an ideal pyrochlore sublattice. Measurement of crystal electric field excitations with inelastic neutron scattering confirms that the moments have perfect Ising character, and further identifies the ground state Kramers doublet as having dipolar-octupolar form with a significant multipolar character. Heat capacity and magnetic neutron diffuse scattering have ice-like features, but are inconsistent with Monte Carlo simulations of the nearest-neighbor and next-nearest-neighbor dipolar spin-ice (DSI) models. A significant remnant entropy is observed as $T\rightarrow0$~K, but again falls short of the full Pauling expectation for DSI, unless significant disorder is added. We show that these observations are fully in-line with what is recently reported for CdEr$_2$Se$_4$, and point to the importance of quantum fluctuations in these materials.

Figures

Tables

Table 1. Table 1: Structural parameters obtained from the XRD and NPD refinements of MgEr 2 Se 4 data. The structural parameters in the lower part of the table are from the 38K NPD refinement.

MgEr2Se4 lattice parameters
space group Fd$\overline{3}$m
	XRD 300K	NPD 38K	NPD 470mK
a	11.5207(14)	11.4999(42)	11.5048(81)
${\chi }^2$	10.41	6.39	8.63
${\chi }^2$ Lebail	11.40	6.83	9.53
Se deficiency (%)	0.00(70)	0.00(98)	0.0(1.2)
Site inversion (%)	0.00(47)	0.0(3.7)	0.0(4.5)

Table 2. Table 2: The full CEF scheme of MgEr 2 Se 4 as calculated from the best fit to the data. The energy levels, relative neutron scattering intensity at 0 K, and wavefunctions are presented for the 8 CEF doublets associated with the ground state manifold.

MgEr2Se4 crystal field levels
$n$	$\delta E$	$I_{relative}$	${\Psi }_n$
0	0	-	$\pm 0.9165\left(7\right)|\pm 15/2⟩+0.360\left(1\right)|\pm 9/2⟩\pm 0.158\left(2\right)|\pm 3/2⟩-0.073\left(2\right)|\mp 3/2⟩\pm 0.0036\left(7\right)|\mp 9/2⟩+0.0035\left(14\right)|\mp 15/2⟩$
1	4.155	1	$0.734\left(5\right)|\pm 13/2⟩\mp 0.488\left(1\right)|\pm 7/2⟩+0.435\left(10\right)|\pm 1/2⟩\pm 0.177\left(2\right)|\mp 5/2⟩+0.0504\left(6\right)|\mp 11/2⟩$
2	6.279	0.0517	$\pm 0.480\left(9\right)|\pm 13/2⟩-0.071\left(5\right)|\pm 7/2⟩\mp 0.868\left(5\right)|\pm 1/2⟩-0.072\left(3\right)|\mp 5/2⟩\pm 0.066\left(4\right)|\mp 11/2⟩$
3	9.193	0.0862	$\pm 0.267\left(8\right)|\pm 15/2⟩+0.252\left(5\right)|\pm 9/2⟩\mp 0.918\left(16\right)|\pm 3/2⟩-0.13\left(11\right)|\mp 3/2⟩\pm 0.040\left(27\right)|\mp 9/2⟩-0.068\left(30\right)|\mp 15/2⟩$
4	10.133	0.2846	$\pm 0.6651\left(9\right)|\pm 11/2⟩-0.7238\left(8\right)|\pm 5/2⟩\mp 0.097\left(4\right)|\mp 1/2⟩-0.0107\left(7\right)|\mp 7/2⟩\mp 0.156\left(4\right)|\mp 13/2⟩$
5	27.273	0.0329	$\mp 0.692\left(1\right)|\pm 11/2⟩-0.571\left(2\right)|\pm 5/2⟩\mp 0.187\left(2\right)|\mp 1/2⟩+0.342\left(3\right)|\mp 7/2⟩\mp 0.207\left(2\right)|\mp 13/2⟩$
6	29.91	0.0188	$+0.290\left(1\right)|\pm 15/2⟩\pm 0.8977\left(4\right)|\pm 9/2⟩+0.3155\left(15\right)|\pm 3/2⟩\pm 0.102\left(2\right)|\mp 3/2⟩+0.011\left(3\right)|\mp 9/2⟩\mp 0.0014\left(9\right)|\mp 15/2⟩$
7	29.945	0.0055	$\pm 0.4035\left(14\right)|\pm 13/2⟩+0.7995\left(14\right)|\pm 7/2⟩\pm 0.1098\left(10\right)|\pm 1/2⟩+0.3437\left(3\right)|\mp 5/2⟩\mp 0.269\left(2\right)|\mp 11/2⟩$

Equations

$$H=\sum _{nm}{B}_n^m{O}_n^m,$$

$$\begin{gathered} |{\psi }_0^{\pm }⟩=\pm 0.9165(7)|\pm 15/2⟩+0.3600(11)|\pm 9/2⟩ \\ \pm 0.1581(16)|\pm 3/2⟩-0.0731(15)|\mp 3/2⟩ \\ \pm 0.0036(7)|\mp 9/2⟩+0.0035(14)|\mp 15/2⟩, \end{gathered}$$

$$⟨{\psi }_0^{+}|J_x|{\psi }_0^{-}⟩=⟨{\psi }_0^{+}|J_y|{\psi }_0^{-}⟩=0$$

$$H=H_{\text{CEF}}+H_Z.$$

$$\begin{gathered} H=-3J_{nn}\sum _{⟨i,j⟩}{\mathbf{S}}_i^{z_i}\cdot {\mathbf{S}}_j^{z_j}-3J_{nnn}\sum _{⟨⟨i,j⟩⟩}{\mathbf{S}}_i^{z_i}\cdot {\mathbf{S}}_j^{z_j} \\ +\frac{3D_{nn}}{5}{r}_{nn}^3\sum _{i<j}\left(\frac{{\mathbf{S}}_i^{z_i}\cdot {\mathbf{S}}_j^{z_j}}{|{\mathbf{r}}_{ij}{|}^3}-\frac{3({\mathbf{S}}_i^{z_i}\cdot {\mathbf{r}}_{ij})({\mathbf{S}}_j^{z_j}\cdot {\mathbf{r}}_{ij})}{|{\mathbf{r}}_{ij}{|}^5}\right), \end{gathered}$$

$$C_V(T)=\frac{dE(T)}{dT}$$

$$H=\sum _{nm}{B}_n^m{O}_n^m,$$

$$H=H_{\text{CEF}}+H_Z.$$

$$M_n=\frac{⟨{\psi }_n|\mathbf{H}\cdot \hat{J}|{\psi }_n⟩}{\left|\mathbf{H}\right|}.$$

$${\int }_0^{\frac{\pi }{2}}{g}_l\sum _n{p}_n(\mathbf{H},\theta )M_n(\mathbf{H},\theta )\sin (\theta )d\theta ,$$

Full text

I Introduction

One of the most notable aspects of the spin-ice class of compounds is the wide range of interesting behaviors they exhibit. Originally of interest as a magnetic analogue of water ice and a playground for thermodynamic modelsHarris et al. (1997, 1998); Bramwell and Harris (1998); Moessner (1998); Bramwell and Gingras (2001); Ramirez et al. (1999), the field has expanded dramatically to include an assortment of interesting variations, including Kagomè iceTabata et al. (2006); Fennell et al. (2007), ordered spin iceMirebeau et al. (2005), dynamic spin iceZhou et al. (2008) and quantum spin ice (QSI)Gingras and McClarty (2014); Savary and Balents (2017); Sibille Romain et al. (2018). The QSI materials are prime candidates for a U(1) quantum spin liquid (QSL)Huse et al. (2003); Hermele et al. (2004); Savary and Balents (2012); Lee et al. (2012); Gingras and McClarty (2014); Hao et al. (2014); Petrova et al. (2015); Shannon et al. (2012); Chen (2016); Benton et al. (2012), which itself may have multiple variationsTaillefumier et al. (2017) including distinct, symmetry-enriched phasesHuang et al. (2014); Li and Chen (2017).

The unifying feature in this panoply is the underlying “classical” spin-ice model, wherein spins on a pyrochlore lattice with local $\langle 111\rangle$ Ising anisotropy and ferromagnetic interactions freeze into an extensively degenerate “ice” phase, characterized by a 2-in-2-out (TITO) constraint on the constituent tetrahedra. This TITO constraint famously maps onto the Bernal-Fowler ice rules for proton-oxygen bond lengths in frozen waterBernal and Fowler (1933), and the associated remnant ‘Pauling entropy’ as $T^{*}~{}\rightarrow~{}0~{}\text{K}$Pauling Linus (1935); Ramirez et al. (1999) remains the primary experimental signature of a classical spin-ice state. The TITO constraint can further be mapped to a divergence-free flux, allowing one to reinterpret the ice as a “magnetic Coulomb” phaseHenley (2010) wherein thermodynamic properties can be calculated by considering a gas of magnetic monopolesCastelnovo et al. (2008); Morris et al. (2009); Kadowaki et al. (2009).

In real spin ice materials, the local Ising condition is a result of trigonal crystal electric fields (CEF), and effective ferromagnetic interactions emerge from summing nearest-neighbor exchange and dipole termsBramwell and Gingras (2001); Gardner et al. (2010). This dipolar spin-ice (DSI) modelden Hertog and Gingras (2000); Melko and Gingras (2004); Isakov et al. (2005) is sufficient to explain the origin of the ice state, and has been successful in reproducing measured heat capacityZhou et al. (2012) (HC) and basic features of neutron diffuse scattering patternsFennell et al. (2007); Morris et al. (2009); Chang et al. (2010); Sen et al. (2013) in known classical spin ices; this includes Ho2Ti2O7Harris et al. (1997); Fennell et al. (2009), Dy2Ti2O7Ramirez et al. (1999), and associated stannates ($R_{2}$Sn2O7)Ke et al. (2007, 2008) and germanates ($R_{2}$Ge2O7)Zhou et al. (2011); Hallas et al. (2012). The breadth of behaviors described above, however, is a testament to the importance of further degeneracy breaking terms. Further neighbor exchange is needed to explain details of diffuse scattering and reproduce measured critical fieldsRuff et al. (2005); Melko and Gingras (2004). Quantum fluctuations result from either transverse molecular fieldsTomasello et al. (2015) or from multipolar superexchange interactionsOnoda and Tanaka (2011); Rau and Gingras (2015); Iwahara and Chibotaru (2015), with the latter invoked to explain experimental data in Pr2Sn2O7 and Pr2Zr2O7Onoda and Tanaka (2010); Kimura et al. (2013); Petit et al. (2016); Sibille et al. (2016). The unique dipolar-octupolar (DO) character of moments in materials such as Nd2Zr2O7Xu et al. (2016); Huang et al. (2014); Lhotel et al. (2015) and Ce2Sn2O7Sibille et al. (2015) is linked to the possibility of symmetry-enriched QSL phasesHuang et al. (2014); Li and Chen (2017).

There is thus clear motivation to extend the study of spin-ice physics to materials beyond the 227 oxides, with different variations in local structure and interactions. The cubic spinels (AB2X4) are prime candidates, as they share the same Fd$\bar{3}$m space group and pyrochlore sublattice as the 227 compounds, but differ in the octahedral arrangement of local chalcogen anions about B-site spin positionsLee et al. (2010); Takagi and Niitaka (2011); Radaelli (2004). Sizable trigonal CEFs create $\langle 111\rangle$ easy axes on some B-site ions, which play a defining role for material propertiesMun et al. (2014); MacDougall et al. (2012); Ehlers et al. (2012). Ferromagnetically coupled $\langle 111\rangle$ easy axis spins reminiscent of spin ices have been reported in several spinels leading to frustration observed through diffuse scattering in single crystals Tomiyasu et al. (2011), or leading to two-in-two-out ordered states of the B-site sublattice MacDougall et al. (2012, 2014); Ma et al. (2015); Lee et al. (2017). In the singular system, CdEr2Se4, remnant Pauling entropy has been reportedLago et al. (2010), and a very recent study on the same material has claimed DSI-like spin-correlations and an anomalously fast monopole hopping rateGao et al. (2018).

Here, we present data on a new spinel, MgEr2Se4, which provides another interesting counterpart to known spin-ice materials. We provide x-ray (XRD) and neutron powder diffraction (NPD) data which confirm an ideal pyrochlore sublattice of Er3+ moments, but with a cubic lattice parameter $\approx 10\%$ larger than Dy2Ti2O7. Inelastic neutron scattering (INS) data reveal that the moments have ideal Ising anisotropy, and further show that they have a significant multipolar character, in fact demonstrating the characteristic DO symmetryHuang et al. (2014). Both heat capacity and magnetic diffuse scattering data exhibit qualitative features of classical spin ice correlations. Follow-up Monte Carlo (MC) simulations, however, show that the collective data are inconsistent with nearest-neighbor and next-nearest neighbor DSI models. This may be a natural consequence of the multipolar character of the Er3+ moments, which seed significant quantum fluctuations.

II Sample preparation and crystal structure

Polycrystalline samples of MgEr2Se4 were prepared via a two-step solid state reaction, following the method described by Flahaut Flahaut et al. (1965). The precursors MgSe and Er2Se3 were prepared by the direct reaction of stoichiometric amounts of the elements at $650^{\circ}$C. Stoichiometric quantities of the two precursors were then combined, pelletized and reacted in vacuum at $1000^{\circ}$C for two days, this step was repeated at least one more time for the precursors to fully react. Structure and purity were confirmed using a PANalytical X’Pert3 X-ray powder diffractometer at the Center for Nanophase Materials Sciences at Oak Ridge National Laboratory (ORNL). NPD measurements were performed with the HB-2A powder diffractometer at ORNL’s High Flux Isotope Reactor, using a 3.6 g sample and neutron wavelengths of $\lambda=1.54~{}\text{\AA}$ and $\lambda=2.41~{}\text{\AA}$, with collimators open-21′-12′ and open-open-12′, respectively. Structural refinements were performed using the FULLPROF Rodríguez-Carvajal (1993) software suite. Additional measurements in a magnetic field used the CTAX instrument with $\lambda=5~{}\text{\AA}$ neutrons. INS was performed with the SEQUOIA Granroth et al. (2006) fine-resolution Fermi chopper spectrometer at ORNL’s Spallation Neutron Source (SNS). Measurements were collected with incident energies $E_{i}=30~{}\text{meV}$ and $E_{i}=50~{}\text{meV}$ with the fine Fermi chopper spinning at frequencies of 300 Hz and 360 Hz respectively. Magnetization and specific heat measurements were performed in the Seitz Materials Research Laboratory at Illinois using a Quantum Design MPMS3 and PPMS, respectively.

A large volume sample of MgEr2Se4 was prepared for exploration with neutron scattering and, unless stated otherwise, was used to obtain all data presented in the figures below. Purity and structure were studied with both x-ray diffraction (XRD) and neutron powder diffraction (NPD), and diffraction patterns on our primary sample are shown in Fig. 1, along with the results of FULLPROF refinements. Impurity peaks in both patterns are denoted by crosses, and were largely accounted for by the orthorhombic phase of Er2Se3 (1.8–5.0 %), which is consistent with a small amount of Mg evaporating during synthesis. In addition to Er2Se3, refinements also showed small amounts of Er2O2Se (2.3–2.5 %) and elemental Er (1.2 %), as well as some small unindexed impurity peaks which are not consistent with any known compound. The best fit refinement implies that the sample had purity of 94.66(62) % and 91.4(1.5) % by mass from the XRD and NPD fits respectively. A model independent estimate of the impurity fraction of 7.8(1.2) % by mass was obtained by comparing the integrated intensity of XRD Bragg peaks associated with the majority phase and the sum intensity of everything else above background. The weighted average of the three estimates gives a value of 6.2(1.2) % impurity phase in the sample.

Both NPD and XRD show that the MgEr2Se4 phase is of high quality and shows no observable defects in the structure. To test for any structural defects, we allowed Se occupancy, Mg occupancy as well as Er and Mg site inversion to vary. Results of best fits are shown in Table 1 and show no such defects with bounds of $<1~{}\%$.

As an additional model independent check for point defects, we performed Le Bail refinements Le Bail et al. (1988) of our data and compared the $\chi^{2}$ of those fits to our best Rietveld refinement. In Lebail fits, the structure factor is not calculated, and instead every peak height is allowed to vary and are fit independently – effectively identifying the ideal description for peaks associated with a single phase in a mixed powderToby (2006). In the current case, we see that the $\chi^{2}$ achieved through a Le Bail peak-by-peak fitting of the majority phase is no smaller than that achieved via the above Rietveld refinement. This is a powerful result, which effectively eliminates the existence of cation inversion, off-stoichiometry on the Se-sublattice, or any other point defect which has the capacity to change the height in a neutron scattering pattern.

III Inelastic neutron scattering

The local CEF environment of Er3+ was determined with INS. In Fig. 2(a), we show a representative INS data set collected using $E_{i}=30$ meV at $T=5$ K, wherein CEF excitations out of the ground state appear as bright dispersionless modes near $E\approx 4$ meV and $E\approx 10$ meV. To access transitions form higher energy levels additional measurements were performed at $T=40$ K and $T=150$ K, the same temperatures were also measured at $E_{i}=50$ meV. In Figs. 2(b–d) and (f–h), we show cuts obtained by integrating the INS data over momentum interval $Q=[2,2.5]\text{\AA}^{-1}$, chosen to maximize the available energy range. A background contribution (shown as a red line) is interpolated from hand picked points at energy transfers away from CEF peak positions. The scattering intensities of the six data sets were fit simultaneously to expectations from the model crystal field Hamiltonian

[TABLE]

where $O_{n}^{m}$ are the Stevens’ operatorsStevens (1952). Fits were performed using a mix of random walk grid search and gradient search methods explained in detail in Appendix C and began with an initial guess calculated using a point charge modelHutchings (1964) and the known structure. The CEF levels and electron wavefunctions were calculated with the quantization axis along the local $\langle 111\rangle$ directions, with only coefficients $B_{2}^{0},B_{4}^{0},B_{4}^{3},B_{6}^{0},B_{6}^{3},B_{6}^{6}$ which can be non-zero by symmetry. Predicted peaks were convoluted with a Voigt profile to account for instrument resolution. Through simultaneous consideration of nearly two dozen observed peaks, we determined the six most likely CEF parameters for MgEr2Se4 to be (in meV): $B_{2}^{0}=-4.214(63)\times 10^{-2},B_{4}^{0}=-6.036(30)\times 10^{-4},B_{4}^{3}=-1.3565(67)\times 10^{-2},B_{6}^{0}=3.264(16)\times 10^{-6},B_{6}^{3}=-3.791(75)\times 10^{-5}\text{and},B_{6}^{6}=2.194(65)\times 10^{-5}$. The scattering pattern associated with these parameters is denoted by solid lines in Figs. 2(b–d) and (f–h), and successfully reproduces both intensity and position of all considered modes and predicts no errant peaks. The data further reveals the existence of two small peaks at energies $E\approx 2$ meV and $E\approx 5$ meV, with spectral weight consistent with the $\approx 7\%$ impurity in the sample.

Figure 2(e) shows the crystal field levels calculated from the above $B_{n}^{m}$, along with the various transitions observed in our scattering data. The associated wavefunctions are given in their entirety in Appendix C, and the ground state Kramers doublet is

[TABLE]

which implies perfectly Ising spins

[TABLE]

with moment $m_{z}=g_{J}\bra{\psi_{0}^{+}}J_{z}\ket{\psi_{0}^{+}}=8.3(1)\mu_{B}$, where the Landé g-factor $g_{J}=\frac{6}{5}$ for Er3+. This wavefunction also facilitates significant non-dipolar superexchange interactions, as discussed below. The two lowest excited doublets are at energies $E_{1}=4.02(2)~{}\text{meV}$ and $E_{2}=6.40(2)~{}\text{meV}$, significantly larger than interaction energies determined below, but still low enough to impact thermodynamic properties at temperatures $T>5$ K.

IV Magnetization

To check the validity of the CEF fits, inferred levels were used to calculate the magnetization in the paramagnetic phase for a range of applied fields. For the case of low lying CEF excitations, the effect of mixing of excited CEF levels must be taken into account for calculating the magnetization Bonville et al. (2013).

Magnetization was obtained using a non-interacting model with total Hamiltonian of the $J=15/2$ Er*+3* multiplet CEF plus the Zeeman energies

[TABLE]

The CEF Hamiltonian is defined as before in Eq. 1 and the Zeeman term $H_{Z}=-g_{J}\mathbf{H}\cdot\hat{J}$. For powder averaging, the magnetization was calculated for over 1500 different local applied field directions for each temperature and field value. Finally Boltzmann statistics were used to find the occupation of each perturbed CEF level and then calculate the associated magnetization. The results are shown as solid lines in Fig. 3. As one can see, this calculation largely reproduces the magnitude, temperature and field dependence of measured magnetization with zero fit parameters. Particularly notable, the calculation was able to reproduce the linear field dependence for $H>2$ T, which we confirmed is the result of the field perturbation of local eigenstates in excited doublets and is not captured in the simplified pseudospin 1/2 models for magnetization that has been used to argue for Ising behavior in past rare earth compoundsBramwell et al. (2000); Xu et al. (2015); Sibille et al. (2015); Hallas et al. (2015); Anand et al. (2016). The overall agreement here provides strong support for the inferred levels above and the Ising character of the Er3+ moments at temperatures $T<5$ K. Although we used a non-interacting model, we believe the agreement at $T=2$ K is to be expected given the small exchange and dipole energies compared to 2 K and the degree of frustration in this material. Additional detail on the calculation method and a comparison to a pseudospin 1/2 models are given in Appendix D.

V Heat Capacity

Measured heat capacity in the range $0.45~{}\text{K}<T<20~{}\text{K}$ is shown in Fig. 4(a), with solid lines representing contributions from phonons, CEFs, impurities and magnetic degrees of freedom, and the linewidth representing the uncertainty. The CEF contribution was calculated exactly from the multi level partition function given by the CEF scheme determined above, and contains non-trivial contributions from both the E1 and E2 doublets. The contribution to the heat capacity from the impurity phases was taken into account by subtracting off the heat capacity expected from a system with transition energies at $E\approx 2$ meV and 5 meV, and appropriately scaled to be between 5 and 9 %; molar mass of the sample was also accordingly scaled. The energies chosen for the impurity phase come from extra modes observed in INS. A small amount of an impurity was found to order from the neutron powder diffraction and was included in the error bars of the impurity contribution – more details can be found in Appendix B. Phonons were modeled by fitting data in the range $7~{}\text{K}<T<25~{}\text{K}$ to the Debye model after CEF and impurity contributions were subtracted. The black line in Fig. 4(a) is the sum of these contributions, and is seen to perfectly describe data above 10 K. The remaining contribution was entirely attributed to the Er3+ moments on the pyrochlore sublattice. This Er*+3* magnetic contribution is dominated by a single broad peak near $T^{*}\approx 1.1~{}\text{K}$, which has a height and position broadly consistent with expectations for the DSI modelden Hertog and Gingras (2000). A similar analysis was done with the heat capacity data of CdEr2Se4 taken from Ref. Lago et al., 2010, after suitably updating the treatment of CEF energy levels using recently measured values taken from Ref. Gao et al., 2018. We plot the resulting magnetic contribution in CdEr2Se4 as red crosses in Fig. 4(a) and the corresponding calculated entropy as a solid red line in Fig. 4(b). We find both the heat capacity and residual entropy of the two erbium selenium spinels to be remarkably similar, adding to confidence in the data and analysis.

To gain further insight into the exact temperature dependence of the heat capacity, we performed MC simulations on a 2048-site cluster with periodic boundary conditions in all directions, as described in Appendix A. This simulation used the Hamiltonian:

[TABLE]

where the parameter $J_{nn}$ ($J_{nnn}$) represents nearest (next-nearest) neighbor exchange interactions, and the strength of the dipole interaction was fixed to $D_{nn}=1.06$ K, as determined from the measured structure. For $J_{nnn}=0$, we confirmed that our calculations match previously published resultsden Hertog and Gingras (2000). We found that the primary effect of low $J_{nnn}$ was to symmetrically broaden and increase the height of the peak in heat capacity, regardless of sign. As discussed in the Appendix A, sufficiently large $J_{nnn}$ are seen to drive the system to either a $Q=0$ or $Q=X$ ordered state for ferromagnetic or antiferromagnetic interactions, respectively.

The best fit of our magnetic heat capacity data gave values $J_{nn}=0.06$ K and $J_{nnn}=-0.1$ K, producing the curve shown in the inset of Fig. 4(a). This fit described the data adequately (reduced $\chi^{2}=16.7$), and significantly better ($\chi^{2}=68.8$) than the curve expected using exchange parameters estimated from fits of NPD data, $J_{nn}=0.06$ K and $J_{nnn}=0$ K, discussed in more detail below. In fact, though small on the scale of $D_{nn}$, we note that our MC calculations were found to be entirely incompatible with the HC data without assuming a ferromagnetic $J_{nn}$ and setting $\left|J_{nnn}/J_{nn}\right|>1$. Both observations stand in contrast to known 227 classical ice systems, but the former may be consistent with the near $90^{\circ}$ Er-Se-Er superexchange path between nearest neighbors in the spinel structure. As discussed below, the sizable value for Jnnn is inconsistent with our magnetic diffuse scattering data and may point to other relevant physics.

In Fig. 4(b), we plot the change in entropy from 0 K obtained from the integrated magnetic heat capacity from data in fields $H=0$ T and $H=0.25$ T, where MC data was used to extrapolate below the first data point at $T=0.45$ K . The relatively small field of $H=0.25$ T was chosen in order to minimize changes to the low lying CEF levels and the broadening of those peaks in the HC data, allowing consistent analysis of data taken both in and out of field. Shading represents experimental uncertainty, which is dominated by uncertainty in the impurity volume fraction. The $H=0$ T data reveal a sizable residual entropy which is partially relieved with small applied fields, broadly consistent with spin-ice behaviorRamirez et al. (1999), but significantly less than the value of $1/2\ln\left(3/2\right)$ predicted for the DSI modelPauling Linus (1935); Ramirez et al. (1999). This is true for both the current data, but also the data taken from Ref Lago et al., 2010 on related material CdEr2Se4 after correcting for contributions from subsequently measured CEF excitations measured recentlyGao et al. (2018). Intriguingly, one sees that the remnant entropy of both systems approach the same value, but fall far short of the full Pauling value. This agreement despite the differing level of purity in the two materialsLago et al. (2010); Gao et al. (2018) builds confidence that the data reflect intrinsic physics. The deviation from Pauling entropy implies that some TITO spin configurations are being removed from the degenerate manifold by an interaction term outside the DSI model.

Though consideration of the solid curve in Fig. 4(a) might imply the degeneracy breaking term could be $J_{nnn}$, this conclusion is not supported by our NPD measurements of magnetic diffuse correlations, as discussed below.

VI Neutron powder diffraction

A summary of our NPD results is presented in Fig. 5. In Fig. 5(a), we show that the diffuse background at $T=38$ K is dominated by paramagnetic scattering, which fit well to Er3+ form factor squared and was used to normalize subsequent data. The paramagnetic scattering was then used to plot data at lower temperatures on an absolute intensity scale- an important step in the analysis of the diffuse scattering, since MC simulations reveal that much of the relevant information regarding next-nearest neighbor interactions is encoded in the scattering intensity at particular values of $Q$. In panels (b)-(d), we plot on an absolute scale the scattering data taken with neutron wavelength $\lambda=2.41$ Å  at $T=0.47$ K, 1.5 K and 4 K, with the contribution at $T=38$ K subtracted to isolate the magnetic contribution. Low temperature magnetic correlations are largely short-ranged and consistent with an ice-like state. A handful of Bragg peaks were observed in the $T=0.47$ K pattern only, reliably associated with the small impurity fraction and subtracted from the pattern in Fig. 5(b). A second data set taken with $\lambda=1.54$ Å  neutrons is included in the same panel, which is largely consistent with the $\lambda=2.41$ Å  data except for discrepancy at the lowest angles, which we attribute to increased background from the proximity to the direct beam at $\theta=0$ in the $\lambda=1.54$ Å  data. The inset of Fig. 5(b) shows NPD data under applied field where the short-ranged ice correlations partially give way to the magnetic Bragg peak at the (2 0 0) position with $H=0.2$ T, consistent with the partial recovery of remnant entropy over the same field range.

Most significant, however, are the solid curves in Fig. 5(b)-(d), which represent the predicted scattering pattern associated with the spin configurations from the above MC consideration of heat capacity data. Though the $T=4$ K data is largely consistent with MC predictions, data at the lower two temperatures deviate significantly in the region $Q\approx 0.55$ Å*-1*. This is the (1 0 0) Bragg position, and can be interpreted as an excess of predicted $Q=X$ correlations driven by the sizable $J_{nnn}$ needed to reproduce the width of the heat capacity peak. When using the same Hamiltonian to fit the NPD data instead (dashed curves), we get $J_{nnn}=0$ which minimizes the $Q=X$ correlations. Not only are these fits poor, but as seen above in consideration of Fig. 4, the inferred parameters lead to a systematic underestimation of heat capacity. We thus conclude that the next-nearest-neighbor DSI model is incapable of explaining the collective behavior of MgEr2Se4, leading to consideration of other effects.

VII Discussion and conclusions

In the context of the above data, it is important to consider the potential role of random fields due to local disorder, which have been suggested as a possible route to a QSL state in Pr pyrochloresMartin et al. (2017); Wen et al. (2017). We consider this explanation unlikely here, given the strong agreement between results on MgEr2Se4 and CdEr2Se4Lago et al. (2010); Gao et al. (2018). To address this possibility, however, we have measured magnetization and heat capacity of a second sample of MgEr2Se4, with a scattering pattern which implies significant structural disorder. As shown in Appendix E, magnetization measurements on this sample were also well explained by our CEF calculations, implying local Ising physics. However, heat capacity revealed a peak which is heavily skewed towards higher temperatures and, surprisingly, exhibits full remnant Pauling entropy. Thus, it seems that disorder impedes, rather than encourages, the mechanisms leading to deviations from the DSI model.

We thus consider one last possibility: that significant quantum fluctuations are driven by transverse spin couplings. On its face exotic, this option is in fact the least speculative, as quantum fluctuations have been predicted for Kramers doublets of the form we have observedOnoda and Tanaka (2011); Huang et al. (2014); Iwahara and Chibotaru (2015); Rau and Gingras (2015); Li and Chen (2017). In contrast to the dominant dipolar character of moments in Dy2Ti2O7 and Ho2Ti2O7Rau and Gingras (2015), the current work and Ref. Gao et al., 2018 show that moments in the Er-spinels contain a sizable ($\approx$1/3) leading order multipolar correction. Such corrections should create transverse exchange couplings $J_{\perp}\propto\delta_{i}^{2}/\alpha^{2}$, where $\alpha$ and $\delta_{i}$ are the coefficients of the $\ket{15/2}$ and next leading order $\ket{J_{z}}$ term in the ground state wavefunctionRau and Gingras (2015). Yb2Ti2O7, known for complex quantum behavior from proximity to competing phases Robert et al. (2015); Jaubert et al. (2015) including a quantum spin liquid state Savary and Balents (2013), has large transverse coupling experimentally determined as $J_{\perp}/J_{z}\approx 0.30$Ross et al. (2011). This level of transverse coupling has a large effect on the material’s propertiesChanglani (2017). Comparing to MgEr2Se4, our data implies a multipolar exchange that gives $J_{\perp}/J_{z}\approx 0.15$. This is not a negligible effect, and should have immediately measurable consequences. The anomalously fast monopole hopping rates recently reported for CdEr2Se4Gao et al. (2018) may be one such example. More direct confirmation may come from diffuse scattering measurements of single crystals, which would also facilitate tests of novel predictions for materials with DO doubletsHuang et al. (2014); Li and Chen (2017).

Taken together, the collective data on MgEr2Se4 paint a picture of a spinel-based pyrochlore which provides an interesting counterpart to existing 227 oxides. The diffraction, heat capacity and inelastic neutron scattering results above leave very little doubt that this material contains the lattice, Ising anisotropy and ferromagnetic exchange necessary to drive spin ice behavior, and there is strong circumstantial evidence to infer significant quantum fluctuations. We further note that MgEr2Se4 is just one member of a series of magnesium rare earth selenidesHigo et al. (2017); Reig-i Plessis et al. (2018), some of which have an even larger capacity for quantum effects. These results, coupled with recent work on CdEr2Se4 and CdEr2S4Gao et al. (2018), may portend the opening up of a new class of magnetic spinel chalcogens, which can contribute meaningfully to the current research into pyrochlore materials.

Acknowledgements.

The authors acknowledge useful discussions with P. Schiffer, G. Sala and G. Chen. 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. Synthesis and thermodynamic measurements were carried out in the Materials Research Laboratory Central Research Facilities, University of Illinois. Scattering measurements were conducted at the Center for Nanophase Materials Sciences, at the High Flux Isotope Reactor and at the Spallation Neutron Source, each DOE Office of Science User Facilities operated by the Oak Ridge National Laboratory. This work is part of the Blue Waters sustained petascale computing project, which is supported by the National Science Foundation (award numbers OCI-0725070 and ACI-1238993) and the State of Illinois.

Appendix A Monte Carlo simulation

We performed Monte Carlo simulations on a 2048-site cluster with periodic boundary conditions in all directions. Our Hamiltonian (Eq. 3) includes nearest and next nearest neighbor Ising interactions and long range dipole-dipole interactions:

[TABLE]

To deal with the dipole term, one wants to sum over an infinite number of images. Because the sum is conditionally convergent, the order in which this sum is taken affects the answer. One approach to choosing the order of this sum is to use the Ewald technique. An alternative approach (see section 4 of Ref. de Leeuw et al., 1980), which we utilize, is to regularize the sum spherically by multiplying the contribution of the image which is in box $n=(n_{x},n_{y},n_{z})$ by $\exp\left(-s\left|n\right|^{2}\right)$. This regularization forces the sum to be absolutely convergent for all $s>0$. We then numerically extrapolate to $s\rightarrow 0$ by evaluating the real space component from several values of finite $s$. A sample extrapolation plot for a fixed $(i,j)$ is shown in Figure 6. Therefore both Ewald and this approach account for the long-range part of the interaction beyond simple truncation. We benchmarked this method against the results obtained via Ewald summation den Hertog and Gingras (2000) obtaining the same results. All data presented in the main text is obtained using the extrapolated parameters.

The specific heat,

[TABLE]

is computed by taking the derivative of this spline. Spin configurations were fed into the program SPINVERT Paddison et al. (2013) to obtain predicted powder-averaged diffuse neutron scattering patterns, $I(Q)$.

For a given choice of $J_{nn}$, it was found that sufficiently large $|J_{nnn}|$ drove a transition into a long-range ordered state. Figure 7 shows the neutron $I(Q)$ and real space pattern of spins for the case of both ferromagnetic and antiferromagnetic $J_{nnn}$. Both ordered states preserve the two-in-two-out constraint of the spin-ice ground state. The ferromagnetic state prefers a state which preserves the symmetries of the Fd$\bar{3}$m space-group, and thus demonstrated spin-spin correlations at locations consistent with a Q=0 state. Real space spin configurations indicate a similar state as preferred by Ho2Ti2O7 and Dy2Ti2O7 for applied fields $\mathbf{H}\parallel[001]$Fennell et al. (2005), or by ferrimagnetic spinelsMacDougall et al. (2012). The expected scattering pattern $I(Q)$ for this state is shown in Fig. 7(a) and the corresponding magnetic structure is depicted in Fig. 7(c).

The antiferromagnetic interaction case preferred a state which broke Fd$\bar{3}$m symmetry, and demonstrated distinct anti-correlations between chains of spins along the $[110]$ direction which are antiparallel to neighboring chains. Neutron intensity indicates a $Q=X$ phase, in that it shows correlations at the cubic $(100)$ and equivalent Bragg positionsHarris et al. (1997). Although similar to the $Q=X$ phase favored by $\mathbf{H}\parallel[110]$ fields in that both have antiparallel chains of spins along the $[110]$ directionFennell et al. (2005), the current phase is actually distinct, in that there is no net moment. These plots are shown in Fig. 7(b) and (d) for $I(Q)$ and real space respectively.

Appendix B Magnetic Properties of the Impurity Phase

In Fig. 8 we show the low temperature NPD patterns at $T=1.5$ K and 0.47 K, without the $T=38$ K pattern subtracted. In addition to broad features associated with spin-ice correlations and Bragg peaks associated with the lattice, the data shows a series of weak Bragg magnetic peaks which appear only at the lowest temperature. These peaks were not indexable in the Fd$\bar{3}$m space group of the spinel structure and had total weight of $0.55(5)~{}\%$ of the diffuse correlations. We thus associate them with the same impurity phase discussed in the main text and we didn’t include data points at those peaks in Fig. 5(b), for cosmetic purposes only. Neither the position nor the weight of the magnetic impurity peaks are capable of accounting for the large peak predicted by MC simulations using best fit parameters for heat capacity data, and the presence or absence of these peaks in the NPD data do not change the conclusion of this work in any way. We estimated the potential contribution to the heat capacity due to the onset of spin order in a $0.55(5)~{}\%$ impurity phase, and incorporated this value into the error bars when determining magnetic heat capacity in the main text.

Appendix C INS fitting method and results

The resulting wavefunctions for all CEF levels are shown in Table 2. The code for fitting the INS data to the CEF model was written in MATLAB. The assumptions for the refinement are as follows: the CEF levels are thermally populated according to the partition function $Z=\sum_{n}\text{exp}\left(-\beta E_{n}\right)$, where $E_{n}$ is the energy of the $n^{\text{th}}$ CEF level. Excitation energies are determined by diagonalizing the Hamiltonian

[TABLE]

where $B_{n}^{m}$ are the crystal field parameters and $O_{n}^{m}$ are Steven’s operators. In order to isolate the unique solution where the wavefunctions are maximally parallel/antiparallel to the $\langle 111\rangle$ directions, a small field of $10^{-7}$ T along the $\langle 111\rangle$ is added to the potential; we note that this field is too small to change the energy of the calculated CEF levels. Peak intensities are given by $p_{n}\bra{\psi_{n}}J_{\alpha}\ket{\psi_{m}}^{2}$ where $p_{n}$ is the probability of an Er3+ ion being in the $\psi_{n}$ state, and $J_{\alpha}=J_{-}+J_{+}+J_{z}$. The fitting was done using a random grid search method that explores the six dimensional phase space of the six non-zero $B_{n}^{m}$ coefficients. The variation within this phase space was done by projecting random vectors in the phase space and then finding the least squares minima of the simulated pattern along those vectors. The overall minimum is then taken as the next starting point, and the process is repeated until convergence. The program was initially run with the lowest-lying CEF excitation fixed to $E=4.1$ meV, thereby limiting the search to vectors in the five dimensional manifold that satisfied this condition. After initial convergence, this condition was relaxed. To check against false minima, the program was run eight separate times, and it was confirmed to converge to the same values.

In order to find errors of the $B_{n}^{m}$ CEF parameters we use a gradient method which is considerably faster than the grid search method, although it is less robust against false minima. For each $B_{n}^{m}$ we fix the parameter’s value and minimize $\chi^{2}$ by moving along the gradient in the remaining five non-fixed CEF parameters. This is repeated for $B_{n}^{m}$ fixed at a value progressively farther from the minima until $\chi^{2}$ has increased by one, which we define as the upper and lower bounds of the error. By doing this for all of the parameters we get the error bars for each CEF parameter.

Appendix D Magnetization calculation

Magnetization curves in the main text were obtained through calculations which took into account the full CEF Hamiltonian for the $J=\frac{15}{2}$ multiplet of Er*+3*. This method allowed us to describe the moment of the material at both higher temperatures, where multiple CEF levels are occupied, and at higher fields, where mixing of the states leads to an increased moment.

In order to calculate the moment, we ignored interactions between moments and treated the problem in the single ion picture. The total Hamiltonian is thus the CEF plus the Zeeman energies

[TABLE]

The CEF Hamiltonian is defined as before in Eq. 1 and the Zeeman term $H_{Z}=-g_{J}\mathbf{H}\cdot\hat{J}$ where $g_{J}$ is the Landé g factor for Er3+. The combined Hamiltonian was diagonalized to give the energies $E_{n}$ and wavefunctions $\psi_{n}$ which were used to find the partition function $Z=\sum_{n}\exp(-\frac{E_{n}}{k_{b}T})$ and the contribution to the moment

[TABLE]

We then powder averaged by integrating over the polar angle $\theta$

[TABLE]

where $p_{n}$ is simply $\frac{\exp(-\frac{E_{n}}{k_{b}T})}{Z}$. Due to symmetry we did not need to average over the azimuthal angle, and only needed to integrate to $\frac{\pi}{2}$. This gave the full powder averaged magnetization per Er atom including the effects of the mixing of higher energy CEF doublets.

The effects of higher energy wavefunctions were found to be particularly important to the calculation of magnetization. To highlight this fact, Fig. 9 shows a comparison of the data for this material calculated with and without these effects considered.

At any fields higher than 1 T, the difference between the calculated magnetization and the magnetization from a simple two level model is considerable. In Fig. 9(b) the same calculation is repeated but with $+15.845$ meV artificially added to all of the excited energy levels, in order to bring the overall first excited energy to 20 meV. This value is much closer to the energy of the 227 rare earth pyrochlores studied, and we can see that the effect now becomes far less important; the resulting curves are quite similar to the moment calculated considering only the ground state doublet.

Appendix E The Effects of Disorder

In order to explore the effects of disorder on the material properties of MgEr2Se4, we performed a series of thermodynamic measurements on a second, less pure powder sample with a demonstrably higher level of local disorder. In Fig. 10(a) we show the XRD pattern from this sample (which we call “sample-B”), fit using a standard Rietveld refinement ($\chi^{2}=5.7$), while Fig. 10(b) presents the results of Le Bail analysis on the same dataset ($\chi^{2}=10.6$). As compared to our primary sample (henceforth referred to as “sample-A”), the XRD data on sample-B revealed a marginally higher fraction of impurities ($\approx 10\%$ total), and a Le Bail fit which improved $\chi^{2}$ considerably over refinement values. Though the peak positions are consistent with a cubic Fd$\bar{3}$m space group, the inability of the standard refinement to describe peak heights within the spinel model, even allowing for variations in site occupancy, reveals the presence of a significant level of structural disorder.

Figure 11 shows the magnetization of sample-B over a range of fields and temperatures, which are interesting to compare to measurements on sample-A presented in the main text. Although INS measurements were not performed on sample-B, we assume a similar CEF to calculate expected magnetization curves shown as solid lines in Fig. 11. Despite not having a separate INS study of the CEF levels of sample-B we find that this agreement in magnetization data shows that the CEF is not significantly modified by disorder, as would be expected for a Kramers ion.

In contrast, heat capacity is modified significantly by disorder effects, as revealed by Fig. 12. In Fig. 12(a), we show a comparison of the heat capacity of the two samples, again with the best fit MC curve for sample-A. The peak in the heat capacity for sample-B is reduced by almost a factor of two, and we have confirmed that there is no spin configuration in the next-nearest-neighbor dipole-ice model capable of reproducing this behavior. Further inspection reveals that the heat capacity of sample-B is not uniformly reduced, but rather that the peak is skewed to higher temperature. This leads to a long high-temperature tail wherein the curve for sample-B lies above sample-A. Integrating the area under these curves leads to the associated entropy curves in Fig. 12(b). Quite surprisingly, we find that sample-B recovers full Pauling residual entropy as $T\rightarrow 0$ K, in direct contrast to the conclusions of the main text for sample-A. This suggests that disorder acts to hinder the mechanisms which are leading to the reduction of residual entropy in relatively pure samples of MgEr2Se4.