Structural and magnetic properties of Fe-Co-C alloys with tetragonal deformation: a first-principle study

Fe-Co alloys with induced tetragonal strain are promising materials for rare-earth-free permanent magnets. However, as ultrathin-film studies have shown, tetragonal Fe-Co structures tend to a rapid relaxation toward a cubic structure as the thickness of the deposited film increases. One of the main methods of inducing the stable strain in the bulk material is interstitial doping with small atoms, like B, C, or N. In this work, we present a full configuration space analysis in density functional theory approach for (Fe$_{1-x}$Co$_x$)$_{16}$C supercells with a single C impurity in one of the octahedral interstitial positions and for the full range of Co concentrations $x$. We discuss all assumptions and considerations leading to calculated lattice parameters, mixing enthalpies, magnetic moments, and averaged magnetocrystalline anisotropy energies (MAE). We present a comprehensive qualitative analysis of the structural and magnetic properties' dependence on short- and long-range ordering parameters. We analyzed all unique Fe/Co atoms occupancies at all stoichiometric concentrations possible in 2x2x2 supercell based on 2-atom tetragonal representation. We rely on the thermodynamic averaging method and large sample count to obtain accurate MAE values. We reevaluate several chemical disorder approximation methods, including effective medium methods (virtual crystal approximation and coherent potential approximation) and special quasirandom structures method applied to Fe-Co-based alloys. We observe a structural phase transition from the body-centered tetragonal structure above 70% Co concentration and confirm the structural stability of Fe-Co-C alloys in the tetragonal range. We show the presence of a broad MAE maximum of around 50% Co concentration and notably high MAE values for Co content $x$ as low as 25%. In addition, we show a positive correlation between MAE and mixing enthalpy.


I. INTRODUCTION
Permanent magnets are an indispensable part of modern technology.Among their main characteristic parameters are the energy product (BH) max and coercive field H C .(BH) max determines the efficiency of a permanent magnet and mainly depends on the saturation magnetization M S and coercive field.
Most of the current high-end magnets, with outstanding performance, contain rare-earth elements, such as samarium in SmCo 5 and neodymium in Nd 2 Fe 14 B. However, rare-earthbased magnets have limitations, such as the relatively low Curie temperature of neodymium magnets, which is insufficient for many applications.Moreover, concerns have risen recently about the rare-earth market fragility, which manifested in the so-called rare-earth crisis in 2011 [1,2].Hence, intense research for rare-earth-free permanent magnets has been conducted in the following years.Many potential candidates have been discovered, including MnBi, MnAl, and FeNi magnets [3,4].Currently, rare-earths prices tend towards levels similar to those during the crisis period, encouraging further efforts towards developing efficient rare-earth-free permanent magnets.
Neise et al. [14] showed that the discrepancies between the theoretically predicted MAE and the measured values could be attributed to the virtual crystal approximation (VCA) utilized by Burkert et al.Using 2 × 2 × 2 supercell approach with atoms arrangements modeled according to the most random nearest neighbors patterns, they showed that ordered phases of Fe 1−x Co x have larger MAE than disordered ones, which was confirmed later by Turek et al. [15].They also proposed the preparation of the Fe 1−x Co x epitaxial films along the Bain path [6], which has since been realized by Reichel et al. [16][17][18] on the Au x Cu 1−x buffer, offering a possibility to tailor the lattice parameter in a wide range [19].
Turek et al. further improved the theoretical prediction, ascribing again the calculated versus experimental MAE differ-ence (of the order of 3 -4) to the VCA.Utilizing a more sophisticated method of the chemical disorder approximation, namely coherent potential approximation (CPA) [20], they obtained MAE of much lower and a less sharp maximum of 183 µeV atom −1 spanning a wider range between about 0.5 and 0.65 Co concentration for c/a ≈ 1.22 [15].They also showed that ordering of the Fe 1−x Co x alloys towards L1 0 phase (derived from B2 CsCl structure elongated along the z-axis) could significantly increase the MAE (by a factor between 2 and 3) to 450 µeV atom −1 for Fe 0.4 Co 0.6 and 580 µeV atom −1 for L1 0 Fe 0.5 Co 0.5 -corresponding well with theoretical K u of 520 µeV atom −1 from Ref. [14].Experiments and further calculations have shown that Fe 1−x Co x bct thin films are prone to a rapid relaxation towards the bodycentered cubic (bcc) structure above the critical thickness of about 15 monolayers (about 2 nm) [16,21].
Additions of small interstitial atoms such as B, C, and N were proposed to stabilize the necessary tetragonal distortion by the formation of Fe 1−x Co x martensite phase.Using special quasirandom structures (SQS) method [22] in (Fe 1−x Co x ) 16 C supercells, multiple authors obtained a bct structure with c/a lattice parameters ratio as high as 1.12 -1.17 [23,24].Several experimentally obtained systems have confirmed these predictions [16-18, 25, 26], although there is still plenty of room for further improvements.Two above-mentioned MAE enhancement methods, namely (i) strain induced by a lattice mismatch between two epitaxially grown layers and (ii) spontaneous lattice distortion due to impurities, are summarized in the recent review by Hasegawa [3].
Steiner et al. performed an Fe 1−x Co x case study by averaging over completely random structures in a 2 × 2 × 2 supercell [27].They suggested that proper caution has to be placed on the averaging method since CPA and VCA are effective medium methods that do not describe local structure relaxation and reduced symmetry.Despite their concerns, they obtained MAE values similar to the CPA results reported previously by Turek et al. [15].Since then, many articles have focused on a supercell approach applied to selected cases of Fe 1−x Co x doped with boron [28], carbon [29], and nitrogen [29][30][31], mostly regarding (i) the L1 0 phase derived from B2 (CsCl) structure strained along the z-axis, or (ii) the Fe 0.4 Co 0.6 disordered alloy.For (Fe 1−x Co x ) 2 B, Däne et al. performed a sampling of the full configuration space of the 12atom supercell, again using the argument that VCA and CPA do not correctly describe the distribution of possible values of MAE and the influence of chemical neighborhood and local geometry optimization.They observed a significant spread of the MAE values with an overall average in good agreement with the experiment.They argue that treating a "true" disorder is certainly beneficial.They also noted that it is necessary to average over sufficiently large supercells, as the supercell size can significantly affect the MAE values obtained [32].
The discussion about configuration space analysis is connected with symmetry and ordering in the supercell.Given the vast data set regarding multiple structures in a single crystal system, analysis of ordering towards specific structures is straightforward to implement; it provides more insight into physical phenomena occurring.Works on energy states of closely related structures reach the '30s-'60s of the 20 th century, including contributions from Bethe, Bragg, Williams, Warren, and Cowley in short-range and long-range order analysis methods of that period [33][34][35][36][37]. Recently, a notable example of ordering effects analysis closely related to our work includes research on the FeNi ordering towards the L1 0 phase performed by Izardar, Ederer, and Si [38][39][40].
Here, we present a complete analysis of all stoichiometric (Fe 1−x Co x ) 16 C compositions modeled in a 2 × 2 × 2 supercell.We consider all possible symmetrically inequivalent arrangements of Fe and Co atoms.The aim of the study is to predict the phase stability and intrinsic magnetic properties for the full range of concentrations of the (Fe 1−x Co x ) 16 C system and place it in the frame of works on F-Co, Fe-Co-B, Fe-Co-N, and Fe-Co-C alloys.To achieve it, we study the full configuration space of the 17-atom representation of the Fe-Co-C system and explore this approach to crystallize the most effective method of similar analyses for future applications.

