Modulated non-collinear magnetic structure of (Co_0.97⁵⁷Fe_0.03)₄Nb₂O₉ as revealed
by Mössbauer spectroscopy

Recently, the corundum-type compounds, M₄A₂O₉ (M = Mn, Co, Fe, A = Nb, Ta) Panja et al. (2021, 2018); Maignan and Martin (2018); Khanh et al. (2016); Ding et al. (2020); Zhang et al. (2021), have drawn great interests due to their rich physics such as magnetoelectric (ME) effect Panja et al. (2021, 2018), magnetic anisotropy Panja et al. (2021); Ding et al. (2020) and magnetism induced spontaneous electric polarization in Fe₄Ta₂O₉ Panja et al. (2018); Maignan and Martin (2018). Among these compounds, Co₄Nb₂O₉, which is the most studied one, has been found to show large ME effect below the Néel temperature of $T_{N}\sim 27\,K$ Kolodiazhnyi et al. (2011); Fang et al. (2015). Cross coupling between the electric field induced magnetization and magnetic field controlled electric polarization have been experimentally observed on powder and single crystalline samples Fang et al. (2015). The microscopic mechanisms underlining these interesting properties remain to be settled. It is believed that the detailed magnetic structures play very important roles in understanding these interesting properties Zhang et al. (2021); Panja et al. (2021); Solovyev and Kolodiazhnyi (2016); Yanagi et al. (2018); Deng et al. (2018); Matsumoto and Koga (2019).

The crystal structure of Co₄Nb₂O₉ is rather simple and have been investigated both by polycrystalline and single crystal X-ray and neutron diffraction methods Bertaut et al. (1961); Khanh et al. (2016); Deng et al. (2018); Ding et al. (2020). The magnetic structure was first reported by Bertaut et al. Bertaut et al. (1961) to be antiferromagnetically coupled ferromagnetic Co²⁺ chains with their magnetic moments along the crystal $c$-axis. It has been shown later by Khanh et al. Khanh et al. (2016) that this magnetic structure is incompatible with the experimentally measured ME effect on single crystal samples. A different magnetic structure with all the magnetic moments nearly parallel to the [11̄0] direction has been proposed based on single-crystal neutron diffraction. Within this model, the magnetic moments exhibit a small canting towards the $c$-axis but the projected moments on the $ab$-plane are all parallel with each other. Later on, both powder and single crystal neutron diffraction data have revealed another distinct non-collinear magnetic structure without any moment canting to the $c$-axis Deng et al. (2018); Ding et al. (2020). Moreover, different spin-flop behaviors have been reported in Co₄Nb₂O₉ and magnetic structures associated with these magnetic anomalies have been suggested to explain the observed ME effect Kolodiazhnyi et al. (2011); Fang et al. (2015); Khanh et al. (2016).

Clearly, one has to know the precise magnetic structure to understand the observed interesting ME effects. To reconcile the above mentioned diverse magnetic structures Bertaut et al. (1961); Khanh et al. (2016); Deng et al. (2018); Ding et al. (2020), we present the results of a local probe of the magnetic structure of (Co_0.97⁵⁷Fe_0.03)₄Nb₂O₉ compound with Mössbauer spectroscopy. Electric field gradient (EFG) tensor has been calculated by density functional theory and served as the coordinate system to determine the directions of the magnetic moments. Tentative fittings of our Mössbauer spectra with previously published magnetic structures have all failed. Therefore, we proposed a more complex helicoidally modulated magnetic structure to satisfactorily describe the whole series of our Mössbauer spectra. The proposed non-collinear magnetic structure may be promising in the understanding of the complex ME observed in related systems.

Polycrystalline samples of (Co_0.97⁵⁷Fe_0.03)₄Nb₂O₉ and Co₄Nb₂O₉ were prepared by using the conventional solid state reaction technique Kolodiazhnyi et al. (2011). To prevent reaction of the sample with oxygen, both mixing and reaction procedures were made inside an argon glove box with the oxygen level controlled below 0.1 ppm. A total of $\sim$1 gram Co/⁵⁷Fe, Co₃O₄, and Nb₂O₅, in a proportion to meet the correct oxygen content, were used and the reaction temperature was fixed at 1100 ^oC. Several reactions with thorough intermediate grindings were needed to improve homogeneity of the doped ⁵⁷Fe and to reduce foreign phases such as CoNb₂O₆.

