Above-room-temperature strain-stable half-metallic phase owing to a high thermoelectric response in CaCu₃Cr₂Re₂O₁₂

Abstract

We present the electronic, magnetic, and thermoelectric properties of the CaCu₃Cr₂Re₂O₁₂ quadruple double perovskite oxide (QDPO) using ab initio calculations. Strong antiferromagnetic interactions between Cu²⁺↑Cr³⁺↑Re⁵⁺↓ result in a ferrimagnetic (FiM) ground state. A half-metallic (HM) state is evident with a finite energy gap (E_g) of 1.62 eV in the spin-majority channel (N^↑), which is large enough to ensure a long mean free path for spin and prevent spin-flipping with a colossal spin-filtering ability. Due to a significant E_g in N^↑, the system displays a high figure of merit (0.76) at 300 K. The calculated spin moments of +0.48/+2.55/−0.92 μ_B and spin-magnetization densities on the Cu/Cr/Re further verify the FiM state. The magnetic phase transition yields a Curie temperature (T_C) of 340 K, which is below the experimental value (360 K). Because thermal energy near T_C disrupts the magnetic ordering, magnetization is consequently reduced, which is also reflected in the susceptibility curve. Additionally, the FiM state of the structure was confirmed under modulated applied magnetic fields. Finally, the HM FiM state is also verified under a moderate biaxial ([110]) strain of ±5%. Hence, this work provides deep insights into this QDPO, highlighting its stable HM FiM behavior, above-room-temperature T_C and high thermoelectric response, which promises potential applications in spintronics.

---

1 Introduction

In the realm of magnetic materials fabricated at the ultra-small scale, special focus is placed on those composed of transition metals (TMs) such as Fe, Co, and Ni, along with their alloys.¹ One of the most promising types of magnetic materials is half-metallic (HM) compounds, such as ferromagnetic (FM) and ferrimagnetic (FiM) materials. These compounds exhibit metallic conductivity in one spin channel while behaving as semiconductors or insulators in the opposite spin channel.² This means that the electron carriers are entirely (100%) spin-polarized (SP) at the Fermi level (E_F), producing HM magnetic motifs that are highly essential for spintronics applications, essentially enabling the manipulation and storage of information based on electron spins, such as in tunneling devices, spin transistors, and magnetic logic gates.^3–9 It has been observed that tunneling conduction in a spin-based trilayer essentially depends on the SP of the electrodes,^10,11 consequently resulting in a giant tunneling magnetoresistance (TMR), which can be utilized in high-density data storage and extremely efficient magnetic sensors.¹² However, an enormous TMR ratio in trilayer devices is susceptible to temperature (temp.), and it sharply decreases with rising temp. due to a substantial reduction in SP magnitude near the magnetic transition temp. (MTT). Hence, materials that maintain maximum (100%) SP charge carriers at room temp. with an MTT > 300 K are highly beneficial for developing future high-temp. spintronic devices. The first predicted HM FM material was the NiMnSb half-Heusler alloy¹³ followed by CrO₂ ^14,15 and Fe₃O₄ ^16,17, which were also proposed as HM and investigated comprehensively. In this regard, perovskite oxides (POs) have become very important candidates, especially after the discovery of sizable magnetoresistance in La_1−xSr_xMnO₃ ¹⁸ and Tl₂Mn₂O₇ ¹⁹, further emphasizing the importance of HM structures. Similarly, double perovskite oxides (DPOs) with an ordered A₂BB′O₆ form, where the B and B′ ions are usually 3d and 4/5d TMs,^20,21 are highly beneficial for spintronic applications.²² Excitingly, an HM FiM state is anticipated in the Sr₂FeMoO₆ DPO, which has a metal-to-insulator transition (MIT) at 410 K and also demonstrates a large TMR ratio at room temperature. In this context, the Fe³⁺ 3d⁵ and Mo⁵⁺ 5d¹ ions with are coupled antiferromagnetically, resulting in a FiM spin-ordering (SO).³ Similarly, an HM FiM behavior is observed in Sr₂CrReO₆ owing to it having the highest Curie temp. (T_C) of 635 K among the DPO family.^23,24 However, it is theoretically established that SP in HM DPOs such as Sr₂FeMoO₆ substantially depends on temp. and cation disorder, as well as being proportional to the spin-magnetization of the cores.²⁵ In the past decade, numerous efforts have been dedicated to develop robust 2D FM/FiM materials with high T_C. One effective approach involves modulating the magnetism in existing 2D layered materials through post-treatment techniques, such as doping, ion or charge implantation,^26,27 strain engineering,^28,29 and defect engineering.^30,31 The research discovered that the stress-induced transformation of two-dimensional materials from an indirect to a direct energy gap (E_g) significantly enhances photoluminous (PL) efficiency, making them more suitable for optoelectronic devices.^32–41 Thus, novel nano-electronic and optoelectronic devices with stress-adjustable properties hold significant potential for spintronic applications. Moreover, to optimize material properties, the doping/ion implantation strategy has been explored, in which additional magnetic cations are introduced at the A-site in the DPO with a ratio of 1 to 3. This results in a new A- and B-site ordered quadruple DPO (QDPO) keeping an form, which possesses a very high MTT and magnetization (M). For instance, the antiferromagnetic (AFM) interaction between the magnetic cations induces a FiM state in the CaCu₃Fe₂Re₂O₁₂ QDPO, which features a very high MIT of 560 K and 100% SP. In this state, the spin-majority and spin-minority channels (N^↑ and N^↓) exhibit nonconductive and conductive traits, respectively, along with a large total magnetic moment (m_t) of 8.7 μ_B.⁴² In this continuum, a newly ordered NaCu₃Fe₂Os₂O₁₂ QDPO was synthesized, which exhibits an HM FiM nature with a wide energy gap (E_g) in N^↑ and spin ordering of Cu²⁺↑Fe³⁺↑Os^5.5+↓, showing a T_C above 380 K.⁴³ Combined experimental and theoretical studies demonstrate that the FiM LaCu₃²⁺Co₂²⁺Re₂^5.5+O₁₂ QDPO is nearly 100% ordered with Cu²⁺↑Fe³⁺↑Os^5.5+↓ SO and a T_C of 150 K.⁴⁴ Wang et al. theoretically predicted a series of HM FiM ACu₃Fe₂Re₂O₁₂ (A = Ca, Sr, Ba, Pb, Sc, Y, La) QDPOs, where the AFM coupling between CuFe and Re is dominant.⁴⁵ It was established that compounds with A²⁺ = Ca/Sr/Ba/Pb and A³⁺ = Sc/Y/La ions exhibit T_C values greater than 405 K and 502 K, respectively. Similarly, LaCu₃Fe₂Re₂O₁₂ is a HM material with a very large T_C of 710 K, where a strong AFM coupling of Cu²⁺↑Fe³⁺↑Os^4.5+↓ results in a FiM lowest energy state.⁴⁶ It was also observed that the saturated total m_t of 7.0 μ_B at 710 K increases to 8.0 μ_B at 2 K. Very recently, Zhang et al. experimentally verified the HM FiM nature in the CaCu₃Cr₂Re₂O₁₂ (CCCRO) system, maintaining a spin ordering of Cu²⁺↑ Cr³⁺↑ Re⁵⁺↓ above room temp. (360 K).⁴⁷ However, in this stable charge configuration, the calculated m_t of 5.0 μ_B is significantly greater than the experimentally measured value of 1.6 μ_B per f.u. Therefore, the intriguing features of this new class of materials, as outlined above, motivate us to explore the different characteristics of the recently synthesized CCCRO QDPO through density functional theory (DFT) calculations. These calculations demonstrate a FiM spin ordering, accompanied by a m_t of 5.0 μ_B per formula unit (f.u.), and a half-metallic nature, featuring a 1.62 eV gap in the N^↑ channel. This results in a high T_C of 340 K. The material shows metallic behavior in the N^↓ channel, with a p-type thermoelectric (TE) response. Moreover, the HM phase of the system remains stable under biaxial strain, and the material exhibits promising TE and catalytic properties, making it a strong candidate for future device applications. Finally, this paper aims to develop a theory of the magnetization process for CCCRO QDPO-based ultra-thin spin valves, which have proven to be effective as high-temperature storage devices.

