Introduction

Frustrated magnetism is a fruitful frontier in the contemporary field of quantum materials wherein mutually incompatible interactions, driven by the arrangement of the magnetic ions and their exchange couplings, can lead to a variety of unusual phenomena1,2. Due to the dependence of the frustrated interactions on lattice geometry, materials with under-studied lattice structures may offer untapped opportunities for the exploration of interesting and exotic physics.

A long-serving example of the rich physics afforded by frustrated systems is illustrated by the A\({}_{2}^{3+}\)B\({}_{2}^{4+}\)O7 magnetic pyrochlore oxides. In these compounds, A is a rare-earth ion (or Y3+) and B is typically a transition metal, with both ions potentially carrying a magnetic moment and residing on the vertices of two distinct but interpenetrating lattices of corner-sharing tetrahedra)3,4,5. For example, these materials have attracted interest for hosting classical spin ice physics6,7 in Ho2Ti2O76,8,9 and Dy2Ti2O710, and potentially their quantum analogs5,11, unconventional ordered states12,13,14, order-by-disorder as in Er2Ti2O715,16,17,18 and possibly other, more exotic, quantum spin liquids19,20,21.

Over the past few years, a recently reported related class of systems, the “breathing” pyrochlores, has garnered excitement. In these, the alternating corner-sharing tetrahedra of the lattice are distinguished by two different bond lengths, breaking inversion symmetry and leading to competition between inter- and intra-tetrahedron exchange interactions. The difference in size between the two types of tetrahedra can be characterized by the “breathing ratio”, defined as \(d^{\prime} /d\) where d and \(d^{\prime}\) are the nearest-neighbor bond distances between magnetic ions in small and large tetrahedra, respectively. Unfortunately, studies of breathing pyrochlores have been limited to a small number of systems due to the scarcity of successfully synthesized compounds, with single-crystal samples being even more elusive.

One family of breathing pyrochlore materials that have been studied in detail is the Li(In,Ga)Cr4(O,S,Se)8 spinels22,23,24,25,26. These compounds have breathing ratios of \(d^{\prime} /d\approx 1\) and are thus in a regime of effective (competing) interactions close to the regular, non-breathing, pyrochlore lattice. LiInCr4S8, LiGaCr4S8, and CuInCr4S8 (initially characterized in the 1970s27) are few members of this family which have been synthesized successfully. While these materials have provided interesting opportunities to study inter-tetrahedron coupling effects for the case of \(d^{\prime} /d\approx 1\) where the up and down tetrahedra have similar sizes, they leave unanswered questions about what may be possible with a more exaggerated breathing ratio. The large S = 3/2 spin for the Cr3+ ions in these materials may also suppress quantum effects28 while an S = 1/2 breathing pyrochlore magnet is less likely to be classical.

Study of the Li(In,Ga)Cr4(O,S,Se)8 has also been hindered by the lack of single-crystal samples, which complicates the determination of the magnetic Hamiltonians of these systems28. Ba3Yb2Zn5O11, however, lies at the opposite extreme where the breathing ratio is large, \(d^{\prime} /d\approx 2\)29. This large breathing ratio leads to very different sizes for the up and down tetrahedra, resulting in rather isolated small tetrahedra with weak inter-tetrahedron interactions. Characterization of Ba3Yb2Zn5O1129,30,31,32,33 has confirmed this picture of an interesting kind of “molecular magnet” of largely, it appears, independent tetrahedral spin clusters, each exhibiting a pair of degenerate ground states. The presence of these emerging low energy doublet degrees of freedom necessitates the consideration of inter-tetrahedron interactions, even in limit of large \(d^{\prime} /d\) breathing ratio. For example, how this degeneracy is resolved by the collective physics of the tetrahedra and the nature of the ultimate ground state of Ba3Yb2Zn5O11 remain open questions.

Turning now our focus to Ba3Yb2Zn5O11, we begin with the single-ion physics. The magnetism of Ba3Yb2Zn5O11 originates from the Yb3+ ions where the J = 7/2 manifold of the isolated ion is split by the C3v symmetric crystal field. This crystal field leads to a Kramers doublet separated by a gap of 38.2 meV from the remaining levels32. Given this splitting is large relative to the scale of the ion-ion interactions34, a description in terms of an effective S = 1/2 degree of freedom for each Yb3+ ion is possible5,34. Analysis of inelastic neutron scattering data on powder samples found that antiferromagnetic Heisenberg exchange is dominant29 with a surprisingly strong, but sub-dominant, Dzyaloshinskii-Moriya interaction30,31,33 within each small tetrahedron. The absence of magnetic ordering down to ~0.1 K31,35, despite a Curie-Weiss temperature of −6.1 K (as measured in the present paper), the observation of a residual entropy of \({k}_{{{{\rm{B}}}}}\ln (2)/4\) per Yb29 and a dispersionless excitation spectrum30,31,33 is consistent with these smaller tetrahedra remaining mostly decoupled (i.e., paramagnetically) down to T ~ 0.1 K.

The unexpected observation of a robust Heisenberg plus Dzyaloshinskii–Moriya interactions, with negligible symmetric anisotropies between the effective S = 1/2 Yb3+ spins in Ba3Yb2Zn5O1134, has led to insights in understanding related ytterbium compounds such as the chalcogenide spinels36,37,38 AYb2X4 (with A = Cd, Mg and X = S, Se) and the nearly perfect S = 1/2 honeycomb Heisenberg antiferromagnet YbCl339,40,41.