Phase purity was checked by room temperature X-ray powder diffraction (XRPD) and the refinements were done by using the FullProf suite ful . Magnetic properties were measured using a dc superconducting quantum interference device (SQUID) magnetometer (Quantum Design). Mössbauer measurements were performed in transmission geometry with a conventional spectrometer working in constant acceleration mode. A 50 mCi $\gamma$-ray source of ⁵⁷Co embedded in Rh matrix and vibrating at room temperature was used. The drive velocity was calibrated using $\alpha$-Fe foil for high velocity measurements and sodium nitroprusside (SNP) for low velocity measurements. The isomer shift quoted in this work are relative to that of the $\alpha$-Fe at room temperature.

The computational work was carried out within the ELK code elk , which is based on the full potential linearized augmented plane waves (FP-LAPW) method. The Perdew-Wang/Ceperley-Alder local spin density approximation (LSDA) exchange-correlation functional Perdew and Wang (1992) was used. LSDA+U calculation was done in the fully localized limit (FLL) and by means of the Yukawa potential method Liechtenstein et al. (1995) with a screening length of $\lambda=3.0$. Slater integrals are calculated according to $\lambda$ and the resulting Coulomb interaction parameters are $U=3.0$ eV and $J=0.88$ eV which are similar to previously used values Solovyev and Kolodiazhnyi (2016). The muffin-tin radii $R_{MT}$ were set to 2.33 a.u., 2.07 a.u., and 1.40 a.u. for Co, Nb, and O atoms, respectively. The plane-wave cutoff was set to $R_{MT}\times|\textbf{G}+k|_{max}=7.0$ and the maximum G-vector for the potential and density was set to $|\textbf{G}|_{max}=12.0$. k-point meshes of $4\times 4\times 2$ a total of 32 k-points were used due to very time consuming calculations when spin orbital coupling (SOC) were included. Experimental lattice parameters of Co₄Nb₂O₉ at 50 K ($a=b=5.180$ Å and $c=14.163$ Å) were used for our calculation Ding et al. (2020).

In Fig. 1, the Rietveld refinements of the room temperature XRPD data for both the pure Co₄Nb₂O₉ compound for comparison and the ⁵⁷Fe doped (Co_0.97⁵⁷Fe_0.03)₄Nb₂O₉ are shown. P3̄c1 space group characteristic of the desired 429-phase Bertaut et al. (1961); Kolodiazhnyi et al. (2011) with a small amount of CoNb₂O₆ impurity 126-phase Zheng et al. (2021), which is usually found in the title compound, were used to refine our XRPD data. The amount of the impurity 126-phase were refined to be about 1.8 wt% and 1.4 wt% for Co₄Nb₂O₉ and (Co_0.97⁵⁷Fe_0.03)₄Nb₂O₉, respectively. In agreement with earlier reports Maignan et al. (2021); Kolodiazhnyi et al. (2011), the determined lattice parameters are $a=5.1702(1)$ and $c=14.1306(4)$ for (Co_0.97⁵⁷Fe_0.03)₄Nb₂O₉, which are a little larger than the values of $a=5.1680(1)$ and $c=14.1268(4)$ for Co₄Nb₂O₉ due to bigger ionic radius of Fe²⁺ than Co²⁺, indicating the successful doping of the ⁵⁷Fe atom into the lattice of the parent 429-phase.

The temperature dependence of the magnetic susceptibility, $\chi(T)$, measured in zero-field-cooling mode shown in Fig. 2 were fitted with the Curie-Weiss law $\chi(T)=\chi_{0}+C/(T+\theta)$ in the paramagnetic regime. The fitted Weiss temperatures are $-72\,K$ and $-71\,K$ and the effective magnetic moments ($\mu_{eff}$) are $5.10\,\mu_{B}$ and $5.15\,\mu_{B}$ for Co₄Nb₂O₉ and (Co_0.97⁵⁷Fe_0.03)₄Nb₂O₉, respectively. The obtained $\mu_{eff}$ values are close to previously reported values and indicate the high spin state of the Co²⁺ ions with $S=3/2$ Khanh et al. (2016). However, these values are significantly larger than the spin only effective magnetic moments of $\sim$3.87 $\mu_{B}$ for $S=3/2$, indicating unquenched contribution from the angular momentum to the effective magnetic moments via strong spin orbital coupling Khanh et al. (2016); Solovyev and Kolodiazhnyi (2016). As shown in the inset of Fig. 2, the Néel transition temperatures were found to be $T_{N}\sim 27$ K and $T_{N}\sim 30$ K for Co₄Nb₂O₉ and (Co_0.97⁵⁷Fe_0.03)₄Nb₂O₉, respectively. The slightly enhanced $T_{N}$ suggests that the relatively small doping of $\sim 3\%$ ⁵⁷Fe is not totally effect-less to the parent compound Maignan et al. (2021) but we expect that the impact should be small with such a small doping level for insulators. The upturn at low temperatures can be ascribed to a weak amount of paramagnetic impurities usually observed in polycrystalline samples that can not be seen by our XRPD and Mössbauer spectroscopy.

