Ab initio Modeling of Phonons in the Family of Quasi-One-Dimensional Antiferromagnets AFeX2

Abstract

A review of the recent development in ab initio calculations of the vibrational properties of quasi-one-dimensional (1D) antiferromagnets with AFeX₂ structure is presented. Density functional theory (DFT + U) was applied to calculate the phonon modes specific to each element in the structure and the corresponding partial density of states (PDOS). The calculations revealed a strongly non-Debye phonon DOS. Using these results, the nuclear inelastic scattering spectra, temperature dependence of the Lamb-Mössbauer factor, infrared (IR) absorption strength, and phonon-specific heat were derived by direct summation over the phonon modes. The calculations demonstrate good agreement with the experiments and pave the way for understanding the anomalous magnetic properties of AFeX₂ quasi-one-dimensional antiferromagnets at a quantitative level.

1. Introduction

Compounds and artificial materials with reduced dimensionality attract unflagging attention because of a variety of unusual physical phenomena studied in these systems during the past several decades [1, 2, 3, 4, 5, 6]. The greatest mystery is of course the high-T_c superconductivity (HTSC)—a phenomenon inextricably linked to the reduced dimensionality of the structures where it is observed [7]. An additional, but a strongly indicative feature of copper-based HTSCs, is a strong interrelation between the superconductivity and magnetism in these compounds. Indeed, HTSC can be considered as a state emerging from the parent low-dimensional antiferromagnetic order upon doping by charge carriers. Even more striking relevance and interplay between magnetism and superconductivity are observed in so-called “iron-based superconductors” —iron-containing pnictides, oxypnictides, and chalcogenides (see Refs. [8, 9, 10, 11, 12] and references therein). Depending on composition, stoichiometry, and substitutions, they can undergo a transition into an antiferromagnetic or superconducting ground state (if any). In the lattice structure, iron forms mostly planes with pnictogen or chalcogen ligands, so the reduced dimensionality always accompanies the iron-based superconductivity. The microscopic origin of superconductivity and the role of the crystal structure in these compounds are still under debate. This also concerns the important question about localization or delocalization of the magnetic moments of iron. Therefore, a deeper insight into the nature of superconductivity can be achieved by means of an extensive study of the materials containing building blocks similar to these systems. In this respect, we note that superconductivity is induced by pressure, even in the quasi-one-dimensional (1D) spin-ladder systems BaFe₂X₃ (X = S, Se) composed of FeX₄ tetrahedra building blocks as well [13, 14, 15].

In the past decades, crystallographic studies of ternary metal chalcogenides A_xFe_yX_z, (A = alkali metal, Tl; X = S, Se chalcogen) have revealed various distinct compositions containing tetrahedral [FeX₄] structural units, as there are A₅FeX₄, AFe₂X₂, AFe₂X₃, A₃FeX₃, AFeX₂, and A₃Fe₂S₄ [16, 17, 18, 19, 20, 21, 22, 23]. In particular, the two-dimensional (2D) superconductors K_xFe_2 − ySe₂ and (Tl, K)Fe_2 − xSe₂ exhibit layers of FePn or FeCh (Pn = pnictogen, Ch = chalcogen) tetrahedra as common structural units [24, 25, 26, 27, 28].

The discovery of the superconductivity under high pressure in the ladder compounds BaFe₂S₃ and BaFe₂Se₃ [13, 14, 15] opens a new area of research among the iron-based superconductors’ family. Iron-based chalcogenides AFe₂X₃ (A = K, Rb, Cs or Ba, and X = chalcogen) have been attracting attention due to their unique quasi-one-dimensional crystal structure and magnetism [13, 14, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44]. In this family, FeX₄ tetrahedra share their edges and form two-leg ladders of Fe sites. In other words, the dominant crystal structure comprises pairs of chains, “legs,” with the interchain bonds of a strength similar to that along the legs, dubbed the “rungs.” Note that the magnetic structure of BaFe₂Se₃ [31, 33] is a 1D analog to the block magnetism in A₂Fe₄Se₅ [45]: four Fe spins in the two-leg ladder form a ferromagnetic block while neighboring blocks are coupled antiferromagnetically. In contrast to BaFe₂Se₃, magnetic structures of BaFe₂S₃ and AFe₂Se₃ (A = K, Rb, Cs) [13, 32, 34, 35] are of the stripe type, in which the magnetic moments couple ferromagnetically along the rung and antiferromagnetically along the leg. It is remarkable that this family of materials is not layered, like all other iron-based superconductors, but instead has the geometry of a two-leg ladder. It has been shown experimentally that BaFe₂S₃ becomes superconducting at pressures above 10 GPa and with the highest critical temperature of T_c = 24 K [14].

In addition to the already mentioned AFe₂X₃ (A = Ba, Rb, K, Cs; X = S, Se) compounds, there is a group of iron-chalcogenide materials AFeX₂ and A₃Fe₂S₄ (A = K, Rb, Cs, Tl; X = S, Se) which also possess 1D structures, though in this case simple chains as opposed to ladders. For example, Na₅FeS₄, A₃FeX₃ (A = Na, Cs; X = S, Se), AFeX₂ (A = K, Cs; X = S, Se), and A₃Fe₂S₄ are built from discrete tetrahedral [FeS₄]^5-complexes [17], edge-sharing double-tetrahedral [Fe₂X₆]^6-complexes [20, 46, 47], or 1D linear and zigzag _∞¹[FeX_4/2]^-chains [20, 48, 49], respectively. The evaluation of magnetic susceptibility data yielded a systematic reduction in the localized iron-spin value from S = 5/2 of discrete tetrahedra via S = 3/2 in double-tetrahedral complexes and down to S = 1/2 in linear-chain compounds [17, 20, 46, 47, 48, 49, 50].

In the present chapter, we focus on the quasi-1D, linear-chain, ternary iron chalcogenides AFeX₂ (A = K, Rb, Tl; X = S, Se) [51, 52, 53, 54]. The crystal structure is shown in Figure 1. The common motif for all of these compounds is the presence of tetrahedral [FeX₄] structural units, arranged edge-sharing into linear chains along the c direction and separated by the A atoms. Although superconductivity was not observed in these linear-chain compounds so far, they can be regarded as paradigm model systems, because the short Fe-Fe distance within the chains drives strong covalence effects. As a consequence, spin reduction and delocalization of the charge carriers are expected at the verge of 1D metallic conductivity. All the experimental results presented below have been measured on high-quality single crystals, the synthesis of which has been mastered in Augsburg and Kishinev by Dr. Tsurkan’s research group.

Figure 1.

Crystal structure of AFeX₂. The FeX₄ tetrahedra, with Fe (red sphere) in the center and X (yellow sphere) at the corners, are highlighted in a transparent yellow color. Blue large spheres denote A (adapted from Ref. [55]).

