Iron in Hydroxyapatite: Interstitial or Substitution Sites?

Iron-doped hydroxyapatite (Fe-HAp) is regarded as a promising magnetic material with innate biocompatibility. Despite the many studies reported in the literature, a detailed theoretical description of Fe inclusions is still missing. There is even no consensual view on what kind of Fe defects take place in Fe-HAp—iron interstitial or calcium substitutions? In order to address these questions, we employ modern first-principles methodologies, including hybrid density functional theory, to find the geometry, electronic, magnetic and thermodynamic properties of iron impurities in Fe-HAp. We consider a total of 26 defect configurations, including substitutional (phosphorus and calcium sites) and interstitial defects. Formation energies are estimated considering the boundaries of chemical potentials in stable hydroxyapatite. We show that the most probable defect configurations are: Fe3+ and Fe2+ substitutions of Ca(I) and Ca(II) sites under Ca-poor conditions. Conversely, Fe interstitials near the edge of the hydroxyl channel are favored in Ca-rich material. Substitutional Fe on the P site is also a probable defect, and unlike the other forms of Fe, it adopts a low-spin state. The analysis of Fe K-XANES spectra available in the literature shows that Fe-HAp usually contains iron in different configurations.


Introduction
Ceramic materials based on calcium apatites are important substances in medicine, biology, ecology and catalysis [1][2][3][4]. Among the apatites, hydroxyapatite (Ca 10 (PO 4 ) 6 (OH) 2 = HAp) is singled out as the main inorganic component of vertebrate bones and teeth. Biogenic HAp differs from the synthetic compound due to the presence of foreign atomic (Mg 2+ , Na + , Cl − ) and molecular (CO 2− 3 , SiO 4− 4 ) ions. Although the role of most impurities is largely unknown, some studies suggest that the presence of a small concentration of iron (below 157 ppm) prevents the loss of calcium from bones and desorption of enamel from teeth [5,6]. Iron-doped HAp (Fe-HAp) has also attracted interest due to its intrinsic magnetic properties [7]. It is a promising material for magneto-resonance imaging agents [8], heat centers for magnetic hyperthermia [4,9], sunscreen filter creams [10], antimicrobial coatings [11], and to be incorporated in drug delivery systems [12,13].
Despite nearly five decades of research on Fe-HAp and related systems [6,[14][15][16][17][18][19], our knowledge regarding the atomistic and electronic structure and iron defects is still very limited. The main difficulty arises from the quality of the samples-single-phase products are hard to synthesize, iron metal and iron oxide often segregate at the surface, and the admixture of stable phosphates, in particular tricalcium phosphate (Ca 3 (PO 4 ) 2 = TCP), is usually observed [16,18]. Discrepancies between the results from different authors are often attributed to different preparation conditions, including sintering temperature and atmosphere, the duration and temperature of doping treatments, the source of precursors, etc. A comprehensive review of experimental methods and results can be found in Ref. [6].
The HAp structure contains two symmetrically distinct positions of calcium, Ca(I) and Ca(II) (Figure 1), both of which are candidates for iron replacement. The pioneering study of Fe-HAp using Mössbauer (MB) spectroscopy by Ok [14] indicated the preference of Fe Ca(I) over Fe Ca(II) , while subsequent MB studies [15,16] concluded on comparable amounts of Fe Ca(I) and Fe Ca(II) defects. In both cases, electron paramagnetic resonance (EPR) data indicated the presence of isolated high-spin Fe 3+ ions with 4 or 6 neighbors, in coexistence with superparamagnetic iron oxide on the surface. The MB experiments of Boda et al. [6] suggest that Fe 3+ Ca(I) defects are more probable when conventional sintering is used during sample preparation, whereas Fe 3+ Ca(II) substitutions appear as a result of hot pressing. Interstitial iron (Fe i ) defects have been reported as well. Analysis of X-ray diffraction (XRD) data shows that iron can occupy 2b or 12i Wyckoff positions [5,[18][19][20]. The first site corresponds to the formation of a linear O-Fe-O structure in the center of the OH channel (site A at Figure 1), while the latter involves a three-fold coordinated Fe atom displacement towards the wall of the channel (site C at Figure 1). However, from XRD and Raman scattering data, Kato et al. [19] argues on the formation of four-fold coordinated Fe. The absence of Fe Ca in these experiments was postulated based on the lack of CO 2− groups [6,16]. These are commonly present in biogenic HAp, and lead to the formation of several defects, including: [CO 2 Besides the lack of agreement between different MB studies, XRD data also holds controversy due to a low electronic contrast between Fe 3+ and Ca 2+ ions [18]. Hence, atomic-level simulations are regarded as a source of clarification. Using empiric interatomic potentials, Jiang et al. [15] concluded that the Fe 2+ ion finds its favorite site at the Ca(I) position, while Fe 3+ favors the Ca(II) site. More recently, Zilm et al. [17] investigated this problem using semi-local density functional theory, but the results were limited to the 2+ oxidation state. Despite enlightening, both studies focused on Fe Ca substitutions only, with alternative substitutional and interstitial forms of iron being overlooked. Hence, we may conclude that the starting point of the present work is a largely unknown picture of Fe impurities in HAp regarding their local structure, oxidation and magnetic states, and optoelectronic properties.
Recently we demonstrated [21] that the use of hybrid density functional theory (hybrid-DFT) can describe the electronic properties of HAp with an accuracy that rivals with highly accurate Green's function (G) method with screened interaction (W), so called GW approximation [22]. The GW and hybrid-DFT methods provided band gaps of HAp close to 7 eV, but generalized gradient functional resulted in a 30 % smaller gap [23][24][25][26][27][28]. Errors of the same magnitude are therefore expected for the calculated optoelectronic properties of defects in HAp, as shown for the case of the OH vacancy [21].
In view of the above, we present a hybrid-DFT study of the most probable Fe-related defects in hydroxyapatite (substitutional and interstitial). The results are interpreted and scrutinized from the perspective of the available experimental data. The paper is organized as follows: Section 2 outlines the methodologies and computational details. Section 3 discusses the results, subdivided in Subsections pertaining the chemical phase diagram of HAp (Section 3.1), the structure and thermodynamics of neutral and charged defects (Sections 3.3 and 3.4, respectively), and a simulation of the Fe K-XANES spectra of Fe-HAp available in the literature (Section 3.5). Finally, Section 4 summarizes the results.

Computation Methods
HAp crystallizes in the form of an molecular ionic crystal of symmetry P6 3 /m (#176 in the crystallographic tables), whose hexagonal unit cell encloses two Ca 5 (PO 4 ) 3 OH units [29,30]. Figure 1 illustrates the supercell used in this study, comprising 8 unit cells (16 HAp formula units). HAp primitive cell contains two and four inequivalent calcium and oxygen sites , respectively. The Ca(I) cation columns are surrounded by O(I) and O(II) from PO 4 anion groups, while mirror-symmetric O(III) sites and Ca(II) ions form a hexagonal channel enclosing the hydroxyl anions O(IV)H. Variation of alignment of OH dipoles in these channels lead to different HAp phases. Namely: (i) a hexagonal disordered phase, with random orientations of OH dipoles; (ii) a hexagonal ordered phases, where OH dipoles are all oriented along the same direction, or (iii) a monoclinic (P2 1 /b) phase, made of two cells repeated along a basal lattice vector (a or b), with the first possessing a ...-OH-OH-... column, while the second one showing an opposite ...-HO-HO-... ordering [31]. The last phase shows anti-ferroelectric properties [32]. Given that the electronic band structure of the above polytypes is rather similar [27], the impact of OH-flipping on the problems being investigated is expected to be minor. Henceforth, we will consider the introduction of iron defects in the ordered hexagonal phase.

First-Principles Defect Energetics
The calculations were performed using density functional theory using the QUAN-TUM ESPRESSO package [33,34]. The many-body electronic potential was evaluated using the hybrid density functional of Heyd, Scuseria and Ernzerhof (HSE) [35,36]. The Kohn-Sham eigenstates obtained using PBE [37] exchange-correlation functional were used to initiate the self-consistent cycle. Due to severe limitations of the semi-local functionals (such as PBE) in describing the electronic structure of HAp in the band gap vicinity [21], the application of hybrid-DFT turns out to be critical in the evaluation of the stability of charged defects.
By opting for a local treatment of the electronic correlation, we expect some error in the description of the d-shell in Fe. Using a GGA+U approach would lead to a narrow gap material, and hence to an artificial hybridization between the d-levels of Fe and those from the conduction and valence band edges of HAp. Unfortunately, it was not possible to combine the non-local HSE functional, with the Hubbard-like local correction to the electronic correlation on the Fe atom.
Core electron states were described by means of optimized norm-conserving Vanderbilt pseudopotentials [38,39], while the Kohn-Sham problem was solved within the plane wave formalism with cut-off energy up to 60 and 240 Ry to expand the wave functions and semi-local potential. The exact exchange operator was evaluated on a grid which corresponded to a plane wave cut-off of 120 Ry. The occupation of states in the vicinity of the Fermi energy were smeared out with a Gaussian function of width 0.002 Ry to improve convergence.
The equilibrium lattice parameters of the unit cell were obtained a = 9.481 Å and c = 6.859 Å using HSE exchange-correlation functional. These values are particularly close to a = 9.417 Å and c = 6.875 Å from an experimental report [30].
Defects were introduced in orthorhombic supercells made up of 8 HAp primitive cells (spacegroup P6 3 /m), containing a total of 352 atoms. The smooth dispersion of the band structure of the large HAp supercell allowed us to use a single point (k = Γ) for sampling the Brillouin zone (BZ). Convergence tests showed that the total energy of such supercells obtained with Γ only sampling, differs by less than 0.1 eV (0.3 meV/atom) from a calculation where the zone was sampled over a 2×2×2 grid of k-points.
Defect containing supercells were first relaxed on a PBE level, since claculation of forces on a hybrid-DFT level is prohibitively expensive for large systems. During the relaxations, coordinates of all atoms were varied until the maximum force became less than 0.4 mRy/Bohr ≈ 10 meV/Å. On a second step the total energy was calculated on hybrid-DFT level for obtained structures. It was shown [40][41][42] that relative errors in energies obtained within this methodology usually have an order of 10 meV.
The formation energy E f of a defect in a crystalline sample can be expressed as [43]: where the main contribution is the energy difference between defective (E d ) and pristine (E HAp ) crystals. R and q denotes the defect configuration and charge state, correspondingly. Last two terms account for stoichiometric and charge differences between the defect and pristine cells. In particular, µ i is the chemical potential of species i which must be added (∆n i > 0) or removed (∆n i < 0) to or from the ideal crystal to create the defect, respectively. The last term in Equation (1) accounts for the exchange of electrons between the defect with charge q and an electron reservoir with chemical potential µ e = E v + E F . Here E v and E F denotes the energy of the valence band top and Fermi level, respectively. The Fermi energy can vary within the band gap (0 ≤ E F ≤ E g ), where E g is the band gap width), depending on the doping of the crystal. The variation of chemical potentials µ i is also limited, with the upper limit taking place if HAp becomes unstable with respect to formation of a compound φ made of n φ i elements of species i and with heat of formation ∆H φ f . The chemical potentials of O, H, Ca and P in standard phases, µ 0 i , were calculated from the energy per atom in molecular oxygen, molecular hydrogen, calcium metal and black phosphorous. The molecules (O 2 in the spin-triplet state or H 2 ) were placed in a simulation cubic box of 20 Å size. The chemical potential of iron was found from the bcc-Fe phase with the BZ sampled over a 16×16×16 k-point mesh. The resulting ground state spin density corresponded to a magnetic moment of 2.80 µ B per cell. The deviation from the experimental value of 2.22 µB results from a known over-localization of the exchange interactions in bulk Fe when using hybrid functionals [44]. An identical result was obtained using projected augmented-wave potentials. Nevertheless, it is noted that hybrid functionals have successfully been used for the study of iron oxides and defects in oxides (see for instance [45,46]).
The heat of formation of the compounds φ used in Equation (2) was estimated from their hybrid-DFT total energies (E φ ) with respect to energies of their constituents in their standard phases (µ 0 i ), Due to periodic boundary conditions, the calculation of charged defects (q = 0) is accompanied by a compensating uniform charge density of opposite sign [47]. The artificial interactions between the periodically repeated charged defects and the background lead to the deviation of calculated periodic total energy E d from E d , in Equation (1). The energy correction E corr , so that E d = E d + E corr , is obtain following the recipe of Lany and Zunger [48], where α M is the Madelung constant of the HAp supercell with edge length L, and is the static dielectric constant ( ≈ 11 [21]). For a singly charged defect the correction amounts to about 0.04 eV. However, this figure can grow considerably in the case of multi-ionized defects (E corr ≈ 0.33 eV for q = +3). The magnetization of the Fe defects was calculated from spin-polarized electron densities n ↑ (r) and n ↓ (r), while formal oxidation states of iron ions were deduced from the occupations of d-orbitals as proposed by Sit et al. [49].

X-ray Absorption Near Edge Structure (XANES)
In order to compare the defect structures found from the first-principles total energy calculations, with those reported in the literature, we simulated the Fe K-XANES spectra for each structure. The calculations were performed using a full-multiple scattering method and the "muffin-tin" approximation for the interatomic potential as implemented in the FDMNES code [50]. The size of the atomic cluster and spectral convolution parameters were adjusted using the spectra of magnetite Fe 3 O 4 as reference. The radius of the cluster was estimated of 7 Å which gives cluster of Fe 1 Ca ∼30 P ∼15 O ∼70 H ∼5 , with exact composition dependent on the structure. The used approach cannot reproduce pre-edge features, which is a well known limitation of a description based on dipolar transitions and the singleparticle approximation [51,52].

Chemical Stability Diagram
The formation of the most probable iron-related defects in HAp depends most notably on the chemical potentials of the several elements involved. Each chemical potential can vary within a certain range, limited by thermodynamic stability conditions of the HAp crystal itself. We estimate the ranges using Equation (2) with respect to a set of boundary phases {φ}. The methodology used to find the structure and energy of each phase was identical to that used for the HAp supercell. This includes the exchange-correlation functional and energy cutoffs. The first candidates for the bordering phases are those involved in the HAp synthesis. There is a variety of production routes, but we will only consider reactants whose species are common to those found in HAp, namely CaO [6], Ca(OH) 2 [53,54] (including the energy, cell volume and bulk modulus) is in line with the usual accuracy of hybrid-DFT. The phosphoric acid and water were calculated in a gas phase (single molecules in large periodic box), hence the lack of cell volume and bulk modulus.  Table 1 and Equation (2), we can estimate the chemical potential ranges within which HAp is thermodynamically stable. We can also reduce the number of independent chemical potentials, for instance to µ H , µ Ca and µ P , by expressing the chemical potential of oxygen as a function of µ i of the remaining elements (c.f. Equation (2)), Let us look first at the range of hydrogen chemical potentials. In this case, phases with give lower (upper) bounds for µ H . In particular, the planes corresponding to CaO, TCP and pure O 2 may form lower bounds, while planes for hydrogen-containing phases may be responsible for upper bounds.
According to our calculations, the HAp stability domain is limited by CaO, Ca(OH) 2 , P, H 2 , DCPD and DCPA phase planes. Figure 2 illustrates the domain of chemical stability of HAp as a convex faceted hull with vertices P i (with i = 1, . . . , 11). These correspond to regions in phase space where three different phases coexists with HAp. Table S1 in Supplementary Materials contains the calculated values of chemical potentials at these intersection points. We note that there is no boundary with Ca metal, with the upper limit for µ Ca − µ 0 Ca being −1.8 eV . This means that if µ Ca = µ 0 Ca was assumed, formation energies of defects involving Ca substitutions would be affected by an error of ∼ 2 eV. In the analysis regarding substitutional iron defects, we will consider the following chemical phase space conditions:

Notation for Defect Structures
Before describing the structure of the iron defects let us first establish a few notation rules. We denote the defect structure as M n S , where M is a substituting species, most notably iron. However, this can also be a vacancy V or group of atoms (e.g., OH). S denotes the position of the defect, which can be a crystalline atomic site (Ca(I), Ca(II) or P) or an interstitial site (i(A), i(B), etc). The flipping of an hydroxyl unit is denoted as OH HO . The number n is the formal oxidation state of iron, estimated from the occupation of the d-orbitals [49]. Table (2) shows the correspondence between the used notation, Kröger-Vink notation and the net charge of the defective supercell. We found a systematic correspondence between the formal oxidation state of iron atom and the defect charge state. This indicates that the gap states of the Fe defects under investigation (which can trap electrons or holes) belong to the 3d shell. Table 2. Correspondence between the defect charge state q (in units of e), Kröger-Vink notation and that used in this work (right-most column).

Notation q
Kröger-Vink This Work

Neutral Iron Defects
We start by considering neutral iron defects based on structures previously proposed in the literature. Table 3 summarizes the relevant data collected. The first seven rows correspond to the Fe Ca substitutions, while the next five to iron interstitials. There is no apparent correlation between the defect structure and sample preparation conditions (e.g., Fe concentration or synthesis temperature).
We considered Ca(I), Ca(II), as well as P substitutions, the later being unexplored in the literature. Figure 3a-c illustrate the substitutional defect structures obtained after relaxation, while Figure 3d-g depict the interstitial structures. The relaxation of the structure where iron was initially set up with the Fe i(A) configuration (Wyckoff position 2b), resulted in another structure, where iron is displaced away from the center of the channel. This new structure is denoted as Fe i(C) (Figure 3f) and approximately corresponds to Wyckoff position 12i. In order to achieve a stable Fe i(A) configuration, one of the nearest OH groups has to be flipped (Figure 3d). Another interstitial position with the iron ion displaced from the center of the channel is labeled as Fe i(B) (Figure 3e) and corresponds to Wyckoff position 2c. The Fe i(C) structure may be slightly modified by flipping of the neighboring OH group, leading to a minute (∼ 0.1 eV) benefit in the total energy. Such small difference is lower than the flipping of isolated OH, estimated as 0.22 eV. Our result suggest that OH flipping may be stimulated by the presence of interstitial iron in the channel.   Finally, we considered iron inserted in the region between PO 4 groups (Wyckoff position 6g), marked as Fe i(D) (Figure 3g). At this location, the negatively charged PO 3− 4 groups are expected to screen the positive charge of the iron ion. Although we scanned the relative stability of other structures, those whose formation energy was above 7 eV were discarded and not investigated further (e.g., Fe OH and Fe i in the vicinity of the Ca(I) column). Table 4 shows the formation energies of the most stable defects, estimated according to Equation (1) using chemical potentials at extreme points of the HAp stability diagram, namely for material grown under Ca-and P-rich (P 7 and P 8 ), Ca-poor (P 3 ) and P-poor (P 1 ) conditions. At Ca-and P-rich conditions, the Fe 0 i(A) defect is likely to be the most probable as it shows the lowest formation energy. Other stable defects (within less that 1 eV above the ground state) are Fe 0 i(B) and Fe 2+ Ca substitutions. Other interstitial sites (Fe 0 i(C) , Fe 0 i(D) ) have higher formation energies. However, they still should be considered since hole or electron trapping may stabilize them.  Figure 2) leads to an easier depletion of phosphorus and to the stabilization of Fe 5+ P as well. Subsequently, at P-poor conditions this effect is further enhanced and Fe 5+ P becomes nearly 4 eV more stable than Fe 2+ Ca , and about 7.5 eV than the most favorable interstitial defect, Fe 0 i(A) . Table 4  The structure of Fe i(A) with two O(IV) neighbors shows Fe-O distances ∼ 0.2 Å longer than those observed by XRD [18]. This may result from the XRD analysis, which gives an ordered structure where oxygen atoms are fixed to crystallographic sites, or may reflect the distances of a positively charged state, where the iron cation is closer to the oxygen anions. Fe i(C) configurations have four-fold coordinated iron. This is at variance with the XRD study of Gomes et al. [18], where three-fold coordinated iron was found, but is in line with XRD results of Kato et al. [19]. The only three-fold coordination of Fe that we find is in the Fe i(B) structure ( Figure 4e). However, the geometry differs from that proposed in Ref. [18], where Fe connects to one oxygen atom from a neighboring PO 4 group and two O atoms from OH groups. In general, the calculated Fe-O distances look slightly overestimated. Of course, the picture could improve when considering charged cells.

Charged Iron Defects
Iron (Fe:4s 2 3d 6 ) is not isoelectronic with respect to the species being replaced (Ca:4s 2 and P:3s 2 3p 3 ). Hence, Fe impurities are expected to create states within the band gap of HAp. These states are localized and may act as hole or electron traps. The description of the Fe-HAp requires the consideration these cases, which can be accounted for by changing the occupation of the highest occupied gap states of the defective supercell. In a real crystal, the capture of a hole (creation of a local positive charge) can be compensated by numerous possibilities: cation vacancies V − Ca , V − H , foreign anion interstitials like [CO 2− 3 ] i , etc., thus leading to a lowering of the Fermi energy (electron chemical potential), We leave the exact mechanisms of charge compensation out of the scope of this work.
An qualitative picture of the charge states allowed for each defect can be found from the respective Kohn-Sham levels in the band gap [61]. Positive charge states (hole trapping) requires the presence of filled states in the gap, while the negatively charged defects (electron trapping) requires the presence of empty states. Figure 4 illustrates the energy of the one-electron states in the band gap obtained from spin-polarized calculations of substitutional and interstitial defects in the neutral charge state. Substitutional iron on Ca sites, Fe 2+ Ca , have six electrons in the 3d shell, but after the first ionization all filled states move below the valence band top. Hence, further ionization would require an energy equivalent to the band gap width (∼7 eV), which essentially tell us that Fe 2+ Ca defects can be single donors, but not double donors. Unoccupied states of Fe 2+ Ca are rather close to the conduction band, suggesting that they are not acceptors either (can not trap electrons).
Iron on the phosphorous site in the Fe 5+ P state (neutral charge state) has three filled one-electron states and 7 empty states in the gap (3d 3 configuration). We will show below that Fe 5+ P can trap one hole or up to two electrons, becoming Fe 6+ P or Fe 4+ P and Fe 3+ P states, respectively.
Neutral interstitial defects show only filled states in the gap. These can trap up to three holes, thus leading to Fe + i , Fe 2+ i and Fe 3+ i , respectively. We found that the structure of the defects depend strongly on the local charge. Figure 5 illustrates some of the most remarkable structural changes of the defects induced by the capture of electrons and holes. Additional reconfigurations are depicted in Figure S1 in Supplementary Information. Table 5    In general, an increase of the positive charge leads to a decrease of Fe-O bond lengths. This effect is strikingly illustrated in Figure 5h,i, which showcase the Fe 3+ P -Fe 6+ P sequence of defects. This rule is understandably violated when there is a change in the coordination of the Fe ion, and therefore, a significant modification of the local electrostatics. Examples are the increases of R Fe-O in Fe 0 i(B) , Fe + i(C) , or Fe 2+ i(D) upon hole capture, (Figure 5f,g), (Figure S1p,q), (Figure S1y,z). The above processes are reversible, i. e., electron trapping at (or hole emission from) the positively charged defects result in the recovery of the longer bond lengths. As expected, HAp cations repel the iron ion. In the case of Fe 3+ i(C) without OH flipping, this effect leads to migration of a proton from OH to a close PO 4 group (see Figures 5d,e). We found that the high-symmetry defect configuration Fe 2+ i(A) is not stable and spontaneously transforms to Fe 2+ i(C) . Hence, upon hole capture by (or electron emission from) Fe + i(A) , the iron moves toward the edge of the OH channel (see Figures 5a,b), The reverse process involves overcoming a barrier (not calculated), and Fe + i(A) is not recovered spontaneously (see Figure 5b,c).
The formation energy of a charged defect is a function of the Fermi energy (c.f. Equation (1)). This dependence is clearly illustrated in Figure 6a,b, which show the results for Ca-and P-rich material and for Ca-poor HAp, respectively. The red shadow area on both diagrams indicates the whole range for the formation energy of Fe P . The upper bound corresponds to P-rich conditions, while the lower bound to P-poor material. The solid red line in Figure 6b shows the formation energy of Fe P under Ca-poor conditions. Thick dashed-dotted lines correspond to formation energies of Fe Ca defects, and thick solid line and other non-solid lines correspond to interstitial defects (see legend).  (Table 4), and slanted -to non-neutral with charge marked by numbers from −2 to +3.
According to our results, the phosphorous substitutions, especially Fe 3+ P (charge state q = −2) is a rather stable species when there is abundance of electrons in the material (n-type HAp). This is even more evidenced in Ca-poor and P-poor conditions, where in p-type material we expect Fe 4+ P , Fe 5+ P and even Fe 6+ P to become more stable and compete with other species, namely Fe i and Fe Ca .
However, the phosphorous substitutions were not previously considered in the literature, since the replacing of phosphorus by iron cation requires the breaking of pretty stable P-O bonds of PO 4 group. Alternatively, the replacement of whole PO 4 group by FeO 4 one will provide the same structural result. The obtained most probable Fe P configuration is tetrahedral coordinated ferric cation ([FeO 4 ] 5− group) is not unusual and can be found in magnetite [52] or in iron-phosphate glass [62].
Regarding substitutional Fe at calcium sites, we find that Fe 3+ Ca(I) is the most probable form under Ca-poor and p-type conditions (Figure 6b). We note that despite Fe 2+ Ca(II) having lower formation energy than Fe 2+ Ca(I) in n-type HAp, the Fe 3+ P species is even more stable, thus making the phosphorus substitution more probable than a replacement of calcium.
In contrast, under Ca-and P-rich conditions ( Figure 6a) the interstitial defects are expected to prevail, especially in p-type and intrinsic HAp. Depending on the Fermi level location, the most probable states are Fe 3+ i(C) , Fe 2+ i(C) and Fe + i(A) defects. The two-fold coordinated Fe 0 i(A) defect, is more favorable in n-type HAp, yet again, Fe 3+ P is more stable and more likely to form.
We can compare the obtained local atomic structure of iron defects in HAp (see Tables 4  and 5) with those reported in the literature (Table 3). Iron substitutions with long distances R Fe-O ≥ 2.2 Å and 6 oxygen neighbors [15,16] are best described by Fe 2+ Ca(I) . The slightly more favorable Fe 2+ Ca(II) substitution has 5 oxygen neighbors, which could explain the result of Boda et al. [6]. The iron with only two neighbors in the study of Gomes et al. [18] could be described by Fe + i(A) , although the small Fe-O distances of 1.7 Å are only reproduced by Fe 6+ P , which is rather unstable. Iron with three [18] or four [19] oxygen neighbors at distances R Fe-O ≥ 1.8 Å can be accounted for by Fe 3+ i(C) or Fe 3+ P defects. In summary, we find that Fe-HAp can contain both substitutional and interstitial defects depending on the preparation conditions. The phosphorus substitutions have iron in low-spin states (M z ≤ µ B ,), making them less useful for many applications envisaged for Fe-HAp. To avoid those defects the synthesis should be performed closer to P-rich and p-type conditions. In that case we expect the formation of high-spin defects with M z = 5µ B (Fe 3+ Ca(I) and Fe 3+ i(C) ) and M z = 4µ B (Fe 2+ i(C) ).

Fe K-XANES of Fe-HAp
The near-edge structure of X-ray absorption (XANES) spectra is particularly sensitive to the details of local atomic structure of the absorbing element [51]. X-ray absorption spectroscopy has been applied to materials without long range order, and that includes the Fe-HAp. We consider the experimental spectra of Fe-HAp published by Gomes et al. [18]. We also keep the notation of the original study regarding the experimental conditions, i.e., 15Fe-500 and 15Fe-1100 which correspond to 15 mol % of Fe per HAp unit cell of samples sintered at 500°C and 1100°C, respectively.
The top three (black colored) curves in Figure 7 shows the experimental data of Gomes et al. [18] for Fe-HAp and for magnetite (Fe 3 O 4 ). The latter was used as a reference for the alignment of the theoretical energy scale to the experimental one. Greek letters mark the main spectral features: pre-edge (α), the main peak (β), its satellite peaks (β and β ), and a more distant peak (γ). The vertical dashed line in Figure 7 provides guidance for the relative positions of the minimum between features β and γ.
In order to determine the types of Fe defects in the 15Fe-500 and 15Fe-1100 samples we simulated the Fe K-XANES spectra for each structure considered in Sections 3.3 and 3.4. Figure 7 shows the simulated spectra of the most probable defect structures. Figure S2 in Supplementary Information shows all calculated spectra. The comparison of experimental and simulated spectra for magnetite Fe 3 O 4 (dotted curves in Figure 7) reveals the main insufficiencies of the simulation method, including (i) the lack of pre-edge features (α) and (ii) a poor reproduction of the β satellite. However, the intensities and energy positions of the main features (β, β and γ) of the magnetite spectrum are correctly reproduced. All spectral features of another reference iron oxide, hematite Fe 2 O 3 , could be reproduced (not shown) using exact the same calculation scheme. In this case, the differences between experimental and simulated spectra show qualitatively the same insufficiencies as for magnetite.  Figure 7. Comparison between experimental Fe K-XANES spectra of Fe-HAp published in Ref. [18] (top three curves) and simulated spectra (all other curves) as described in the text.
The simulated spectra for Fe-HAp show good sensitivity with respect to their atomistic structures and charge states. In most cases the increase of the iron oxidation state shifts the position of the main peak (β) to higher energies. This follows from an increase of the local positive charge and binding energy of K-electrons. An exception is seen for the Fe 3+ i(D) spectrum ( Figure S2d), which is explained by significant changes in its local atomic structure (coordination change from 4 to 6) upon ionization.
The relative intensities I β and I β , corresponding to features β and β from the simulated spectra, are sensitive to the defect type: phosphorus substitutions have I β < I β , calcium substitutions show I β > I β , and most interstitial defects have I β ∼ I β . Additionally, the spectra show considerable differences concerning the position and shape of the γ feature.
The spectrum of the 15Fe-1100 sample shows comparable β and β intensities (I β ∼ I β ), suggesting that most Fe defects are of interstitial character. Conversely, the higher intensity of peak β in the spectrum of sample 15Fe-500 may indicate the presence of substitutional Fe on calcium sites or Fe n i(A) defect (n = 0, 1+). The latter scenario is in line with the conclusions of Gomes et al. [18], where the formation of the high-symmetry Fe i(A) structure was proposed in the 15Fe-500 and 15Fe-800 samples, whereas Fe i(C) -type defects were suggested for the high-temperature treated 15Fe-1100 sample.

Conclusions
We performed a hybrid-DFT study of iron defects in hydroxyapatite. The most favorable structures (and respective charge states) were used to simulate Fe K-XANES spectra of Fe-HAp. These were compared to experimental data available in the literature. The results allow us to draw the following conclusions:

2.
Under P-poor conditions, phosphorous substitutions are the most favorable, resulting in the formation of low-spin Fe defect states. However, p-type HAp may contain high spin iron interstitials Fe 3+ i(C) .

3.
Under Ca-poor conditions, the high-spin calcium substitution, Fe 3+ Ca(I) , is the most favorable species. Fe 2+ Ca defects have higher formation energy comparing to Fe 3+ P and Fe 3+ Ca(I) .

4.
Under Ca-and P-rich conditions, interstitial iron atoms in the OH channel are the most prominent. Depending on the position of Fermi level the most favorable are Fe 3+ i(C) , Fe 2+ i(C) and Fe + i(A) . When compared to Fe i(C) , the last structure involves the flipping of a nearby hydroxyl unit. The OH flipping does not introduce a significant change to the formation energy of Fe 3+ i(C) and Fe 2+ i(C) defects.

5.
High-spin iron defects are Fe 3+ Ca(I) and Fe 3+ i(C) . These are both expected in p-type HAp. Such configurations are expected to be most useful for materials targeting magnetic hyperthermia or magnetic resonance imaging applications. 6.
The comparison of Fe K-XANES spectra of theoretically predicted defect structures with experimental data [18] confirms the interstitial character of iron defects in samples sintered at high (1100°C) temperature, but does not exclude the substitution defects for samples sintered at lower temperatures.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/ 10.3390/nano11112978/s1, Table S1: Chemical potentials in eV at selected points from the phase stability diagram, Figure S1: Illustration of the local atomic structure of iron atom (>4 Å) for all considered Fe-HAp configurations, Figure S2: Simulated Fe K-XANES spectra for the considered Fe-HAp structures. Pseudopotential files, example input script and quantum espresso output file. Atomic coordinates of Fe-HAp structures in POSCAR format.