Proximity to this regime of weak symmetric exchange anisotropy in breathing pyrochlores with a large breathing ratio (\(d^{\prime} /d\gg 1\)) has also been predicted to lead to rich collective physics when these tetrahedra are weakly coupled. Theoretical proposals include realizations of bosonic analogs of topological band structures such as Weyl touching points42, enhanced stability of quantum spin ice11,43, and possibly the emergence of a rank-2 U(1) spin liquid44. Realization of these kinds of interesting collective physics of the small tetrahedra requires inter-tetrahedron interactions to be at play, and thus their presence and magnitude to be exposed and characterized experimentally.

For the bulk of the experimental data on Ba3Yb2Zn5O11 thus far collected, these have been found to be too small to be observable, with a model of isolated tetrahedra proving sufficient to reproduce quantitatively all of the experimental data above ~1 K30,31,33,35. However, notable discrepancies appear in the heat capacity30,31,35 below ~0.3 K, suggesting the single-tetrahedron model is incomplete at very low energies. Inter-tetrahedron interactions could provide a natural explanation30, though the effects of structural disorder31 must also be considered carefully, given the non-Kramers E doublet ground state of each small tetrahedron30,35 can be split by the non-magnetic disorder. Resolving the origin of these discrepancies through the investigation of powder samples is difficult, not only due to the small energy scales involved, but also due to the directional averaging inherent for powders35. In particular, in neutron scattering, this averaging can render dispersive excitations difficult to resolve beyond a linewidth. Thus, to explore carefully the low-temperature physics of Ba3Yb2Zn5O11 and understand the nature of any inter-tetrahedron interactions, single-crystal studies are indeed required. Moreover, single-crystal studies allow access to the details of field-driven level crossings even at the single-tetrahedron level which is unresolvable using powder samples.

In this work, we report the synthesis and characterization results obtained for the single-crystal samples of Ba3Yb2Zn5O11. We present comprehensive studies of the single-crystal Ba3Yb2Zn5O11 using a wide range of experimental techniques including inelastic neutron scattering (INS) and bulk magnetic susceptibility, as well as ultra-sensitive tunnel diode oscillator (TDO) magnetic susceptibility measurements. Our single-crystal measurements allow us to directly probe the finite wavevector dependence of the magnetic excitations and anisotropic response at low temperatures and high magnetic fields. We find that the single tetrahedron model maintains its qualitative agreement with inelastic neutron measurements for most magnetic fields. On the other hand, our diffuse neutron scattering data collected at 70 mK shows evidence of deviations from the single tetrahedron model near the critical magnetic field attributed to the first level crossing and falls short of explaining the very low energy behavior. All of these results suggest that the inter-tetrahedron interactions could be responsible for the low-energy and low-temperature physics, as depicted, for example, by the diffuse neutron scattering above the critical field of the first level crossing in this system. These results are also in agreement with the previously reported low-temperature specific heat measurements where the simple single tetrahedron model fails to reproduce the observed data at low temperatures35. This work paves the way for further theoretical and experimental investigations aimed at shedding light on the unsolved ground state of this interesting breathing pyrochlore system.

Results

Magnetic property measurements

The magnetic susceptibility, χ, was measured on our single-crystal sample down to 1.8 K under an applied magnetic field μ0H of 0.01 T along the crystallographic [111] direction using the aforementioned Cryogenic S700X SQUID system, as shown in Fig. 1c. A broad maximum is observed at ~4 K, similar to that in the powder sample29. The Curie-Weiss temperature θCW = −6.85 K obtained from the Ba3Yb2Zn5O11 single-crystal susceptibility data by fitting χ = C/(T − θCW) in the temperature range of 10 K < T < 30 K, agrees well with the previously reported value θCW = −6.7 K using a powder sample29.

Fig. 1: Structural details and characterization of the single crystal.
figure 1

a The crystal structure of Ba3Yb2Zn5O11. b Single-crystal sample of Ba3Yb2Zn5O11 used for single-crystal neutron scattering experiment and the Laue back-scattering diffraction pattern of plane perpendicular to the [111] direction. c Magnetic susceptibility χ vs temperature, T, down to 1.8 K. (Upper Inset) Field-dependent magnetization M(H) at different temperatures T = 1.8 K (blue curve) and T = 2.5 K (red curve), with μ0H sweep rate of 0.03 T/min. The inflexion point in the M vs H data corresponds to a transition from an approximate total spin Stot = 0 state to an approximate Stot = 1 state. (Lower Inset) First derivative of magnetization M(H) as function of field H at temperatures of T = 1.8 K (blue curve) and T = 2.5 K (red curve). d The Rietveld refinement of the neutron powder diffraction pattern of Ba3Yb2Zn5O11 at 2 K collected using NOMAD at the Spallation Neutron Source of ORNL. The red and blue solid lines represent the calculated intensity and the difference between observed and calculated intensities (weighted profile R-factor51, Rwp = 4.79%), while the short vertical green marks represent the expected Bragg peak positions. e The pair density function (PDF) data at 2 K along with the refinement is shown. Black open circles, red lines, and blue lines represent the PDF data, the refinement using the PDFGUI software50 with assumed space group \(F\bar{4}3m\) and the difference between PDF data and refinement, respectively.

Applying larger magnetic fields, Ba3Yb2Zn5O11 can be tuned through two transitions. In powder samples, these appear near the critical fields μ0Hc1 ≈ 3.5T and μ0Hc2 ≈ 8.8 T31,33. We show in Fig. 2a the energy level results calculated using the single tetrahedron model as a function of field. The theoretical values of μ0Hc1 and μ0Hc2 are indicated by vertical dash lines. As alluded to earlier, the TDO frequency f is an extremely sensitive measurement of the magnetic susceptibility45,46,47 and reveals fine details of changes in the magnetic response of the sample. We show in Fig. 2b–d the measured TDO frequency in Ba3Yb2Zn5O11 over a range of magnetic fields along the [111] direction, with strength 0 T ≤ μ0H 12 T. Our SQUID and TDO results collected on the single-crystal reveal that the critical fields for μ0H[111] [Fig. 2c, d] are nearly identical to those reported previously on powder samples33.

