Ferroelectric polarization of hydroxyapatite from density functional theory

Abstract

The theoretical ferroelectric polarization of the low-temperature (monoclinic, P2₁) phase and the high-temperature (hexagonal, P6₃) phase of hydroxyapatite Ca₁₀(PO₄)₆(OH)₂ is calculated based on the density functional theory (DFT). In the monoclinic structure, the value of ferroelectric polarization is found to be 9.87 μC cm⁻² along the [001] direction. In the hexagonal structure, the ferroelectric polarization is 7.05 μC cm⁻² along the [001] direction. The main contribution to the electric polarization comes from ordered hydroxyl OH⁻ anions for both phases, although the inorganic Ca₅(PO₄)₃ apatite framework also gives a non-negligible contribution. A detailed analysis of ferroelectric polarization and structural change of the hydroxyapatite is presented for a better understanding of this important biomaterial.

---

Introduction

Hydroxyapatite (HAp) Ca₁₀(PO₄)₆(OH)₂, the main constituent of human bone¹ and teeth enamel,² has been the subject of a considerable number of both experimental and theoretical studies.^3–24 HAp's applications span from dental to orthopedic implants, including bone tissue engineering scaffolds, coatings, fillers and many others.²² Pasteris et al.²⁵ used laser Raman microprobe microscopy to study the degree of hydroxylation and the state of atomic order of several natural and synthetic calcium phosphates including apatite of biological, geological and synthetic origins. However, a weak OH-band occurs in dentin whereas a strong OH-band was suggested in tooth enamel.²⁵ In addition, recent studies have also shown that HAp particles inhibit the growth of many kinds of cancer cells.¹³ Perhaps as important as water for nature, HAp has also been chosen as one of the key ingredients of the molecular machinery of vertebrate life.¹⁶

Although being one of the prevalently studied biomaterials for orthopedic, dental, protein purification and stem cell applications, there is less research done on the electrical properties of hydroxyapatite.¹⁹ Since the prediction of the possibility of piezo- and pyro-electricity in hydroxyapatite, several experimental and theoretical reports have led to the new level of understanding the electrical behaviors of calcified tissues in vertebrates.¹⁹ Both bone and dentine show interesting electrical properties foremost of which is piezoelectricity, i.e. the possibility of generating electrical charge at the surface when it is subjected to a stress.²⁶ This surface charge induced by pressure has been considered to be connected with the generation of the bone functional shape in accordance with the Wolff's law.²⁷ Bone also possesses a spontaneous dipolar electrical polarization due to its pyroelectric property,²⁸ which could have functional distribution and physiological effects. Also, the ability to create discrete electrostatic domains on nanocrystalline films of hydroxyapatite will open the possibility of understanding how surface charge influences biological interactions. Kumar et al. presented some results of in vivo experiments in the literature.²⁹ Positive and negative polarization of HAp plates implanted in calvarial bones of rats resulted in an enhanced osteoblast activity. Also, protein adsorption onto the HAp surface implanted in canine femora and rat tibiae, respectively, resulting in improved osteoconduction on both positively and negatively polarized HAp surfaces in vivo was reported. Though there is evidence that the polarized surface can affect the biological response, further researches are needed to understand the in vivo behaviors of charged HAp surfaces.

The two principal crystalline components of bone are the mineral HAp and the protein collagen, both of which may present piezoelectricity.¹⁶ In the case of HAp, piezoelectricity is due to the polarization of OH⁻ groups confined to the channel formed by the overlapping hexagonal calcium atoms. According to Rietveld analysis carried out on the synthetic sample, HAp has been found to be predominantly hexagonal at room temperature,⁷ but the diffraction pattern of this material might be interpreted as a mixed phase: 23% monoclinic P2₁/b and 77% monoclinic P2₁.

