Breathing chromium spinels: a showcase for a variety of pyrochlore Heisenberg Hamiltonians

TL;DR

This paper reveals that long-range Heisenberg interactions govern the magnetic phenomena in chromium breathing pyrochlore materials, uncovering diverse spin liquid behaviors influenced by chalcogen choice, which could guide future studies on spin liquids.

Contribution

It demonstrates that the magnetic interactions in these materials are primarily long-range Heisenberg types, a previously unrecognized aspect, and classifies their complex spin liquid behaviors based on chalcogen variation.

Findings

Long-range Heisenberg interactions dominate the magnetic behavior.

Different chalcogens lead to various classical spin liquid states.

The materials exhibit diverse degeneracy manifolds, including Coulomb and spiral spin liquids.

Abstract

We address the long-standing problem of the microscopic origin of the richly diverse phenomena in the chromium breathing pyrochlore material family. Combining electronic structure and renormalization group techniques we resolve the magnetic interactions and analyze their reciprocal-space susceptibility. We show that the physics of these materials is principally governed by long-range Heisenberg Hamiltonian interactions, a hitherto unappreciated fact. Our calculations uncover that in these isostructural compounds, the choice of chalcogen triggers a proximity of the materials to classical spin liquids featuring degenerate manifolds of wave-vectors of different dimensions: A Coulomb phase with three-dimensional degeneracy for LiInCr4O8 and LiGaCr4O8, a spiral spin liquid with two-dimensional degeneracy for CuInCr4Se8 and one-dimensional line degeneracies characteristic of the face-centered cubic antiferromagnet for LiInCr4S8, LiGaCr4S8 and CuInCr4S8. The surprisingly complex array of prototypical pyrochlore behaviors we discovered in chromium spinels may inspire studies of transition paths between different semi-classical spin liquids by doping or pressure.

Tables

Table 1. Table 1: Exchange couplings of six chromium breathing pyrochlore compounds . The first six lines are calculated for for the room temperature structures, while the last two lines are for low temperature structures. All Hamiltonian parameters are calculated within GGA+ U 𝑈 U with J H = 0.72 subscript 𝐽 H 0.72 J_{\rm H}=0.72  eV. The values of the J i subscript 𝐽 𝑖 J_{i} are given with respect to classical spins of length S = 3 / 2 𝑆 3 2 S=3/2 and without double counting of bonds [see Eq. ( 1 )].

Material	T(K)	$U$ (eV)	$J$ (K)	$J^{\prime }$ (K)	$J_2$ (K)	$J_{3a}$ (K)	$J_{3b}$ (K)	${\theta }_{\mathrm{CW}}$ (K)
\ceLiInCr4O8	RT	1.91	59.8(2)	22.0(2)	0.3(1)	1.9(1)	0.9(1)	-332
\ceLiGaCr4O8	RT	1.50	66.2(3)	100.0(2)	0.7(1)	2.0(1)	1.3(1)	-658
\ceLiInCr4S8	RT	1.77	-0.3(1)	-28.0(1)	0.7(1)	5.3(1)	2.4(1)	30
\ceLiGaCr4S8	RT	1.85	-7.7(1)	-12.2(1)	1.2(1)	6.1(1)	3.0(1)	-20
\ceCuInCr4S8	RT	1.43	14.7(1)	-26.0(1)	1.1(1)	6.4(1)	4.5(1)	-70
\ceCuInCr4Se8	RT	1.68	-25.4(2)	-31.0(1)	0.3(1)	4.8(1)	3.9(1)	135
\ceLiInCr4O8	20	1.77	40.6(2)	40.8(2)	0.3(1)	2.0(1)	0.9(1)	-332
\ceLiGaCr4S8	10	1.80	-10.8(1)	-9.5(1)	1.2(1)	6.2(1)	3.1(1)	30

Equations

$$H=\sum _{i<j}{J}_{ij}{\mathbf{S}}_i\cdot {\mathbf{S}}_j.$$

$$\begin{split}F(Q)&=-0.3094e^{-0.0274\big{(}\frac{Q}{4\pi}\big{)}^{2}}+0.36804e^{-17.0355\big{(}\frac{Q}{4\pi}\big{)}^{2}}\\ &\mspace{20.0mu}+0.6559e^{-6.5236\big{(}\frac{Q}{4\pi}\big{)}^{2}}+0.2856\,,\end{split}$$

$$\chi (\mathbf{k})\sim \frac{1}{\frac{2}{3T}-\frac{1}{J(\mathbf{k})}}.$$

$${\mathbf{h}}_i\equiv \frac{\delta H}{\delta {\mathbf{S}}_i},$$

Full text

Introduction

Over the past thirty years, materials with magnetic moments on the vertices of networks of corner-shared triangular or tetrahedral units have played center stage in the experimental search for exotic magnetic states driven by competing, i.e. frustrated, interactions Lacroix2011 . In this context, the pyrochlore lattice of corner-shared tetrahedra has emerged as a quintessential example of high frustration in three dimensions with a multitude of new phenomena having been uncovered Gardner2010 .

In particular, for rare-earth pyrochlore oxides with a trivalent $4f$ magnetic ion, there is a wealth of stoichiometric compounds available in large single-crystal form, necessary for detailed neutron scattering studies. This, along with the synergetic dialogue between theory and experiments, has led to the experimental discovery and the theoretical rationalization of numerous phenomena like spin ice physics, magnetic moment fragmentation, order-by-disorder, fragile splayed-ferromagnetism and candidate spin liquid phenomenology Hallas2018 ; Rau2019 .