2 Computational details

To study the physical aspects of the considered material, DFT calculations were performed within the framework of the full-potential linearized augmented plane-wave (FP-LAPW) method as implemented in the WIEN2k code.⁴⁸ The generalized gradient approximation (GGA)⁴⁹ was used for the exchange–correlation function along with the Hubbard parameter (U) of 4.0, 3.5, and 2.4 eV for the Cu 3d, Cr 3d, and Re 5d orbitals, respectively, to account for correlation effects.⁵⁰ Additionally, spin–orbit coupling (SOC) effects were included in the scalar relativistic form due to the presence of a heavy Re element. Inside the atomic sphere, l_max = 12 was set as an upper limit, and the plane-wave cutoff was chosen as R_mt × K_max = 7 and G_max = 24. The K_max was chosen to ensure no charge leakage occurred and the total energy converged. A 10 × 10 × 10 k-mesh was used, generating 76 points within the full Brillouin zone, sufficient for energy convergence. Moreover, the energy convergence as a function of the k-grid has been checked carefully, as shown in Fig. 1S of the ESI.† This clearly displays that the chosen k-point grid is large enough for the energy convergence. The atomic positions were fully relaxed until the atomic forces reached ∼5 mRy per a.u. Self-consistency was ensured by total energy (E_t) and charge convergence criteria of 10⁻⁵ Ry and 10⁻⁴ C, respectively. Furthermore, T_C calculations were carried out using the Vampire code.⁵¹ To perform the calculations, first, the experimental atomic sites are fully relaxed until the charge and forces converged to a minimum criterion mentioned above. After that, the Self-Consistent Field method is used to converge the E_t up to five decimal places, which is a very high metric to accurately described the electronic structure of the compound.

3 Results and discussion

The crystal diagrams of the A–B site ordered CCCRO QDPO are shown in Fig. 1, with different magnetic orderings (unconstrained system), exhibiting the space group Pn3̄ (No. 201).⁴⁷ The experimentally observed lattice constants of a = b = c = 7.4152 Å ⁴⁷ were used for the present study and our DFT calculated value of 7.4097 Å is in excellent agreement with the experiment. The Wyckoff positions for Ca, Cu, Cr, Re, and O ions are 2a(0.25, 0.25, 0.25), 6d(0.25, 0.75, 0.75), 4b(0, 0, 0), 4c(0.5, 0.5, 0.5), and 24h(0.5583, 0.7559, 0.0662), respectively. First, we investigated the magnetically stable state of the system by considering four possible spin orderings (SOs): FM (Cu²⁺↑, Cr³⁺↑, Re⁵⁺↑), FiM-1 (Cu²⁺↑, Cr³⁺↑, Re⁵⁺↓), FiM-II (Cu²⁺↓, Cr³⁺↑, Re⁵⁺↓), and FiM-III (Cu²⁺↓, Cr³⁺↑, Re⁵⁺↑), as displayed in Fig. 1(a–d), respectively. In the case of the FM SO, the spins of all magnetic ions are aligned in the same direction (see Fig. 1(a)). However, in FiM-1 SO, the Cu and Cr spins are parallel to each other, but anti-parallel to that of Re (see Fig. 1(b)). For FiM-II SO, the Cu and Re spins are aligned with each other but anti-aligned to that of Cr (see Fig. 1(c)). Finally, in FiM-III SO, the spins of Cr and Re are aligned in the same direction as each other, whereas they are anti-aligned with the Cu spin (see Fig. 1(d)). Henceforth, the Cu²⁺↑Cr³⁺↑Re⁵⁺↓ (FiM-1) configuration is the lowest energy state with energy differences of −105/−78/−55 meV to the FM/FiM-II/FiM-III SO. This state has a m_t of 5.0 μ_B per f.u., which is in excellent agreement with recent experimental observations.⁴⁷ Notably, the lengths of the arrows in the figures represent the magnitude of each atom’s spin magnetic moment (m_s). Next, we studied the structural stability of the CCCRO material in its stable FiM-I phase by calculating the formation and cohesive energies (E_for. and E_coh.) as:where is the DFT-computed E_t of the CCCRO structure, while E^bulk_Ca, E^bulk_Cu, E^bulk_Cr, E^bulk_Re, and are the DFT-determined E_t for the Ca, Cu, Cr and Re ions and O molecules in their bulk stable phases, respectively. Similarly, E^iso._Ca, E^iso._Cu, E^iso._Cr, E^iso._Re, and are the DFT-measured E_t for the isolated Ca, Cu, Cr and Re ions, and oxygen molecule, accordingly. The calculated E_for./E_coh. per atom is −2.31 eV/−4.89 eV. The negative signs for both energies confirm that the structure is thermodynamically stable and supports an exothermic reaction. Moreover, the mechanical stability of the structure is utilized by calculating the elastic tensors (C_ij). The computed three independent elastic parameters, namely C₁₁, C₁₂, and C₄₄, satisfied the basic demands and Born criteria as (C₁₁ − C₁₂) > 0, (C₁₁ + 2C₁₂) > 0, C₁₁ > 0, C₄₄ > 0, and C₁₂ < B < C₁₁ (ref. 52) for the mechanical stability. Hence, the determined C₁₁, and C₁₂, and C₁₂ are 289, 211, and 140, correspondingly, which satisfies the above-mentioned mechanical stability standard.