Zahn et al.⁸ carried out molecular dynamics simulation and verified that the disordering of hydroxide ions orientations occurs in such a way that collective reorientation of OH⁻ ion rows is not observed. It appears to be a property solely due to the inherently disordered, nonstoichiometric P6₃/m HAp and not of its pure, monoclinic, P2₁/b counterpart.³⁰ Furthermore, the crystal structure of stoichiometric HAp is monoclinic (P2₁/b),⁴ but due to foreign ion inclusion, vacancy formation, and Ca²⁺ depletion, it becomes hexagonal (P6₃/m)³ for more disordered biological HAp. Horiuchi et al.²¹ reported experimentally that an antiferroelectric (AFE) to paraelectric (PE) phase transition accompanied by a critical slowing down and a sharp increase in dielectric relaxation strength was not observed at the phase transition. This implies a weak interaction between the dipoles of the neighboring OH⁻ columns. Tofail et al.⁶ employed calculations with empirical interatomic potentials to show that nonpolar–polar transition could be induced by temperature. In 2013, Lang et al.¹² demonstrated that hydroxyapatite exhibits ferroelectricity by using piezoresponse force microscopy. Recently it was shown that polarized HAp (with an electrically charged surface) is more active, which is more attractive for interactions with living cells, attaching more living cells (osteoblasts) on its surface.^10,12,17 Such effect is significant for medical applications, particularly for implant coatings. Therefore, a more detailed theoretical study about the ferroelectric (FE) polarization is mandatory. To the best of our knowledge, the theoretical analysis of ferroelectricity from first principles is missing in the literature.

In the present paper, we report on the density functional theory (DFT) study of monoclinic hydroxyapatite and its hexagonal phase. Our calculation allows the direct computation of the ferroelectric polarization of HAp, an important information otherwise difficult to access from experimental data.

The paper is organized as follows. First, we give a brief account of the computational details. Then, we discuss our results by considering two structures, namely the low-temperature (monoclinic, P2₁) and the high-temperature (hexagonal, P6₃) phases, and their calculated ferroelectric polarization. The total ferroelectric polarization is further decomposed into the hydroxyl and calcium apatite contributions.

Computational methods

All investigations were performed using density functional theory³¹ as implemented in the Vienna ab-initio Simulation Package (VASP).^32,33 Kohn-Sham equations were solved using the projector-augmented-wave (PAW) method with the PBEsol exchange–correlation functional.³⁴ The energy cutoff for the plane waves expansion was set to 600 eV, 4 × 2 × 6 and 4 × 4 × 6 Monkhorst–Pack grids of k points were used for the monoclinic (P2₁/b or P2₁) and hexagonal (P6₃/m or P6₃) structures respectively. The convergence in total energy and Hellmann–Feynman forces were set to 1 μeV and 0.01 eV Å⁻¹, respectively. For all structures, we have chosen the experimentally determined lattice constant parameters and performed relaxation of internal atomic positions. The ferroelectric polarization has been calculated using the Berry phase approach by building an appropriate PE reference state and considering the suitable path in the configurational space connecting PE to the FE state. For the calculation of the ferroelectric polarization, it is convenient to introduce a reference paraelectric structure, i.e. centrosymmetric. We can choose for simplicity an antiferroelectric structure. By comparing the reference structure and the ferroelectric one and connecting them through an appropriate linear interpolation of the atomic positions, one can avoid the spurious presence of polarization quanta, which would lead to erroneous estimates of polarization. Furthermore, we can also analyze the atomic displacements responsible for ferroelectricity. In our work, both reference structures are obtained by displacement of the atomic positions. In particular, we divide the structures into two parts: one made of distinct hydroxyl units (OH⁻), labeled as H, and the other consisting of the inorganic framework [Ca₅(PO₄)₃]⁺, called Ap. In PE structures, for the H part, the hydrogen and the oxygen atoms are equally distributed along [001] direction, forming a paraelectric arrangement of O and H atoms. At the same time, we also transform the inorganic framework Ap into a centrosymmetric structure. Therefore, the FE state (λ = 1) is obtained from the centric (λ = 0) state by shifting of the different functional groups.

The ferroelectric polarization is defined as ΔP = P^FE − P^PE = ΔP_ion + ΔP_el, where ΔP_ion is the sum of the products of the deviations in the position of each ion in the unit cell (with respect to the centric reference phase) with the nominal charge of its rigid core; the electronic contribution ΔP_el has been calculated by using the Berry phase approach^35,36 which correctly accounts for quantum-electronic effects. Graphical rendering of the crystal structures was performed using VESTA.³⁷

Results