In contrast with the above rare-earth systems characterized by strongly anisotropic interactions at the $O(1)$ K scale, there is a paucity of magnetic compounds involving transition metal ions on a pyrochlore network and coupled via a primarily isotropic Heisenberg exchange of a high ($O(10^{2})$ K) energy scale. The availability of such compounds in single-crystal form would empower researchers to explore and possibly discover novel collective behaviors existing over a wider and more experimentally accessible temperature window. Specifically, this would allow exposing new ways in which perturbations, which might have appeared inconsequential from a cursory assessment, ultimately emerge as the dictating forces behind the low-temperature and low-energy scale physics of the material. In that regard, the recent successful synthesis of magnetic pyrochlore fluorides $ABM_{2}$F7 ($A=$Na; $B=$Ca, Sr; $M=$ Ni, Co, Fe, Mn) in large single crystals is indeed an exciting and most welcome development in the field of highly-frustrated magnetism Plumb2019 . There is, however, a short-coming with the latter compounds whose implications remain to be fully ascertained: The cation disorder on the non-magnetic $A$ and $B$ sites is expected to randomize the $M$-$M$ superexchange and, most likely, ultimately drive the low-temperature state of these systems to a semi-classical spin glass phase, albeit with some nontrivial spin dynamics Plumb2019 . In this context, the disorder-free so-called breathing chromium spinels, which define a fairly broad range of materials (e.g. \ceLiInCr4O8, \ceLiGaCr4O8, \ceLiInCr4S8, \ceLiGaCr4S8, \ceCuInCr4S8 and \ceCuInCr4Se8) Okamoto2013 ; Duda2008 ; Tanaka2014 ; Nilsen2015 ; Okamoto2015 ; Lee2016 ; Saha2016 ; Okamoto2017 ; Wawrzynczak2017 ; Pokharel2018 ; Okamoto2018 , constitute a significant opportunity for carrying out the above research program in and below the $O(10^{2})$ K temperature scale.

The introduction of a breathing degree of freedom to the regular pyrochlore lattice, characterized by the alternation of small and large tetrahedra, results in two inequivalent nearest-neighbor exchange couplings, with $J$ for the small and $J^{\prime}$ for the large tetrahedra, and adds a new level of complexity and thus richness to the regular pyrochlore Heisenberg Hamiltonian Lacroix2011 . In the breathing chromium spinels, of chemical formula $AA^{\prime}$Cr4X8 ($A$=Li, Cu; $A^{\prime}$=Ga, In; $X$=O, S, Se), the magnetic Cr3+ spin $S=3/2$ resides on a breathing pyrochlore lattice (see Figure 1). Of particular interest, thanks to their significantly different size, the ions at the $A$ and $A^{\prime}$ sites do not chemically admix (i.e. they order), and the breathing chromium spinels are not subject to the $A/B$ site disorder and consequential exchange randomness of the aforementioned $ABM_{2}$F7 pyrochlores. While an ordered arrangement of the $A$ and $A^{\prime}$ cations was first reported over half-a-century ago Joubert1966 , it is only recently that breathing chromium spinels have become a subject of focused attention Okamoto2013 . A variety of experimental results have been reported, including a complex sequence of structural and magnetic transitions in \ceLiInCr4O8 and \ceLiGaCr4O8 Nilsen2015 ; Lee2016 ; Tanaka2014 ; Saha2016 as well as a half-magnetization plateau in high magnetic field measurements on \ceLiInCr4O8 Okamoto2017 .

On the theoretical front, the properties of the classical nearest-neighbor breathing pyrochlore Hamiltonian have been investigated Benton2015 , and spin-spin correlations in the ground state Tsunetsugu2017 and the effects of spin-lattice coupling Aoyama2019 have been explored. Possible implications of the topological aspects of magnetic excitations in the breathing pyrochlore lattice have been pointed out Li2016 ; Ezawa2018 . However, to the best of our knowledge, there appears to have been no attempt to (i) provide a microscopic determination of the Cr-Cr exchange interactions and to (ii) investigate using an unbiased theoretical framework the implications of the concurrent thermal and quantum fluctuations operating in these compounds. Such an endeavor is necessary to gain a deeper microscopic understanding of the already available experimental results beyond the simplest interpretation currently provided as well as identifying new and exciting avenues of investigation in these materials – this is the purpose of the present work. As we shall discuss later on, different signs (ferromagnetic or antiferromagnetic) exchange couplings $J$ and $J^{\prime}$ are realized in the oxides ($X=$O), sulfide ($X=$S) and selenides ($X=$Se) breathing chromium spinels. This provide an uncharted territory for the investigation of the nontrivial role played by long-range interactions beyond $J$ and $J^{\prime}$ in these materials. This leads us to propose that breathing chromium spinels constitute a promising platform to explore the role of strong competing long-range interactions in a semi-classical ($S=3/2$) regime. In particular, we predict these to engender new forms of nontrivial correlations in the cooperative paramagnetic regime characterized by different classically degenerate ground state manifolds of different dimension and intricate magnetic orders at low temperatures. Examples include emergent effective Heisenberg antiferromagnetism on the face-centered cubic lattice Benton2015 , a Coulomb phase Moessner1998 and a spiral spin liquid Bergman2007 .

Results

Breathing pyrochlore Hamiltonians

Room temperature structures.- We first perform electronic structure calculations for all presently known room temperature structures of chromium breathing spinels. We use the structures of Okamoto et al. (Ref. Okamoto2013, ) for \ceLiInCr4O8 and \ceLiGaCr4O8, the structures of Okamoto et al. (Ref. Okamoto2018, ) for \ceLiInCr4S8, \ceLiGaCr4S8 and \ceCuInCr4S8 and that of Duda et al. (Ref. Duda2008, ) for \ceCuInCr4Se8. The main exchange connectivity of all structures is illustrated in Figure 1. We follow the experimental literature Okamoto2013 ; Tanaka2014 ; Nilsen2015 ; Lee2016 in referring to the exchange couplings within the small and large tetrahedra as $J$ and $J^{\prime}$, respectively. As in the isotropic pyrochlore lattice Iqbal2017 , there are twelve second nearest neighbors $J_{2}$ and twelve third nearest neighbors which split into two symmetry inequivalent classes; six of these, named $J_{3a}$, have an in-between Cr3+ ion, and the six others, named $J_{3b}$, do not.

Figure 2 presents the first main result of our work. The exchange couplings as defined in the Heisenberg Hamiltonian of the form

[TABLE]