Fig. 1 Crystal structures of the CaCu₃Cr₂Re₂O₁₂ quadruple structure in a (a) ferromagnetic (FM), (b) ferrimagnetic (FiM)-I, (c) FiM-II, and (d) FiM-III spin ordering.

To describe the electronic structure of the CCCRO, we calculated and presented the SP total/partial density of states (TDOS/PDOS) within the GGA+U and TDOS for GGA+U+SOC methods in Fig. 2(a, a′ and b), respectively. Fig. 2(a) suggests that the system is HM within the GGA+U method with fully SP conduction electrons, which is consistent with the experimental findings.⁴⁷ A large E_g of 1.62 eV remains in the N^↑ channel, and a few TDOS states appear at the Fermi level (E_F) in N^↓, leading to a metallic behavior of N^↓. Consequently, the entire system illustrates an HM phase. The broad E_g in N^↑ is sufficiently large to prevent spin-flipping, which helps to avoid thermal excitation of the electrons and ensures that the HM nature of the system is stable at or above room temp. to some extent. From Fig. 2(a′), we observe that only Re 5d and O 2p states (where hybridization occurs between Re 5d and O 2p states near the E_F) cross the E_F in N^↓, transitioning from the valence band (VB) to the conduction band (CB), which is responsible for the metallicity. In contrast, no ion states cross or lie at the E_F in the N^↑, confirming the 100% SP at the E_F. Additionally, the GGA+U+SOC calculated TDOS (see Fig. 2(b)) exhibits the metallic behavior of the system. Furthermore, the calculated SP band structures within the GGA+U and GGA+U+SOC methods (shown in Fig. 3(a/a′ and b), respectively), validate the HM and metallic nature of the structure as demonstrated by the TDOS in Fig. 2(a and b). To further confirm the HM state of the structure, we calculated the spin-polarized TDOS of the system using the hybrid functional PBE0, as shown in Fig. 2S of the ESI,† which is implemented in Wien2K and gives almost the same high accuracy in results as found in the case of HSE06, which is implemented in the VASP code. One can see that the calculated TDOS exhibits a HM state where N^↑ contains a large E_g, while N^↓ is showing a metallic behavior within PBE0 (see Fig. 2S of the ESI†), as predicted in the case of the GGA+U method (see Fig. 2(b)).

Fig. 2 Calculated spin-polarized density of states (DOS) in the CaCu₃Cr₂Re₂O₁₂ quadruple structure: (a/a′) total/partial DOS for the GGA+U and (b) total DOS for the GGA+U+SOC method,.

Fig. 3 Calculated spin-polarized band structures in the CaCu₃Cr₂Re₂O₁₂ quadruple structure, for the (a/a′) spin-majority/spin-minority channel within the GGA+U (1st-row) and (b) GGA+U+SOC (2nd-row) methods.

Next, we discuss the system’s magnetism, where the energetically stable FiM-I configuration Cu²⁺↑Cr³⁺↑Re⁵⁺↓ produces an m_t of 5.0 μ_B per f.u., which is higher than that of other spin configurations. This value is reasonable, as also mentioned in ref. 47. However, it remains considerably higher than the experimentally measured value of 1.6 μ_B per f.u.⁴⁷ The calculated m_s for the Cu, Cr, and Re ions within the muffin-tin sphere are 0.48/2.54/−0.92 μ_B using the GGA+U method, which changes slightly to 0.47/2.53/−0.89 μ_B upon including spin–orbit coupling (SOC). The “−” sign in the Re m_s indicates that its spin is oppositely aligned with the Cu and Cr spins, confirming the FiM state. The Cu, Cr, and Re ions are in the +2 (3d⁹), +3 (3d³), and +5 (5d²) states, respectively, holding electronic configurations of t³_2g↑t³_2g↓e²_g↑e¹_g↓, t³_2g↑t⁰_2g↓e⁰_g↑e⁰_g↓, and t²_2g↑t⁰_2g↓e⁰_g↑e⁰_g↓ with spin quantum numbers , respectively. Hence, the computed m_s on each ion also verifies their respective valence states. The induced orbital magnetic moment (m_orb.) on the Cu/Cr/Re ion is 0.06/−0.05/0.23 μ_B. Compared to that on the Cu and Cr, the sizable magnitude of m_orb. on the Re ion highlights the strong SOC influence. Furthermore, as Cu lies in a more than half-occupied d state, its m_s and m_orb. are aligned with each other (i.e., both have a “+” sign), as discussed in ref. 53. Conversely, Cr and Re ions lie in less than half-filled d states, so their m_s and m_orb. are anti-aligned. For Cr, the signs are “+/−” for m_s/m_orb., while for Re, they are “−/+”, as described in ref. 54. Thus, the above discussions support the AFM interactions between the Cu/Cr and Re ions.

To further directly confirm the m_s magnitude on the magnetic ions and the FiM SO in the system, we visualized the three-dimensional spin-magnetization density iso-surfaces with an iso-value of ±0.005 e Å⁻³ in Fig. 4. Firstly, the high densities around the Cr ion, compared to those around the Cu and Re ones, establish that Cr is the primary contributor to the system magnetism and confirm the calculated m_s values on each ion. Secondly, the Cu and Cr density colors are the same (red), in contrast to that of Re (light green). This indicates that the Cu and Cr spins are aligned in the same direction (parallel to each other), but they are anti-aligned to the Re-spin, which further supports the FiM-I (Cu²⁺↑Cr³⁺↑Re⁵⁺↓) SO. Additionally, the combined e_g and t_2g orbital features are evident in the Cu and Cr/Re ion densities as they occupy +2 (e²_g↑e¹_g↓) and +3 (t³_2g↑)/+5 (t²_2g↑) states, respectively. Therefore, the AFM coupling between Cu/Cr and Re ions can be understood using the super-exchange interaction concept, where the possible electron transfer between them is illustrated in Fig. 5. Since Cr and Re maintain +3 and +5 valence states, as t_2g states are less than half-filled in both cases (see Fig. 5), virtual hopping between these t_2g orbitals leads to AFM coupling, resulting in a FiM state.