The HAp unit cell lattice parameters are fixed at experimental values³⁸ a = b = 9.423 Å and c = 6.883 Å for the high-temperature (hexagonal, P6₃/m) phase, and a = 9.426 Å, b = 18.856 Å and c = 6.887 Å for the low-temperature (monoclinic, P2₁/b) phase as shown in Table 1 and Fig. 1. The hexagonal HAp phase consists of 44 atoms in the unit cell, while the monoclinic HAp structure has 88 atoms because the unit cell is doubled along the b axis. For the low-temperature (monoclinic, P2₁) and the high-temperature (hexagonal, P6₃) structures, the atomic positions have been relaxed under fixed experimental lattice parameters. We introduce in the modeling two new structures, P2₁/b, and P6₃/m space group as reference centrosymmetric lattices, shown in Fig. 2. Here, a schematic representation of the OH electric dipole displacement from paraelectric to the ferroelectric configuration is illustrated. We caution the reader that only the polarization for λ = 1 has a real physical meaning, while the intermediate values are only computational states which allow monitoring the continuous evolution of the polarization as a function of λ, thus excluding the inclusion of quanta of polarization. We have also done additional calculations for the centrosymmetric configuration while keeping the O–H bond as what it is. We relaxed the structure and compared its energy with non-centrosymmetric phase, confirming that the energy difference between them is about 0.02 eV per unit cell and the ferroelectric phase is slightly higher in energy. The Fig. 3 shows a variation of the electric polarization of P6₃ (left) and P2₁ (right) structures as a function of the normalized amplitude of the polar distortion between centric (λ = 0) and polar (λ = 1) configurations.

Table 1 Experimental lattice parameters³⁸ of the hydroxyapatite crystal used in our present DFT calculation

Ca10(PO4)6(OH)2	Lattice parameters
Crystal symmetry	a/Å	b/Å	c/Å	α/°	β/°	γ/°
P63/m; polar P63	9.423	9.423	6.883	90	90	120
P21/b; polar P21	9.426	18.856	6.887	90	90	119.97

---

Fig. 1 Crystal structures of the low-temperature monoclinic P2₁ phase (a) and (b), and the high-temperature hexagonal P6₃ phase (c) and (d) hydroxyapatite Ca₁₀(PO₄)₆(OH)₂.

Fig. 2 Schematic representation of the OH⁻ electric dipole evolution from the paraelectric to ferroelectric configuration.

Fig. 3 Variation of the electric polarization of (a) P6₃/m to P6₃, and of (b) P2₁/b to P2₁ structures as a function of the amplitude of the polar distortion between centric (λ = 0) and polar (λ = 1) configurations.

It is useful to separate the polarization contributions coming from the two parts: the floppy component (OH⁻), (the H sublattice), and the inorganic framework [Ca₅(PO₄)₃]⁺, namely Ap. Then, we evaluate the FE polarization. Here HAp means the total value of ferroelectric polarization by considering the correlated distortions of the H and Ap together; H means the FE value of OH⁻ by considering the centrosymmetric inorganic framework; Ap means the FE value of [Ca₅(PO₄)₃]⁺ under centrosymmetric OH⁻ lattice; H + Ap means the simple algebraic summation of H and Ap. These results are displayed in the Fig. 4.

Fig. 4 The contribution of different parts of the ferroelectric polarization. HAp is the total value of ferroelectric polarization; H gives the FE value of OH⁻ within the centrosymmetric inorganic framework, and Ap gives the FE value of [Ca₅(PO₄)₃]⁺ within centrosymmetric OH⁻ lattice; H + Ap stands for the simple algebraic summation of H and Ap values.

Discussion

The hydroxyapatite structure consists of a pseudo-hexagonal network of apatite's tetrahedra with Ca²⁺ ions in the interstitial sites and columns of anions oriented along the c axis. Based on previous refinements of diffraction data, two phases have been suggested: a disordered hexagonal structure with P6₃/m symmetry and a monoclinic structure with P2₁/b symmetry. The monoclinic unit cell is derived from two hexagonal unit cells with the b axis parameter being twice that of the hexagonal value and involves an ordered arrangement of the anion columns.⁷ It has been proposed that the monoclinic phase is associated with stoichiometric apatites.

In order to study the ferroelectric order, we consider a paraelectric (centric) structure in terms of which one can present the polar crystal structure as due to symmetry-lowering structural distortion from the high symmetry centric group. It is useful to introduce a parameter λ which is the normalized amplitude of the distortion connecting the centric to the polar structure: λ = 1 represents the compound in the P6₃ or P2₁ space group, λ = 0 represents the compound in the P6₃/m or P2₁/b space group, respectively. We found that the systems have FE polarization of P = −7.05 μC cm⁻² and −9.87 μC cm⁻² along c axis for the hexagonal and monoclinic cases, respectively. Experimentally¹² the remnant polarization in hydroxyapatite film is ca. 25% of a PZT film, which possesses a polarization value of about 20–70 μC cm⁻². Our calculated results (7–10 μC cm⁻²) are very close with such experimental report.