Crystallographic and magnetic data of AFeX₂ (A = K, Rb, Tl; X = S, Se) are summarized in Table 1. With respect to their magnetic properties, these iron chalcogenides can be divided into the following two groups:

1. monoclinic TlFeS₂, TlFeSe₂, KFeSe₂, and RbFeSe₂ with the magnetic moments ordered perpendicular to the chains; and

2. monoclinic KFeS₂, RbFeS₂ with the magnetic moments slightly tilted from the chain axis in the ordered state.

Compound	Space group	d (Å) (Fe-Fe)	TN (K)	μord (μB)	Moment orientation	Reference
KFeS2	C2/c	2.70 uniform	250	2.43	13° from chain	[20, 48, 50, 56, 57]
RbFeS2	C2/c	2.71 uniform	188	1.83	║ chain and slightly tilted	[20, 48]
TlFeS2	C2/m	2.65 dimerized	196	1.85	┴ chain	[21, 57, 58, 59]
KFeSe2	C2/c	2.81 uniform	310	3	┴ chain	[20, 48]
RbFeSe2	C2/c	2.83 uniform	250	2.66	┴ chain	[20, 48, 55]
TlFeSe2	C2/m	2.74 uniform	290	2.1	┴ chain	[58, 60]

Table 1.

Crystallographic and magnetic data for the AFeX₂ (A = K, Rb, Tl; X = S, Se) (adapted from Ref. [55]).

In all these compounds, iron is formally trivalent with 3d⁵ electronic configuration. However, the ordered magnetic moment is found to be far below the formal ionic high-spin value of 5μ_B, which indicates significant 3d delocalization. This has been explained by the intimate Fe-Fe contact due to the small intra-chain Fe-Fe distance, which is not much larger than the Fe-Fe distance (2.48 Å) in metallic iron. Therefore, along the chains a certain degree of itinerancy and even 1D conductivity can be expected.

Concerning the spin susceptibilities, Faraday magnetic balance and SQUID (superconducting quantum interference device) measurements of K, Rb, and Tl compounds AFeX₂ (A = K, Rb, Tl; X = S, Se) revealed the magnetic susceptibility that increases with raising the temperature (to a magnitude between 4 × 10⁻⁴ emu/mol and 2 × 10⁻⁴ emu/mol at room temperature) for rubidium (Rb) and thallium (Tl) compounds [55, 58] and passing through the maximum for the potassium compound [50]. The behavior of the susceptibility with temperature maximum is typical for one-dimensional (1D) antiferromagnetic (AFM) Heisenberg chains [61, 62, 63, 64], with T_max = 565 K for KFeS₂, whereas RbFeSe₂ and TlFeX₂ (X = S, Se) [55, 58] compounds exhibit a continuous linear increase (up to 720 K, 550 K, and 500 K for RbFeSe₂, TlFeS₂, and TlFeSe₂, respectively) and do not show any tendency for saturation. This kind of a linear increase of the susceptibility is rather unusual for 1D antiferromagnetic (AFM) Heisenberg chains of localized spins, which typically exhibit a susceptibility maximum at a temperature comparable to the intra-chain exchange [61, 62, 63, 64]. Such susceptibility maxima had been approximated using the model of an AFM spin S = 1/2 chain [50]. On the other hand, for RbFeSe₂, TlFeS₂, and TlFeSe₂, the quasi-linear susceptibility increase above T_N without any tendency to saturate [55, 58] can hardly be described within the spin-chain model. It seems that in these compounds with the strong direct Fe-Fe exchange, a fraction of d-electrons is close to delocalization.

The specific heat C_P(T) has been measured in several works [33, 34, 37, 65, 66, 67, 68] in the temperature range from liquid helium (L_He) up to room temperature (RT). Measurements reveal a weak lambda anomaly at T_N, sometimes very tiny if visible at all. As a first approach, the temperature dependence of C_P(T) was fitted within classical phenomenological models including Debye and Einstein terms. The difference between the measured and the fitted C_P(T) was considered as the magnetic contribution to the specific heat C_M(T). The entropy changes associated with the magnetic ordering transition were estimated from the area under the lambda anomaly in C_P/T but have shown unrealistically small values, typically below 1 JK⁻¹ mol⁻¹. At the antiferromagnet-paramagnet (full order-disorder) transition, the expected entropy change is close to Rln2 = 5.76 JK⁻¹ mol⁻¹ if the iron ion is in the low-spin S = 1/2 state, which is an order of magnitude larger than the one obtained by the routine described above. Some authors, anyway, conclude that the small entropy change at the magnetic transitions is an indication of the low S = 1/2 spin state of the iron ion in AFeX₂ ternary chalcogenides.

The close connection between magnetic and superconducting properties in widely studied iron-containing superconductors does not allow one to study these phenomena separately [69], therefore, the study of the magnetic subsystem of the low-dimensional iron compounds is an actual problem. Two of the important quantities characterizing magnetic systems are the magnetic heat capacity and magnetic entropy. Determination of these quantities and their temperature dependences is necessary for constructing theoretical models of interactions, magnetic excitations, and superconductivity in the system.

Summarizing the introduction into the subject matter of AFeX₂ compounds, we note below a list of contradictions that motivated our studies:

• the magnetic moment per iron ion in the range of 1.83–3.0 μ_B is considerably reduced with respect to the expected one from ionic Fe³⁺ spin S = 5/2 configuration (∼5.0 μ_B). At the same time, the low-spin S = 1/2 state by no means fits the observed magnitude of the magnetic moment;

• the values of the hyperfine field at the Fe nuclei obtained by Mössbauer spectroscopy measurements are also reduced and do not fit either the high-spin or the low-spin state of the iron ion;

• the increase of the magnetic susceptibility with raising the temperature is unusual (the Curie-Weiss descending behavior could be expected), however, for quasi-1D spin chains the Bonner-Fischer model predicts an increase of the susceptibility in a certain temperature range starting from the origin. At the moment, it is the S = 1/2 spin-chain model that can be applied to the particular case with great care, maybe as a hint to a possible approach to resolve the problem;

• the entropy change at T_N, estimated by phenomenological approaches, is below the lowest limit for S = 1/2. This indicates the insufficiency of the standard treatment of the temperature dependence of specific heat and the presence of strong fluctuations in the spin system within the entire temperature range because of the quasi-1D character of these compounds;

• ab initio calculations were not applied to AFeX₂ ternary chalcogenides before. Nowadays, it is a powerful instrument to study all kinds of physical properties of matter with reduced dimensionality.

2. Density functional theory approach for the phonon density of states

Density Functional Theory (DFT) calculations are a powerful tool to calculate the phonon density of states. This is achieved by combining DFT with phonon calculations, which are typically performed using the so-called supercell approach. This method allows for the theoretical prediction of phonon frequencies at various points in the Brillouin zone, from which the phonon density of states (PDOS) can be derived. Moreover, the PDOS provides comprehensive important information about the vibrational states of a material, which is crucial for understanding various physical properties, such as absorption spectra or specific heat.