Fig. 4 Calculated three-dimensional spin-magnetization density iso-surfaces with an iso-value of ±0.005 e Å⁻³ in the CaCu₃Cr₂Re₂O₁₂ quadruple structure. For simplicity, Ca and oxygen ions are omitted.

Fig. 5 Schematic representation of the antiferromagnetic superexchange coupling between Cu²⁺ (3d⁹)/Cr³⁺ (3d³) and Re²⁺ (5d²) orbitals via O⁻² (2p²) in the CaCu₃Cr₂Re₂O₁₂ quadruple structure.

The total m_s of an atom is the result of the combined contributions from the m_s values of its electrons, including the m_orb. component. In FiM materials, the magnetism primarily originates from the m_s of the unpaired metal ions within the crystalline lattice. Consequently, these materials display magnetic moments with spontaneous long-range ordering, driven by a quantum mechanical phenomenon called the exchange or interchange interaction. The exchange interaction, which governs SO, can be categorized into short-range (direct), long-range, or indirect exchange interactions. This interaction is characterized by the presence of two distinct sublattices, A and B, and at least three different exchange interactions: J_AA, J_AB, and J_BB.⁵⁵ The ions in sublattice A are oriented in one direction, whereas those in sublattice B are oriented in the opposite direction. For example, the exchange interaction results from the parallel alignment of atomic m_s in FiM materials. Hence, to figure out the origin of the high T_C in the CCCRO system, the magnetic exchange couplings were computed using the Heisenberg spin method:^56,57

Here, J_ij represents the interaction between the magnetic ions i and j due to exchange forces. If J_ij is positive (J_ij > 0), an FM interaction occurs, allowing the spins to align in a parallel configuration (↑↑). In contrast, if J_ij is negative (J_ij < 0), an AFM interaction arises, resulting in an antiparallel spin alignment (↑↓). FM materials exhibit the former behavior, while AFM and FiM are associated with the latter. S_i/S_j is the spin vector of the Cu/Cr/Re ion at the i/j site and represents the spin variables. For simplicity, only the three nearest neighbors J_ij, namely J_Cu–Cr, J_Cu–Re, and J_Cr–Re, are taken into account. To do this, J_ij in four varying SOs, namely FM (Cu²⁺↑Cr³⁺↑Re⁵⁺↑), FiM-I (Cu²⁺↑Cr³⁺↑Re⁵⁺↓), FiM-II (Cu²⁺↓Cr³⁺↑Re⁵⁺↓), and FiM-III (Cu²⁺↓Cr³⁺↑Re⁵⁺↑), as shown in Fig. 1, is calculated as:⁴⁶

E_t(FM) = E₀ + 12J_Cu–CrS_CuS_Cr + 12J_Cu–ReS_CuS_Re + 12J_Cr–ReS_CrS_Re (4)

E_t(FiM-I) = E₀ + 12J_Cu–CrS_CuS_Cr − 12J_Cu–ReS_CuS_Re − 12J_Cr–ReS_CrS_Re (5)

E_t(FiM-II) = E₀ − 12J_Cu–CrS_CuS_Cr + 12J_Cu–ReS_CuS_Re − 12J_Cr–ReS_CrS_Re (6)

E_t(FiM-III) = E₀ − 12J_Cu–CrS_CuS_Cr − 12J_Cu–ReS_CuS_Re+ 12J_Cr–ReS_CrS_Re (7)

Our results show that the three considered couplings are AFM. In the case of Cu²⁺↑–Cr³⁺↑, J_ij = 0.957 meV for the atom interactions closest to Cu–Cr and becomes −0.374 meV for the second-to-nearest ones, hinting at the AFM coupling between Cu²⁺↑ and Cr³⁺↑. However, the AFM Cu²⁺↑–Re⁵⁺↓ and Cr³⁺↑–Re⁵⁺↓ interactions are much stronger than those of Cu²⁺↑–Cr³⁺↓, resulting in a long-ranged Cu²⁺↑–Cr³⁺↑–Re⁵⁺↓ SO (FiM-I) as the lowest energy state of CCCRO. The calculated values of the AFM exchange coupling constants J_Cr–Re = −30.21 meV and J_Cu–Re = −28.57 meV for the interactions S_Cr↑–S_Re↓ and S_Cu↑–S_Re↓ are substantially higher than the 5.85 meV for J_Cu–CrS_Cu↑S_Cr↑. Therefore, the much stronger AFM interactions between the Cu/Cr and Re ions result in a high T_C = 340 K for the system.

Magnetic materials typically form magnetic domains, regions where the magnetic moments are highly aligned, and the boundaries between the magnetic domains (domain walls) shift in the presence of an external magnetic field (H), leading to changes in the material’s M. The domain structure organizes itself to minimize the overall magnetostatic energy of the material. The domains’ size, shape, and orientation are determined by the interactions between the exchange, magnetostatic, and anisotropic energies of the system. Consequently, exposure of a material to an external H induces an M within the material. The magnitude of the M is directly proportional to the applied H, as given by the relation:

M = χH

In this equation, χ represents the material’s magnetic susceptibility, a dimensionless quantity in SI units. The behavior of solid materials under an external H is affected by factors including their atomic structure, electromagnetic excitation, pressure, and temperature. Hence, to visualize the magnetic phase transition at T_C in the CCCRO QDPO, M versus T_C calculations have been performed (presented in Fig. 6) using the Vampire code.^51,58 The computed T_C is 340 K, ∼7% lower than the experimentally observed value of 360 K.⁴⁷ This discrepancy can primarily be attributed to inherent differences in computational modeling and experimental conditions. As the temp. approaches T_C, thermal energy disrupts the ordered alignment of m_s, corresponding to a sharp decrease in M. This results in a corresponding decrease in the m_s of the constituent ions due to the reduced number of unpaired d-electrons. This phase transition is critical as it signifies the loss of spontaneous M and the onset of enhanced thermal fluctuations, which can be seen in the χ curve shape as depicted in Fig. 6. This implies some level of anisotropy within the material and the influence of opposing magnetic interactions.⁵⁹ Importantly, the sharp peak in χ observed at T_C indicates the material’s heightened sensitivity to H due to increased SO.

Fig. 6 Calculated susceptibility and magnetization curves versus Curie temperature of the CaCu₃Cr₂Re₂O₁₂ quadruple structure.