Fig. 2: Energy levels and field-induced transitions.
figure 2

a Field-dependent energy levels calculated by the single-tetrahedron model. Irreducible representations of the zero field (Td) and finite field (C3) symmetry groups are indicated. Critical fields where there is a ground state level crossing (μ0Hc1 = 3.65 T and μ0Hc2 = 9.48 T) indicated by vertical dash lines are close to the experimental values. The low energy regime (e.g., as would be thermally populated at T 1 K) is indicated by a gray shaded box. b Tunnel diode oscillator (TDO) frequency (f) as a function of magnetic field H[111] up to 20 K. Traces are shifted vertically for clarity. Small dips in frequency in addition to two level crossings are marked by arrows. The nature of these small anomalies calls for further investigation. The sharp features at higher temperatures, e.g., near ~6 T, are associated with background signal from our experimental apparatus. c, d f vs H[111] as a function of magnetic field H[111] up to 0.880 K, near the ground state level crossings at c 3.5 T, and at d 8.8 T. Additional features observed ~3.5 T are marked by gray dash lines. Insets in c and d show f vs H[111] for the upsweep (black) and downsweep (red) at T = 0.041 K, where no hysteresis is seen.

As shown in Fig. 2c, the TDO data reveals the presence of two shoulder-like structures for field values just below the first level crossing at 3.5 T and at the lowest temperature, and become more distinguishable as temperature is increased. These features may suggest the presence of additional level crossings near μ0Hc1 that are not captured by the single-tetrahedron model. We leave a detailed exploration of the nature of these features to future work. Similar to our SQUID data, we do not see obvious hysteresis either near μ0Hc1 or μ0Hc2. This is in contrast to ref. 33, where a field-sweep-rate dependent hysteresis was reported on powder samples of Ba3Yb2Zn5O11 and attributed to an avoided level crossing of the ground and first excited state near μ0Hc1.

However, in a [111] field, the crossing at μ0Hc1 is not avoided; symmetry forbids any coupling between these states [see Fig. 2a]. We note that in our data no hysteresis is observed [see insets of Fig. 2c, d], in contrast to that seen in ref. 33. Note that depending on the magnitude of the field sweep rate relative to the energy splitting at the avoid crossing, non-adiabatic transitions between the ground and first excited states are possible and their probability can be estimated using the usual Landau-Zener formula. Significant non-adiabatic probability implies hysteresis in the state of the system (and thus magnetization) as the field is swept up or down. The authors of ref. 33 attribute their observed hysteresis to a non-equilibrium effect due to the finite rate of change in the magnetic field strength. However, estimates of the probability of non-adiabatic transitions (using the standard Laudau–Zener formula48) for any reasonable splitting size [O(meV)] and field sweep rate [\(O({{{\rm{mT/min}}}}]\)] are negligible. With increasing temperature to 20 K [Fig. 2b], the sharp features associated with the two ground state level crossings at μ0Hc1 and μ0Hc2 get broader and diminished, while additional features as indicated by arrows, appear as the temperature is raised.

From measurements on our single crystal, we have exposed signatures of two field-driven transitions, as expected from the theoretical single-tetrahedron model. We next consider how the spin-spin correlations, as probed by neutron scattering measurements on our single crystal, evolve as a function of field and, in particular, across the two transitions at μ0Hc1 and μ0Hc2.

In the following sections, we present the results of three complementary neutron scattering experiments. In order, these address the presence or absence of structural disorder (Sec. “Structural characterization”), the finite-energy magnetic excitations (Sec. “Energy-resolved inelastic scattering”) as well as an indirect measurement of the low-energy excitations via an energy-integrated experiment (Sec. “Energy-integrated diffuse scattering”).

Structural characterization

One of the main concerns that arose in previous studies on Ba3Yb2Zn5O1130,31,35 was whether or not any disorder played a role in the discrepancies observed in their reported low-temperature measurements relative to the theoretical results on the above single tetrahedron model, especially in the heat capacity data31,35. In an attempt to address this concern and investigate the detailed crystal structure at low temperature and any possible disorder, neutron scattering measurements were performed on Ba3Yb2Zn5O11 powder sample using the Nanoscale-Ordered Materials Diffractometer (NOMAD) at Oak Ridge National Laboratory. The data was collected for selected temperatures between 100 K to 2 K. Rietveld refinements were performed using the GSAS II software on the collected data. No structural transition is observed down to 2 K and the system can be well described with the \(F\bar{4}3m\) space group, as shown in Fig. 1d. To understand the local atomic structure of Ba3Yb2Zn5O11, the pair distribution function (PDF) of the total scattering data was also analyzed49. The PDF data, G(r), collected at 2 K is shown in Fig. 1e. Data were refined using the PDFGUI software50 with the \(F\bar{4}3m\) space group. As evident from the refinement, the observed data are well represented by our model with \(F\bar{4}3m\) symmetry with weighted R-factor51 Rw = 5%. No peak splitting or significant mismatch is resolvable; our measurements thus do not provide any evidence for local structural disorder in this system down to 2 K. To further strengthen this statement, we also performed high quality single-crystal x-ray diffraction measurements using our single-crystal sample. The details of these measurements and the results are summarized in the Supplementary Information. The refinement converged to a parameter set with residual factors R1 = 0.0107 (F2, I > 2σ(I)) and wR2 = 0.0205 (F2, all reflections)52. The occupancy factors for each atomic site was refined and found to be ~1 within the error bars, providing further evidence for the absence of site mixing in the Ba3Yb2Zn5O11 single-crystal samples. (Refer to Supplementary Table 1 and Supplementary Table 2 in the Supplemental Information for details on the refinement results).