The ⁵⁷Fe Mössbauer spectra of (Co_0.97⁵⁷Fe_0.03)₄Nb₂O₉ taken in the paramagnetic temperature range, $T>T_{N}$, shown in Fig. 3 were fitted with two doublets corresponding to the two Co/⁵⁷Fe sites in the crystal structure Bertaut et al. (1961), namely the relatively flat Co/⁵⁷FeO₆ layer, Co/⁵⁷Fe(1)/site 1, and the highly buckled Co/⁵⁷FeO₆ layer, Co/⁵⁷Fe(2)/site 2. The fitted values of the isomer shift (IS) are $\delta_{1}(300\,K)=1.092(1)$ mm/s and $\delta_{2}(300\,K)=1.050(1)$ mm/s, respectively. These values are typical for high-spin ions Fe²⁺ ($d^{6},S=2$) located in FeO₆ octahedrons Gütlich et al. . The values of the quadrupole splitting (QS) are $\Delta_{1}(300\,K)=0.192(4)$ mm/s and $\Delta_{2}(300\,K)=0.889(2)$ mm/s, respectively. To identify the ⁵⁷Fe doping site, we made calculation of the QS for the non-magnetic phase. The EFG tensor was calculated by

where $V^{{}^{\prime}}_{C}$ is the Coulomb potential with the $l=m=0$ component removed in each muffin-tin. The calculated $V_{zz}(1)=-1.03\times 10^{21}$ V/m² and $V_{zz}(2)=-3.60\times 10^{21}$ V/m² with $\eta(1)=\eta(2)=0$. The corresponding QS (with $\Delta=eQV_{zz}/2$) values are $\Delta_{c1}=-0.171$ mm/s and $\Delta_{c2}=-0.599$ mm/s for the flat Co/⁵⁷Fe(1)/site 1 and the buckled Co/⁵⁷Fe(2)/site 2, respectively. These values agree reasonably well with our experimental results.

With decreasing temperature, we observed obvious increasing of the separation between the two peaks of both doublets, corresponding to larger QS values of $\Delta_{1}(35\,K)=0.763(3)$ mm/s and $\Delta_{2}(35\,K)=1.807(4)$ mm/s. These high QS values indicate that the ⁵⁷Fe nuclei are located in crystal sites with stronger EFG, which reflects the asphericity of the charge-density distribution near the probing nucleus. We also noticed that the fitting quality decreases with decreasing temperature with only two simple doublets, which may be caused by several reasons. For example, 1) short range magnetic correlations that appear well above $T_{N}$; 2) cluster relaxation effects due to inhomogeneous doping of the ⁵⁷Fe etc., which will be investigated in the future with more measurements in this temperature range.

As shown in Fig. 4, below $T_{N}\approx 30$ K, complex Zeeman splittings appear in our Mössbauer spectra. To fit these spectra, considering that the QSs are large, we used the full Hamiltonian of hyperfine interactions in the coordinate system of the principal axes ($|V_{zz}|\geq|V_{yy}|\geq|V_{xx}|$) of the EFG tensor Gütlich et al.

where $\hat{I}$ and $\hat{I}_{x},\hat{I}_{y},\hat{I}_{z}$ refer to the nuclear spin operator and operators of the nuclear spin projections onto the principal axes and $Q$ denotes the quadrupole moment of the nucleus. $\phi$ and $\theta$ are the azimuthal and polar angles of the hyperfine magnetic field in the EFG coordinate system, respectively. $\eta=(V_{yy}-V_{xx})/V_{zz}$ represents the asymmetry parameter of the EFG at the nucleus.