Recently, a material which is applicable for spintronics application, contain the M process of the spin valve of HM ferrites, has been attracting deep attention. M as a function of the applied H is represented by the M–H curve, which reveals the saturation M that is reached beyond a certain threshold of H. As H decreases to zero after saturation, M does not revert to zero, a behavior referred to as hysteresis. This property is critical for various applications, such as magnetic sensing devices. Therefore, the magnetic properties of CCCRO are evaluated by measuring M in terms of the applied H and temp. of 2 T/5 T and 2/100/200/300 K, respectively, as shown in Fig. 7(a) and (b). At low temp. (2 K), as shown in Fig. 7(a), the material exhibits strong SO and quickly drives the system to saturation, displaying sharp transitions in the hysteresis loop. As the temp. increases to 100 K, thermal fluctuations begin to reduce the M, however, the material still reaches a saturation, as found under 2 T (see Fig. 7(a)), though more gradually. At 5 T (see Fig. 7(b)), the material achieves saturation more effectively, but the saturated M is lower at 2 K. At higher temp. (200 K and 300 K), the material’s M continues to decrease due to an increase in thermal agitation, and the transition to saturation becomes more gradual. However, the high field (5 T) still drives the material towards saturation, but the overall M is lower.

Fig. 7 Variation of the saturated magnetization (M) with applied magnetic field (H) for (a) 2 T and (b) 5 T at a fixed temperature of 2/100/200/300 K in the CaCu₃Cr₂Re₂O₁₂ quadruple structure.

By analyzing the band structure, the material’s metallic, semiconducting, or insulating nature can be inferred from its transport properties. A decrease in the Hall coefficient (R_H) with increasing temp. suggests metallic behavior, while a Seebeck coefficient (S) greater than zero signifies p-type semiconducting characteristics. Thus, the transport properties of the CCCRO QDPO were computed using the BoltzTraP2 software⁶⁰ under the constant relaxation time approximation (τ = 1 × 10⁻¹⁹ s). The electrical conductivity per unit relaxation time exhibits a distinct temp. dependence due to the SP electronic structure. There is a significant increase in σ for non-metallic N^↑, from 3.18 × 10¹⁹ S m⁻¹ at 300 K to 4.0 × 10¹⁹ S m⁻¹ at 800 K, as thermal excitation promotes electrons across the gap. In contrast, N^↓ shows a weaker temp. dependence, from 2.58 × 10¹⁹ S m⁻¹ at 300 K to 2.80 × 10¹⁹ S m⁻¹ at 800 K, as conductivity is primarily limited by electron-phonon scattering (see Fig. 8(a)). Initially, N^↑ has a higher S (1.92 × 10⁻⁵ V K⁻¹) at 300 K due to limited thermally excited carriers, while N^↓ exhibits a lower S (−2.94 × 10⁻⁶ V K⁻¹) at 300 K. As temp. increases (up to 800 K), the S for N^↑ decreases (−3.73 × 10⁻⁶ V K⁻¹) as thermal excitation saturates available states, while the S for N^↓ increases (8.66 × 10⁻⁶ V K⁻¹) due to enhanced carrier diffusion (Fig. 8(b)). After a certain temp., S stabilizes for both channels, indicating saturation in the TE potential.

Fig. 8 Calculated (a) electrical conductivity per unit relaxation time , (b) Seebeck coefficient (S), (c) thermal conductivity per unit relaxation time , (d) molar magnetic susceptibility (χ), (e) Hall-coefficient (R_H), and (f) figure of merit (ZT) of the CaCu₃Cr₂Re₂O₁₂ quadruple structure.

The thermal conductivity per unit relaxation time increases linearly with temp. primarily due to phonon–electron scattering, peaking around 800 K for both channels. For N^↑, the conductivity is slightly higher, with the value increasing from 2.88 × 10¹⁴ W m⁻¹ K at 300 K to 1.02 × 10¹⁵ W m⁻¹ K at 800 K (Fig. 8(c)). This difference arises from the increased availability of high-energy states in N^↑ due to thermal excitation across the E_g. N^↓, being metallic, shows a less pronounced increase, from 2.01 × 10¹⁴ W m⁻¹ K at 300 K to 6.27 × 10¹⁴ W m⁻¹ K at 800 K, due to the relatively constant DOS at the E_F. χ is higher for N^↓ than N^↑ at lower temp. as N^↓ is metallic and has a higher DOS at the E_F. At 300 K, χ for N^↓ is 1.59 × 10⁻¹¹ and decreases to 1.57 × 10⁻¹¹ at 800 K due to thermal excitation disrupting spin alignment. For N^↑, χ increases from 1.35 × 10⁻¹¹ at 300 K to 1.47 × 10⁻¹¹ at 800 K as thermal excitation promotes electrons across the E_g, allowing for stronger spin alignment. However, χ for N^↑ saturates around the T_C due to phase transition effects (Fig. 8(d)). The Hall coefficient (R_H) for N^↑ remains positive across this temp. range, indicating hole conduction, while for N^↓, it decreases with increasing temp., signifying metallic behavior (shown in Fig. 8(e)). In N^↓, R_H decreases with temp. as thermally activated carriers dominate, signaling the metallic behavior. For N^↑, R_H increases as thermal excitation generates additional carriers, enhancing the R_H. Further, the molar specific heat capacity (C_V) increases with temp. following the Debye Law, reflecting stronger atomic binding and activation of additional vibrational modes. Despite differing electronic properties, C_V remains nearly identical for both spin channels because it is governed primarily by phonon vibrations in the lattice, which are independent of the SP electronic structure. Finally, the thermal performance of the systems as quantified by the figure of merit,

, is significantly higher for the N^↑ channel at room temp. This is due to its enhanced S and lower κ_e, which collectively results in improved TE performance. At ambient temp., the ZT values for N^↑ and N^↓ are 0.76 and 0.001, respectively, highlighting the superior TE efficiency of the N^↑ channel. This behavior is depicted in Fig. 8(f), which illustrates the temp. dependence and contrast in ZT values for the two spin channels.

4 Strained system