We should note that, as a general rule, DFT calculations underestimate and inaccurately predict the forces between atoms due to certain limitations inherent in the approximation methods used for the exchange-correlation functional. The exchange-correlation functional is a crucial component in DFT that describes the quantum mechanical interactions among electrons, including exchange and correlation effects. Especially, in the case of compounds with electronic d-states, this is important. However, finding an exact representation of this functional is a significant challenge, and approximations must be made.

For instance, commonly used approximations, such as the Local Density Approximation (LDA) or the Generalized Gradient Approximation (GGA), can lead to inaccuracies in describing a system where electronic correlations play a significant role, such as in transition metal complexes or systems with strong covalent bonding. So, both of these complicating factors are native to the compounds under consideration because at least they contain iron ions.

Typically, in order to improve the accuracy of force predictions, researchers use more advanced functionals or add specifically designed corrections, such as the DFT + U methods or the use of hybrid functionals that mix DFT with Hartree-Fock exchange. These methods can significantly improve the accuracy of DFT predictions, including the forces between atoms.

DFT + U is an extension of DFT that is designed to improve the description of strongly correlated systems, particularly those containing transition metals, rare earth elements, or systems where d or f electrons play a significant role. In the cases of such materials, DFT + U can provide a more accurate description of electronic structures, magnetic properties, and the geometry of molecules and solids than standard DFT. It captures better the on-site electron-electron repulsion in these systems, which simple DFT underestimates. The U term in DFT + U refers to an added on-site Coulomb interaction term that corrects for the self-interaction error often seen in standard DFT calculation. The additional parameter U needs to be chosen appropriately for the system under study.

While it can lead to more accurate predictions for certain properties in correlated systems, the need for these empirical parameters can be seen as a drawback compared to parameter-free DFT methods. While DFT + U is a valuable tool for improving the accuracy of DFT calculations in systems with strong electron correlation effects, the choice of U is crucial and can significantly influence the results; thus, it requires careful calibration or selection based on experimental data.

So, we can conclude that the exact estimation of the interatomic forces and perfect representation of the crystal structure are two main necessary conditions to correct prediction of the phonon density of states of a compound.

Doubtless, the modern X-ray diffraction measurement techniques provide crystallographic structure data of high accuracy, even in the case of typical, generally available commercial equipment.

Conversely, the accuracy of DFT calculated forces between atoms might need meaningful efforts and application of front-line experimental methods. In our studies, we have used three ways for calibration of interatomic forces by using infrared (IR) spectroscopy, Mössbauer spectroscopy, and the method of nuclear inelastic scattering of synchrotron radiation. The detailed description of how to use such experimental data to calibrate the calculated phonon density of states is given below. Leaping ahead, we note that, as it is shown below, such high-end experimental approaches by turns need to be and might be improved by their combination with advanced DFT calculation methods.

All the ab initio DFT calculations were performed by means of the widely used Vienna ab initio simulation package (VASP 5.3) [70, 71, 72, 73]. The electron-ion interactions were considered using the projector-augmented wave (PAW) method. This is a frozen-core method that employs the exact wave-functions of the valence electrons, instead of pseudo-wave functions [74]. Exchange and correlation corrections were taken into account by means of the Perdew-Burke-Ernzerhof (PBE) generalized gradient approximation (GGA) [75]. The detailed description of calculated parameters, where used, is given in the corresponding cited publications. Here, we just note that the Medea-Phonon software was used to obtain the phonon dispersion and density of states in harmonic approximation [76]. In the so-called direct approach to the lattice dynamics, the forces acting on all atoms are evaluated by a set of finite displacements of a few atoms within an otherwise perfect crystal. The phonon dispersion was calculated using the totally optimized equilibrium structure. In all calculations, the spin polarization due to the antiferromagnetic order of the compounds was taken into account, using the antiferromagnetic spin pattern previously obtained by neutron diffraction [48].

3. Experimental methods

Single crystals of all compounds under consideration were grown by the Bridgman method. The structural details of the crystals were investigated by conventional X-ray diffraction on powdered single crystals at room temperature using a STADI-P (STOE & Cie GmbH, Germany) diffractometer with CuK_α radiation. The data were analyzed by standard Rietveld refinement using the program FULLPROF [77]. For each compound, we could not detect any impurity phases above the background. The structural analysis confirmed the monoclinic structure for all systems and gave lattice parameters similar to those presented in previous studies [48].

The specific heat was measured by a relaxation method using a commercial physical properties measurement system PPMS-9 (Quantum Design, USA) in the temperature range 1.8 ≤ T ≤ 300 K and in magnetic fields up to 90 kOe.

All the Mössbauer spectra presented below were measured with an accuracy of ±0.1 K over the whole temperature range on a conventional constant-acceleration spectrometer (WissEl, Germany) equipped by a room-temperature rhodium-matrix cobalt-57 gamma-radiation source. Preliminarily, the spectrometer was calibrated at room temperature with an α-iron foil. For the measurements, the absorber was prepared by rubbing the crystals under consideration on a glue tape and packing them into a holder closed with thin aluminum foil. The preparation procedure was done in an inert argon-atmosphere dry box.

The nuclear inelastic scattering [78, 79] experiments were performed at the Dynamics Beamline P01 of PETRA III synchrotron (DESY, Germany) [80]. Utilizing the nuclear gamma resonance of ⁵⁷Fe at 14.413 keV, the measurements were carried out with an inline high-resolution monochromator providing an energy bandwidth of 0.9 meV full width at half maximum (FWHM). The samples with natural ⁵⁷Fe isotope content were measured at a stabilized temperature of 295 K.

Infrared (IR)-absorption spectra were obtained in the wavenumber range of 100–600 cm, which corresponds to the frequency range of 3–18 THz, by using an IFS 113v (Bruker, Germany) FTIR spectrometer.

4. Calculated and experimental results and their interlinked synergy

4.1 Element-specific phonons in RbFeSe₂ and their density of states. DFT approximation

The recital of our results will be given by starting with RbFeSe₂ and continued by the related compound KFeS₂. This kind of order is mainly conditioned just by the chronology of the studies and does not try to depict any sophisticated rule.

Figure 2 presents the DFT computed phonon density of states of RbFeSe₂. Each vibrational mode is simultaneously Raman and infrared active and belongs to one of the two irreducible representations of the C₂ point-group symmetry corresponding to the RbFeSe₂ structure. The vibrational modes of the selenium atoms exhibit relatively low frequencies and provide the major contribution to the frequency range of 1–3.5 THz. Iron atoms possess the highest oscillation frequencies, and their vibrations contribute a substantial part to the high-frequency range of 6.5–9.5 THz. The vibrational modes of the selenium atoms are distributed over the full frequency range up to 9.5 GHz, with maxima in the PDOS both in low- and high-frequency regimes. The Debye-type acoustic PDOS of the rubidium atoms without any optical contributions indicates that the Rb atoms oscillate almost decoupled from the iron atoms in the chain. Such behavior is typical for weakly bound Rb atoms.