We now discuss in more detail the separate role of the H and Ap on the polarization of the different HAp structures by fixing the centrosymmetric inorganic framework or by fixing the centrosymmetric OH⁻ lattice, respectively. For P6₃ structure, we obtain the total polarization P(HAp) = −7.05 μC cm⁻² from the fully relaxed structure. Assuming that the inorganic framework is centrosymmetric, we find the OH⁻ dipoles to be P(H) = −10.41 μC cm⁻². On the other hand, when fixing the OH⁻ coordinates to be centrosymmetric, the P(Ap) = 2.48 μC cm⁻² value is obtained, provided by the [Ca₅(PO₄)₃]⁺ sublattice. Finally, we have a simple algebraic summation of P(Ap) + P(H) = −7.93 μC cm⁻² which is close to the total polarization. For P2₁ structure, similarly we have total polarization P(HAp) = −9.87 μC cm⁻², while the separate P(H) = −8.74 μC cm⁻² and P(Ap) = −0.80 μC cm⁻², and P(Ap) + P(H) = −9.54 μC cm⁻² as well.

The origin of the different polarization values of these two structures can be understood. We note that the difference of the polarization value of the two polar structures is not too large, for the P6₃ and P2₁ structures. From calculated polarization of H and Ap components for the two structures, we found that the Ap contributes to the polarization difference more than the H part, even though Ap has only a small contribution to the total polarization.

The main mechanism for the transformation of P6₃ into P2₁ is that the [PO₄]³⁻ tetrahedron structure of Ap has a rotation of about 5° in monoclinic phase.⁷ The difference of H is mainly due to the variation of bond length between oxygen and hydrogen, the bond length of OH⁻ is 1.05359 Å in P6₃ and 1.01883 Å in P2₁, which results in different polarization values. For P6₃, the polarization of H should be larger than that of P2₁, but the polarization of Ap is opposite in sign, therefore the total polarization of P6₃ is smaller than that of P2₁.

Conclusions

We have analyzed the ferroelectric polarization of the low-temperature (monoclinic, P2₁) phase and the high-temperature (hexagonal, P6₃) phase of hydroxyapatites, on the basis of density functional theory calculations and the modern theory of ferroelectric polarization. We decomposed the total dipole moment contributions into hydroxyl groups and calcium phosphate apatite framework. The major polarization is from ordered hydroxyl anions for both phases. However, the inorganic [Ca₅(PO₄)₃]⁺ may have a sign change for polarization between these two polar structures. Detailed analysis of ferroelectric polarization and structural change of the HAp is provided. Our DFT results are useful for the better understanding of this important biomaterial in many actual and prospective applications.

Acknowledgements

This work was supported by the National Natural Science Foundation of China (51672171, 11274222), the National Key Basic Research Program of China (2015CB921600), and the Eastern Scholar Program from the Shanghai Municipal Education Commission. The Special Program for Applied Research on Super Computation of the NSFC-Guangdong Joint Fund (the second phase), the supercomputing services from AM-HPC, and Shanghai Supercomputer Center are also acknowledged.

References