It is important to note that, the eigenvalues of such Hamiltonian depend on a number of parameters, such as ($eQV_{zz}$, $B$, $\eta$, $\theta$, and $\phi$), and some of them are correlated. Thus, these parameters cannot be determined independently from Mössbauer spectra, but only in certain combinations. Moreover, the analysis of the Mössbauer spectra can only determine the direction of the magnetic moments in the coordinate system of the principal axes of the EFG tensor. It is necessary to know the relative direction of the EFG coordinate system with respect to the crystal lattice coordinate before we can determine the magnetic structure by fitting of our Mössbauer spectra.

Therefore, before any detailed analysis of our low temperature Mössbauer spectra, we made theoretical analysis based on first-principle electronic structure calculations to study the EFG tensor and other hyperfine interactions. First, we did calculations without considering SOC. To establish the collinear in-plane antiferromagnetic structure as a starting point, an internal magnetic field was applied within the muffin-tin of Co atoms in the [$\pm 1,0,0$] direction which was used to break the symmetries and was reduced to almost zero in the end by setting a reducing factor of 0.95 for each self-consistent field calculation circle. The converged magnetic moments are spin only in this case and amounts to $\mu_{S1}=2.63$ $\mu_{B}$ and $\mu_{S2}=2.59$ $\mu_{B}$ for the two Co atoms. These values are actually close to the reported values from single-crystal neutron diffraction Ding et al. (2020). The calculated $V_{zz}(1)=-8.15\times 10^{21}$ V/m² and $V_{zz}(2)=-8.29\times 10^{21}$ V/m² with $\eta(1)=\eta(2)=0$. Then, the QS was calculated by $\Delta=eQV_{zz}/2(1+\eta^{2}/3)^{1/2}$ to be $\Delta_{1}=-1.36$ mm/s and $\Delta_{2}=-1.38$ mm/s. These values are inconsistent with the values determined at 35 K from above both in magnitude and relative ratio of the two sites which may be an indication of redistribution of the electron charges at the probing sites that happens when going through the antiferromagnetic transition. The magnetic hyperfine field for each atom was calculated according to ref. Blügel et al. (1987) implemented within the ELK code. The fermi contact field were calculated to be $B_{F}(1)=32.0$ T and $B_{F}(2)=30.6$ T. These values are much too large when compared with the magnetic splittings of our Mössbauer spectrum at 4 K as shown in Fig. 4 (a). When the spin dipole part was considered, these values get even larger to be $B_{FD}(1)=36.2$ T and $B_{FD}(2)=33.6$ T. Since the hyperfine field is mainly contributed by three parts as $B_{hf}=B_{F}+B_{D}+B_{L}$, we have to consider the orbital part $B_{L}$ in the present case. In fact, for high-spin Fe²⁺, $B_{L}$ can be as large as $\sim 20$ T and opposite to $B_{F}$ Gütlich et al. .

Therefore, we made calculations with SOC included. In this case, an internal magnetic field was applied within the muffin-tin of Co in the [$\pm 1,0,\delta_{l}=0.057$] direction to break the symmetries, which reduces to nearly zero as in our earlier calculations without SOC. Free rotation of the local moment towards $c$-axis and within the $ab$-plane were both allowed due to Dzyaloshinskii-Moriya like interactions as a result of SOC. Due to this moment rotation, our calculation converges extremely slow, indicating that the SOC plays an important role in determining the ground state of Co₄Nb₂O₉. However, after several thousands of iterations, we found that the total energy changes of the system are very small, $\sim 0.1$ meV. Most importantly, we found that the changes of the EFG tensor is negligibly small, e.g. the change of $V_{zz}$ is smaller than $\sim$1 $\%$ between the last 2000 iterations. From our calculation, the main component ($V_{zz}$) of the EFG tensor is roughly directed along the $c$-axis of the crystal structure, with a maximum deviation angle of $\sim$3.3 ^o for all the 8 Co atoms in the unit cell. Therefore, we fixed the direction of the magnetic moments to the direction of the experimental easy axes of the magnetic moments determined by a tentative fit of the 4.2 K Mössbauer spectrum. In this case, our calculation was easily converged to a level of $\sim 0.27$ $\mu eV$.