Energy-resolved inelastic scattering

To date, only powder samples of Ba3Yb2Zn5O11 have been measured with inelastic neutron scattering30,31,35. Within the constraints of powder averaged inelastic neutron scattering data, the experimentally observed low energy excitations have been found to be in agreement with the single tetrahedron model30,31,33. In order to study the detailed physics of this system beyond the single tetrahedron model, such as the details of the energy levels crossing depending on the direction of applied field, single-crystal neutron scattering experiments performed at zero and finite fields (applied along different lattice directions) are necessary. Single-crystal samples also allow us to investigate in detail the anisotropic behavior of the excitations compared to previous powder neutron scattering studies.

To this end, we have performed inelastic neutron scattering experiments using our single-crystal sample of Ba3Yb2Zn5O11 as explained in Sec. “Neutron scattering” with applied field H[111] (at DCS) and \({{{\bf{H}}}}\parallel [1\bar{1}0]\) (at CNCS). The field evolution of the low energy excitations for a few selected fields between μ0Hc1 and μ0Hc2 are shown in Fig. 3 with applied field H[111]. Figure 3 shows the wavevector dependence for each energy band as a function of field, in which each cut was made over a narrow energy window.

Fig. 3: Inelastic neutron scattering results.
figure 3

Inelastic neutron scattering from a single-crystal sample of Ba3Yb2Zn5O11 measured at T = 70 mK using the Disk Chopper Spectrometer (DCS) for a series of fields H[111] a 0 T b 5 T, and c 10 T. For each field, the left panel shows the dynamical structure factor as a function of energy transfer and momentum along [h, −h, 0] direction. An integration over the full K range in the [k, k, −2k] direction was performed, i.e., covering the Q range from −2 to 2 along that direction. On the right, constant energy intensity maps in the [h+k, −h+k, −2k] plane are shown for selected dispersionless bands which are integrated along energy with an energy width specified to fully capture each band for each magnetic field.

Similar to powder samples30,31,33, dispersionless (flat) excitations are observed at zero and finite fields in our single-crystal samples as well. These flat bands are single-tetrahedron excitations (discrete energy levels) arising from the interactions between the effective spin-1/2 states of the four Yb3+ ions on each of the individual (small) tetrahedra. A slight broadening of the excitations greater than the resolution of the instruments is observed, although there is no clear dispersion visible. We note that some intrinsic broadening of the levels is to be expected due to any dispersion-induced by inter-tetrahedron interactions; the slightness of this broadening and inability to resolve any dispersion of those excitations is consistent with the smallness of these couplings.

Having access to single-crystal samples, we can explore the wavevector dependence of the intensity of each dispersionless mode in three-dimensional reciprocal space and check for the theoretically expected symmetry of the neutron scattering intensity in the scattering plane in zero and finite fields. In zero field, the presence of full Td symmetry and time-reversal symmetry imposes a sixfold symmetry in the [h + k, −h + k, −2k] plane. Thus, each band exhibits a pinwheel-type pattern which is sixfold symmetric in this scattering plane, as shown in the three right panels of Fig. 3a. See also Fig. 4a, upper two panels for energy windows [0.45, 0.60] meV and [0.65, 0.85] meV vs theoretical calculations (lower two panels of Fig. 4a for corresponding energy windows).

Fig. 4: Dynamical structure factor.
figure 4

Comparison of the experimental (top row of each panel) and theoretical (bottom row of each panel) dynamical structure factor for the model of Eq. (1) with the parameters of Eq. (4) for several energy windows at different magnetic fields, including the atomic form factor of Yb3+. The energy windows are convoluted with a finite Gaussian to emulate experimental broadening and resolution. In the following, experimental plots symmetry operations were applied to improve the statistics. H[111] at 0 T: sixfold rotation and H > 0 T: threefold rotation. \({{{\bf{H}}}}\parallel [1\bar{1}0]\) at 0 T: mirror symmetry operations centered at [h,−h, 0] and [0, 0, k] directions and H > 0 T mirror symmetry operations centered at [0, 0, k]. H[111] data were measured at DCS with incident neutron energy Ei = 5 meV and only the detector background subtracted (i.e., not including the high-temperature background). \({{{\bf{H}}}}\parallel [1\bar{1}0]\) data was measured at CNCS and 100 K data were used as the background to eliminate unavoidable nearly temperature-independent background scattering observed due to scattering from the magnet.

Upon application of a [111] field, we observe an evolution of the intensity as a function of wavevector which changes qualitatively across each critical field [refer to Fig. 2a]. The details of the redistribution of intensity depend on the field direction: for [111] fields, only a threefold symmetry remains, while for fields along \([1\bar{1}0]\), only a single twofold symmetry is present. The zero field Hamiltonian \({{{{\mathcal{H}}}}}_{{{{\rm{eff}}}}}\) in Eq. (3) has a mirror symmetry, σ = −C2, that is broken by finite field. However, this mirror symmetry operation reverses a [111] magnetic field and therefore, combining it with time-reversal, yields an additional symmetry

$$I({{{\bf{Q}}}},E;{{{\bf{B}}}})=I({{{{\bf{C}}}}}_{2}{{{\bf{Q}}}},E;{{{\bf{B}}}}),$$
(1)

for the [111] field direction.

We thus effectively have a set of twofold axes in the plane perpendicular to [111], giving a symmetry group of D3 rather than only a C3 symmetry in zero field. Under a finite magnetic field, this effective inversion symmetry [Eq. (1)] is broken and the intensity pattern becomes threefold symmetric. This can be seen in Fig. 3b, c, which shows the E-Q plots and the energy-integrated cuts for each energy band for a [111] field (5 T) above the first level crossing transition at Hc1 and one (10 T) beyond the second level crossing transition at Hc2.