A. System preparation
We used the full-potential local-orbital (FPLO18.00)code [41,42] with the generalized gradient approximation (GGA) exchange-correlation functional in the Perdew, Burke, and Ernzerhof (PBE) [43] parametrization for all calculations.The use of FPLO was dictated by, inter alia, the inherent implementation of the full-potential approach (i.e., omitting the crystalline potential shape approximation), and the expansion of the extended states in terms of localized atomic-like numerical orbitals basis [41,42].The full-potential approach is particularly essential for accurately determining a subtle quantity such as MAE.Another important factor in choosing FPLO is the very high performance of the code, at the expense of the lack of multithreading.In our approach, scaling multiple single-thread calculations up in an embarrassingly parallel manner is the optimal solution.
Initially, we built a 2 × 2 × 2 supercell of the 2atom Fe 1−x Co x body-centered system representation in the P4/mmm space group (s.g.123).The result is a computational cell containing a total of 16 Fe/Co atoms.Initial atomic positions were assumed to be perfect (0, 0, 0) and (1/2, 1/2, 1/2) in each unit cell, and a single C atom was introduced as an octahedral interstitial dopant on the (0, 0, 1/4) site in the supercell.The resultant structure is shown in Fig. 1(a).Structures visualizations were prepared in VESTA software [44].The carbon concentration in the prepared models is about 6 at% and 1.25 wt% (1 C atom per 16 TM atoms).Initial atomic positions were optimized for Co concentrations equivalent to all stoichiometric cases in the 17-atom supercell (Fe 16 C, Fe 15 CoC, Fe 14 Co 2 C, ..., Co 16 C).At this stage, we used VCA for the disorder treatment, 6 × 6 × 6 k-point mesh, 10 −5 density and 10 −7 Ha (∼2.72 10 −5 eV) energy convergence criteria and 10 −3 eV Å −1 force tolerance for initial optimization.Cell volume and c/a optimization were performed based on a thirdorder surface fit to energy versus compuational cell volume in the 160 -208 Å 3 range, incremented by 1 Å 3 and c/a ratios in the 1.05 -1.16 range, incremented by 0.01.Uniaxial elongation of the cell was assumed after Reichel et al. [24].The preparation of the VCA system ended with a full optimization of atomic positions for the minimum of the mentioned fit.We used a scalar-relativistic approach with the same parameters as before.An exemplary resultant structure for the Fe 8 Co 8 C system is shown in Fig. 1

(b).
In the final step of structures' preparation, atomic sites were populated with all possible discrete, stoichiometric, geometrically inequivalent Fe/Co occupations.The equivalency was determined based on the initial, perfect body-centered tetragonal geometry.4 195 unique combinations were obtained out of 65 534 total combinations without repetitions, including 748 unique combinations out of 12 870 for the Fe 8 Co 8 C case alone.The criterion of identity between the combinations was the equity of all interatomic distances between all atom types, i.e., Fe-Fe, Co-Co, Fe-Co, Fe-C, and Co-C in the initial, perfect supercell.It can be proven that it is unambiguous and directly couples each combination with the distribution of minority atoms in the supercell, such as the short-range ordering parameter described later.This approach provided us with a relatively simple method for preliminary analysis.Electron density was then converged in the scalar-relativistic mode, using 9 × 9 × 9 k-points over the entire Brillouin zone, following five additional force optimization steps for every structure to prevent numerical artifacts.For this step of the calculations, convergence criteria were set at 10 −6 density and 10 −8 Ha (∼2.72 10 −6 eV).One of the final Fe 8 Co 8 C structures is presented in Fig. 1(c).
Relevant magnetic parameters were derived based on the converged electron density and systems' energies, as described later.Those include magnetocrystalline anisotropy energies (MAE), mixing enthalpies (∆H mix ), magnetic hardness parameter (κ), Bethe short-range order parameter (σ), Warren-Cowley short-range order parameter (α XY ) for first coordination shell, and long-range ordering parameter towards B2 phase (S ).Specific equations and methods relevant to detailed parts of the presented work are introduced further along-side the results.

B. Assumptions and ensemble averaging methods
We estimate our MAE results for each data point to be within 15% relative error due to relatively low k-point mesh.Obviously, obtaining accuracy within 1% for each considered structure would be highly valuable.However, raising the accuracy would greatly increase the computational cost beyond current capabilities.Obtained system energies and the mixing enthalpies are much more accurate.Bound by this limitation, we focus on qualitative trends and averages in more subtle values, such as MAE.We assume the error imposed by the low k-point mesh for each data point is random and noncumulative.
We utilize thermal averaging after Däne et al. [32] to include influence of non-optimal ground level energy states: where E ν denotes the total energy of a unique atomic arrangement combination ν, MAE ν represents its magnetocrystalline anisotropy energy, and n ν is the number of geometrically equivalent configurations.An important part of the discussion is whether the averaging assumed in Eq. 1 is proper.Foremost, we acknowledge the fact that at room temperature, a vast part of the system does not occupy the ground state, which is calculated in plain DFT.It results, e.g., in the real magnetic moments being lower than predicted.A fact more important for us is that Eq. 1 does not count factors such as the energy barrier height between atom arrangements in the cell.In fact, if the energy barrier is high enough, simple arithmetic averaging should be more appropriate.The height of the energy barrier between the conformations could be obtained by, for example, the nudged elastic band (NEB) method [45,46].However, it would be computationally not yet feasible to obtain heights of all possible transitions [47].Obtaining at least a few values of the barriers in the near future could be beneficial.The solution is, however, not compatible with our methods.Less accurate but less costly linear scaling DFT methods could be utilized to obtain rough values of the barriers.Moreover, this thermodynamic approach results in the configurations' statistical distribution corresponding to slow cooling.We can further assume that despite the obtained result do not rely solely on the most optimal atomic arrangements, the lowest energy structures vastly contribute to the overall MAE.Overall, Eq. 1 certainly works for situations corresponding to slow cooling of the alloy.Hence, it is another assumption in our work that applies to thermal averages.
Apart from the assumptions, an important factor to note is the notation we use to describe various C impurity nearestneighbor patterns.Those designations (Fe-C-Fe, Co-C-Co, and especially Fe-C-Co) should not be mistaken with the common Fe-Co-C system designation, which we also utilize in this work.