From our DFT results, it is reasonable to assume that the direction of $V_{zz}$ is along the $c$-axis in the fitting procedure of our Mössbauer spectra. According to previously reported magnetic structures Khanh et al. (2016); Deng et al. (2018); Ding et al. (2020), we tried to fit our Mössbauer spectra measured below $T_{N}$ with three different models. Model I): canting of the magnetic moments out of the $ab$-plane and free rotation between the two sites are both allowed. Model II): canting of the magnetic moments out of the $ab$-plane is allowed but with a single projected direction in the $ab$-plane for the two sites Khanh et al. (2016). Model III): the magnetic moments are confined within the $ab$-plane but free rotation between the two sites are allowed Deng et al. (2018); Ding et al. (2020). According to the above three models, two sextets were fitted to the experimental data by considering different directions of the hyperfine magnetic field within the coordinate system of the principal axes of the EFG tensor. This is in accordance with the commensurate propagation vector k=0 found in the above three models. However, as shown in Fig. S1 - S3 in the supplemental material, all the above mentioned three models exhibit obvious discrepancies from the experimental data at low temperatures. Actually, only the spectrum taken at 30 K can be reasonably described with a magnetic structure similar to the one proposed by Khanh et el. Khanh et al. (2016). With decreasing temperature, the quality of the fits get worse gradually as indicated also by the $\chi^{2}$ values shown in Fig. S1 - S3. Since our investigated sample was actually doped with ⁵⁷Fe, there might be local randomness that introduced by possible inhomogeneous doping effects. This may give a distribution of the measured hyperfine magnetic field. Indeed, we also tried to fit our Mössbauer spectra using the above models by including hyperfine field distribution effects Rancourt and Ping (1991). But, one can see from Fig. S4 that there is no improvement when compared with the simple two sextets model. Actually, the temperature evolution of the Mössbauer spectra also rules out the possibility that local randomness to be the reason why our fitting got worse and worse with decreasing temperatures. These results suggest that the magnetic structure of (Co_0.97⁵⁷Fe_0.03)₄Nb₂O₉ or even Co₄Nb₂O₉ changes from a nearly collinear simple structure Khanh et al. (2016) at higher temperature to a more complex possibly non-collinear one at lower temperatures.

Therefore, to describe the low temperature spectra satisfactorily, we used a more sophisticated model by taking into account the features associated with an incommensurate helicoidal structure as were used to described the complex Mössbauer spectra of BiFeO₃ Sando et al. (2013), 3$R$-AgFeO₂ Sobolev et al. (2017), Fe₃PO₇ Sobolev et al. (2018), FeP Sobolev et al. (2016) and many others with modulated magnetic structures. In this case, the experimental spectrum is approximated as a superposition of two sets of Zeeman patterns corresponding to the two crystal sites. $\phi$ and $\theta$ angles in equation III for each Zeeman pattern can be expressed by the rotation angle $\omega$ within the magnetic moment rotation plane, which varies continuously in the $0\leq\omega_{i}\leq 2\pi$ interval for an incommensurate helicoidal structure. The angle, $\Theta_{n}$, between the normal direction of the rotation plane and the $c$-axis/$V_{zz}$ was fitted as a free parameter. To take into account of possible anharmonicity (bunching) of the spatial distribution of the magnetic moments of Co²⁺/⁵⁷Fe²⁺, a Jacobian elliptic function was used Sobolev et al. (2017, 2016)

where $K(m)$ is the complete elliptic integral of the first kind, $\lambda$ the helicoidal period, and $m$ is the anharmonicity parameter related to the distortion of the helicoidal structure.

To take account the local field anisotropy that distorts from circular helicoid, we consider the anisotropic hyperfine coupling tensor A. We assume that A is diagonal with respect to the principal axes of the EFG tensor and the hyperfine field can be written as Sobolev et al. (2017)