Figure 3b shows the results obtained for a 5 T applied field, where energy cuts obtained for the band centered at E = 0.25 meV reveal broad features centered around \([20\bar{2}]\) that have the expected threefold symmetry in the [h + k, −h + k, −2k] plane. It is clear that the obtained pattern for 5 T is different from the results collected in zero field. For the next higher excitation band at E = 0.65 meV, the pattern is rotated by 30 degrees compared to the band at E = 0.25 meV. Refer also to Fig. 4c (upper two panels) for energy windows [0.20, 0.40] meV and [0.55, 0.70] meV as well as theoretical calculations (lower two panels of Fig. 4c for corresponding energy windows). Similarly, the intensity pattern in Q space in this scattering plane changes as we move up in energy with each flat excitation showing a distinct pattern.

Similar to the 5 T results, the 10 T field data also demonstrates a subsequent change in the flat band intensity scattering following the second level crossing at Hc2. For the energy transfer range determined by our incident energy, we could only observe two flat excitation bands at 10 T. The band centered at E = 0.25 meV shows a similar threefold pattern for this field as observed for 5 T. In contrast, the second band at E = 1.4 meV shows an intensity that is more concentrated at small wavevector; likely a transition between approximate Stot = 2 states of mostly ferromagnetic character. More explicitly, referring to Fig. 2a, one can roughly track the approximate Stot = 2, \({S}_{{{{\rm{tot}}}}}^{z}=+2\) level (emerging from the E + T2 manifold at μ0H = 0) to the ground state at 10 T (through several avoided level crossings), as well as the Stot = 2, \({S}_{{{{\rm{tot}}}}}^{z}=+1\) state to the level at 1.4meV.

In order to explore the anisotropy of the neutron scattering intensity as a function of the field direction, we also performed inelastic neutron scattering experiments with \({{{\bf{H}}}}\parallel [1\bar{1}0]\), perpendicular to the [111] direction. As with the [111] direction, we consider individual energy cuts at several fields ranging from 0 T to 10 T, comparing the applied field along \([1\bar{1}0]\) direction to the analogous cut for [111] field. The results are shown in Fig. 4. For each panel, the experimental data and theoretical calculation using single tetrahedron model are shown for both field directions, H[111] and \({{{\bf{H}}}}\parallel [1\bar{1}0]\). We emphasize that Fig. 4 shows good agreement between the experimental data and the single tetrahedron model calculations, for all energy cuts and for both field directions, and within the constraints from experimental noise, background and resolution. All qualitative features for each energy cut are well captured by the theoretical calculations. This indicates that the nature of the dispersionless bands at energies E > 0.3 meV are well explained by the single tetrahedron model.

To further assess the agreement of the inelastic scattering data against the single tetrahedron model, we show in the Supplementary Information the constant energy cuts in Fig. S3 for H[111] and Fig. S4 for \({{{\bf{H}}}}\parallel [1\bar{1}0]\) along [k, k, −2k] and [0, 0, k], respectively. For both field directions, the experimental data agrees well with the calculation using a single tetrahedron model.

The consistency of the single-crystal data and the theoretical model fit to data from polycrystalline samples is non-trivial, given the different sample growth processes (powder vs. single crystal). Further, given the information lost due to powder averaging, one may worry that the theoretical model may not have been fully constrained by the polycrystalline data30,35. The agreement reported here is evidence that the single-crystal and powder samples are exhibiting the same low-energy physics, and the model determined from powder samples remains valid, at least for energies O (0.1 meV).

Energy-integrated diffuse scattering

While, the inelastic scattering data reported in Sec. “Energy-resolved inelastic scattering” provides detailed information about the finite-energy excitations, due to the challenges arising from experimental background, instrument resolution and the limited flux, it offers limited information about the physics of the system at low energies, where one might expect to see evidence of physics beyond the single tetrahedron model. Therefore, in order to capture this very low-energy physics, a diffuse neutron scattering experiment, which integrates over all energies, was performed using the Wide Angle Neutron Diffractometer (WAND2) instrument at Oak Ridge National Laboratory. WAND2 features a highly efficient high resolution 3He 2-dimensional position sensitive detector which enables it to map a large portion of reciprocal space for single crystals. One of the specialized data collection purposes of WAND2 is the measurement of diffuse scattering in single crystals. Figure 5 shows the neutron scattering data measured at WAND2. Data was collected at 70 mK and fields up to 4.8 T. High-temperature data collected at 50 K and zero field, was used as background and subtracted from the 70 mK data for each field. The incident neutron wavelength at WAND2 is λ = 1.486 Å, with corresponding incident energy Ei = 37.04 meV. Due to the inability to discriminate the energy and considering the energy resolution (>3 meV), one can expect the data obtained at WAND2 to amount to energy-integrated scattering up to ~2 meV. The diffuse scattering intensity (static structure factor) is given by integrating Eq. (7) over energy

$$I({{{\bf{Q}}}},B)\propto F{({{{\bf{Q}}}})}^{2}\mathop{\sum}\limits_{\mu \nu }({\delta }_{\mu \nu }-{\hat{Q}}_{\mu }{\hat{Q}}_{\nu })\left\langle {\mu }_{-Q}^{\mu }{\mu }_{Q}^{\nu }\right\rangle$$
(2)

where, again, F(Q) is the atomic form factor of Yb3+53. Per its definition, this is symmetric under Q → −Q, even for nonzero magnetic fields and unpolarized neutron scattering intensity. Combining this with the threefold rotational symmetry, the integrated intensity has an enlarged sixfold symmetry for a [111] field. We note that this additional symmetry arises due to the symmetrization of the unpolarized intensity and the integration of over-energy transfer that reduces the structure factor to a static moment-moment correlation function.