1. P. Fratzl, H. S. Gupta, E. P. Paschalis and P. Roschger, J. Mater. Chem., 2004, 14, 2115–2123.
2. S. Busch, U. Schwarz and R. Kniep, Chem. Mater., 2001, 13, 3260–3271.
3. M. I. Kay, R. A. Young and A. S. Posner, Nature, 1964, 204, 1050–1052.
4. J. C. Elliott, P. E. Mackie and R. A. Young, Science, 1973, 180, 1055–1057.
5. S. Nakamura, H. Takeda and K. Yamashita, J. Appl. Phys., 2001, 89, 5386.
6. S. A. M. Tofail, D. Haverty, K. T. Stanton and J. B. McMonagle, Ferroelectrics, 2005, 319, 117–123.
7. D. Haverty, S. A. M. Tofail, K. T. Stanton and J. B. McMonagle, Phys. Rev. B: Condens. Matter Mater. Phys., 2005, 71, 094103.
8. O. Hochrein, R. Kniep and D. Zahn, Chem. Mater., 2005, 17, 1978–1981.
9. S. B. Lang, S. A. M. Tofail, A. A. Gandhi, M. Gregor, C. Wolf-Brandstetter, J. Kost, S. Bauer and M. Krause, Appl. Phys. Lett., 2011, 98, 123703.
10. V. S. Bystrov, E. Paramonova, D. Yu, A. Katashev, A. Karlov, N. Polyaka, A. V. Bystrova, A. Patmalnieks and A. L. Kholkin, J. Phys.: Condens. Matter, 2011, 23, 065302.
11. F. Chiatti, M. Corno and P. Ugliengo, J. Phys. Chem. C, 2012, 116, 6108–6114.
12. S. B. Lang, S. A. Tofail, A. L. Kholkin, M. Wojtas, M. Gregor, A. A. Gandhi, Y. Wang, S. Bauer, M. Krause and A. Plecenik, Sci. Rep., 2013, 3, 2215.
13. S. Ezhaveni, R. Yuvakkumar, M. Rajkumar, N. M. Sundaram and V. Rajendran, J. Nanosci. Nanotechnol., 2013, 13, 1631–1638.
14. A. A. Gandhi, M. Wojtas, S. B. Lang, A. L. Kholkin, S. A. M. Tofail and A. Bandyopadhyay, J. Am. Ceram. Soc., 2014, 97, 2867–2872.
15. M. Gregor, T. Plecenik, S. A. M. Tofail, M. Zahoran, M. Truchly, M. Vargova, F. Laffir, G. Plesch, P. Kus and A. Plecenik, Appl. Surf. Sci., 2015, 337, 249–253.
16. V. Uskokovic, RSC Adv., 2015, 5, 36614–36633.
17. V. S. Bystrov, J. Coutinho, A. V. Bystrova, Y. D. Dekhtyar, R. C. Pullar, A. Poronin, E. Palcevskis, A. Dindune, B. Alkan, C. Durucan and E. V. Paramonova, J. Phys. D: Appl. Phys., 2015, 48, 195302.
18. V. S. Bystrov, Ferroelectrics, 2015, 475, 148–153.
19. S. A. M. Tofail, A. A. Gandhi, M. Gregor and J. Bauer, Pure Appl. Chem., 2015, 87, 221–229.
20. M. Šupová, Ceram. Int., 2015, 41, 9203–9231.
21. N. Horiuchi, N. Wada, K. Nozaki, M. Nakamura, A. Nagai and K. Yamashita, J. Appl. Phys., 2016, 119, 084903.
22. F. Chiatti, M. Delle Piane, P. Ugliengo and M. Corno, Theor. Chem. Acc., 2016, 135, 15.
23. S. B. Lang, Phase Transitions, 2016, 89, 678–694.
24. J. Dhal, S. Bose and A. Bandyopadhyay, Mater. Sci. Eng., C, 2013, 33, 3061–3068.
25. J. D. Pasteris, B. Wopenka, J. J. Freeman, K. Rogers, E. Valsami-Jones, J. A. M. van der Houwen and M. J. Silva, Biomaterials, 2004, 25, 229–238.
26. E. Fukada and I. Yasuda, J. Phys. Soc. Jpn., 1957, 12, 1158–1162.
27. H. M. Frost, Angle Orthod., 1994, 64, 175–188.
28. S. B. Lang, Nature, 1966, 212, 704–705.
29. D. Kumar, J. P. Gittings, I. G. Turner, C. R. Bowen, L. A. Hidalgo-Bastida and S. H. Cartmell, Acta Biomater., 2010, 6, 1549–1554.
30. A. Cockbain, Mineral. Mag., 1968, 36, 654–660.
31. P. Hohenberg and W. Kohn, Phys. Rev. [Sect.] B, 1964, 136, 864–876.
32. G. Kresse and J. Hafner, Phys. Rev. B: Condens. Matter Mater. Phys., 1993, 47, 558–561.
33. G. Kresse and J. Furthmüller, Comput. Mater. Sci., 1996, 6, 15–50.
34. A. E. Mattsson, R. Armiento and T. R. Mattsson, Phys. Rev. Lett., 2008, 101, 239701.
35. R. Resta, Rev. Mod. Phys., 1994, 66, 899–915.
36. R. D. King-Smith and D. Vanderbilt, Phys. Rev. B: Condens. Matter Mater. Phys., 1993, 47, 1651–1654.
37. K. Momma and F. Izumi, J. Appl. Crystallogr., 2011, 44, 1272–1276.
38. M. Espanol, J. Portillo, J.-M. Manero and M.-P. Ginebra, CrystEngComm, 2010, 12, 3318.

---