where $\mu_{x,y,z}$ are projections of the ⁵⁷Fe spin moment on the principal axes of the EFG tensor, and $A_{\perp}\equiv A_{xx}=A_{yy}$, $A_{\parallel}\equiv A_{zz}$ were used since the direction of $V_{zz}$ coincides with the crystal $c$-axis. The hyperfine field component projected onto the field rotation plane can be written as $B_{\omega_{i}}=\mu(A^{2}_{\parallel}cos^{2}\omega_{i}+A^{2}_{\perp}sin^{2}\omega_{i})^{1/2}$. However, in first order, only the component $\textbf{A}\cdot\mathbf{\mu}$ parallel to the spin direction will affect the magnitude of the hyperfine field. Then we can rewrite the expression used in our fitting procedure as $B_{\omega_{i}}(\parallel\mu)=B_{\parallel}cos^{2}\omega_{i}+B_{\perp}sin^{2}\omega_{i}$ Sobolev et al. (2017). In our model, one of the anisotropy direction was assumed to be perpendicular to the plane formed by the easy magnetization direction and the $c$-axis/$V_{zz}$. The easy direction of the magnetic moments for the two sites are described by four angles, namely $\Theta_{ea}(1)$/$\Theta_{ea}(2)$ and $\Phi_{ea}(1)$/$\Phi_{ea}(2)$. During our fitting procedure, we found that using two $\Theta_{n}(1)$/$\Theta_{n}(2)$ angles for the two rotation planes does not improve much of the fitting quality than using only one for both sites. Thus, in the final fits, one single $\Theta_{n}$ for the two rotation planes was used in the fitting model. The asymmetry parameter, $\eta$, were also fixed to the calculated values of $\eta(1)=\eta(2)=0.15$ for both sites. Moreover, due to the small value of $\eta$, we can not accurately determine the value of $\Phi_{ea}$, e.g. the fitted standard deviation is very large. Therefore, we have fixed $\Phi_{ea}$ to the 90^o for both sites.