are shown as a function of interaction strength $U$ in the GGA+$U$ functional. For all six compounds, a unique $U$ value can be determined at which the couplings yield the experimental Curie-Weiss temperature (vertical lines). These fixed parameter values lead to the sets of Hamiltonian couplings given in Table 1. The $U$ values are very reasonable Tapp2017 for Cr3+ and are all located in a narrow interval $1.43\,\text{eV}\leq U\leq 1.91\,\text{eV}$, indicating that the great diversity among the chromium breathing pyrochlores, varying from antiferromagnetic to ferromagnetic Heisenberg models, can all be described with roughly the same exchange correlation functional. We find that both oxides, \ceLiInCr4O8 and \ceLiGaCr4O8, are dominated by antiferromagnetic $J$ and $J^{\prime}$ couplings (Figures 2A and B), with longer range couplings being negligibly small. Interestingly, it is found that in \ceLiGaCr4O8 the large tetrahedron is associated with the larger exchange coupling. The calculated breathing exchange anisotropy for \ceLiGaCr4O8 is $B_{f}\equiv J/J^{\prime}=0.66$, very close to the experimental estimate $B_{f}=0.6$ Okamoto2013 . However, according to our calculation, \ceLiInCr4O8 at $B_{f}\equiv J^{\prime}/J=0.37$ is far less anisotropic than assumed so far ($B_{f}=0.1$) Okamoto2013 . For the sulfides \ceLiInCr4S8, \ceLiGaCr4S8 and \ceCuInCr4S8 (Figures 2D-F), the large tetrahedra is characterized by ferromagnetic exchange $J^{\prime}$. However, the biggest surprise here is the presence of substantial second- and third-neighbor couplings. As each Cr3+ ion has a total of 24 such bonds, compared to six $J$ or $J^{\prime}$ bonds, they can have a substantial influence on the behavior of the materials. Interestingly, the fact that the experimentally determined Curie-Weiss temperature of \ceLiGaCr4S8 is small is thus readily explained by a cancellation of ferromagnetic $J$ and $J^{\prime}$ by antiferromagnetic $J_{2}$ and $J_{3a}$, $J_{3b}$ couplings rather than by opposite signs of $J$ and $J^{\prime}$ as hypothesized in Ref. Okamoto2018, . Opposite signs of $J$ and $J^{\prime}$ are only found for \ceCuInCr4S8. Finally, the selenide \ceCuInCr4Se8 (Figure 2C) is dominated by ferromagnetic $J$ and $J^{\prime}$ couplings; however, even for this compound, the third neighbor couplings $J_{3a}$ and $J_{3b}$ are antiferromagnetic and, most importantly, non-negligible.

Temperature dependence.- For most materials, only a room temperature structure is on record. However, Ref. Nilsen2015, gives a $T=20$ K structure for \ceLiInCr4O8, while Pokharel et al. Pokharel2018 provide a $T=10$ K structure for \ceLiGaCr4S8. The two Cr-Cr distances within small and large tetrahedra for these two structures are tabulated along with the measurements of the six room temperature structures in SI Table S2. While the breathing anisotropy of the oxide \ceLiInCr4O8 decreases significantly at low temperature, it is essentially temperature independent for the sulfide \ceLiGaCr4S8. Studying the two cubic low temperature structures can give us an indication of how temperature dependent the electronic behavior of the breathing pyrochlore materials is within the cubic $F\bar{4}3m$ structure. Clearly, a systematic study of the temperature dependence of the exchange couplings would require more detailed structural information for the other materials. The last two lines of Table 1 show the exchange interactions calculated for the two low temperature structures. Even though the low temperature structure of \ceLiInCr4O8 has a reduced breathing anisotropy of $r^{\prime}/r=1.029$, compared to $r^{\prime}/r=1.051$ at room temperature, the exchange couplings $J$ and $J^{\prime}$ become nearly identical (see SI Fig. S1a), with the room temperature exchange anisotropy of $J^{\prime}/J=0.37$ having been eliminated. For \ceLiGaCr4S8, while the low temperature structural anisotropy is still significant at $r^{\prime}/r=1.070$, the two ferromagnetic couplings $J$ and $J^{\prime}$ become very similar near the $U$ value determined by the experimental $\uptheta_{\rm CW}$ (see SI Fig. S1B). Thus, we find that the nearest-neighbor exchanges $J$ and $J^{\prime}$ remain antiferromagnetic for \ceLiInCr4O8 and ferromagnetic for \ceLiGaCr4S8, but in both cases they become very similar at low temperatures. At the same time, the further neighbor couplings remain very small for the oxide and substantial for the sulfide.

Energy mapping clearly identifies three different types of breathing pyrochlore Hamiltonians: Mostly antiferromagnetic oxides, sulfides with ferromagnetic large tetrahedron and antiferromagnetic longer range $(J_{2},J_{3\rm a},J_{3\rm b})$ exchange, and the mostly ferromagnetic selenide perturbed by crucial antiferromagnetic longer range interactions. Thus, we organize our pseudofermion functional renormalization group analysis of the Hamiltonians into oxide, sulfide and selenide sections.

Pseudofermion functional renormalization group calculations

Oxides