DOS and E_g are crucial descriptors of the electronic properties of materials. The DOS represents the number of electronic states at each energy level, while the E_g refers to the energy difference between the CB and VB. Thus, the orbital hybridization between the O 2p orbitals and the Cr 3d and Cu 3d orbitals plays a central role in determining these properties. Changes in atomic distance, induced by strain, affect the overlap between these orbitals, thereby influencing the E_g and DOS peak height. Interatomic distances are reduced in materials under compressive (comp.) strain, leading to more substantial orbital overlap and hybridization. This typically results in a narrowing of the E_g and a reduction in the sharpness of the DOS peak. Conversely, under tensile (tens.) strain, interatomic distances are increased, which weakens orbital overlap and can lead to a widening of the E_g and a sharper DOS peak. Therefore, we perform DFT calculations to simulate the effects of comp. and tens. strain on the CCCRO QDPO. A comp. strain is applied by reducing the lattice constant, bringing the atoms closer together, and tens. strain is applied by increasing the lattice constant, causing the atoms to move farther apart. To assess the stability of the HM phase in the CCCRO QDPO, we examine the effects of biaxial ([110]) strain on the electronic and magnetic properties within a reasonable range of ±5%. First, the structural robustness of the aforementioned systems is evaluated by plotting the E_for./E_coh. in blue/red with respect to the strain, as shown in Fig. 9(a). Both quantities remain negative throughout the entire applied strain range. However, structures under comp. strain exhibit lower stability compared to those under tensile strain, as the values of E_for./E_coh. shift toward higher (more positive) energies with increasing strain from 0% to −5%. Next, to quantitatively analyze the variations in the electronic structure of the system resulting from strain, we display only the calculated E_g of N^↑ in Fig. 9(b), since N^↓ remains metallic across the entire strain range considered. It is observed that E_g decreases solely from 1.62 to 1.24 eV with the varying compressive strain from 0% to −5%. The change in E_g is mainly due to the charge re-distribution (Cu 3d, Cr 3d and O 2p) state hybridization. It is observed that the Re–Cr bond lengths are least affected for the entire strain ranges, but the Cr–Cu bond lengths increase significantly (from 3.21 Å to 3.70 Å). This also changes the interaction forces among the ions linked with the stretched bond lengths. As a result, the energies around the E_F change (the ground state energies of the CB/VB show a decreasing/increasing trend, respectively). This causes the decrease in E_g. However, small changes were found for the cases of tensile strain compared to the ones without stress (0%).

Fig. 9 Calculated (a) formation/cohesive energy (E_for./E_coh.) in blue/red color and (b) energy gap (E_g) in the spin-majority channel with respect to ±5% biaxial ([110]) strain in the CaCu₃Cr₂Re₂O₁₂ quadruple structure.

Under comp. strain, the interatomic distances between O, Cr, and Cu atoms decrease, bringing their atomic orbitals into closer proximity. This leads to stronger orbital hybridization between the O 2p and the Cr 3d and Cu 3d orbitals. The closer proximity of the orbitals causes the electronic states to become more delocalized, reducing the band gap E_g as the conduction and valence bands move closer together. This narrowing of the E_g indicates a transition toward a more metallic or semiconducting state. Also, the DOS peak height decreases because the increased overlap between the atomic orbitals leads to a broader distribution of electronic states near the E_F. As a result, the DOS peak becomes less sharp due to the more delocalized electronic states. As shown in Fig. 3S of the ESI,† this effect is observed for the Cr and Cu 3d states. Fig. 4S of the ESI† shows similar behavior for the O 2p states under comp. strain. In contrast, under tens. strain, the interatomic distances between the atoms increase, which causes the O 2p orbitals and Cr/Cu 3d orbitals to become more localized. This increases the distance of the O 2p orbitals from the Fermi level E_F, while the Cr and Cu 3d orbitals are pushed even further away due to decreased orbital overlap. The reason for this decrease is the increasing localization of these orbitals, which leads to a reduction in the number of available electronic states near the E_F, thus decreasing the DOS peak height. For the O 2p orbitals, the peak height increases under tens. strain. This is because as the strain increases, the O 2p orbitals become more localized near the E_F, which increases the DOS in this region. Despite the increase in peak height for the O 2p orbitals, the total DOS does not show a significant increase under tens. strain. This is because, although the O 2p orbitals become more localized and their peak height increases, the overall distribution of electronic states remains narrow, and the states of the Cr and Cu 3d orbitals are less delocalized. Therefore, the overall DOS profile does not exhibit a significant broadening, and the peak height for the TDOS remains nearly unchanged.

For qualitative evaluation, we plotted the TDOS for various comp. and tens. strains of the motif, as shown in Fig. 10. This established the fact that N^↑ maintained a definite E_g value, while in N^↓, states lie at the E_F for −1/−3/−5% (comp.) and +1/+3/+5% (tens.) strained systems, as shown in (Fig. 10(a–c) and 10(a′)–(c′)), respectively. Hence, this ensures the robustness of the HM nature of the system under the considered ±5% biaxial strain and maybe above it. Finally, we plotted the calculated m_s values on the Cu/Cr/Re ion as a function of strain in Fig. 11. In the case of the Cr/Re ion, a small reduction/increment in m_s magnitude is observed when comp./tens. strain is adjusted from −1% to −5%/+1% to +5%. However, it almost remains the same for the whole strain range on the Cu ion. So, an increase/decrease in m_s value for the comp./tens. strain can be readily explained by taking the hybridization factor into account. As the comp. strain amplitude increases, the ions move closer together, resulting in stronger hybridization. This leads to greater charge transfer among the ions, reducing the m_s magnitude of the individual magnetic ions. In contrast, the ions move further apart for tensile strained systems, which weakens the hybridization between them. This causes less charge transfer among ions and enhances the m_s value, as illustrated in this article.⁶¹ Additionally, both Cu and Cr ions have “+” signs on their m_s, whereas it’s a “−” for that of the Re ion over the whole strain range. This further endorses the AFM coupling between Cu/Cr and Re, which turns out to have a FiM SO.

Fig. 10 GGA+U calculated total density of states (TDOS) for (a/a′) −/+1%, (b/b′) −/+3%, and (c/c′) −/+5% compressive/tensile biaxial ([110]) strains in the stable ferrimagnetic-I spin ordering of the CaCu₃Cr₂Re₂O₁₂ quadruple structure.

Fig. 11 Calculated partial spin magnetic moments (m_s) on the Cu, Cr, and Re ions for ±5% biaxial ([110]) strain in the stable ferrimagnetic-I ordering of the CaCu₃Cr₂Re₂O₁₂ quadruple structure.

5 Conclusion