III. RESULTS AND DISCUSSION
A. Structural properties  We will first discuss the structural parameters of the alloys under consideration.During the VCA geometry optimization, we observed a structural phase transition from body-centered tetragonal (bct) to face-centered cubic (fcc) structure, which occurs between 11 and 12 Co atoms in the supercell (between 69 and 75% Co concentration), see Fig. 2. It corresponds to the well-known phase transition towards hexagonal closepacked structure for high Co concentration in Fe 1−x Co x .The fcc structure is the closest to the hcp structure we can obtain under the assumed constraints.Although unstable at the standard conditions, the fcc structure for pure Co has been obtained in the high-pressure regime by Yoo et al. [48].
Unit cell volume decreases monotonically with Co concentration after a weak peak for a single Co atom in the supercell, with a significant drop with the transition from bct to fcc structure.Distinct maximum in unit cell volume has been argued by Pauling and other authors, as brought recently by Díaz-Ortiz et al., to be of the same nature as a peak in magnetization (Slater-Pauling curve) [49,50].The weak maximum we obtained stays in contradiction with the expected, Slater-Pauling-like shape of the curve brought to attention by Prinz [51] and successfully reproduced in calculations, e.g., by Díaz-Ortiz et al. [50] and Steiner et al. [27], with a distinct maximum at around 20-30% Co in Fe 1−x Co x .We ascribe this discrepancy to the presence of the dopant atom in the unit cell.Nevertheless, a noticeable positive deviation from Vegad's law is apparent.A similar influence of the small interstitial dopant on the structural (and magnetic) parameters of the system has been observed by Chandran et al. for the (Fe 1−x Co x ) 16 N 2 system [30].
The exact lattice parameters obtained using the VCA in the bct regime are a ranging from 2.75 Å for Fe  [16].The result for equiatomic (Fe 0.5 Co 0.5 ) 16 C (188 Å 3 ) is close to values obtained by Khan and Hong in equiatomic (Fe 0.5 Co 0.5 ) 32 C (about 187 Å 3 ) [29] and (Fe 0.5 Co 0.5 ) 32 N (about 188 Å 3 ) [28].It is also close to the result by Odkhuu and Hong for (Fe 0.5 Co 0.5 ) 16 N (about 190 Å 3 ) [31].Similar values have also been presented for B-doped Fe 1−x Co x alloys by Reichel et al. [24].This slight overestimation of the transition metal alloy lattice parameter is an expected behavior of the applied PBE exchange-correlation functional.Diaz-Ortíz et al. provided an excellent review of structural parameters, magnetic moments, and stabilities of Fe 1−x Co x alloys calculated from first principles.They listed several other results of unit cell volume for Fe 1−x Co x , ranging from 180 to 190 Å 3 per 16-atom cell [50].Most importantly, Delczeg-Czirjak et al. showed that lattice parameters do not exhibit any significant dependency on the atomic configuration exemplified by the C impurity nearest neighbors [23].We followed the assumption of not optimizing lattice parameters for every configuration, as it would be too computationally demanding.
Derived lattice parameters lead to the c/a ratio in the bct regime rising from 1.07 in the case of Fe 16 C to 1.12 for Fe 5 Co 11 C.It is in agreement with the initial assumption of Burkert et al. [5] and following theoretical estimations of uniaxial strain induction by interstitial impurities [23,24].Reichel et al. presented experimental c/a value of 1.05 for B-doped Fe 0.38 Co 0.62 , and c/a for (Fe 0.4 Co 0.6 ) 16 C equal 1.03-1.04,which is lower than the value of approximately 1.10 close to earlier calculations results present in the literature, and also predicted by us.They provided several possible reasons for the observed difference in their work [24].The phase transition from bct to fcc for (Fe 1−x Co x ) 16 C has been also previously reproduced computationally by Delczeg-Czirjak et al. for Co concentration around 65 at% [23].Uniaxial strain in the order of a few percent has been numerously shown to lead to reasonable MAE values [5,14,15,27], which can be further improved, e.g., by buffer-induced effects in thin-film applications [16,18,24,47].