Figure 2.

Phonon density of states in RbFeSe₂, calculated by DFT as a function of frequency: element specific for Rb, Fe, and Se atoms (from the top toward the bottom) and in total (bottom), adapted from Ref. [51].

As we noted above, it is known that the DFT approach usually underestimates vibrational frequencies. The suitable correction factor might be estimated by comparison of calculated vibrational frequencies and exact periodicity parameters of the oscillation modes obtained by an experimental technique. In our particular case, the frequency-scale renormalization for the phonon modes has been done in three different ways:

1. by comparison of the experimentally measured temperature dependence of the Lamb-Mössbauer factor with the one calculated from the iron atoms’ partial PDOS in accordance with the known relation between the Lamb-Mössbauer factor and the mean-square displacement of the resonant ion, while the last one is directly expressed via the PDOS. To accomplish that, the approximation procedure was done with the only fitting parameter, which is the oscillation frequency scale recalibration factor;

2. by combined assay of the experimentally measured IR-absorption spectrum of the compound and the corresponding PDOS since high-frequency phonons are optically active according to the calculations;

3. by direct comparison of calculated and experimentally measured partial PDOS of iron atoms obtained with nuclear inelastic scattering.

4.1.1 Calibration of RbFeSe₂ phonon frequency scale

To use the first option, we start from the initial definition of the Lamb-Mössbauer factor as the temperature dependence of the ratio of recoil-free to total nuclear resonant absorption given by [81, 82]:

$$f_M\left(T\right)=exp\left\{-\frac{⟨x^2⟩E_{\gamma }^2}{{\left(\mathrm{\hslash }c\right)}^2}\right\},$$E1

where E_γ denotes the gamma-quantum energy, c the velocity of light, and ħ the Planck constant. The value of the Lamb-Mössbauer factor directly establishes the Mössbauer spectral area for the spectrum of the corresponding iron atoms. And both of them can be considered like related by multiplier quantities.

The mean-square displacement <x²> is given by the expectation value of the squared vibrational amplitude of the iron ions along the propagation direction of the gamma radiation. Its temperature dependence can be expressed via the phonon density of states [82] as:

$$⟨x^2⟩=\frac{\mathrm{\hslash }}{2M}\int \frac{1}{\omega }coth\left(\frac{\mathrm{\hslash }\omega }{2k_{B}T}\right)g_{\mathit{Fe}}\left(\omega \right)\mathit{d\omega },$$E2

with the phonon frequency ω, the partial PDOS g_Fe of the Mössbauer nuclei, the (effective) mass M of the Fe nucleus, and the Boltzmann constant k_B. Then, we calculate the temperature dependence of the Mössbauer spectral area, Eq. (1), by numerical integration of Eq. (2) inserting the PDOS of iron g_Fe(ω), given in Figure 2.

Figure 3 shows the adjustment of the calculated to the experimental temperature dependence of the Mössbauer spectral area. In the analysis, the oscillations of the iron ions along X, Y, and Z directions have been taken into account without weighting, because the X-ray diffraction data did not reveal any significant texture of the samples. The bare atom mass of iron has been inserted for M. The only fitting parameter was the scale-correction factor for the frequency as discussed above. The best agreement between the calculation and experiment has been achieved by adjusting the frequency scale by a factor of 1.08.

Figure 3.

Temperature dependence of the Mössbauer spectral area of RbFeSe₂: red solid symbols – experiment; blue solid line – calculation utilizing the iron partial PDOS (adapted from Ref. [51]).

To do the second option we compared, as presented in Figure 4, the experimental IR-absorption spectrum with the calculated PDOS, where the frequency scale was corrected by the factor of 1.08, and we also performed the fitting procedure to find the best agreement between the DFT and experimental data. The IR-absorption spectrum shows the best agreement with the calculated PDOS in the case of a scaling correction factor of 1.06 (solid bar chart). The 2% difference in the scaling-factor values obtained with Mössbauer and optical experiments we qualitatively explain by possible slight crystal structural distortions that result from surface oxidation of the compound during IR-absorption spectroscopy measurements, giving rise to a corresponding phonon-frequency shift. As the IR-absorption measurements were performed in ambient air atmosphere, they should be more sensitive to modifications of the surface than the Mössbauer measurements, which were conducted under protective argon atmosphere, both the preparation of the absorber and the measurements. Nevertheless, such a proximity of these two scaling-factor values anyway appears to make these results quite persuading.

Figure 4.

Infrared-absorption spectrum of RbFeSe₂ compared with the calculated phonon density of states of the compound. The frequency scale of the PDOS is corrected in accordance with the Mössbauer data (dashed bar charts) and to get the best coincidence of the peaks of PDOS maxima with the IR absorption maxima (solid bar charts) (adapted from Ref. [51]).

Now, we begin to consider the results obtained with using the high-end experimental technique based on the nuclear inelastic scattering (NIS) of synchrotron radiation on resonant iron nuclei. The Fe nuclear inelastic scattering spectrum of RbFeSe₂ is shown in Figure 5 (blue solid hexagons) along with the instrumental function measured simultaneously with the spectrum (green solid hexagons in Figure 5).

Figure 5.

⁵⁷Fe nuclear inelastic scattering spectrum of RbFeSe₂ (blue hexagons) and instrumental function (green hexagons) (adapted from Ref. [52]).

The partial iron PDOS of RbFeSe₂ was determined from the NIS spectrum following the procedure described in Ref. [83] and is depicted in Figure 6 (blue open hexagons). The magenta solid line in Figure 6 shows the DFT calculated partial PDOS for the iron atoms in RbFeSe₂ recalibrated by using the correction factor of 1.08. The PDOS calculated by DFT consists of a discrete set of contributions from phonon modes, while the NIS-derived PDOS is a quasi-continuous sequence due to the settings of the beamline setup, but mainly due to the energy resolution of the monochromator, which amounts to 0.9 meV. To compare both with each other, the DFT calculated partial PDOS of the iron atoms was convoluted with a Gaussian profile of 0.9 meV FWHM, simulating the resolution of the monochromator. The result of the convolution was finally recalibrated as mentioned above yielding the solid magenta line presented in Figure 6.

Figure 6.

Partial phonon density of states for iron atoms of RbFeSe₂ obtained from the nuclear inelastic scattering spectrum (blue open hexagons) and from ab initio calculations (magenta line, see description in the body text) (adapted from Ref. [52]).