These ab initio calculations lay a foundational understanding for developing a quadruple-double perovskite oxide, CaCu₃Cr₂Re₂O₁₂, with potential applications in spintronics. We demonstrated that magnetic properties can be engineered through doping (charge/ion implantation) to create ferromagnetic (FM) and various ferrimagnetic (FiM) spin orderings (SOs) to figure out the magnetic ground state of the system, which is Cu²⁺↑Cr³⁺↑Re⁵⁺↓, while keeping a reasonable total magnetic moment of 5.0 μ_B per f.u. along with partial spin moments of 0.48/2.54/−0.92 μ_B on the Cu/Cr/Re ions. This leads the system into a FiM phase, which is further confirmed by plotting the three-dimensional spin-magnetization density isosurfaces. Interestingly, the system exhibits a half-metallic (HM) nature with a sufficiently large energy gap (E_g) of 1.62 eV in the spin majority channel (N^↑), which is an essential parameter to avoid spin-flipping above room temperature (temp.), with a Curie temperature (T_C) of 340 K. Moreover, the transport properties of the system reveal its metallic nature in the spin minority channel (N^↓) and an insulating trait in the spin majority channel (N^↑), as evidenced by the increasing/decreasing Hall coefficient (R_H) for N^↑/N^↓ and increasing electrical conductivity with temperature. The positive Seebeck coefficient confirms p-type behavior, while thermal conductivity trends reflect the influence of phonon–electron interactions. These findings underscore the material’s potential for thermoelectric (TE) applications, driven by its favorable charge carrier dynamics and heat transport traits. Moreover, the computed T_C is further explained regarding magnetization (M) and magnetic susceptibility parameters. We have explicitly discussed the behavior of the spin-polarized (SP) electron population within the material by applying an external magnetic field (H) of varying strength (2 T and 5 T), enabling the tuning of magnetic properties. Further, it is predicted that the system demonstrates a robustness of its HM state against a reasonable range of ±5% biaxial ([110]) strains. Hence, the combination of a robust HM FiM state, large E_g in N^↑, high T_C (above room temp.), and reasonable TE response of the CaCu₃Cr₂Re₂O₁₂ structure opens pathways for further research in spintronics, offering insight into how materials can be tailored for applications such as memory devices, quantum computing, and spin-based transistors at high temperature.

Data availability

The datasets used and/or analyzed during the current study are available from the corresponding author on reasonable request.

Author contributions

Shehla Yasmeen: writing – original draft, investigations, formal analysis, data curation. S. Nazir: writing – review and editing, validation, supervision, project administration, conceptualization.

Conflicts of interest

The authors declare no competing interests.

Acknowledgements

The computational work was supported by University of Sargodha, Sargodha, Pakistan.

References