B. Mixing enthalpy and basic magnetic properties versus Co concentration
A basic parameter describing the system is the mixing enthalpy.It provides information about the tendency towards the formation of respective structures instead of separation into their constituent phases (in this case, pure Fe-and Cobased phases).For each structure, we calculated mixing enthalpy ∆H mix between bct Fe 16 C and fcc Co 16 C using equation analogous to the one used by Díaz-Ortiz et al., for convenient comparison with their results [50]: as it, in fact, is the same quantity they calculated for ordered Fe 1−x Co x structures in 2 × 2 × 2 supercells.The results, presented in Fig. 3 Overall, the magnitude of mixing enthalpies suggests good mixing potential, comparable to both TM alloys and steels.Moreover, the shape of the curve suggests the stability of each of the structures relative to neighboring ones, up to 11 Co atoms in the system, or up to the calculated bct-fcc transition.Furthermore, a slight asymmetry in the dependence of mixing enthalpy on x can be observed.On average, the systems closer to the Co-side have lower energies, especially for Co-C-Co systems.However, the absolute minimum for Co-C-Co systems occurs for Fe 8 Co 8 C. For Fe-C-Co, and especially Fe-C-Fe systems, the minimum is moved to the left.The effect of ordering on the mixing enthalpy will be discussed in the following sections.On average, for the region around the equiatomic  Such behavior contradicts the negligence of the direct chemical neighborhood of the impurity atom in earlier works of Khan and Hong [28,29,52].However, we will try to show that despite notable influence on exact quantitative results, neglecting the direct C neighborhood does not alter the qualitative trends in the (Fe In Fig. 3(b), we see a decrease in average spin magnetic moments per TM atom with increasing Co concentration.The average magnetic moment on an Fe atom in Fe 16 C is 2.38 µ B , and the average magnetic moment on a Co atom in Co 16 C is 1.53 µ B .There is a positive deviation from a linear change with x, similar to the Slater-Pauling-like characteristics of unit cell volume versus x dependency.As seen in partial Fe and Co contributions to the average spin magnetic moment, this deviation from a linear trend stems from the Fe contribution.The partial contribution from Co magnetic moments increases linearly.However, as opposed to pure Fe 1−x Co x results reported by Bardos [53], we do not observe a characteristic, sharp maximum related to Slater-Pauling behavior.There is a considerably low deviation in average Fe, Co, and total TM magnetic moments across different configurations.The structural phase transition, between 11 and 12 Co atoms, affects magnetic moments on both Fe and Co atoms, but the change is minimal.
Giannopoulous et al. found magnetization in thin films of (Fe 0.45 Co 0.55 )-C with 20 at% C to be in range of 1600 emu cc −1 [54], which translates to about 2.05 µ B atom −1 .In the literature review performed by Diaz-Ortíz et al., as well as in their own results, we can find average magnetic moments in bcc Fe and bcc Co ranging from 2.13 to 2.35 µ B on Fe atoms and from 1.53 to 1.77 µ B on Co atoms.Their MBPP/PBE (mixed-basis pseudopotential code) calculations for ordered Fe-Co phases yield a total magnetic moment of 2.36 µ B atom −1 for Fe 3 Co DO 3 phase, 2.29 µ B atom −1 for Fe-Co B2 phase, and 2.00 µ B atom −1 for FeCo 3 DO 3 phase [50].Similarly, Chandran et al. reported from VASP/GGA that Fe bcc has a magnetic moment of 2.22 µ B atom −1 , and Co bcc has a magnetic moment of 1.59 µ B atom −1 , not counting for the orbital moment contribution, which for both systems should be around 0.10-0.15µ B atom −1 [30].For C-doped systems, Delczeg-Czirjak et al. found in SPR-KKR/PBE (spin-polarized relativistic Korringa-Kohn-Rostoker) with CPA that the average magnetic moment drops from 2.2 µ B atom −1 in systems with the composition close to Fe 0.4 Co 0.6 to around 1.8 µ B atom −1 in systems with the compositions close to (Fe 0.4 Co 0.6 ) 16 C [23].
Possible giant MAE values are the property, which initially brought attention to the Fe 1−x Co x system.Hence, MAE is among the first characteristics of the system to consider.We calculated MAE according to the formula: where E 100 and E 001 denote the system's energies in the [1 0 0] and [0 0 1] magnetization axis directions (hard and easy axis in the bct structure, respectively).More precisely, we performed a single-step of fully-relativistic calculations in two orthogonal directions, [1 0 0] and [0 0 1], over a charge den-sity self-consistently converged in the scalar-relativistic approach [55], a method proven previously to be both accurate and effective [56,57].Figure 3(c) presents MAE versus x for all configurations, as well as thermodynamical averages according to Eq. 1 and assuming T = 300 K for each Co concentration.To provide an approximate scale in MJ m −3 , we assume a uniform, average cell volume of 186 Å 3 across all Co concentrations and TM atoms configurations.Vertical histograms are scaled to fit the width between points and represent the data spread.There is apparently unimodal distribution of all MAE results for the whole x range among all configurations.A bimodal distribution can be observed closer for the lowest energy configurations results, with MAE values being either very high or near zero.We can observe that MAE varies hugely between configurations, with the absolute maximum for 7 Co atoms in the 16 TM-atom supercell.With more than 11 Co atoms in the system, we observe a rapid decrease and change in the sign of MAE, associated with the phase transition.The high difference in MAE between individual configurations is consistent with similar results for ordering towards L1 0 phase in equiatomic FeNi obtained by Izardar, Ederer, and Si.Though we focus on qualitative trends with low convergence criteria for each data point, they conducted a full convergence for several dozen structures [38,40].
Focusing on qualitative trends, in Fig. 3 [27,58].Moreover, we can observe a few high-MAE configurations among the 5% most preferable ones, see green histograms in Fig. 3(c).The thermodynamically averaged MAE values over 5% of the lowest energy configurations overestimate averages of all symmetrically non-equivalent configurations.It suggests the nonnegligible influence of high-energy (and hence low probability) structures stemming from their quantity.
Our quantitative MAE results can be placed in the context of numerous works describing selected atomic configurations in pure Fe-Co, as well as B-, C-, and N-doped systems, realized both experimentally and by DFT calculations to date.Giannopoulos et al. found experimentally K u for C-doped Fe 0.45 Co 0.55 thin films to be in order of 0.8 MJ m −3 [54], exact same value as obtained by Re-  ichel et al. for (Fe 0.4 Co 0.6 ) 0.98 C 0.02 thin films [17].Reichel et al. have also shown from combined DFT and experimental analysis that the (Fe 0.4 Co 0.6 ) 32 C system possesses slightly lower MAE of the order of 0.5 MJ m −3 and much higher stability for relatively thick films [16].They also reported B-doped Fe 1−x Co x alloys to behave similarly, with a little higher MAE than C-doped system [24].Odkhuu and Hong provide similar results for (Fe 1−x Co x ) 16 N 2 [31].Delczeg-Czirjak et al. found MAE for Fe 6 Co 10 C to be in the order of 51 µeV atom −1 or 0.75 MJ m −3 as calculated in WIEN2k/SQS, higher than their SPR-KKR/CPA calculations (41.6 µeV atom −1 ) [23].For B2 Fe-Co-C and Fe-Co-N systems, Khan and Hong reported MAE values of 0.65 and 0.58 MJ m −3 , respectively [29].Overall, our results agree well with previous calculations and experiments wherever direct comparison is possible.Qualitative trends among major magnetic properties are similar and quantitative results lie close to previous DFT data.However, the dataset we provide is vastly greater than anything currently available in the literature.

C. Magnetocrystalline anisotropy energy and magnetic hardness in relation to the mixing enthalpy
To systematize the dataset, we first analyze the dependency of MAE on the mixing enthalpy.This dependency for all configurations is shown in Fig. 4(a).We see an increase of MAE with lowering the system enthalpy, indicating the preference towards high-MAE structures.There is a significant scatter of values for separate systems around the average.Systems with the dopant atom neighbored by two Co atoms have noticeably larger MAE and lower mixing enthalpy relative to the systems with Fe-C-Fe and Fe-C-Co nearest neighbors (NN) sequence.
To further explore the usability of investigated structures, we calculate magnetic hardness.It is a parameter describing the system resistance towards spontaneous selfdemagnetization and can be defined as [4]: where K 1 is the magnetic anisotropy constant, M S is the saturation magnetization, and µ 0 is the vacuum permeability.A simple empirical rule is that a permanent magnet candidate needs κ greater than 1 to resist self-demagnetization. κ is a useful technical value, as plenty of magnets with relatively low MAE values are manufactured widely due to their high magnetic hardness and low materials cost.
In the case of the Fe-Co-C system, numerous experimental realizations showed a possibility of further amendment of the system to at least double its MAE by tuning the c/a ratio, where interstitial doping can be combined with growth on specifically tailored substrates [17,23,24,26].We also previously showed the positive effect of 5d doping of a similar system [57].Hence, we are interested in promising compositions showing at least semi-hard magnetic properties due to C-doping alone.Skomski and Coey described systems with κ around 0.5 as semi-hard [4].We mark the κ = 0.5 value in Fig. 4(b).In our estimation, we assume K 1 equals MAE, as defined before.Saturation magnetization is derived from the calculated total magnetic moment and cell volume.Thus, we can expand Eq. 4 to the form: where i is the atomic site in the computational cell, M i is the total magnetic moment of the atom occupying site i, and V is the computational cell volume.Figure 4(b) presents the resultant magnetic hardness versus mixing enthalpy relation.It is similar to the MAE dependency on the mixing enthalpy, presented in Fig. 4(a).The magnetic hardness of many configurations exceeds the conventional limit of 0.5 for semi-hard magnetic materials but does not exceed 0.9, remaining below the limit for hard magnetic materials.Odkhuu and Hong reported similar values of κ, ranging from 0.5 to 1 for the (Fe 1−x Co x ) 16 N 2 system [31].From Eq. 5, we can see that there are two main ways to improve the magnetic hardness of the sample.We can either improve MAE or reduce magnetic moment.For permanent magnet applications, we are at the same time interested in as high saturation magnetization as possible.It implies that improving magnetic anisotropy while maintaining relatively high magnetic moments is of interest.Alternatively, achieving high magnetic hardness at the cost of magnetic moment can be beneficial in case of sufficient economic advantage.Relatively negligible changes in total magnetic moment across configurations with the same Co content suggest that, in our case, the MAE changes are a decisive factor in the magnetic hardness variations for different configurations.Either way, both pathways for MAE improvements are feasible in the Fe-Co-C system.

D. Magnetic moments
Looking into the dataset, we focus on average magnetic moments per TM atom in the system, along with the spread of the values in different atomic configurations.Figure 5 summarizes results for exemplary Co concentrations x, 25%, 50%, and 75%.Presented trends in average Fe, Co, and total spin magnetic moments -dependencies on mixing enthalpy and short-range ordering, and their distribution -are representative.Similar results in the literature are scarce, in contrast to analyses of TM magnetic moments on different impurity atom coordination shells, performed by, e.g., Delczeg-Czirjak et al. and Khan et al. [23,28,29,52].
As presented in Fig. 5(a), for low Co concentration, particularly the low enthalpy configurations are associated with high average magnetic moment on Co atoms.It can be explained by the preferred Fe-C-Fe neighborhood, as the dopant atom tends to lower magnetic moments on neighboring atoms.Delczeg-Czirjak et al. shown that TM atoms adjacent to the C impurity in the (Fe 1−x Co x ) 16 C system have significantly reduced magnetic moments [23].For intermediate Co content (exemplified by the Fe 8 Co 8 C system), there is no significant correlation between the average total magnetic moment and mixing enthalpy neither on Fe nor on Co atoms in the bct range.For x = 0.75 (in the fcc range), a preference towards higher Fe and lower Co magnetic moments emerges.We can observe that despite the average total spin magnetic moment on Fe and Co atoms varying considerably between configurations, the average total spin magnetic moment per atom remains almost constant.Spin magnetic moment on Co atoms remain close to 1.5 µ B atom −1 , as predicted by linearity in its linear partial contribution to the total average spin magnetic moment in the supercell.
The trend can be seen more clearly in Fig. 5(d-f) where we present histograms of the average Fe, Co, and total magnetic moments in the structures.In general, the magnetic moment on Fe atoms depends much more on their chemical neighborhood than the magnetic moment on Co atoms.In Fig. 5(d), we see that on the Fe-rich side of the concentration range, for Fe 12 Co 4 C, the total magnetic moment in the system, 2.233 µ B , remains almost constant across all configurations with a triple standard deviation of 0.03 µ B .A similar trend can be observed for the average Fe magnetic moment (2.48±0.08 µ B atom −1 ).However, for average Co magnetic moments (1.57µ B atom −1 ), we can see that the triple standard deviation is relatively high and equals 0.31 µ B atom −1 .On the Co-rich side, for Fe 4 Co 12 C alloy, see Fig. 5(f), we notice that the total magnetic moment in the system also remains almost constant (1.70±0.02µ B atom −1 ).Still, we observe a noticeable variation of 0.16 µ B atom −1 around the average value of Co magnetic moments (1.45 µ B atom −1 ).However, average magnetic moments on Fe atoms, 2.48 µ B atom −1 , vary considerably across different configurations, in the range of ±0.42 µ B atom −1 , which yields almost 34% relative variability between lowest and highest Fe magnetic moment value.In Fig. 5(e), presenting results for Fe 8 Co 8 C, we observe moderate variation in average magnetic moments on both Fe and Co atoms, in the range of 2.53±0.20 µ B atom −1 on Fe, 1.56±0.22µ B atom −1 on Co, and a total magnetic moment in the system of 2.03±0.05µ B atom −1 .
Again, a major driving factor in the spread of magnetic moments across all structures can be the magnetic moment lowering by the neighboring C impurity, which is most prominent on Co atoms, as presented for numerous Fe-Co-based systems by Khan et al. [28,29,52].Moreover, a similar result for N-doped B2 Fe-Co was obtained by Chandran et al..They obtained magnetic moments being reduced from 2.78 µ B to 2.09 µ B between next-nearest and nearest neighbors of the dopant for Fe atoms and from 1.76 µ B to 1.12 µ B for Co atoms, with magnetic moment fluctuations propagating into next-nearest neighbors [30].
To explore other factors influencing magnetic moments in the system, we can use a local neighborhood-based order parameter σ of Bethe, which can be defined for a binary alloy as [35]: where p XY denotes the probability of finding an XY nearest neighbor pair.
Though developed for equiatomic systems, σ derived from Eq. 6 also provides useful information for non-equiatomic binary systems, as it depicts changes in the system with increasing content of NN pairs of non-similar atoms.In that case, σ generally takes values between -1 and 1, with positive values indicating preference towards dislike (in our case, Fe-Co) atomic pairs in the structure and negative values indicating preference towards same-atom type pair (Fe-Fe and Co-Co).However, both minimum and maximum achievable σ changes with the system composition and supercell size, σ min being in <-1,0> range (likewise atom pair affinity) and σ max in <0,1> range (dislike atom pair affinity).
Considering different atomic configurations for particular Co concentrations makes it possible to determine the effect of the former on the values of magnetic moments on individual atoms.Díaz-Ortiz et al. shown for Fe 1−x Co x that the average magnetic moment does not change significantly with ordering [50].Similarly, for special quasirandom structures (SQS), comparing C impurity local neighborhood, Delczeg-Czirjak et al. did not find any relevant change of total magnetic moment in the Fe-Co-C system.The magnetic moment for the (Fe 0.5 Co 0.5 ) 4 C in their work remained at around 1.8 µ B atom −1 [23].Indeed, for (Fe 1−x Co x ) 16 C, we do not see any significant change in the average spin magnetic moment with the local chemical neighborhood, as shown in Fig. 5(gi).Only a slight increase in the average spin magnetic moment with short-range ordering can be observed for the Fe 8 Co 8 C system, presented in Fig. 5(h).It validates effective medium approaches, such as VCA and CPA, to work for disordered (Fe 1−x Co x ) 16 C, similarly to Fe 1−x Co x , the latter pointed by Díaz-Ortiz et al. [50].As for the average Fe and Co magnetic moments, we can see the variation across different structures drops with short-range ordering, indicating a strong contribution from Fe-Co NN interaction.It is consistent with a known strong Fe-Co d orbital hybridization and exchange interaction [31].For any specific minority atoms concentration in our computational cell, the σ range is restricted due to limitations induced by the composition and system size, as described above.

E. Ordering and its influence on magnetic parameters
Apart from average magnetic moment dependence on the short-range ordering, we can explore the ordering effect on other important system properties, including mixing enthalpy, magnetocrystalline anisotropy energy, and magnetic hardness.Figure 6 presents aggregated results for Co content between 3 and 11 atoms in the system, in the bct region.We do not present results for lower Co concentrations because they cover only a small number of configurations and do not have reasonable statistics.Figure 6(a) shows mixing enthalpy decrease with the increase of short-range Fe-Co ordering, i.e., the fraction of Fe-Co pairs among all NN pairs.It might indicate system stabilization by Fe-Co nearest neighbor and Co-Co or Fe-Fe nextnearest neighbor interaction.As previous studies have shown, in the case of the N-doped B2 phase, nearest neighbors Fe-Co exchange integral and next-nearest neighbors Co-Co integral calculated by Odkhuu and Hong contribute the most to magnetic ordering [31].Hence, we ascribe the system stabilization to the same interactions.
Figure 6(b) shows the distribution of MAE in structures with different atomic configurations.Both the highest and lowest MAE for single configuration can be observed for σ equal to 0. For the highest σ values, MAE converges to around 85 µeV atom −1 for the Co-C-Co NN sequence and to around 10 µeV atom −1 for the Fe-C-Fe NN sequence.For negative σ values, which can be associated with low Co concentrations, MAE drops to 0 µeV atom −1 .It can be deduced that the NN ordering influences the MAE by strong Fe-Co interplay.Nevertheless, the factor that contributes most to the overall behavior of the MAE relative to order is the direct immediate chemical neighborhood of the impurity atom.In Fig. 6(c), we present thermodynamic averages, according to Eq. 1. Bars represent the range of MAE values obtained in calculations.
We observe no significant correlation between average MAE and atoms distribution for σ > 0. The most probable MAE for σ equal to 0 is quite high regardless of the dopant neighborhood.The changes in MAE described above are clear, though the scatter of MAE values for various individual structures is substantial.
Taking all the above into account, the configuration space of Fe-Co-C alloys can be somewhat effectively reduced to random nearest-neighbor patterns.Still, it should be done cautiously and can lead to substantial errors, though any anomalies should be evident in the results.Along with the low average magnetic moment dependence discussed above, the lack of strong MAE dependence on the short-range ordering implies that Fe-Co-C retains the properties of a random alloy, similarly to pure Fe-Co.Thus, methods relying on conformational space reduction by neighbor patterns analysis, such as SQS, yield a non-negligible error, similar to effective medium methods, as noted before by Díaz-Ortiz et al. [50].In future studies, it should be decided on a case-by-case basis whether the trade-off between the significant reduction in computation time in approximate (SQS-type) methods and the accuracy and ability to obtain a complete picture of the system in methods that allow order-dependence analysis is justified.
Figure 6(d) presents a similar picture for magnetic hard-ness.We can see that usable magnetic hardness can be obtained for systems around and above σ = 0.For highly ordered systems, the first coordination shell of the dopant plays a key part.Above σ = 0.4, only Co-C-Co and part of Fe-C-Co systems retain magnetic hardness in the semi-hard region.The interesting part is the negative-σ side of Fig. 6(b-d).We observe that where Fe-Fe and Co-Co interactions dominate, MAE and hence the magnetic hardness drops.Although the σ is convenient and effective parameter in analyzing the aggregated results, especially showing the linear decrease of mixing enthalpy with increasing dislike atom pairs content in the supercells, it lacks one property necessary to conduct a complete analysis.It conveys a strict order parameter definition only for equiatomic binary alloy.Namely, its expected value, the same as the value for a completely disordered alloy, is not always equal to zero and depends on minority atoms concentration c m as 4(c m − c 2 m ).For equiatomic alloy, σ equals 0 for completely disordered alloy and takes values up to 1 (or -1) for completely ordered alloys.
To investigate the properties of disordered alloys in a broad concentration range, we use Warren-Cowley short-range order parameter α [36,37], which for the first coordination shell (α Fe,Co I -shortened further to α) can be simplified as: where c A denotes the concentration of atom type A, p AB /2c B = P AB equals the conditional probability of finding an atom of type B at the first coordination shell of the randomly selected atom of type A, and when substituted, gives the exact Warren-Cowley formulation.Structures with all α parameters (for different coordination shells) equal to 0 are disordered, and structures with α i equal to 1 (or -1) are perfectly ordered on coordination shell i.For simplicity, in Fig. 7, we present only MAE versus α dependency.Generally, in an infinite crystal, α takes values between 2c A c B −1 2c A c B and 1 [59].We get only zero to negative α values due to the small computational cell size.Overall, the plot is similar to the positive σ part of Fig. 6(b) and (c) taking into account that preferred dislike atom type coordination is associated with positive σ, but negative α.The most probable MAE value is proportional to the ordering for Co-C-Co systems and, to some extent, for others.Apart from that, we want to highlight three main observations.Firstly, there is a considerable spread in values for random alloys (for α = 0).It is further indicator implying that certain methods of configurational space reduction, like SQS, are inherently predestined to fail in proper Fe 1−x Co x -based alloys MAE predictions, and the uncertainty of such results can be, in fact, substantial.Secondly, same as for σ and similarly to order parameters in recent works by Izardar and Ederer for L1 0 FeNi [38,39], the MAE value converges towards a reasonably high MAE value for perfectly ordered systems.Lastly, in all (Fe-C-Fe, Fe-C-Co, and Co-C-Co) systems, there is a group of configurations that possess high MAE, increasing with ordering.We remind here that α = −1 structures are ordered.For Fe-C-Fe and Fe-C-Co systems, the average MAE value diverges and eventually suddenly drops for high-order structures -a behavior described above for Bethe σ dependencies.
From the comparison of high-order structures calculations to the random Fe 1−x Co x alloy, Díaz-Ortiz et al. deduced that ordered structures are stable, with B2 phase among them [50].Structures predicted by them, namely DO 3 , L6 0 , and B2, as well as similar phases such as L1 2 exhibit a high degree of short-range σ and α ordering, as calculated according to Eqs. 6 and 7. Wu et al. have, similarly, reported stability of Fe-rich DO 3 , and equiatomic B2 phases [58], and Odkhuu and Hong postulated B2 Fe-Co to be a good matrix for low-energy high-MAE N-doped phases [31].One of the very first works on the topic of strained Fe 1−x Co x system treated with CPA effective medium approximation by Turek et al. researched L1 0 ordering influence on MAE in the system [15].L1 0 and B2 phases differ only by lattice parameters c/a ratio, where L1 0 is an fcc-like structure and B2 is close to bcc.As such, we also checked specifically the B2 ordering in the low c/a regime for the C-doped Fe 1−x Co x alloy.
For this purpose, we use the long-range order parameter S of a binary alloy, which is defined in relation to a specific structure, in our case -B2-like Fe 8 Co 8 C. Ordering towards B2 and its equivalent L1 0 phase has been studied in VCA and CPA approaches in several works to date, including one by Turek et al. [15].The parameter S value equal to 1 is associated with a perfect ordering towards the chosen structure (in our case -an ideal crystal in the B2 type), and S equal to 0 represents an absolute lack of the ordering of the given type.Though, a system without ordering towards one structure can be perfectly ordered towards another structure, such as L1 2 structure having a zero S towards L1 0 , both being highly-ordered fcc-like structures and having a high degree of nearest-neighbor ordering.Long-range ordering parameter S can be represented in general as follows [33][34][35]60]: where p denotes the probability of finding an atom of a given type on the expected atomic site.For two-atom type 2 × 2 × 2 supercell and B2 ordering we expand it as: where N I denotes the number of minority atoms close to z = 0 or z = 0.5c plane, N II denotes the number of minority atoms close to z = 0.25c or z = 0.75c plane, and N is the sum of minority atoms in the system.The sites are visualized in Fig. 8.An effectively similar approach has been used recently by Izardar et al. studying equiatomic FeNi L1 0 binary phase [38,39].Parameter S provides a linear scale, similar to one applied by Turek et al. [15].In Fig. 9, we show ordering towards B2 structure dependencies analogous to Fig. 6, presenting results for short-range ordering parameter σ.As the S parameter towards B2 considers only equiatomic systems, the results aggregated are for Fe 8 Co 8 C only.Similarly to σ-dependency, Fig. 9(a) presents a monotonic decrease in mixing enthalpy with B2 ordering in Fe 8 Co 8 C. The energy of configurations with the Co-C-Co NN sequence is, on average, significantly lower than the energy of configurations with the Fe-C-Co NN sequence, which is, in turn, lower than the energy of Fe-C-Fe systems.This fact is independent of the ordering.Perfectly ordered B2 structure with C dopant between two Co atoms possesses the lowest energy.
In Fig. 9(b), we see multiple atomic configurations deviating vastly from the average.In fact, the single highest MAE value, which is twice the average, can be observed for S = 0.5.The associated structure is presented in Fig. 10(c).The qualitative agreement of MAE averages, presented in Fig. 9(c), with the work of Turek et al. is good.We can see that MAE does not follow any specific trend with B2 ordering.For low ordering towards the B2 phase, we can see both very high and very low MAE values.MAE value converges towards a reasonably high 85 µeV atom −1 for perfect B2 order and Co-C-Co configuration.Conversely, for C impurity in the Co plane (neighbored by two Fe atoms), MAE converges towards a low value of approximately 10 µeV atom −1 .These are exactly the same MAE values as for most positive sigma and most negative alpha parameters, see Figs. 6 and 7.It is, in fact, the same structure, visualized further in Fig. 10.Magnetic hardness versus B2 ordering, shown in Fig. 9(d), have to be similar to MAE since the system magnetization has been shown above to not depend on the ordering.The main conclusion is that for higher ordering, only systems with Co-C-Co and Fe-C-Co NN sequences possess usable magnetic hardness.Similarly to σ, for low B2 ordering, we can still observe many individual atomic arrangements with hardness above 0.5.
It might be tempting to dive more deeply into the evaluated atomic occupation configurations individually, with particular emphasis on the high-symmetry structures.However, such analysis is beyond the scope of this work, as we rely on error cancellation due to the high sample count.A detailed look at the specific structures would require a much finer k-point mesh and fine atomic positions optimization of such atomic arrangements.Nevertheless, to emphasize possible further paths of Fe-Co-C system investigation, we present in Fig. 10 four selected low-energy, high-MAE structures (a, b, c, and e), as well as a high-energy, low-MAE, perfectly ordered B2 structure (d).We found that high-order structures for as low as 25% Co concentration can indicate usable magnetic properties.Interestingly, the lowest energy structure for Fe 12 Co 4 C is the Co interlayer in the plane farthest away from the C impurity.Since the price of Fe is negligible in the overall price of an Fe-Co alloy, those are promising candidates for future permanent magnets.As for qualitative trends, we observe the L1 2 structure among the lowest energy systems for high Co concentrations in the fcc regime.Despite the structure changes towards bct with lowering of the Co content, the atomic occupations for low-energy Fe 12 Co 4 C remain the same as in the high-Co L1 2 phase, presented in Fig. 10(e).