For \ceLiInCr4O8 and \ceLiGaCr4O8, the $J$ and $J^{\prime}$ couplings forming the tetrahedra are both antiferromagnetic. With only these two couplings, one would expect the typical bow-tie features with the associated pinch points in the reciprocal space susceptibility at sufficiently low temperatures (see Ref. Iqbal2019, ). In this previous study, we showed that the inclusion of even very weak second nearest-neighbor couplings can dramatically alter the susceptibility profile of the nearest-neighbor ground state phase and destroy the pinch points and bow-ties Conlon2010 . As stated above, for both breathing pyrochlore oxides, our DFT analysis reveals the presence of additional antiferromagnetic $J_{2}$, $J_{3a}$, and $J_{3b}$ couplings (see Table 1). At the classical level, it is known that when the dominant interactions $J$ and $J^{\prime}$ are large and antiferromagnetic such that the four spins at the vertices of each tetrahedra can be assumed to sum to zero approximately, then a third neighbor coupling $J_{3a}$ in an ideal pyrochlore lattice is equivalent at temperature $T\ll{\rm min}(J/k_{\rm B},J^{\prime}/k_{\rm B})$ to a $J_{2}$ coupling of opposite sign Chern2008 ; Conlon2010 . Thus, effectively, both chromium oxides are approximately described by a pyrochlore Hamiltonian with antiferromagnetic $J$, $J^{\prime}$ and ferromagnetic $J_{2}$: $(J_{2}-J_{3a})/\bar{J}=-0.025$ for \ceLiInCr4O8 and $-0.011$ for \ceLiGaCr4O8 ($\bar{J}=\sqrt{J^{2}+J^{\prime 2}}$ is the overall nearest-neighbor energy scale, listed in SI Table S1. According to Ref. Iqbal2019, , we expect that the spectral weight moves away from the pinch point and forms two symmetrical maxima, resulting in a hexagonal cluster pattern of intensities in the $[hhl]$ plane. Indeed, using the exchange couplings obtained for the room temperature structure (see Table 1), we employ PFFRG to calculate the magnetic susceptibility profiles which are shown in Figures 3A-F. They are found to closely resemble the susceptibility of the isotropic $J_{1}$-$J_{2}$ pyrochlore magnet for antiferromagnetic $J_{1}$ and ferromagnetic $J_{2}$ Iqbal2019 . It is worth noting that, due to the relatively weaker effective second nearest-neighbor coupling in \ceLiGaCr4O8, the pinch points and bow-ties are better preserved at low temperature for \ceLiGaCr4O8 (see Figure 3D) in comparison to \ceLiInCr4O8 (Figure 3A). However, from our PFFRG calculations we find that for both oxides, and at all temperatures down till the RG flow breakdown, the maxima of the susceptibility are always located at the high-symmetry $W$ point in the extended Brillouin zone, i.e., at a $\mathbf{q}=\frac{2\pi}{a}(2,1,0)$-type ordering vector (see SI section S2), as found in Ref. Iqbal2019, for the nearest-neighbor isotropic $S=3/2$ Heisenberg antiferromagnet. A direct space schematic illustration for the spin configuration corresponding to this order is shown in Figure 4A. In direct space, for this spin configuration, the pyrochlore lattice decomposes into two pairs of fcc sublattices. Within each pair, the two fcc sublattices are separated from each other by a nearest neighbor vector which is perpendicular to the $\frac{2\pi}{a}$ component of the ordering vector, e.g., for a $\frac{2\pi}{a}(2,1,0)$-ordering vector, one pair of fcc sublattices is separated by the vector $\frac{a}{4}(1,0,1)$ and the other pair is separated by the vector $\frac{a}{4}(1,0,-1)$. All fcc sublattices show a stripe antiferromagnetic order, i.e., spins are antiferromagnetically aligned in the $\frac{2\pi}{a}$-direction and ferromagnetically in the perpendicular directions. Within each pair, the neighboring spins of different sublattices are aligned antiparallel in spin space. The two pairs can be rotated freely with respect to each other. In reciprocal space, this configuration leads to main Bragg peaks at all $\frac{2\pi}{a}(2,1,0)$-type vectors which share the same $\frac{2\pi}{a}$ component, e.g., $\frac{2\pi}{a}(2,1,0)$, $\frac{2\pi}{a}(-2,1,0)$, $\frac{2\pi}{a}(0,1,2)$, $\frac{2\pi}{a}(0,1,-2)$, $\frac{2\pi}{a}(2,-1,0)$, $\frac{2\pi}{a}(-2,-1,0)$, $\frac{2\pi}{a}(0,-1,2)$, and $\frac{2\pi}{a}(0,-1,-2)$. In the reciprocal plane perpendicular to the $\frac{2\pi}{a}$-direction, the relative orientation of the sublattices results in subdominant Bragg peaks with half the intensity of the dominant ones at the $\frac{2\pi}{a}(1,0,1)$-type vectors in this plane, i.e., at $\frac{2\pi}{a}(1,0,1)$, $\frac{2\pi}{a}(-1,0,1)$, $\frac{2\pi}{a}(1,0,-1)$, $\frac{2\pi}{a}(-1,0,-1)$. For the classical model of \ceLiGaCr4O8, this state is the exact ground state given by the Luttinger-Tisza method and confirmed with the iterative minimization procedure. However, we find that this state is strongly dependent on the interplay between the two kinds of third nearest-neighbor couplings. This leads to the finding that for \ceLiInCr4O8, due to the larger difference between $J_{3a}$ and $J_{3b}$, an incommensurate order on the sublattices is stabilized as the ground state found by iterative minimization. We would like to emphasize that the stabilization of the orders discussed above rests crucially on the presence and the relative magnitude of the third nearest-neighbor couplings. Indeed, classically, the pure $J_{1}$-$J_{2}$-model on a regular $J_{1}\equiv J=J^{\prime}$ pyrochlore lattice would feature a $\mathbf{q}=\mathbf{0}$ state, completely different from the $\frac{2\pi}{a}(2,1,0)$-order found here.

To make contact with experimental results, it is useful to assess the temperature dependence of the structure factor as calculated within PFFRG. The relation $T=(2\pi/3)S(S+1)\Lambda$ can be used to relate the infrared cut-off $\Lambda$ employed in the PFFRG framework to temperature $T$ Iqbal2019 ; Iqbal2016 . The relation is obtained by comparing the classical limit ($S\to\infty$) of PFFRG, where only the RPA diagrams contribute, i.e. a mean-field description, with the conventional spin mean-field theory formulated in terms of temperature $T$ instead of $\Lambda$. The resulting estimated ordering temperatures are given in SI section S2. As explained in detail in the Appendix A of Ref. Iqbal2019, , in the $S\to\infty$ limit, the absence of higher diagrammatic orders in $1/S$ in the one-loop PFFRG used here implies that, in particular for the nearest-neighbor pyrochlore antiferromagnetic Heisenberg model, a spurious divergence of the susceptibility at finite temperature occurs which is not expected from what is well established for the classical problem Moessner1998 ; Iqbal2019 . Thus, our estimates of the ordering temperatures in all of the pyrochlore magnets considered here are systematically too high Iqbal2018 . Our comparisons to an experimental temperature $T$ are shown for a $\Lambda$ obeying $\Lambda/\Lambda_{\rm c}=T^{\rm exp}/T^{\rm exp}_{\rm c}$. It is expected that multi-loop implementations of PFFRG will lead to significant improvements for the calculated temperature scales Rueck2018 ; Kugler2018 ; however, these advanced schemes are still under development.

We have also performed the PFFRG calculation for the \ceLiInCr4O8 Hamiltonian corresponding to its low temperature structure (see Table 1, last two lines). The calculated susceptibility profile (see SI Fig. S4A-C) is seen to be virtually indistinguishable from that obtained from the room temperature structure shown in Figures 3A-C even though the breathing anisotropy of the Hamiltonian is significantly different, having evolved from $J^{\prime}/J=0.37$ at room temperature to $J^{\prime}/J\cong 1$ at $T=20$ K. The general observation of the magnetic response being largely independent of the breathing anisotropy can be understood from a classical picture. As discussed in SI section S3, the eigenvalues of the interaction matrix $\mathcal{J}({\bf k})$ for a $J$-$J^{\prime}$-only model exhibit flat bands at the lowest energies which determine the typical bow-tie features of pyrochlore antiferromagnets. Interestingly, these flat bands persist for any $J>0$ and $J^{\prime}>0$ Benton2015 (see SI Fig. S7A-D in SI section S3), rendering the momentum space profile of the magnetic susceptibility independent of $J^{\prime}/J$. As a consequence, longer range $J_{2}$, $J_{3a}$, $J_{3b}$ have the opportunity to play a crucial role in determining the magnetic behavior of the system. In the present case, however, we find that the effective second-neighbor coupling of \ceLiInCr4O8 at 20 K given by $(J_{2}-J_{3a})/\bar{J}=-0.030$ is nearly unchanged compared to room temperature (see SI Table S1) thus leading to almost identical magnetic responses.

It is now interesting to compare and verify the accuracy of the Hamiltonian determined for \ceLiInCr4O8 against polarized neutron scattering data from powder samples of Okamoto et al. Okamoto2015 . It was shown by Benton et al. Benton2015 that within the self-consistent Gaussian approximation (SCGA), the experimental estimate $J^{\prime}/J=0.1$, or even a range $0.05\lesssim J^{\prime}/J\lesssim 0.15$, is consistent with the neutron data. We also calculate the product $|F(Q)|^{2}S(Q)$ of the magnetic form factor $|F(Q)|^{2}$ and the powder averaged structure factor $S(Q)$ using

[TABLE]

for all the materials considered here Brown2004 . The $S(Q)$ is calculated using both the quantum $S=3/2$ PFFRG and classical Monte Carlo (CMC) methods. We show in Figures 5A-B the comparison to the $T=30$ K and the $T=150$ K neutron scattering data. The PFFRG calculations for the DFT determined Hamiltonians for \ceLiInCr4O8 at room temperature ($J^{\prime}/J=0.37$) as well as at $T=20$ K ($J^{\prime}/J=1$), upon inclusion of the second- and third-neighbor couplings, give an excellent agreement with the $T=30$ K neutron data. In contrast, the CMC results, both for the room temperature as well as $T=20$ K couplings, progressively lose agreement as $Q$ is lowered. This observation may hint towards the non-negligible role of quantum fluctuations in accurately capturing the correlations especially at low $Q$ values. The excellent agreement of the earlier SCGA calculations of Ref. Benton2015, is also likely rooted in the fact that in this approach the “hard constraint”, namely, that the spin vector on each site has magnitude $S$, is relaxed and only implemented on average. This may tacitly incorporate some effects of quantum fluctuations by allowing for fluctuations of the classical moments Kimchi2014 . At $T=150$ K, the neutron data essentially follows the $Q$ dependence of the magnetic form factor Okamoto2015 ; while the computed scattering intensities for the room temperature and $T=20$ K Hamiltonians are consistent, the disagreement with the experimental $T=150$ K neutron data is neither satisfactory for the PFFRG nor for the CMC results. A similar observation has already been made by Benton et al. Benton2015 and thus remains unexplained (see also SI section S2).

Sulfides

We have performed PFFRG calculations for all three chromium sulfides, \ceLiInCr4S8, \ceLiGaCr4S8 and \ceCuInCr4S8. As we find very similar susceptibility profiles for the three compounds, we show the one for \ceCuInCr4S8 in Figures 3G-I and the other two in SI section S2, SI Figs. S4D-I. \ceCuInCr4S8 has a ferromagnetic $J^{\prime}$ which is the strongest among all the interactions present in this material. At the low temperatures that we are considering here, we can therefore assume that the four spins on the large tetrahedra point in approximately the same direction, thereby behaving as effective spin-6 entities. Since the large tetrahedra are arranged in a face centered cubic (fcc) magnetic lattice (see SI Figs. S5A-B), the system can effectively be mapped onto a spin-6 Heisenberg Hamiltonian on a fcc lattice with renormalized nearest-neighbor couplings $J_{1}^{\rm fcc}=(J+4J_{2}+2J_{3a}+2J_{3b})/16$. These couplings are all antiferromagnetic with values $J_{1}^{\rm fcc}=1.11$ K for \ceLiInCr4S8, $J_{1}^{\rm fcc}=0.95$ K for \ceLiGaCr4S8 and $J_{1}^{\rm fcc}=2.56$ K for \ceCuInCr4S8. Interestingly, the nearest-neighbor antiferromagnetic classical Heisenberg model on the fcc lattice features a subextensive ($\mathcal{O}[L]$) ground state degeneracy in the classical limit ($S\rightarrow\infty$), with associated wave vectors of the form ${\bf q}=\frac{2\pi}{a}(1,\delta,0)$ and symmetry-related ${\bf q}$ Henley1987 . This degeneracy appears as lines of strong intensity in the magnetic response. However, as further explained in SI section S5, the susceptibility profile of our pyrochlore model differs from that of the regular fcc antiferromagnet since it is modulated by a form factor arising from the presence of a tetrahedral basis of our effective fcc lattice. For the structure factors in the $[hk0]$ plane, this leads to a square ring type structure formed by lines of strong and almost constant intensity Benton2015 , see Figure 3H (and SI Figs. S4E and H). The point ${\bf q}=\frac{2\pi}{a}(1,0,0)$ (and symmetry-related ${\bf q}$-space positions) is special because it resides at the junctions of the ground state manifold ${\bf q}=\frac{2\pi}{a}(1,\delta,0)$ and ${\bf q}=\frac{2\pi}{a}(1,0,\delta)$. Because of this, it is found that collinear ordered states with the ordering wave vector of the type $\mathbf{q}=\frac{2\pi}{a}(1,0,0)$ are selected by thermal Henley1987 as well as quantum Oguchi1985 order-by-disorder mechanism. As the PFFRG method incorporates the combined effects of thermal and quantum fluctuations, and thus of their selection effect, it is interesting to note that our analysis reveals maxima at ${\bf q}=\frac{2\pi}{a}(1,0,0)$ type positions for all three sulfides (see SI section S2). The corresponding direct space spin configuration is illustrated in Figure 4B. The astonishing property that all sulfides show almost identical susceptibility profiles even though the ratio $J_{1}^{\rm eff}/J^{\prime}$ varies significantly between the three compounds can be understood by considering the eigenvalues of the interaction matrix $\mathcal{J}({\bf k})$ for a pyrochlore system with $J>0$ and $J^{\prime}<0$. As shown in SI Figs. S7E-H of SI section S3, the lowest band exhibits line-like degeneracies, regardless of the size of $J/J^{\prime}$, leading to a magnetic response which is largely independent of this ratio Benton2015 .

Plumier et al. Plumier1971 ; Plumier1977 studied \ceCuInCr4S8 using unpolarized neutrons. In order to compare our calculations to this experiment, we subtract the room temperature neutron spectra from the low temperature spectra in order to remove the background. In Figs. 5C-D, this data is compared to the calculated structure factor $S(Q)$ weighted with the magnetic form factor $|F(Q)|^{2}$ of Cr3+ ions, i.e., $|F(Q)|^{2}S(Q)$. For the 40 K neutron data, the experiment matches well with our numerical data from PFFRG up to an arbitrary scale factor. At 85 K, while the peak positions match very well between theory and the experimental intensity profile, the intensity has some deviation at low $Q$. This could be partly due to the uncertainty in the mapping of the flow parameter $\Lambda$ to temperature. Since, 85 K is about 2.4$T_{\rm c}$ for \ceCuInCr4S8, the calculated structure factor is compared at 2.4$\Lambda_{\rm c}$. Nevertheless, the comparison of the spin-spin correlation profile to the available experimental information lends strong support to the validity of the Hamiltonian we determine for \ceCuInCr4S8. On the other hand, CMC calculations are unable to capture the dominant correlations present in the $40$ K experimental data, possibly hinting at the role of quantum fluctuations in accurately capturing the important correlations at low temperature, similar to what we observe in \ceLiInCr4O8.

Selenide

The material \ceCuInCr4Se8 is the only selenide breathing chromium spinel with detailed structure reported for which we can perform calculations. At a first glance, the Hamiltonian parameters (Table 1) appear similar to those for \ceLiGaCr4S8, with a ferromagnetic large tetrahedron. However, the PFFRG result shown in Figure 3J-L displays a completely different spin susceptibility. In contrast to a repetition of the broken high intensity lines representative of the one-dimensional sub-extensive degeneracy that we discussed in the previous subsection, we find the two-dimensional degeneracy of a sphere in reciprocal space. The ferromagnetic large tetrahedron characterized by $J^{\prime}$ coupling form a large ($S=6$) moment fcc lattice with an effective nearest-neighbor interaction given by $J_{1}^{\rm fcc}=(J+4J_{2}+2J_{3a}+2J_{3b})/16$ which evaluates to $J_{1}^{\rm fcc}=-0.55$ K. The ferromagnetic nature of $J_{1}^{\rm fcc}$ in the selenide is indicative of the crucial difference between the selenide and the sulfides (for which the $J_{1}^{\rm fcc}$ are antiferromagnetic): in contrast to the sulfides, nearly negligible $J_{2}$ and smaller $J_{3a}$ and $J_{3b}$ couplings do not fully compensate the substantial ferromagnetic $J$ anymore. Nonetheless, clearly no simple ferromagnetic order is realized in \ceCuInCr4Se8. Rather, the two third-neighbor couplings (taken together) substantially perturb the FM large tetrahedra which leads to the appearance of a two-dimensional degeneracy of classical ground states. At the ordering temperature (SI section S2), the system develops an incommensurate spiral magnetic order with a helix pitch vector $\mathbf{q}\approx\frac{2\pi}{a}(0.521,0,0)$ (and symmetry-related points $(0,q,0)$ and $(0,0,q)$). The corresponding magnetic susceptibility profile obtained from PFFRG at the ordering temperature is shown in Figures 3J-K wherein one observes peaks in $S(\mathbf{q})$ at $\mathbf{q}\neq\mathbf{0}$. The corresponding direct space spin configuration is illustrated in Figure 4C. It is worth emphasizing that the shift of the spectral weight away from $\mathbf{q}=\mathbf{0}$ is the combined effect of both $J_{3a}$ and $J_{3b}$ couplings of a comparable magnitude Tymoshenko2017 . Above the ordering temperature, and within the cooperative paramagnetic regime, we observe that the spectral weight is distributed rather uniformly over the spiral surface suggesting the presence of an approximate spiral spin liquid on the pyrochlore lattice stabilized concomitantly by thermal and quantum fluctuations. While an incommensurate long-range ordered spiral phase of $(q,0,0)$ type has been observed for the material \ceZnCr2Se4 with pyrochlore magnetic lattice Tymoshenko2017 , it is only very recently that a spiral spin liquid has been observed in the chromium spinel MgCr2O4 Bai2019 .

The spin spiral surface found in the cooperative paramagnetic regime in a PFFRG calculation, is also found to be present in the corresponding classical model at $T=0$, as revealed by a Luttinger-Tisza analysis. However, the Luttinger-Tisza eigenstates on this spiral surface generically do not represent real normalizable spin configurations, i.e., the spin length on different fcc sublattices are not equal in magnitude. It turns out that, it is only for the $(q,0,0)$-type ordering vectors (with $q\approx\frac{2\pi}{a}(0.45,0,0)$) that the minimal energy can be achieved by a configuration of normalized spins. This ground state configuration then breaks the symmetry of the lattice as it is only governed by one of the three $(q,0,0)$-type vectors. Indeed, our classical Monte Carlo calculations give an ordering vector of $\mathbf{q}=\frac{2\pi}{a}(0.40\pm 0.04)$ at the transition temperature. Although the spiral surface is inaccessible at $T=0$ due to the violation of the spin length constraint, thanks to thermal and/or quantum fluctuations it can in principle be made accessible. Indeed, for \ceCuInCr4Se8, our PFFRG analysis shows that thermal and quantum fluctuations taken together are able to restore a well-defined spin spiral surface (see Figs. 3J-K) right above the ordering temperature. In contrast, our classical Monte Carlo simulations, right above the ordering temperature, find a highly nonuniform distribution of spectral weight, and thus an absence of spiral spin liquid, however, as the temperature is increased to $\sim 1.5$ times the ordering temperature, a spiral surface appears (compare SI Fig. S8 of SI section S4 with Figs. 3J-K above). These findings lend support to the role played by quantum fluctuations in aiding the stabilization of a spiral spin liquid. It is worth noting that within the ordered phase, this single-$\mathbf{q}$ state does not however represent a single planar spiral throughout the entire lattice. Instead, similar to what is found in \ceCuInCr4S8, the lattice can be divided into two pairs of fcc sublattices. Within each pair, the sublattices are connected by a nearest neighbor vector perpendicular to the ordering vector, and the planar spiral orders within each pair are of equal pitch and phase. However, the two pairs of fcc sublattices are found to be out of phase with respect to each other by approximately $5^{\circ}$. This offset in phase is a direct and unique consequence of the breathing anisotropy, and is found to vanish in the isotropic pyrochlore lattice.

Discussion

One of the important findings of this study is the significant and heretofore unappreciated role played by long-range exchange interactions on the breathing pyrochlore lattice. While we have previously pointed out the sensitivity of the spin structure factor to the sign and size of the subleading couplings in the $J_{1}$-$J_{2}$ quantum pyrochlore Heisenberg antiferromagnet Iqbal2019 , we have found in the present work for the chromium breathing pyrochlores an impressive demonstration of the decisive role played by longer-range couplings, i.e., the two kinds of symmetry inequivalent third-nearest neighbor couplings $J_{3a}$ and $J_{3b}$. On the other hand, the effects of breathing anisotropy are shown to be surprisingly minor. Indeed, we observe that in \ceLiInCr4O8, the effective ferromagnetic second-neighbor coupling of 2.5% of the average nearest-neighbor coupling far outweighs the substantial room temperature breathing anisotropy of $J^{\prime}/J=0.37$ in its effect on the structure factor profile. In the case of the sulfides, although the $J^{\prime}/J$ ratios vary between $93$ and $-1.8$ among the three compounds, the susceptibility profiles appear all very similar. In this case, this is rooted in the effective mapping of the system to an emerging effective nearest-neighbor fcc Heisenberg antiferromagnet with the fcc lattice sites being occupied by ferromagnetic tetrahedra featuring a large $S=6$ magnetic moment. In this case, we show that the spin structure factor displays line-like degeneracies which are, once more, largely independent of the breathing anisotropy ratio $J^{\prime}/J$. Finally, for the dominantly ferromagnetic selenide \ceCuInCr4Se8, we demonstrate that the combined effect of the two third-neighbor couplings $J_{3a}$ and $J_{3b}$ drastically perturbs the ground state away from the simple ferromagnetic order. Interestingly, we find that the perturbed state corresponds to an approximate spiral spin liquid above the ordering temperature and within the cooperative paramagnetic regime, where the individual wave vectors form a sphere-like manifold in reciprocal space. Our PFFRG analysis indicates that the combined effect of quantum and thermal fluctuations only leads to a weak order-by-disorder selection into an incommensurate spiral state. Approximate spiral spin liquids on three-dimensional lattices are so far known only on the diamond lattice in MnSc2S4 Bergman2007 ; Gao2016 ; Iqbal2018 , with a very recent occurrence having also been reported for the pyrochlore material MgCr2O4 Bai2019 .

Conclusion

We have theoretically investigated six chromium spinels featuring crystalline, and thus magnetic exchange breathing anisotropy. The Hamiltonians which were determined from density functional theory calculations and investigated using the pseudofermion functional renormalization group method showcase a colorful selection of magnetic properties arising from frustrated interactions. We find that the oxide compounds are in a perturbed Coulomb phase Moessner1998 , reminiscent of the famous spin ice materials. The sulfides display an effective fcc lattice Heisenberg antiferromagnetic type behavior, and the selenide features incommensurate magnetic correlations. As a unifying feature, all materials are found to be close to a classical degeneracy which shows up as different variants: The oxides are near a phase featuring an extensive number of classical ground states scaling exponentially in the volume of the system. The selenide is approximately degenerate on a two-dimensional surface with a manifold of states scaling exponentially in $L^{2}$ (where $L$ is the linear dimension of the system). Finally, the sulfides are characterized by approximate line-like degeneracies with the number of ground states scaling exponentially in the linear dimension $L$. This variety of different behaviors in a family of related materials represents an attractive feature which promises a wealth of opportunity for future investigations; on the one hand, doping series interpolating between different types of ground states could be very interesting. On the other, it is an intriguing question how other known but only sketchily characterized breathing chromium spinels like \ceCuGaCr4S8, \ceAgInCr4S8, \ceCuGaCr4Se8, \ceAgInCr4Se8 Haeuseler1977 fit the picture outlined in this study. Of particular importance are experimental investigations with single crystals which would allow to assess our theoretical predictions.

Materials and Methods

Energy mapping method

We determine the electronic structure and total energies of the breathing pyrochlores using the full potential local orbital (FPLO) basis set Koepernik1999 and the generalized gradient approximation (GGA) functional Perdew1996 . The GGA+$U$ Liechtenstein1995 correction to the exchange and correlation functional is used to deal with the strong electronic correlations on the Cr3+ $3d$ orbitals. Two parameters enter these calculations, the on-site interaction strength $U$ and the Hund’s rule coupling $J_{\rm H}$. As the Hund’s rule coupling is an intra-atomic interaction, we do not expect it to vary much from one compound to the other, and thus henceforth fix its value to $J_{\rm H}=0.72$ eV as commonly employed for Cr3+ Mizokawa1996 . Meanwhile, the on-site interaction $U$ was fitted from the experimentally determined Curie-Weiss temperatures as explained below. Heisenberg Hamiltonian parameters are extracted using the energy mapping approach Jeschke2011 ; Iqbal2018 : For this purpose, a $3\times 1\times 1$ supercell of the primitive unit of the cubic $F\bar{4}3m$ structure with $Cm$ space group is constructed. In this structure, there are 9 independent spins (out of 12 in total), allowing 75 spin configurations with different energies. Total energies are converged using $6\times 6\times 6$ $\mathbf{k}$ points. An arbitrary subset of 19 spin configurations allows us to fit six exchange interactions, using a Heisenberg Hamiltonian of the form of Eq. (1). Note that we count every bond only once. The quality of the fit is seen to be excellent (see SI Fig. S2 in SI section S1 for an example), indicating that the method works extremely well for the chromium breathing pyrochlores, may they be dominantly antiferromagnetic or ferromagnetic. Now, the $U$ parameter of the GGA+$U$ functional can be fixed: We calculate the Hamiltonian parameters for a small set of $U$ values, which we choose such that the Curie-Weiss temperatures $\uptheta_{\rm CW}=-\frac{1}{3}S(S+1)(3J+3J^{\prime}+12J_{2}+6J_{3a}+6J_{3b})$ (see Figure 1) calculated from these Hamiltonians cover a range of temperatures which includes the experimental $\uptheta_{\rm CW}$. We find the $U$ value from the condition $\uptheta^{\rm theor}_{\rm CW}=\uptheta^{\rm exp}_{\rm CW}$ by interpolation (see Figure 2). For the breathing pyrochlores, we find the $\uptheta^{\rm theor}_{\rm CW}$ to be monotonous functions of $U$ and thus to allow for a unique solution. Some more discussion on $\uptheta_{\rm CW}$ and $U$ can be found in SI section S1.

Pseudofermion functional renormalization group method

In the PFFRG scheme Reuther2010 , a general spin-$S$ Heisenberg Hamiltonian, e.g. Eq. (1), is treated by first re-expressing the spin operator at each lattice site by $2S$ spin-1/2 degrees of freedom using Abrikosov pseudofermions. Aside from the desired Hilbert space of spin-$S$ states this approach also introduces unphysical states with lower spin magnitudes. Since these states contribute less to the energy of the system compared to the physical ones, they act as excitations and, thus, naturally get excluded at low temperatures in the RG process Iqbal2019 . Next, the pseudo-fermionic Hamiltonian is treated using the Functional Renormalization Group (FRG) formalism Reuther2010 ; Metzner2012 ; Baez2017 ; Buessen2018 . Here, the bare fermionic Green’s function in Matsubara space is equipped with an infrared frequency cut-off $\Lambda$ suppressing the fermionic propagation at low frequencies. The PFFRG flow equations then follow from the $\Lambda$ dependence of the $m$-particle irreducible vertex functions which are given by an exact and infinite hierarchy of coupled integro-differential flow equations Metzner2012 . For a numerical implementation, this hierarchy however needs to be truncated, which is done by neglecting three-particle and higher vertex functions. Still, and most importantly, the truncated version of the method exactly sums up the full set of fermionic Feynman diagrams in two limits taken separately, namely, the large-$S$ limit Baez2017 and the large-$N$ limit Buessen2018 . In doing so, it provides an unbiased analysis of the competition between magnetic ordering and disordering tendencies Iqbal2019 ; Reuther2010 ; Metzner2012 ; Baez2017 .

The flow equations are numerically solved in direct space and, after the Fourier transform of the spin-spin correlations, one obtains the static susceptibility in reciprocal space $S(\mathbf{q})$ as a function of $\Lambda$. Note that the RG parameter $\Lambda$ can be interpreted as the temperature $T$ Iqbal2019 and is therefore used as a proxy to investigate the thermal evolution of the magnetic response. In the present study, we evaluate spin correlators within a cube of edge length of $5a$, where $a$ is the cubic lattice constant, which incorporates a total of 2315 correlated sites, producing well-converged results with a high $\mathbf{q}$-space resolution. Furthermore, we approximate the frequency dependence of the vertex functions by discrete grids containing 64 points for each frequency variable. An ordered system is identified by cusps or kinks in the susceptibility flow signaling a magnetic instability towards either magnetic or valence bond order Reuther2010 . In contrast, a smooth flow of the susceptibility down to $\Lambda\rightarrow 0$ indicates a magnetically disordered state, that is, a putative (quantum) spin liquid Reuther2010 .

Luttinger-Tisza method

In the Luttinger-Tisza method Kaplan2007 ; Lapa2012 the ground state and energy spectrum of a classical spin system are approximately calculated. The classical limit of Eq. (1) consists of replacing the spin operators by classical vectors, which are normalized (the so called strong constraint). To find an approximate ground state, this constraint is relaxed to be fulfilled only on average over the whole system. This allows to Fourier transform the interaction matrix $J_{ij}$ with respect to the underlying Bravais lattice of the actual spin lattice, leading to a matrix $J_{\alpha\beta}(\mathbf{k})$, where $\alpha,\beta$ label the basis points of the lattice.

The eigenvalues $\lambda(\mathbf{k})$ of this matrix subsequently give the energies of spin spiral states with wave vector $\mathbf{k}$, the aforementioned energy spectrum, and the absolute value of the corresponding eigenvector components give the relative length of the vector spins on the different basis points. This means, that for such a state to fulfill the strong constraint, all components have to have equal absolute value.

If the lattice is of Bravais type, there is only one basis point and therefore the Fourier transformed interaction matrix is a scalar, which means this condition is trivially fulfilled, rendering the Luttinger-Tisza method exact on such lattices. In this case, we can also, through the equivalence of the Luttinger-Tisza method and PFFRG for $S\to\infty$ Baez2017 , calculate the classical spin susceptibility

[TABLE]

Iterative minimization method

By employing the iterative minimization method Lapa2012 ; Iqbal2019 , we find the ground state of the classical version of Eq. (1) obtained by replacing the spin operators by normalized vectors. In this scheme, we typically start from a random spin configuration, randomly select one vector-spin $\mathbf{S}_{i}$ at a time and rotate it to make it antiparallel to its local field,

[TABLE]

thus minimizing the energy. In a given sweep, we update all the spins contained in the configuration so that on average each spin is updated once. For our calculations, we use 32 cubic unit cells of the pyrochlore lattice in each direction with periodic boundary conditions, totaling to $16\times 32^{3}=$524,288$$ spins in the system. Convergence is reached once the total energy difference per sweep is less than $10^{-10}$ K. To reduce the likelihood of the system getting stuck at a local energy minimum in configuration space, we carry out the minimization for at least 20 different random initial configurations. We also use specific non-random configurations, e.g. spin spirals, as a starting point, to test exact spin parameterizations and, possibly, find lower minimal energies than those achieved by starting with arbitrary random spin configurations.

Classical Monte Carlo method

We have performed classical Monte Carlo simulations for the model Hamiltonian parameters of the compounds as obtained from DFT given in Table 1, employing the standard single spin-flip technique. For each Monte Carlo update, we perform the Metropolis moves as many times as the number of spins in the configuration such that, on average, each spin is updated once. Starting from a random spin configuration and after allowing for the system to reach thermal equilibrium, we evaluate the Fourier transform of the equal-time spin-spin correlator, i.e., the static structure factor $S(\mathbf{q})$ for 50 different configurations where each configuration is separated by $100$ Monte Carlo updates. At lower temperatures, in order to increase the number of accepted moves, we restrict the newly generated spins to a small cone around its predecessor. The simulations are done for a system of $32\times 32\times 32$ cubic unit cells of the pyrochlore lattice with periodic boundary conditions in each direction, amounting to $16\times 32^{3}$ spins. The entire simulations are independently carried out 15 times and the average of the results is taken for suppression of numerical noise.