Thus, Figure 6 proves excellent agreement between the ab initio calculated PDOS and the NIS-derived PDOS. In the entire frequency range of lattice vibrations, all features of the iron PDOS in RbFeSe₂ are quantitatively described by the calculation. The optimal frequency correction factor of 1.077 ± 0.002 is practically identical to that obtained previously from the analysis of the temperature-dependent Lamb-Mössbauer factor. Such accordance corroborates our previous results and ensures that the temperature dependence of the Lamb-Mössbauer factor enables an accurate calibration of the DFT calculated PDOS of solids.

So, three independent experimental techniques done on different sample pieces, performed with various equipment, including both the commercial one and the unique high-end measuring setup in DESY Centre, gave almost identical results. We also note that, the last comparison of DFT results and data obtained by nuclear inelastic scattering experiments prove the excellent quantitative agreement between them. It is clearly seen that all the features of the iron PDOS are well reproduced within the entire frequency range of vibrations in the RbFeSe₂ lattice. So, for further quantitative analysis of the thermodynamic properties, we use the frequency scaling factor of 1.08.

As an optional side-conclusion, we need to pay attention to the fact that the approximation of the temperature dependence of the Lamb-Mössbauer factor by using the calculated PDOS and the frequency correction factor as the only fitting parameter gave the same frequency scale calibration correction as was obtained by using the high-cost and barely available nuclear inelastic scattering method based on synchrotron radiation. Despite the own difficulties of using the conventional Mössbauer spectroscopy related to operating gamma-radiation sources, this method, beyond doubt, is more accessible than nuclear inelastic scattering measurements depending on synchrotron radiation sources. Thus, both techniques can be equally applied for evaluation of the frequency correction factor for the DFT calculated phonon frequencies.

4.1.2 Temperature dependence of RbFeSe₂ specific heat and its analysis

In the case of RbFeSe₂, the total specific heat is expected to originate from two contributions, (i) the lattice contribution due to acoustic and optical phonons as discussed above and (ii) an additional magnetic contribution determined by the thermal population of excited magnetic states. The phonon density of states enables to directly calculate the lattice contribution to the specific heat by using the harmonic approximation [84]. In the harmonic approximation, the lattice heat capacity can be determined as follows [85]:

$$C_V\left(T\right)={\mathit{DNk}}_B\int {\left[\frac{\mathrm{\hslash }\omega /2k_{B}T}{sinh\left(\mathrm{\hslash }\omega /2k_{B}T\right)}\right]}^2{g}_T\left(\omega \right)\mathit{d\omega },$$E3

where D is the number of degrees of freedom in the unit cell (3 in our case) and g_T(ω) is the total PDOS (Figure 2, bottom panel). Traditionally experimental works [33, 34, 37, 65, 66, 67, 68] present the temperature dependence of the specific heat at constant pressure C_P(T) because the specific heat is normally measured at ambient conditions in a wide temperature range so that the thermal expansion is not negligible.

Since the general thermodynamic approach gives the following relation:

$$C_P\left(T\right)-C_V\left(T\right)=-T{\left[{\left(\partial V/\partial T\right)}_P\right]}^2/{\left(\partial V/\partial P\right)}_T>0,$$E4

therefore, the specific heat at constant pressure C_P(T) is always larger than the one measured at constant volume. In the case of solids, the subtraction result between these two specific heat expressions is given by a simple relation:

$$C_P\left(T\right)-C_V\left(T\right)=-T\frac{{\alpha }^2{V}_0^2}{{\mathit{dV}}_0/\mathit{dP}}={\alpha }^2{\mathit{BV}}_{0}T,$$E5

where α ≡ α(P) denotes the thermal expansion coefficient, V₀ represents the molar volume, and B is the bulk modulus. The bulk modulus of RbFeSe₂ has been estimated directly from the second derivative of the total energy as a function of the unit-cell volume as B = 17.1 GPa. The thermal expansion coefficient α was estimated utilizing the crystal-structure data presented for room temperature and temperature of 14 K in Ref. [48]. This allowed us to express the specific heat at constant pressure as C_P ≈ C_V + 0.008[J/mol K²]·T. Figure 7 compares the calculated lattice heat capacity to the experimentally obtained total specific heat. The difference between the lattice specific heat and the total one is associated with the magnetic contribution to the specific heat and presented in Figure 8. The corresponding magnetic entropy change can be calculated from the experimentally measured magnetic specific heat as

Figure 7.

Temperature dependence of the specific heat C_P of RbFeSe₂: black circles – experimental data, red line – calculated lattice contribution to the specific heat at constant pressure C_P=C_V + 0.008[J/mol K²] T (adapted from Ref. [51]).

Figure 8.

Temperature dependence of the magnetic contribution to the heat capacity C_M(T) (black empty circles) and magnetic entropy change (red empty squares) of RbFeSe₂ obtained as the difference between the experimentally measured specific heat and the calculated lattice contribution (adapted from Ref. [51]).

$$\mathrm{\Delta }{S}_M=\int \left[C_m\left(T\right)/T\right]\mathit{dT}.$$E6

The integration from zero temperature to 290 K gives the value of ΔS_M = 6.03 JK⁻¹ mol⁻¹, which should be considered as an estimation of the lower boundary, because of the temperature limit (T < 290 K) of our equipment.

Note that in quasi-one-dimensional antiferromagnets, magnetic fluctuations are expected to persist up to high temperatures half of the Néel temperature and more above T_N. The experimental entropy change has to be compared to the theoretical values ΔS_M expected at the antiferromagnetic to paramagnetic order-disorder transition for the three possibilities of low spin 1/2, intermediate spin 3/2, and high spin 5/2 for Fe³⁺ (3d⁵) ions, given by Rln2 = 5.76 JK⁻¹ mol⁻¹, Rln4 = 11.52 JK⁻¹ mol⁻¹, and Rln6 = 14.89 JK⁻¹ mol⁻¹, respectively. Thus, for RbFeSe₂, the experimentally obtained value of the magnetic entropy change ΔS_M comes closest to the case of the intermediate-spin state S = 3/2.

In order to probe the local magnetic state of the iron ions in RbFeSe₂, the low-temperature Mössbauer spectroscopy data have been carefully analyzed and the hyperfine field at the ⁵⁷Fe nucleus was determined as H_0hf = 216 kOe. This value conforms well with Mössbauer results for the related alkali-metal and thallium-iron sulfides and selenides. Such a low value of the hyperfine field indicates a strong reduction of the local iron-spin moment. This hyperfine-field value is significantly smaller than that (H_0hf ∼ 500 kOe) typical for the high-spin state of Fe³⁺ ions, but too large for the low-spin state.

Therefore, the experimentally determined hyperfine field supports the intermediate-spin state S = 3/2 of the iron atoms in RbFeSe₂ in agreement with the set of data for ΔS_M. Following Ref. [48] for edge-sharing FeSe₄ tetrahedra in RbFeSe₂, the S = 3/2 intermediate-spin state results from the distortion of the tetrahedra and the related splitting of the Fe³⁺ 3d-orbital triplet and doublet by the corresponding low-symmetry crystal field. Moreover, the partial delocalization of the iron 3d-electrons in RbFeSe₂ is also responsible for reduced values of hyperfine field and magnetic moment, as compared to those values typical for the S = 3/2 case.