IV. SUMMARY AND CONCLUSIONS
We conducted a full configuration space analysis for 2 × 2 × 2 (Fe 1−x Co x ) 16 C supercell based on a 2-atom bodycentered tetragonal unit cell, with a single C impurity at one of the octahedral interstitial positions in the supercell.The calculations were performed using density functional theory (DFT) with the generalized gradient approximation (GGA) using the full-potential local-orbit scheme (FPLO18).
In our tetragonal (Fe 1−x Co x ) 16 C supercells, we observe a structural phase transition from a body-centered tetragonal (bct) to a face-centered cubic (fcc) structure at a Co concentration of about 70 at.%.The lattice parameter c/a ratio in the bct region ranges from 1.07 to 1.12.We calculated relevant  magnetic properties for all non-equivalent Fe/Co atoms arrangements in the computational cell.Since DFT calculations are, by definition, performed for a temperature of 0 K (for the ground state), we used thermodynamic averaging with an assumed temperature of 300 K in determining the average mag-netocrystalline anisotropy energy (MAE) values.Although, as previous experiments have shown, the structure expected above the critical Co concentration (x ≃ 0.7) is hexagonal, the assumed tetragonal geometry of the supercell does not allow this and leads to an fcc structure.
One of the basic features of the supercell geometry we analyzed is the first coordination shell of the C dopant atom.The C atom has two nearest neighboring sites which can be occupied by two Fe atoms, two Co atoms, or one Fe and one Co atom.We found that for low Co concentrations, structures with impurities adjacent to two Fe atoms become more stable.The expected result of the stabilization of the (Fe 0.5 Co 0.5 ) X C alloys by the Co-C-Co nearest neighbor sequence for medium to high Co concentrations is also confirmed in our results.
Although we observe a rather large scattering of magnetic moments for different configurations on both Fe and Co atoms, the total magnetic moment in the supercell remains more or less constant.Average (spin) magnetic moments decrease with increasing Co content, without a clear maximum for intermediate concentrations.
Positive MAE values in the bct region indicate a uniaxial magnetocrystalline anisotropy and show a broad maximum around medium Co concentration (x ≃ 0.5).The calculated course of MAE as a function of Co concentration is in very good quantitative agreement with experimental data, which is a noteworthy improvement over effective medium methods.The magnetic hardness of many configurations exceeds the conventional limit of 0.5 for magnetically semi-hard materials but does not exceed 0.9, remaining below the limit for hard magnetic materials.In addition, for relatively low Co concentrations, on the order of 25%, we have identified a number of energetically stable structures with high MAE values and potential economic significance.
The calculated mixing enthalpy of considered Fe-Co-C alloys is the lowest at around 50% Co concentration.Moreover, the general trends indicate that higher values of MAE (and magnetic hardness) correlate with more negative values of mixing enthalpy, which It shows the better structural stability of high-MAE atomic configurations.
A significant part of the discussion is devoted to determining the effect of ordering on the magnetic properties of the compositions under consideration.We focus on the Bethe and Warren-Cowley short-range ordering parameters and the ordering parameter towards the arbitrarily chosen B2 (CsCl) structure.In the largest range of values of the Bethe shortrange ordering parameter, its increase correlates with an in-crease in MAE, while for the highest values of the parameter (above 0.2) we no longer track correlation.Furthermore, we observe no significant correlation between MAE and the value of the Warren-Cowley short-range ordering parameter and the ordering parameter towards the B2 structure.The direct neighborhood of the impurity dominates MAE value dependencies.On the contrary, we see a clear decrease in the value of the enthalpy of mixing (higher stability) as shortrange and long-range ordering parameters increase.
In summary, we present a relatively simple and effective method for averaging multiple configurations to predict accurate MAE values for the Fe-Co-C system.We show that the method can be made even more efficient by averaging a few percent of the most energetically favorable structures, with little loss in accuracy.In addition, the Fe-Co-C system is a good matrix for further modifications (e.g., induction of additional stresses) stabilized by the Fe-Co nearest neighbor interactions.Considering that B-, C-, and N-doped Fe-Co alloys possess similar structural and magnetic properties, further research of Fe/Co ordering in interstitially-doped Fe-Co can provide much-needed insight towards efficient, rare-earthfree permanent magnet development.