Fig. 5: Diffuse neutron scattering results.
figure 5

Comparison of the experimental (right) and theoretical (left) energy-integrated structure factor at several magnetic fields H[111] (0, 1 T, 3.5 T, and 4.8 T) using the single tetrahedron model. Theoretical plots include the atomic form factor of Yb3+. Experimental data were measured using the Wide Angle Neutron Diffractometer (WAND2). Both experimental and theoretical plots show the data at T = 70 mK after subtracting the background at 50 K at 0 T.

In zero field, the qualitative features of the data are reproduced in the theoretical calculation: a hexagon of high-intensity with maximum intensity along [h, −h, 0] and equivalent directions. In addition, it shows minima in the intensity at zero, as well as at six points equivalent to ~[6, −6, 0]. Six weaker minima are also visible near [2, 2, −4] (and equivalent positions) in both the theoretical calculation and the experimental data. However, the agreement is not perfect; these six minima are more pronounced in the theoretical calculation as compared to the experimental data. The pattern for μ0H = 1 T is similar to the zero field one in both the theoretical calculation and the experimental data. In contrast, the calculated intensity at μ0H = 3.5 T is considerably different from the experimental data. Namely, the high-intensity hexagon seen experimentally is broader than the “spokes”-like pattern of the theoretical one, with an overall more uniform intensity. Further, the minimum near [2, 2, −4] and its equivalents seen in the theoretical calculation are essentially absent in the experimental data. Moreover, the theoretical results for 3.5 T show very similar pattern to the 0 T and 1 T data, with the same set of pronounced minima near [2, 2, −4] and equivalent reciprocal points. At 4.8 T, the theoretical results differ markedly from the ones obtained for the three lower fields (compare the four panels in the left column of Fig. 5). While some of the features observed experimentally in the 4.8 T data, such as the location of the intensity minima near (e.g.) [6, −6, 0], agree with the theoretical results, there are important qualitative differences. These include significantly less intensity modulation in the high-intensity central hexagon in the theoretical calculations, compared to the experimental data, as well as the depth of the intensity minima near [6, −6, 0] and equivalent positions. In comparison, the DCS and CNCS inelastic neutron scattering results – that do not include the low-energy contributions, as discussed in Sec. “Energy-resolved inelastic scattering”, are in overall acceptable agreement with the theoretical calculations using the single tetrahedron model (Fig. 4).

Recalling the comments made at the beginning of this Section “Energy-integrated diffuse scattering” regarding the ability of the WAND2 experiment to probe low-energy excitations, which closer to elastic channel in the inelastic experiments, we interpret the results of Fig. 5 in the following way: as the inter-tetrahedron interactions are expected to be small compared the intra-tetrahedron ones, we can argue that the possible dispersion of the ground state doublet induced by these interactions is also small, with the corresponding excitations lying at very low energies, close to the elastic channel. The differences between the data and the theoretical calculations using the single tetrahedron model may thus be due to low energy diffuse scattering features arising from physics beyond the single tetrahedron model. Such discrepancies at low energy are consistent with the previously reported heat capacity data31,35, which showed a broad peak centered at ~0.1 K, with an entropy change consistent with the release of two degrees of freedom per single tetrahedron. In the single-tetrahedron model, the ground state doublet is degenerate; the complete release of this entropy, indicating a distinctive ground state, thus points again to physics beyond the single-tetrahedron model. The energy scale of this splitting was estimated to be ~0.015 meV at 0 T31—much smaller than the instrumental resolution of our neutron scattering experiments, ~0.05 meV, as in the experiments presented in Sec. “Energy-resolved inelastic scattering”. We indicate this energy scale by a shaded gray box in Fig. 2a. Compared to the single-tetrahedron energy levels, we see that this energy scale is indeed small.

The most natural origin for the low-energy splitting of the doublets may be attributed to inter-tetrahedron interactions or to non-magnetic disorder that lowers the symmetry of each tetrahedral cluster. Recalling the discussion in Sec. “Structural characterization”, our PDF measurements and single-crystal x-ray diffraction results show no observable evidence for disorder in our samples, indicating any non-magnetic structural disorder is small and below the resolution threshold of our experiments. However, given the minute energy scale of the doublet splitting (~0.015 meV), we cannot rule out that correspondingly minute distortions, below the experimental resolution, could potentially account for the entropy release with the disorder splitting the non-Kramers E ground doublet of isolated small tetrahedra.

At higher magnetic fields, away from the level crossings, the ground state is non-degenerate and thus the effect of either small inter-tetrahedron interactions or non-magnetic disorder should be less significant. However, we note that there still exist discrepancies, in the diffuse scattering data from WAND2 as well as in the TDO measurements, with the single tetrahedron model at these higher fields. A similar disagreement was noted in the high-field heat capacity35. The effective strength of the inter-tetrahedron interactions can vary between the different degenerate energy level manifolds of a single tetrahedron and as a function of field as the levels mix. It is therefore not clear whether higher fields will enhance these interactions sufficiently to account for the discrepancies with the data.

Discussion

We reported here the single-crystal synthesis of a breathing pyrochlore. Single-crystal neutron scattering data using high-quality single crystals of the breathing pyrochlore Ba3Yb2Zn5O11 was also presented. Our neutron scattering data show essentially dispersionless bands within instrumental resolution. The single tetrahedron model can overall explain the observed dependence on wavevector for the individual excitation bands of the system. Discrepancies between experimental data and the single tetrahedron model are observed in the diffuse scattering data at higher fields (μ0H > 3T) and T = 0.1 K, indicating that physics not included in the single tetrahedron model is responsible. Our powder diffraction data and pair distribution function analysis show that there is no observable disorder in the samples, which further support the suggestion that inter-tetrahedron interactions are the most likely explanation for the discrepancies between the experimental data and the predictions of the single tetrahedron model. Tuning the strength of the inter-tetrahedron interactions in breathing pyrochlores away from the isolated limit where \(d^{\prime} /d\gg 1\), by changing the breathing ratio, may be a path to realizing exotic states in breathing pyrochlores generally. The array of possibilities for such tuning, e.g. chemical or hydrostatic pressure54,55, or by synthesis of alternative members of this family with different rare-earth ions (other than Yb3+), presents many avenues for future investigations.