One can see from Fig. 4 that using the above model allowed us to satisfactorily describe the entire series of our Mössbauer spectra below $T_{N}$. The fitting quality were found to be much better than previous models, that is, the resulting $\chi^{2}$ were reduced to be smaller than $\sim$1.58 for all these spectra. At 4.2 K, the determined EFGs corresponding to the two sites are $V_{zz}(1)=-6.9(2)$ V/m² and $V_{zz}(2)=-9.6(4)$ V/m², which are close with the calculated values ($V_{zz}(1)=-7.89$ V/m² and $V_{zz}(2)=-8.17$ V/m² from our theoretical calculation with SOC. Moreover, from our DFT calculation, the non-zero $\eta$ values were obtained only after SOC was included, indicating that the strong SOC has a considerable effect on the asymmetry distribution of the charges around the probing nucleus which deserves further theoretical studies. The fitted directions for the easy magnetization are $\Theta_{ea}(1)=74(1)$ ^o, $\Theta_{ea}(2)=69(1)$ ^o in the coordinate system of EFG. The normal direction of the rotation plane was found to be $\Theta_{n}=105(1)$ ^o, resulting in angles of 31 ^o and 36 ^o relative to the corresponding easy magnetization directions of the two sites. The fitted anharmonicity parameter $m$ at 4.2 K are close to $\sim 0$ and $\sim$ 1 for site 1 and site 2, respectively. With increasing temperature, $m$ quickly gets close to $\sim$ 1 for both sites. These values suggest that the modulated component of the magnetic moments are nearly within the $ab$-plane at higher temperatures close to $T_{N}$. This explains why the 30 K spectrum can be well fitted with only two subspectra. The corresponding helicoidal contour for the two sites at 4.2 K, projected on the same atoms, are shown in Fig. 5 (a) and (b). Clearly, at ground state, the helicoidal structure is more isotropic for Co/⁵⁷Fe(1) but more anisotropic for Co/⁵⁷Fe(2). With increasing temperature, the magnetic structure eventually becomes almost collinear near the transition temperature, $\sim 30$ K. The spatial anisotropy of the determined hyperfine magnetic field, $\textbf{B}_{hf}\propto\textbf{A}\bullet\mathbf{\mu_{eff}}$, is largely due to the angular dependency of the hyperfine coupling tensor A and/or anisotropy $\mathbf{\mu_{eff}}$ as were discussed in other systems Sobolev et al. (2018, 2017).

The determined average hyperfine field, $\langle B_{hf}\rangle$, from their distribution $p(B_{hf})\propto(|\partial B_{hf}(\omega)/\partial\omega|)^{-1}$ Sobolev et al. (2016) were shown in Fig. 6 as a function of temperature. Solid lines are power-law, $\langle B_{hf}(T)\rangle$ = B₀(1 - $T/T_{N}$)^β, fits to the experimental data in the temperature range between 15 K and 30 K. The fits lead to $B_{0}(1)=9.3(3)$ T, $\beta(1)=0.29(3)$ and $B_{0}(2)=20(2)$ T, $\beta(2)=0.26(6)$ for Co/⁵⁷Fe(1) and Co/⁵⁷Fe(2) sites, respectively. The values of the critical exponent are between the theoretical value of $\sim 0.23$ expected for two-dimensional XY system and $\sim 0.345$ for three dimensional XY system Taroni et al. (2008); Le Guillou and Zinn-Justin (1980). The fitted values of $B_{0}$ are much smaller than the fermi contact field, as mentioned above, but are much larger than the values from our calculation with SOC included ($B_{FDL}(1)=3.4$ T and $B_{FDL}(2)=7.5$ T), reflecting the delicate balance between the two contributions.

Finally, from the above discussed experimental and calculation results, we should try to understand the magnetic structure of the title compound. As discussed above, our fitting of the magnetic structure to the Mössbauer spectra were done in the coordinate system of the principal axes of the EFG tensor. Therefore, we show our calculated directions of the principal axes of the EFG tensor in the left panel of Fig. 7 as green arrows. The directions of the easy axes are shown together as red arrows. The top view of corresponding layers as labeled in the left panel are shown in the right panel. We can see that there is an inter-layer rotation between neighboring layers for the Co(1) site and an intra-layer rotation for the Co(2) site. The calculated angle amounts to $\sim$28^o as shown in the bottom of the right panel. Clearly, the observed modulated helicoidal magnetic structure arises as a compromise of competing interactions, such as antiferromagnetic $J_{1}$ between Co(1)-Co(1), ferromagnetic $J_{2}$ between Co(1)-Co(2), and single ion anisotropy etc., similar to previously reported systems Sando et al. (2013); Sobolev et al. (2018). If we take the easy axis as the average direction of the magnetic moments as shown in Fig. 7, one immediately finds out that it is a combination of previously reported magnetic structures by different groups with neutron diffraction techniques Khanh et al. (2016); Deng et al. (2018); Ding et al. (2020). The canting angles out of the $ab$-plane amounts to 16 ^o for site 1 and 21 ^o for site 2 , determined from our Mössbauer data, compares well with the magnetic structure reported by Khanh et al. Khanh et al. (2016). On the other hand, the in-plane rotation angle calculated in this work $\sim$28^o, which can not be determined by the current Mössbauer data, is much larger than the values reported by Ding et al. Ding et al. (2020). We also want to emphasize that our proposed modulated helicoidal magnetic structure is in contrast with previous models observed by neutron diffractions Khanh et al. (2016); Deng et al. (2018); Ding et al. (2020). This is not surprising to us since the reported magnetic structures are different from group to group even they were all deduced by neutron diffraction techniques. Moreover, these simple magnetic structures can not fully explain the observed complex ME effects Kolodiazhnyi et al. (2011); Fang et al. (2015); Khanh et al. (2016). There are also accumulating evidence that the magnetic structures are much more complex than previously observed with neutron diffractions in similar systems of Fe₄Nb₂O₉ Zhang et al. (2021) and Co₄Ta₂O₉ Choi et al. (2020). Especially, the observed spontaneous electric polarization in the absence of magnetic field in Fe₄(Nb,Ta)₂O₉ Maignan and Martin (2018); Panja et al. (2018); Chen et al. (2021) systems are still lacking of an explanation. Due to the complexity of the actual ground state magnetic structure and the similarity of the average magnetic structure to the reported models, it could be difficult to be refined by neutron diffraction techniques due to domain and powder averaging effects Ross et al. (2015). Thus, our results call for more precise measurements with specially designed neutron diffraction experiments to reveal the true magnetic ground states of this system.

Concluding, we have carried out ⁵⁷Fe Mössbauer measurements on polycrystalline samples of (Co_0.97⁵⁷Fe_0.03)₄Nb₂O₉ to study its possible magnetic structures. Our results show that this compound exhibits very large electric field gradient (EFG), therefore, the principal axes of the EFG tensor can be used as a good coordinate system to solve its magnetic structure by fitting of the measured Mössbauer spectra. The directions of the principal axes ($|V_{zz}|\geq|V_{yy}|\geq|V_{xx}|$) were calculated by density functional theory with spin orbital coupling effect. From these theoretical calculations, we have proposed a modulated helicoidal magnetic structure to simulate the low temperature Mössbauer spectra since all other previously reported magnetic structures failed to describe these Mössbauer spectra. The average magnetic structure, derived from the fitted easy axis direction of the magnetic moments, can be used to reconcile previously reported conflicting magnetic structures. Our proposed non-collinear magnetic structure may be useful in the explanation of the complex magnetoelectric effects observed in this system.

Project supported by the National Natural Science Foundation of China (Grant Nos. 11704167, 51971221). The authors are grateful to the support provided by the Supercomputing Center of Lanzhou University.