Figure 1 .
Figure 1.Examples of prepared and obtained crystal structures of Fe 8 Co 8 C. Initial cubic supercell -input to virtual crystal approximation (VCA) relaxation (a), structure resultant from VCA geometry optimization (b), and one of the final structures with VCA atoms substituted by Fe and Co atoms (c).Iron, cobalt, and carbon atoms are presented in dark red, light blue, and black, respectively.

)
Number of Co atoms in the supercell

Figure 2 .
Figure 2. Dependency of lattice parameters (c/a) ratio (black) and unit cell volume (red) on Co concentration x in (Fe 1−x Co x ) 16 C system, calculated with FPLO18 in virtual crystal approximation (VCA) with PBE exchange-correlation potential.The dashed line denotes the structural phase transition between body-centered tetragonal (bct) and face-centered cubic (fcc) structures.

5
Co 11 C to 2.82 Å for Fe 15 Co 1 C, and c ranging from 3.01 Å for Fe 5 Co 11 C to 3.07 Å for Fe 12 Co 4 C. Resultant optimized volume of the bct systems ranges from 185.8 Å 3 for Fe 5 Co 11 C to 192.3 Å 3 for Fe 15 Co 1 C. Consistency with Fe 16 C supercell volume obtained by Delczeg-Czirjak et al. [23] in VASP code (about 196 Å 3 ) is good, as well as comparison to experimental value (about 183 Å 3 ) obtained by Reichel et al. for (Fe 0.4 Co 0.6 ) 0.98 C 0.02 (a), correspond well with the aforementioned data for Fe 1−x Co x .The absolute values of ∆H mix (up to 8 mRy atom −1 ) are only slightly lower in comparison with up to 9 mRy atom −1 calculated by Díaz-Ortiz et al. [50].It indicates the stability of both disordered and ordered (Fe 1−x Co x ) 16 C alloys with a minor structure destabilization by the dopant.
Fe 8 Co 8 C, the energy of the systems with C impurity neighbored by two Co atoms is lower compared to systems with the C atom adjacent to two Fe atoms or one Fe and one Co atom.This is consistent with the observation by Delczeg-Czirjak et al. [23] that the energy of Fe-Co-C systems depends mainly on the direct chemical neighborhood of the impurity atom, with a preference towards Co-C-Co nearest neighbors sequence.A similar effect has been calculated by Chandran et al. for N-doped Fe and Fe 1−x Co x [30].