Methods

Sample Synthesis

A powder sample of Ba3Yb2Zn5O11 [structure shown in Fig. 1a] was synthesized, as reported earlier29, by finely grinding mixed powders of Yb2O3 (99.99%, Alpha Aesar), ZnO (99.99%, Sigma Aldrich), and BaCO3 (99.95%, Alpha Aesar) in the stoichiometric ratio of 1:5:3 and reacting at 1140 C in a box furnace for 48 hours with several intermediate grindings at a lower temperature. In preparation for the single-crystal growth, the powder was compressed hydrostatically into a cylindrical rod, and then sintered at 1150 C in a vertical Bridgman furnace. Finally, large single crystals of Ba3Yb2Zn5O11 were grown using the optical floating zone method. A typical crystal was grown in an O2 atmosphere at 0.9 MPa, with an initial growth speed of 10 mm per hour and, upon stabilization of the liquid zone, at ~6 mm per hour until finished. To ensure uniform homogeneity of the liquid zone, the feed and seed rods were rotated in opposite directions with a rotation speed of 20 rpm during growth.

Susceptibility and magnetization measurements

Magnetic susceptibility, χ, and magnetization, M, measurements were performed on a single-crystal sample of Ba3Yb2Zn5O11 up to μ0H = 7 T, using an in-house Cryogenic S700X SQUID magnetometer at temperatures of 1.8 K to 200 K using a Helium-4 probe. To explore the high-field magnetic properties of the Ba3Yb2Zn5O11 crystals, TDO measurements (described below)56 were carried out at the DC Field Facility of the National High Magnetic Field Laboratory in Tallahassee. A dilution refrigerator and a 3He system were used to cover the temperature range from 41 mK to 20 K. The field-dependent measurements up to 18 T were conducted using a superconducting magnet. The field sweep rate was kept low (0.1 to 0.3 T/min) to minimize the magnetocaloric effect. For a typical TDO measurement, bar-shaped Ba3Yb2Zn5O11 crystals of ~2 mm in length and ~1 mm in transverse width were prepared and placed inside a detection coil, with the [111] direction aligned along the coil axis. Together, the coil and sample within form the inductive component of a LC circuit. The LC circuit, powered by a tunnel diode operating in its negative resistance region, was tuned to resonance in a frequency range between 10 and 50 MHz. The shift in the resonance frequency f, which is related to the change in the sample magnetization M (df/dHd2M/dH2)46, was then recorded. The TDO method measures the resonance frequency to a very high precision56 which enables identifying changes in the magnetic moments down to ~10−12 emu compared to ~10−8 emu in SQUID magnetometry measurements.

Neutron scattering

INS measurements were performed at a number of facilities: using the Disk Chopper Spectrometer (DCS) at the National Institute of Standards and Technology, the Cold Neutron Chopper Spectrometer (CNCS) at Oak Ridge National Laboratory and the Wide Angle Neutron Diffractometer (WAND2) also at Oak Ridge National Laboratory. A superconducting magnet was used in each of the above instrument to provide, in all cases, a vertical magnetic field up to 10 T (DCS), 8 T (CNCS), and 5 T (WAND2). In all experiments, a single-crystal sample of mass of ~1g was mounted on a Cu sample holder in a dilution refrigerator. At DCS and WAND2, the sample was mounted with the [h + k, −h + k, −2k] scattering plane being horizontal and with the magnetic field applied along the [111] cubic direction. Neutron-absorbing Cd was used to shield the sample holder to reduce background scattering. The experiment at CNCS was performed with the sample mounted so the horizontal scattering plane is [hhl] with the magnetic field applied vertically, along the \([1\bar{1}0]\) cubic direction. In all experiments, low-temperature measurements were conducted with applied fields at the 70 mK base temperature and the background determination measurements were conducted under zero magnetic field at 50 K, 100 K, and 50 K for DCS, CNCS and WAND2, respectively. The sample stage was rotated close to 180 to cover few Brillouin zones in all the experiments. After carefully investigating the as-collected (unsymmetrized data), we used the crystallographically-allowed symmetry operations to improve statistics. The specific symmetry operations applied for each figure reported above in Sec. “Energy-resolved inelastic scattering” are noted in the pertinent figure captions. In order to investigate any structural disorder in Ba3Yb2Zn5O11, total neutron scattering measurements were performed on a polycrystalline sample using the NOMAD at the Spallation Neutron Source (SNS) at Oak Ridge National Laboratory (ORNL, USA) with a lowest temperature of 2 K considered. Analysis and visualization of the neutron scattering data were performed using DAVE MSlice57, Mantid58, and Python software packages.

Single-tetrahedron model

Given the large inter-tetrahedron distance of \(d^{\prime} \gtrsim 6\) Å compared to the small intra-tetrahedron distance d ~3.3 Å [as shown in Fig. 1a], to model in Ba3Yb2Zn5O11, we begin by assuming that all inter-tetrahedron interactions are negligible and model the system as a set of decoupled “small” tetrahedra, hence focusing on the intra-tetrahedron interactions, which are expected to dominate. Since the first excited crystal field level is ~38.2 meV above the ground doublet31, we can project the Hamiltonian describing the Yb3+–Yb3+ interactions34 into the lowest ground state doublet, ignoring any higher-order perturbative corrections from virtual crystal field excitations5,17,34,59. Symmetry strongly constrains the allowed exchange interactions between these resulting (from the projection) S = 1/2 effective spins, with a bond symmetry group of 2 mm and full tetrahedral symmetry \(\bar{4}3m\) about the center of each tetrahedron60. These interactions are expected to act pair-wise between the effective Yb3+ spins 1/2 in each tetrahedron34, Si, with the model for each individual tetrahedron taking the form30,31,33