The most significant magnetic contribution to the specific heat appears in the temperature range 30–200 K. Indeed, because of the quasi-one-dimensional structure of RbFeSe₂, one could expect that quantum fluctuations, which otherwise prevent magnetic ordering in ideal 1D spin chains according to the Mermin-Wagner theorem [86, 87], persist well below the three-dimensional (3D) ordering temperature T_N = 248 K. From theoretical point of view, both for 1D Heisenberg- [spin 1/2 to 3] or XY- [spin 1/2] antiferromagnetic spin chains with nearest-neighbor intra-chain exchange J, numerous papers report an approximately linear temperature dependence of the magnetic contribution to the specific heat at temperatures T < J/k_B. For quasi-one-dimensional antiferromagnets with interchain exchange J’ ≪ J, 3D antiferromagnetic order typically sets in at about T_N ∼ (J × J’)^1/2/k_B. This means that the temperature corresponding to the intra-chain exchange energy is expected to be much larger than the Néel temperature, and therefore the temperature range covered by our measurements should undoubtedly match the regime where the quasi-linear behavior of the magnetic specific heat is predicted. However, we have to note here that usual spin-chain compounds, where existing theories are applicable, undergo magnetic order at much lower temperatures and hence, to our best knowledge, there is no appropriate theory to treat our present particular case of magnetic specific heat of quasi-one-dimensional spin-3/2 chains with the interchain interaction triggering the transition to the 3D antiferromagnetic state at a Néel temperature as high as T_N = 248 K.

4.2 Element-specific phonons in KFeS₂ and their density of states. DFT + U approximation

Further we continue with the compound KFeS₂, for which the same study strategy was used. Figure 9 depicts the calculated total and element-specific PDOS as a function of frequency for KFeS₂ obtained by means of the same DFT approach we used above for the RbFeSe₂ compound. The PDOS of KFeS₂ is similar to that of RbFeSe₂, it could be divided into two distinct areas: the low- and high-frequency ones. The low-frequency area extends up to 6 THz, which shows Debye-like frequency dependence (approximately up to 3 THz), and the dominant contribution to this area comes from potassium atoms. The high-frequency area shows a significant amount of phonon density of states and consists of only vibrational modes of iron and sulfur atoms. Because of the large number of high-frequency vibrational modes of Fe and S atoms, the Debye model cannot be applied for a fine analysis of the vibrational properties of the system such as lattice specific heat. Nevertheless, this result obtained by means of a simple DFT approach yielded a phonon density of states, which does not perfectly agree with the experimental results we got for KFeS₂. Indeed, a significant disagreement was found.

Figure 9.

The DFT calculated phonon density of states in KFeS₂: element-specific (K, Fe, and S atoms from the top toward the bottom) and the total PDOS (the bottom).

4.2.1 KFeS₂ PDOS and calibration of the phonon frequency scale

The above-mentioned experimental results for KFeS₂ are presented in Figure 10: (i) the partial iron phonon density of states is presented in the mainframe of Figure 9, and the inset gives the experimental infrared-absorption spectrum. As we noted above, the IR-absorption spectrum measures only the optical phonon modes of all three atoms in the compound, while the NIS technique probes both acoustic and optical phonon modes, but only those of the iron atoms. In the low-frequency range, the iron partial PDOS follows the cubic frequency dependence predicted by the classical Debye model for acoustic phonons. Notably, both NIS and IR measurements reveal two maxima in the optical phonon density at about 9 and 11.5 THz. Such accordance of both independent methods corroborates the high accuracy of the experimental data, while the PDOS calculated by means of the simple DFT-approach is not in accordance with them.

Figure 10.

Partial iron PDOS of KFeS₂ evaluated from the NIS spectrum of KFeS₂ and IR-absorption spectrum (presented in the inset) (adapted from Ref. [53]).

We should also note that the magnetic moment of iron ions in KFeS₂ obtained within the simple DFT-approach does not match well the experimental value obtained by neutron diffraction. This probably indicates strong delocalization of the Fe d-electrons and hence a proper description of KFeS₂ demands to apply the DFT + U approach. The DFT + U method strongly improves the description of materials with mixed type of electron distribution. This material class includes mid-to-late first row transition metal oxides and sulfides, because the self-interaction error is largest for such multiply charged (i.e., +2 or higher) open-shell metal ions [88]. DFT + U theory does require selection of the only parameter, U_eff = U – J, the difference between intra-atomic Coulomb and exchange energies [89]. In a simplified DFT + U approach, these two parameters are not taken into account separately, but their difference is used. Sometimes the interatomic exchange energy is omitted from consideration and the only parameter U is used [90]. So, the value of the parameter U_eff should be preliminarily estimated from available experimental data to ensure the correct use of any DFT + U calculations of any properties of KFeS₂. Here, we estimate U_eff based on experimentally determined electronic and magnetic properties. Note that in case of KFeS₂, both parameters, U and J, have to be considered independently.

In case of magnetically ordered systems, the proper choice of U_eff is indispensable for a correct description of the magnetic ground state and calculation of the magnetic moments [91]. Thus, comparing ab initio results with experimental data allows to find the correct U_eff value just by simple selection. The magnetic ground state and ordered magnetic moment of the Fe³⁺ ions are well known from neutron diffraction measurements [55]. Below T_N = 250 K, KFeS₂ shows antiferromagnetic order with an ordered magnetic moment of 1.9 μ_B [48]. The electronic band gap below T_N was estimated from the temperature dependence of the electrical resistivity of quenched KFeS₂ crystals [92]. To find the proper values of the parameters U and J, we have performed a series of ab initio calculations of the magnetic moment and electronic band gap of KFeS₂ with different values of the parameters. The on-site Coulomb terms U and J have been considered independently for the strongly correlated Fe d-electrons. Both parameters have been adjusted to obtain a good agreement between the calculated and the experimentally determined values of the magnetic moment [48] and the band gap [92]. The best agreement was obtained for U = 1.5 eV and J = 2 eV.

The obtained U and J values allowed us to carry out the calculations of the vibrational properties of KFeS₂ within the DFT + U scheme. Figure 11 shows the frequency dependence of the total and the element-specific PDOS for KFeS₂ calculated by this advanced way. It reveals a complex structure familiar to us, covering two distinct frequency ranges. The potassium phonon modes dominate at relatively low frequencies of 1–5 THz, while the iron atoms possess the highest vibrational frequencies and constitute a substantial part of the high-frequency range of 9–12 THz. The sulfur atoms show vibrational modes in the whole frequency range under consideration.

Figure 11.