1. M. D. Nguyen, H. V. Tran, S. Xu and T. R. Lee, Appl. Sci., 2021, 11, 11301.
2. R. J. Celotta and D. T. Pierce, Science, 1986, 234, 333–340.
3. K.-I. Kobayashi, T. Kimura, H. Sawada, K. Terakura and Y. Tokura, Nature, 1998, 395, 677–680.
4. X. Xiong, J. Kang, Q. Hu, C. Gu, T. Gao, X. Li and Y. Wu, Adv. Funct. Mater., 2020, 30, 1909645.
5. I. Zutic, J. Fabian and S. Das Sarma, Rev. Mod. Phys., 2004, 76, 323–410.
6. A. Fert, Rev. Mod. Phys., 2008, 80, 1517–1530.
7. M. I. Katsnelson, V. Y. Irkhin, L. Chioncel, A. I. Lichtenstein and R. A. de Groot, Rev. Mod. Phys., 2008, 80, 315–378.
8. W. A. Borders, A. Z. Pervaiz, S. Fukami, K. Y. Camsari, H. Ohno and S. Datta, Nature, 2019, 573, 390–393.
9. K. Premasiri and X. P. Gao, J. Phys.: Condens. Matter, 2019, 31, 193001.
10. J. Z. Sun, W. J. Gallagher, P. R. Duncombe, L. Krusin-Elbaum, R. A. Altman, A. Gupta, Y. Lu, G. Q. Gong and G. Xiao, Appl. Phys. Lett., 1996, 69, 3266–3268.
11. T. Obata, T. Manako, Y. Shimakawa and Y. Kubo, Appl. Phys. Lett., 1999, 74, 290–292.
12. W. tin Chen, M. Mizumaki, H. Seki, M. S. Senn, T. Saito, D. Kan, J. P. Attfield and Y. Shimakawa, Nat. Commun., 2014, 5, 3909.
13. R. A. de Groot, F. M. Mueller, P. G. V. Engen and K. H. J. Buschow, Phys. Rev. Lett., 1983, 50, 2024–2027.
14. K. Schwarz, J. Phys. F: Metal. Phys., 1986, 16, 211.
15. R. Wiesendanger, H.-J. Guntherodt, G. Guntherodt, R. J. Gambino and R. Ruf, Phys. Rev. Lett., 1990, 65, 247–250.
16. A. Yanase and K. Siratori, J. Phys. Soc. Jpn., 1984, 53, 312–317.
17. I. Shvets, R. Wiesendanger, D. Burgler, G. Tarrach, H.-J. Guntherodt and J. Coey, J. Appl. Phys., 1992, 71, 5489–5499.
18. S. Jin, T. H. Tiefel, M. McCormack, R. A. Fastnacht, R. Ramesh and L. H. Chen, Science, 1994, 264, 413–415.
19. A. P. Ramirez, J. Phys.:Condens. Matter, 1997, 9, 8171.
20. M. T. Anderson, K. B. Greenwood, G. A. Taylor and K. R. Poeppelmeier, Prog. Solid State Chem., 1993, 22, 197–233.
21. G. King and P. M. Woodward, J. Mater. Chem., 2010, 20, 5785–5796.
22. T. Saha-Dasgupta, J. Supercond. Novel Magn., 2013, 26, 1991–1995.
23. H. Kato, T. Okuda, Y. Okimoto, Y. Tomioka, K. Oikawa, T. Kamiyama and Y. Tokura, Phys. Rev. B:Condens. Matter Mater. Phys., 2004, 69, 184412.
24. P. Majewski, S. Geprägs, O. Sanganas, M. Opel, R. Gross, F. Wilhelm, A. Rogalev and L. Alff, Appl. Phys. Lett., 2005, 87, 202503.
25. O. Erten, O. N. Meetei, A. Mukherjee, M. Randeria, N. Trivedi and P. Woodward, Phys. Rev. B:Condens. Matter Mater. Phys., 2013, 87, 165105.
26. J. Hong, C.-J. Kang and J. Kim, Phys. Rev. B, 2022, 106, 195428.
27. S. Jiang, L. Li, Z. Wang, K. F. Mak and J. Shan, Nat. Nanotechnol., 2018, 13, 549–553.
28. X.-J. Dong, J.-Y. You, B. Gu and G. Su, Phys. Rev. Appl., 2019, 12, 014020.
29. X. Hu, Y. Zhao, X. Shen, A. V. Krasheninnikov, Z. Chen and L. Sun, ACS Appl. Mater. Interfaces, 2020, 12, 26367–26373.
30. A. Avsar, A. Ciarrocchi, M. Pizzochero, D. Unuchek, O. V. Yazyev and A. Kis, Nat. Nanotechnol., 2019, 14, 674–678.
31. Y. Zhao, L. Lin, Q. Zhou, Y. Li, S. Yuan, Q. Chen, S. Dong and J. Wang, Nano Lett., 2018, 18, 2943–2949.
32. X. Cong, Y. Liao, Q. Peng, Y. Yang, C. Cheng, W. Zhang, P. Fang, C. Chen, L. Miao and J. Jiang, RSC Adv., 2015, 5, 59344–59348.
33. P. Nayebi and M. Shamshirsaz, Eur. Phys. J. B, 2020, 93, 170.
34. L. Song, H. Cao, X. Chai, X. Chen, Z. Ye, Y. Zhou, X. Huang, H. Zhu and X. Zheng, Phys. E, 2020, 120, 114057.
35. W. Miao, L. Wang, J. Wang and M. Sun, ACS Appl. Nano Mater., 2023, 6, 20074–20081.
36. Q. Zheng, G. Luo, Q. Liu, R. Quhe, J. Zheng, K. Tang, Z. Gao, S. Nagase and J. Lu, Nanoscale, 2012, 4, 3990–3996.
37. M. Serlin, C. Tschirhart, H. Polshyn, Y. Zhang, J. Zhu, K. Watanabe, T. Taniguchi, L. Balents and A. Young, Science, 2020, 367, 900–903.
38. C. Chen, J. Li and X.-L. Sheng, Phys. Lett. A, 2017, 381, 3337–3341.
39. X. Gao, H. Liu, D. Wang and J. Zhang, Chem. Soc. Rev., 2019, 48, 908–936.
40. Z.-Z. Lin, Q. Wei and X. Zhu, Carbon, 2014, 66, 504–510.
41. J. Xi, Y. Nakamura, T. Zhao, D. Wang and Z. Shuai, Acta Phys–Chim. Sin., 2018, 34, 961–976.
42. W. T. Chen, M. Mizumaki, H. Seki, M. S. Senn, T. Saito, D. Kan, J. P. Attfield and Y. Shimakawa, Nat. Commun., 2014, 5, 3909.
43. X. Wang, M. Liu, X. Shen, Z. Liu, Z. Hu, K. Chen, P. Ohresser, L. Nataf, F. Baudelet and H.-J. Lin, et al., Inorg. Chem., 2019, 58, 320–326.
44. Z. Liu, Q. Sun, X. Ye, X. Wang, L. Zhou, X. Shen, K. Chen, L. Nataf, F. Baudelet and S. Agrestini, et al., Appl. Phys. Lett., 2020, 117, 152402.
45. D. Wang, M. Shaikh, S. Ghosh and B. Sanyal, Phys. Rev. Mater., 2021, 5, 054405.
46. Z. Liu, S. Zhang, X. Wang, X. Ye, S. Qin, X. Shen, D. Lu, J. Dai, Y. Cao, K. Chen, F. Radu, W.-B. Wu, C.-T. Chen, S. Francoual, J. R. L. Mardegan, O. Leupold, L. H. Tjeng, Z. Hu, Y.-f. Yang and Y. Long, Adv. Mater., 2022, 34, 2200626.
47. J. Zhang, Z. Liu, X. Ye, X. Wang, D. Lu, H. Zhao, M. Pi, C.-T. Chen, J.-L. Chen, C.-Y. Kuo, Z. Hu, X. Yu, X. Zhang, Z. Pan and Y. Long, Inorg. Chem., 2024, 63, 3499–3505.
48. P. Blaha, K. Schwarz, F. Tran, R. Laskowski, G. K. H. Madsen and L. D. Marks, J. Chem. Phys., 2020, 152, 074101.
49. J. P. Perdew, B. Kieron and E. Matthias, Phys. Rev. Lett., 1996, 77, 3865–3868.
50. E. C. Calderon, J. J. Plata, C. Toher, C. Oses, O. Levy, M. Fornari, A. Natan, J. M. Mehl, G. Hart, B. M. Nardelli and S. Curtarolo, Comput. Mater. Sci., 2015, 108, 233–238.
51. R. F. L. Evans, W. J. Fan, P. Chureemart, T. A. Ostler, M. O. A. Ellis and R. W. Chantrell, J. Phys.: Condens. Matter, 2014, 26, 103202.
52. F. Mouhat and F.-X. Coudert, Phys. Rev. B:Condens. Matter Mater. Phys., 2014, 90, 224104.
53. K. Samanta and S.-D. Tanusri, J. Phys. Soc. Jpn., 2018, 87, 041007.
54. H. L. Feng, M. Arai, Y. Matsushita, Y. Tsujimoto, Y. Guo, C. I. Sathish, X. Wang, Y.-H. Yuan, M. Tanaka and K. Yamaura, J. Am. Chem. Soc., 2014, 136, 3326–3329.
55. Z. Jirák, S. Krupička, Z. Šimša, M. Dlouhá and S. Vratislav, J. Magn. Magn. Mater., 1985, 53, 153–166.
56. M. Pajda, J. Kudrnovsky, I. Turek, V. Drchal and P. Bruno, Phys. Rev. B:Condens. Matter Mater. Phys., 2001, 64, 174402.
57. E. Sasioglu, E. lu, L. M. Sandratskii and P. Bruno, Phys. Rev. B:Condens. Matter Mater. Phys., 2004, 70, 024427.
58. R. F. L. Evans, U. Atxitia and R. W. Chantrell, Phys. Rev. B:Condens. Matter Mater. Phys., 2015, 91, 144425.
59. A. A. Belik, R. Liu, M. Tanaka and K. Yamaura, Molecules, 2024, 29, 5488.
60. G. K. Madsen, J. Carrete and M. J. Verstraete, Comput. Phys. Commun., 2018, 231, 140–145.
61. K. Yang, R. Wu, L. Shen, Y. P. Feng, Y. Dai and B. Huang, Phys. Rev. B:Condens. Matter Mater. Phys., 2010, 81, 125211.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d5ra02609a

---