$$\begin{array}{ll}{{{{\mathcal{H}}}}}_{{{{\rm{eff}}}}}\equiv &\mathop{\sum }\limits_{i=1}^{4}\mathop{\sum}\limits_{j < i}\left[{J}_{zz}{S}_{i}^{z}{S}_{j}^{z}-{J}_{\pm }\left({S}_{i}^{+}{S}_{j}^{-}+{S}_{i}^{-}{S}_{j}^{+}\right)\right.\\&+ \;{J}_{\pm \pm }\left({\gamma }_{ij}{S}_{i}^{+}{S}_{j}^{+}+{{{\rm{h.c.}}}}\right)\\&+ \;\left.{J}_{z\pm }\left({\zeta }_{ij}\left[{S}_{i}^{z}{S}_{j}^{+}+{S}_{i}^{+}{S}_{j}^{z}\right]+{{{\rm{h.c.}}}}\right)\right]-{\mu }_{{{{\rm{B}}}}}\mathop{\sum }\limits_{i=1}^{4}{{{\bf{B}}}}\cdot {{{{\boldsymbol{\mu }}}}}_{i},\end{array}$$
(3)

where we have included a magnetic field B = μ0H with the local Yb3+ magnetic moment operators μi given by

$${{{{\boldsymbol{\mu }}}}}_{i}\equiv {\mu }_{B}\left[{g}_{\pm }\left({{{{\hat{\boldsymbol{x}}}}}}_{i}{S}_{i}^{x}+{{{{\hat{\boldsymbol{y}}}}}}_{i}{S}_{i}^{y}\right)+{g}_{z}{{{{\hat{\boldsymbol{z}}}}}}_{i}{S}_{i}^{z}\right],$$
(4)

and associated g-factors, gz and g±. We borrow definitions from ref. 61 for the bond form factors γij and ζij and the local axes (\({{{{\hat{\boldsymbol{x}}}}}}_{i}\), \({{{{\hat{\boldsymbol{y}}}}}}_{i}\), \({{{{\hat{\boldsymbol{z}}}}}}_{i}\)). Note that there is freedom in choosing the sign of Jz±, as it can be changed by a basis transformation.

The spectrum of this Hamiltonian is partly constrained by tetrahedral symmetry. The sixteen states of the four-spin system break into the irreducible representations of \(\bar{4}3m\)

$${A}_{2}\oplus 3E\oplus {T}_{1}\oplus 2{T}_{2},$$
(5)

under the action of the tetrahedral group35. This gives a level structure that consists of a singlet (A2), three doublets (E), and three triplets (T1 or T2), with the highest lying doublet (E) and triplet (T2) being nearly degenerate; see Fig. 2a. We use the parameters of Rau et al.35, which were fit using data from powder samples in finite magnetic fields, and do not reconsider here these values:

$$\begin{array}{ll}{J}_{zz}&\!\!\!\!=-0.040\,{{{\rm{meV}}}},\quad{J}_{\pm}\,\,=+0.141\,{{{\rm{meV}}}},\\ {J}_{\pm \pm }&\!\!\!\!=+0.160\,{{{\rm{meV}}}},\quad{J}_{z\pm }=+0.302\,{{{\rm{meV}}}},\end{array}$$
(6)

with g-factors gz = 2.726 and g± = 2.301.

Theoretically, we compute the unpolarized inelastic neutron scattering intensity, I(Q, E), using the same methods presented in refs. 30,35, which is given by

$$I({{{\bf{Q}}}},E)\propto \frac{| {{{{\bf{k}}}}}_{{{{\rm{f}}}}}| }{| {{{{\bf{k}}}}}_{{{{\rm{i}}}}}| }F{({{{\bf{Q}}}})}^{2}\mathop{\sum}\limits_{\mu \nu }({\delta }_{\mu \nu }-{\hat{Q}}_{\mu }{\hat{Q}}_{\nu }){S}_{\mu \nu }({{{\bf{Q}}}},E),$$
(7)

where F(Q) is the atomic form factor of Yb3+53. kf and ki are the final (scattered) and incident neutron wave vectors, respectively, and Q ≡ kf − ki. The spin structure factor for each tetrahedron is given by

$${S}_{\mu \nu }({{{\bf{Q}}}},E)\equiv \mathop{\sum}\limits_{nm}{\rho }_{m}\left\langle m| {\mu }_{-{{{\bf{Q}}}}}^{\mu }| n\right\rangle \left\langle n| {\mu }_{{{{\bf{Q}}}}}^{\nu }| m\right\rangle \delta (E-{E}_{n}+{E}_{m}),$$
(8)

where \({\rho }_{m}\equiv {e}^{-\beta {E}_{m}/T}/Z\) is the Boltzmann weight and the Fourier transform of the magnetic moment operator is given by \({{{{\boldsymbol{\mu }}}}}_{{{{\bf{Q}}}}}\equiv \frac{1}{4}\mathop{\sum }\nolimits_{i = 1}^{4}{e}^{i{{{\bf{Q}}}}\cdot {{{\bf{r}}}}}{{{{\boldsymbol{\mu }}}}}_{i}\), with μi given by Eq. (2). Here we note that the measured inelastic neutron scattering data has been divided by factor of kf/ki and thus the comparisons were made between the experimental and theoretical Sμν(Q, E) results.