Frequency dependence of the PDOS in KFeS₂ calculated within DFT + U approach: element-specific (K, Fe, and S atoms from the top toward the bottom) and in total (bottom). (adapted from Ref. [53]).

The comparison of the DFT + U calculated PDOS results with the experimental data of nuclear inelastic scattering measurements and the IR-absorption spectrum shows a good agreement between them. Figure 12 proves satisfactory quantitative accordance between the iron partial PDOS derived from our NIS data and that obtained by DFT + U calculation within the entire frequency range of lattice vibrations in KFeS₂ after the frequency axis has been rescaled by a correction factor of 0.92.

Figure 12.

Frequency dependence of the iron partial PDOS of KFeS₂, as derived from the NIS spectrum (blue circles), and from ab initio calculations (magenta line), as described in the text (adapted from Ref. [53]).

As seen in Figure 10, the experimental infrared-absorption spectrum shows two broad lines. Again, like in case of the NIS data, the best agreement for the resonant frequencies of the two IR absorption peaks with the calculated ones was obtained by rescaling the frequency axis by a correction factor of 0.92.

4.2.2 Temperature dependence of KFeS₂ specific heat and its analysis

We calculated the lattice specific heat of KFeS₂ utilizing the calculation procedure we used above for RbFeSe₂. Thereby, DFT + U calculations were applied to estimate the bulk modulus (B = 10.1 GPa) of KFeS₂ from the second derivative of the total energy as a function of the unit-cell volume. The total energy per unit-cell volume itself was determined from a polynomial fit, truncated after the fourth-order term. The thermal expansion coefficient α was derived from crystal-structure data of two states –at room temperature and at 14 K (Ref. [48]), hence the relation of C_P ≈ C_V + 0.003[J/mol K²]·T was applied to calculate the specific heat depicted in Figure 13 (see the legend). Below 25 K, the calculated lattice contribution reveals the cubic temperature law characteristic of the Debye model for the specific heat in solids. The magnetic specific heat of KFeS₂, obtained after subtraction of the calculated lattice contribution from the experimental data, is shown in Figure 14. The corresponding magnetic entropy change was calculated from the experimental magnetic specific heat by integrating from zero temperature up to 300 K. We got ΔS_M = 2.98 Jmol⁻¹ K⁻¹, which is an estimation of the lower boundary, because of the temperature limit (T ≤ 300 K) of our PPMS-9 setup, compared to the maximum temperature of thermal stability of the crystal of around 900–1000 K [50]. As already mentioned above, for quasi-one-dimensional systems, magnetic fluctuations are expected to persist, even significantly above the magnetic ordering temperature. Therefore, the experimental entropy change, which is still below the value Rln2 = 5.76 Jmol⁻¹ K⁻¹ expected for an iron low-spin S = 1/2 state, suggests S = 1/2 as the most probable spin state, but does not exclude S = 3/2.

Figure 13.

Temperature dependence of the specific heat C_P of KFeS₂ (black circles). The red line indicates the calculated lattice contribution to the specific heat at constant pressure C_P = C_V + 0.003[J/mol K²]T (adapted from Ref. [53]).

Figure 14.

Temperature dependence of the magnetic heat capacity (black empty circles), obtained as the difference between the experimental heat capacity and the calculated lattice contribution, and corresponding change of the magnetic entropy (red empty squares) of KFeS₂ (adapted from Ref. [53]).

Nevertheless, such verification of the accurateness of the model calculations on one hand, and obtaining information on the iron-spin state, yielding the observed magnetic moment, on the other hand, motivates further development of advanced theoretical spin-chain models beyond the well-known Bonner-Fischer model of S = 1/2 linear magnetic chains with anisotropic coupling [61].

As another optional important side-conclusion, we need to pay attention that our approach consistently describes the physical properties of quasi-one-dimensional KFeS₂ within DFT + U approximation. Indeed, it is well known that GGA+U functionals may underestimate band gaps, therefore, hybrid functionals [93, 94, 95, 96] or self-interaction-corrected DFT [97, 98] techniques have been used to obtain more accurate results. We used a strategy of finding values of U and J (U = 1.5 eV and J = 2 eV [99]) parameters within the DFT + U route, which provides the total agreement with the experimentally measured magnetic moment per iron ion (1.9 μ_B) and band gap (0.35 eV). Then, the phonon frequencies were calculated from the eigenvalues of the dynamic matrix. The latter involves the evaluation of interatomic forces in all directions acting on atoms of 18 supercells for the particular KFeS₂ lattice. Thus, for each supercell containing 96 atoms, the dynamic matrix accumulates 96 × 18 × 3 = 5184 force constants. The resulting iron partial PDOS is compared to the experimental PDOS obtained from nuclear inelastic scattering. The quite good agreement opens prospects for the analysis of thermal properties. The overall consistency of the calculated magnetic moment, electronic band gap, and iron PDOS with the experimental results in KFeS₂ makes the DFT + U scheme promising for studying the physical properties of related chain compounds by the low-cost DFT calculation methods. This is of special interest in the light of recent attempts to search for high-temperature superconductivity in such kind of compounds [98].

Now, we should clarify that our studies in the field of quasi-one-dimensional iron chalcogenides presented in this chapter might be applied to other compounds of that family and some results, including complex DFT + U calculations, IR-absorption spectroscopy measurements, Mössbauer spectroscopy, etc., have already been achieved. Nevertheless, in continuation of such studies on further compounds of this family, we inevitably shall have to face with the whole-aspect interpretation challenges, as we have already encountered with two considered compounds. So, despite our best efforts and the application of a wide arsenal of investigation methods and the strategy of synergetic combination of their results, it has become apparent that we are at a crossroads with no proper roadmap to choose the exact way for correct theoretical description of the behavior of iron-spin chains in the systems under consideration. Given the complexity and unique nature of the compounds we had focused on, it was imperative to use new avenues for solving the problems.

So, these innovative, modern, and unconventional approaches we used provided everyone with a wide and strong experimental facts-based fundament for further building a strong and full theoretical description of such compounds. We believe that the fresh breakthrough insights and the application of approaches based on complex theoretical physics methods will be able to overcome the obstacles regarding the full description. The potential benefit of leveraging theoretical physics approaches to devise novel solutions is obvious. We have reached a few milestones during solving the problem of a full-sided description of the family of quasi-one-dimensional chain magnets, but it needs to take a break due to the fact that the absence of key theoretical models became like a barrier to solving the problem. However, the complexity of these compounds allows considering it like model-systems for the verification of the new theoretical models which might be further generalized for the class of such complex magnets.

5. Conclusions

The lattice vibrations of ternary quasi-one-dimensional iron chalcogenides AFeX₂ (A = K, Rb; X = S, Se) were modeled using a unified approach based on density functional theory, taking into account intra-atomic and exchange Coulomb correlations in the d-electron system (DFT + U approach).