Figure 3 .
Figure 3. Mixing enthalpy (∆H mix ), average spin magnetic moments per transition metal atom (M), and magnetocrystalline anisotropy energy (MAE) per transition metal atom versus Co concentration x in (Fe 1−x Co x ) 16 C system, as calculated using FPLO18 and PBE exchange-correlation potential for all nonequivalent Fe/Co site occupancies in a 2 × 2 × 2 supercell.In panel (a), the light blue and dark red colors represent systems with, respectively, two Co and two Fe atoms neighboring the C impurity.Dark grey represents systems with one Fe and one Co atom neighboring the dopant.For readability, the plotted points are slightly shifted for Co-C-Co and Fe-C-Co configurations.The light blue and dark red colors in the (b) panel represent the average contribution of Co and Fe magnetic moments, respectively.Dark gray is the sum of both.Gray histogram (c) represents the aggregation of all results, while the green one represents 5% of the most energetically favorable atomic arrangements.The circles represent respective average values, and the solid lines are averaged splines to guide the eye.Dashed lines indicate the structural phase transition between bct and fcc structures.
(c), we see a broad maximum for x ≃ 0.25 -0.75.According to Eq. 1, we obtained an average MAE of 0.75 MJ m −3 for Fe 8 Co 8 C. MAE decreases by around 20% between x = 0.5 and x ≃ 0.3.It is in contrast to a rapid drop in MAE for low Co concentrations reported by Delczeg-Czirjak et al. (65% drop between x = 0.6 and x = 0.3).We obtained nearly the same MAE values for x ≃ 0.6 and x ≃ 0.3.Intriguingly, we observe several configurations with relatively high MAE values for such low Co concentration as 0.25.Our findings of no- enthalpy (meV atom -1 ) Mixing enthalpy (mRy atom -1 ) Co-C-Co NN sequence Fe-C-Fe NN sequence Co-C-Fe NN sequence