The phonon frequencies were calculated from the eigenvalues of the lattice dynamic matrix. The equilibrium crystal structure optimized with respect to atomic positions, cell shape, and cell volume was used for calculations of the phonon dispersion in harmonic approximation. In the antiferromagnetic phase of the AFeX₂ compounds, the calculations account for the spin polarization in the magnetically ordered state. It is important to note that the MedeA Phonon commercial software was allowed to perform element-specific calculations, thus providing the partial PDOS for every kind of ion in a particular AFeX₂ compound. It appeared that the PDOS, both partial and total, exhibit strongly non-Debye shape: vibration frequencies of cross-linking alkaline ions lie in the low-frequency domain of 2–4 THz and can be more or less satisfactorily approximated by a PDOS of the Debye type. However, the partial PDOS of iron and chalcogen vibrations represent broad bands of low-frequency contributions of strongly non-Debye shape supplemented by high-frequency, double- or triple-peak Einstein-like vibration bands in the frequency range of 6.7–12 THz. The strongly non-Debye shape of the total PDOS is a consequence of the extremely anisotropic lattice of the AFeX₂ family compounds and the complex structure of the elementary cell comprising rigid and soft bonds.

It is not surprising that attempts to build a phenomenological description of the temperature behavior of the heat capacity based on the Debye model led to unsatisfactory results. The addition of Einstein modes, especially their number, for better agreement with experiment was impossible to justify phenomenologically. Good agreement with experiment was found for the temperature behavior of the heat capacity in a model combining one Debye contribution and two Einstein ones, which, however, ran into a contradiction when comparing the change in magnetic entropy during the magnetic disorder/order transition. Estimates from the heat capacity measurements have given approximately 10% of Rln2 = 5.76 JK⁻¹ mol⁻¹—the minimum possible change for spin S = 1/2 at an iron ion. Quite a ridiculous conclusion! Moreover, estimates of the Debye temperature obtained from measurements of the heat capacity and temperature dependence of the Lamb-Mössbauer factor differ by a factor of 2, which is also annoying to the eye. Thus, the application of the results for DFT + U phonon calculations to experiments on heat capacity could resolve the obvious contradictions that the application of phenomenological models leads to.

Before turning to the experiment and C(T) analysis, we noticed that DFT approach usually underestimates frequencies of the lattice vibrations. Fortunately, nuclear inelastic scattering (NIS) technique is a direct tool to obtain the PDOS from measurements, which the calculated PDOS can be compared to, so that, the frequency correction factor could be deduced. It is an advantage of our calculations that they give element-specific PDOS, since NIS is determined solely by iron ion vibrations.

High-quality single-crystalline AFeX₂ samples were grown by means of the Bridgman method. The temperature dependence of the specific heat was measured. Then, to establish the frequency correction factor, nuclear inelastic scattering measurements were carried out and the partial Fe-ion PDOS was derived from the NIS spectrum. To calibrate the frequency scale of the DFT + U calculations, the calculated iron PDOS was directly adjusted to the experimental one. Regarding, for example, RbFeSe₂, a correction factor f = 1.08 was obtained from the best fit for RbFeSe₂ compound. The comparison of the experimental temperature dependence of the Lamb-Mössbauer factor with the calculated one using direct summation over DFT + U phonon modes for iron also gave the correction factor f ≃ 1.08. Moreover, the IR absorption measurements showed the presence of two high-frequency peaks semiquantitatively matching the calculated PDOS after correction of frequency scale by a very similar factor f = 1.06. For DFT + U partial PDOS calculations, the correction of the frequency scale has indisputable implications on rectifying the magnetic contribution to the specific heat and the change of magnetic entropy upon the disorder-order magnetic transition. Consistency of the three experimental techniques in determining the frequency correction factor formed the basis of considering the temperature dependence of the specific heat.

After correction of the frequency scale, the ab initio results for the total PDOS were used to calculate the temperature dependence of the lattice contribution to the specific heat of AFeX₂ without any further model approximations (neither Debye nor Einstein). Note here that also no other fitting parameters were used in the calculation. The magnetic specific heat was obtained by subtracting the calculated lattice contribution from the experimental specific heat in the temperature range of 2–300 K. The result for C_V(T) was converted to C_P(T) with the relationship given in Eq. (5). The resulting magnetic entropy change ΔS_M(RbFeSe₂) ≃ 6.03 J mol⁻¹ K⁻¹ and ΔS_M(KFeS₂) ≃ 2.98 J mol⁻¹ K⁻¹, obtained by integrating the magnetic specific heat divided by temperature C_M/T, is close to 5.76 J mol⁻¹ K⁻¹, expected for the pure low-spin (S = 1/2) state of the iron ion at T_N. Certainly, there is enough space to collect further entropy change up to the crystal stability temperature above 900 К. This suggests that the intermediate S = 3/2 spin state cannot be excluded to be realized for the iron ions in the chains. The measured magnetic moment of 2.66 μ_B and the calculated moment of ∼2.80 μ_B pег iron ion in KFeS₂ are in agreement with the intermediate spin state.

The estimation of the magnetic entropy change provides, therefore, a lower bound for the reduced value of the iron-spin state. This could be important for advancing the theoretical models of spin chains, generalizing the S = 1/2 linear magnetic chains model by Bonner-Fischer.

The approach we proposed and developed dramatically improves the consistency between the magnetic entropy change upon the disorder-order transition in quasi-one-dimensional AFeX₂ compounds, deduced from the specific heat measurements, with the expected entropy change in the spin-chain models. However, it is a bit too early to herald that the accurate quantitative conclusion about the iron ion spin state can be made from the thermal measurements. Indeed, because the magnetic contribution to the entropy change is an integral characteristic over a wide temperature range from absolute zero up to ∼900 K, a limited temperature range of validity of the accepted approximations can introduce sizeable deviations from the true result. Hence, the dependence of the approximations and basic parameters on temperature must be considered, namely: harmonic vibration approximation breaks down; phonon eigenfrequencies depend on temperature, the dependence is expected to be anisotropic; the frequency correction factor can be different for particular constituent ions of the lattice; it can be inhomogeneous in frequency, etc. This tells us that there is still a lot of work to be done further to make the theory quantitative at a convincing level.

Acknowledgments

A.K., M.K., and D.A.T. acknowledge the Kazan Federal University Strategic Academic Leadership Program (PRIORITY-2030). This work was funded as joint project by the Russian Foundation for Basic Research (RFBR) and the German Research Foundation (DFG) under the contract numbers №19-52-12068 and KR2254/3-1, respectively. We are especially grateful to Hans-Christian Wille, Ilya Sergueev, Olaf Leupold, Farit Vagizov, Yury Lysogorskiy, Irina Filippova, and Tatyana Gavrilova for their invaluable contributions to original works.

Conflict of interest

The authors declare no conflict of interest.