Figure 4 .
Figure 4. Magnetocrystalline anisotropy energy (MAE) (a), and magnetic hardness (b) versus Fe 16 C and Co 16 C mixing enthalpy in the (Fe 1−x Co x ) 16 C system (from 3 to 11 Co atoms in the supercell).The light blue color denotes systems with two Co atoms neighboring the C dopant, the dark red color denotes systems with two Fe atoms neighboring the impurity, and the black color denotes systems with the C atom neighbored by one Fe and one Co atom.In panel (b), the dashed line for hardness equal to 0.5 indicates a semi-hard magnetic materials threshold.The results were obtained in FPLO18 with the PBE exchange-correlation potential.

Figure 5 .
Figure 5.The average spin magnetic moment per atom versus the mixing enthalpy of Fe 16 C and Co 16 C (a-c), the same value presented as a histogram (d-f), and versus Bethe short-range ordering (g-i), obtained for all geometrically inequivalent Fe/Co arrangements in a 2 × 2 × 2 supercells of (Fe 1−x Co x ) 16 C, calculated with FPLO18 code and PBE exchange-correlation potential.Results are presented for 4, 8, and 12 Co atoms in the supercell.The Fe 4 Co 12 C fcc structure is metastable.

Figure 6 .
Figure 6.Dependence of mixing enthalpy (a), magnetocrystalline anisotropy energy (b, c), and magnetic hardness (d) on short-range ordering parameter σ in (Fe 1−x Co x ) 16 C structures with from 3 to 11 Co atoms in the supercell.The light blue color denotes systems with two Co atoms neighboring the C dopant, the dark red color indicates systems with two Fe atoms neighboring the impurity, and the black color denotes systems with the C atom neighbored by one Fe and one Co atom.Results were obtained using the FPLO18 code with PBE exchange-correlation potential.Fe-C-Fe and Co-C-Co data points are slightly shifted for better readability.Panel c presents thermodynamic averages according to Eq. 1, and error bars denote the maximum and minimum calculated values.

Figure 7 .
Figure 7. Magnetocrystalline anisotropy energy versus Warren-Cowley short-range order parameter for the first coordination shell of transition metal atoms in (Fe 1−x Co x ) 16 C (from 3 to 11 Co atoms in the supercell).Points represent 300 K thermodynamical averages according to Eq. 1. Results were obtained using the FPLO18 code with PBE exchange-correlation potential.Fe-C-Fe and Co-C-Co data points are slightly shifted for better readability.Bars denote the minimum and maximum values obtained in calculations.

Figure 8 .
Figure 8. Graphical representation of atomic sites relevant in Eq. 9.The presented structure is a 2 × 2 × 2 supercell with a single octahedral C atom (black).Sites I are located close to z = 0 or z = 0.5c plane, and sites II lie close to z = 0.25c or z = 0.75c plane.

Figure 9 .
Figure 9. Dependence of mixing enthalpy (a), magnetocrystalline anisotropy energy (b, c), and magnetic hardness (d) on long-range ordering parameter S in Fe 8 Co 8 C. The light blue color denotes systems with two Co atoms neighboring the C dopant, the dark red color indicates systems with two Fe atoms neighboring the impurity, and the black color denotes systems with the C atom neighbored by one Fe and one Co atom.Results were obtained using the FPLO18 code with PBE exchange-correlation potential.Fe-C-Fe and Co-C-Co data points are slightly shifted for better readability.Error bars on the c panel denote maximum and minimum calculated values.

Figure 10 .
Figure 10.Panels (a-c) present exemplary obtained low-energy, high-MAE, and high-symmetry supercells: Co interlayer separated from C impurity by half of the supercell (a), Co-C-Co B2 (b), and highest MAE Fe 8 Co 8 C (c).Panel (d) presents the high-energy Fe-C-Fe B2 structure, and panel (e) presents low-energy, high-MAE Fe 4 Co 12 C L1 2 .In panel (e), the alternative fcc representation of the L1 2 structure is presented with blue lines.Supercell lattice parameters were optimized in FPLO/PBE with virtual crystal approximation and atomic positions were optimized for a few steps in every atomic position occupancy.
[50]Co x )16C system and possibly in other interstitially doped Fe 1−x Co x systems.Surprisingly, we can observe a tendency towards the energetic preference of systems containing Fe-C-Fe nearest neighbors sequence for low Co concentrations.The rapid increase in mixing enthalpy for Co-rich systems is consistent with mixing enthalpies calculated by Díaz-Ortiz et al. and instability of Co-rich bct alloys observed in experiments[50].
table, positive MAE for low Co concentrations contradict earlier results obtained with effective medium methods.VCA and CPA reported negative MAE for low Co concentrations Fe-Co alloy, as seen on MAE versus c/a versus x maps by Burkert et al. and Turek et al. [5, 15].On the other hand, it is consistent with the findings of Steiner et al., who, from Fe 1−x Co x supercells, reported positive MAE in a much wider Co concentration range, and Wu et al., who reported high MAE for Fe 12 Co 4 C and Fe 11 Co 5 C