*2.2. Planar Density of Unit Cell and Super Cells (2* × *2* × *2) of Hydroxyapatite*

According to the resulting list of planes by X-ray diffraction (Figure 1, Table S1), the planar density values of each diffracted plane in unit cell and super cells (2 × 2 × 2) of hydroxyapatite are calculated in Figure 3, Figures S6 and S7. To calculate the planar density values, diffracted planes were selected from a low angle to a high angle in tandem. It is worth mentioning that the matrix of all super cell lattices was considered to be 2 × 2 × 2. According to the center of atoms, the planar density is calculated by the area of the atoms in the plane divided by the total area of that plane [19].

#### *2.3. Elastic Stiffness Constant and Elastic Compliance of Hydroxyapatite*

Hooke's law is shown in Equation (2); the stress corresponds to the strain for small displacements. It is the basic form, that this symmetry can be converted to the six items of stress (σ) and strain (ε) [20].

$$
\begin{pmatrix}
\sigma\_{\rm{XX}} \\
\sigma\_{\rm{YY}} \\
\sigma\_{\rm{zz}} \\
\sigma\_{\rm{yz}} \\
\sigma\_{\rm{xz}} \\
\sigma\_{\rm{xy}}
\end{pmatrix} = \begin{pmatrix}
\mathsf{c}\_{11} & \mathsf{c}\_{12} & \mathsf{c}\_{13} & \mathsf{c}\_{14} & \mathsf{c}\_{15} & \mathsf{c}\_{16} \\
\mathsf{c}\_{21} & \mathsf{c}\_{22} & \mathsf{c}\_{23} & \mathsf{c}\_{24} & \mathsf{c}\_{25} & \mathsf{c}\_{26} \\
\mathsf{c}\_{31} & \mathsf{c}\_{32} & \mathsf{c}\_{33} & \mathsf{c}\_{34} & \mathsf{c}\_{35} & \mathsf{c}\_{36} \\
\mathsf{c}\_{41} & \mathsf{c}\_{42} & \mathsf{c}\_{43} & \mathsf{c}\_{44} & \mathsf{c}\_{45} & \mathsf{c}\_{46} \\
\mathsf{c}\_{51} & \mathsf{c}\_{52} & \mathsf{c}\_{53} & \mathsf{c}\_{54} & \mathsf{c}\_{55} & \mathsf{c}\_{56} \\
\mathsf{c}\_{61} & \mathsf{c}\_{62} & \mathsf{c}\_{63} & \mathsf{c}\_{64} & \mathsf{c}\_{65} & \mathsf{c}\_{66} \\
\end{pmatrix} = \begin{pmatrix}
\varepsilon\_{\rm{XX}} \\
\varepsilon\_{\rm{yy}} \\
\varepsilon\_{\rm{xx}} \\
\varepsilon\_{\rm{zz}} \\
\varepsilon\_{\rm{xy}} \\
\varepsilon\_{\rm{xy}} \\
\varepsilon\_{\rm{xy}}
\end{pmatrix} \tag{2}
$$

Additionally, Hooke's law can be written (Equation (3)):

σxx = C11εxx+C12εyy+C13εzz+C14εyz+C15εzx+C16εxy σyy = C21εxx+C22εyy+C23εzz+C24εyz+C25εzx+C26εxy σzz = C31εxx+C32εyy+C33εzz+C34εyz+C35εzx+C36εxy σyz = C41εxx+C42εyy+C43εzz+C44εyz+C45εzx+C46εxy σzx = C51εxx+C52εyy+C53εzz+C54εyz+C55εzx+C56εxy σxy = C61εxx+C62εyy+C63εzz+C64εyz+C65εzx+C66εxy. (3)

The elastic stiffness determines the response of crystal to an externally applied stress or strain and provides information about bonding characteristics and mechanical and structural stability [21]. The hydroxyapatite system has five elastic constants (Equation (4)). Therefore, the values of five independent Cij, can be named C11, C12, C13, C33, C44.

**Figure 3.** *Cont.*

**Figure 3.** Array and position of the involved atoms such as (**a**) (020) unit cell, (**b**) (020) super cell, (**c**) (111) unit cell, (**d**) (111) super cell, (**e**) (002) unit cell, (**f**) (002) super cell, (**g**) (102) unit cell and (**h**) (102) super cell (the calculations are in the Supplementary materials).

The crystallographic nature of the hexagonal structure is shown in Figure S3. Furthermore, for conventional hexagonal systems, such as hydroxyapatite, the relationship between Cij and Sij is introduced in Equations (5)–(9) [22,23].

$$\mathcal{S}\_{11} = \frac{1}{2} \left( \frac{\mathcal{C}\_{33}}{\mathcal{C}\_{33}(\mathcal{C}\_{11} + \mathcal{C}\_{12}) - 2(\mathcal{C}\_{13})^2} + \frac{1}{\mathcal{C}\_{11} - \mathcal{C}\_{12}} \right) \tag{5}$$

$$\mathbf{S}\_{12} = \frac{1}{2} \left( \frac{\mathbf{C}\_{33}}{\mathbf{C}\_{33}(\mathbf{C}\_{11} + \mathbf{C}\_{12}) - 2(\mathbf{C}\_{13})^2} - \frac{1}{\mathbf{C}\_{11} - \mathbf{C}\_{12}} \right) \tag{6}$$

$$\mathbf{S}\_{33} = \frac{\mathbf{C}\_{11} + \mathbf{C}\_{12}}{\mathbf{C}\_{33}(\mathbf{C}\_{11} + \mathbf{C}\_{12}) - 2(\mathbf{C}\_{13})^2} \tag{7}$$

$$\mathbf{S}\_{13} = -\frac{\mathbf{C}\_{13}}{\mathbf{C}\_{33}(\mathbf{C}\_{11} + \mathbf{C}\_{12}) - 2(\mathbf{C}\_{13})^2} \tag{8}$$

$$\mathbf{S\_{44}} = \frac{1}{\mathbf{C\_{44}}} \tag{9}$$

According to the Equations (5)–(9), to obtain S values C values are needed. We have used two approaches to obtain C values: First, theoretical calculations were performed via the CASTEP model of materials studio software version 6.0 in the GGA level of theory with a PBE basis set (Figure S4). The second approach is based on ultrasonic measurements (Table S2). The complete set of five elastic stiffness constant values (C11, C12, C13, C33 and C44) of the samples was found from ultrasonic measurements of the phase velocity anisotropy. The stiffness constant values were recorded by utilizing Equations (10)–(15). In these equations, ρ and V are density of sample and velocity, respectively [24–28].

$$\mathbf{C}\_{11} = \rho \mathbf{V}\_{1/\iota \prime}^{2} \mathbf{C}\_{22} = \rho \mathbf{V}\_{2/2}^{2} \tag{10}$$

$$\mathbf{C}\_{66} = \rho \mathbf{V}\_{1/2}^2 = \rho \mathbf{V}\_{2/1'}^2, \mathbf{C}\_{55} = \rho \mathbf{V}\_{1/3}^2 = \rho \mathbf{V}\_{3/1}^2 \tag{11}$$

$$\mathbf{C}\_{12} = \sqrt{\left(\mathbf{C}\_{11} + \mathbf{C}\_{66} - 2\rho \mathbf{V}\_{12/12}^2\right) \left(\mathbf{C}\_{22} + \mathbf{C}\_{66} - 2\rho \mathbf{V}\_{12/12}^2\right)} - \mathbf{C}\_{66} \tag{12}$$

$$\mathbf{C}\_{\mathbf{44}} = \varrho \mathbf{V}\_{\mathbf{2/3}}^2 = \varrho \mathbf{V}\_{\mathbf{3/2}}^2 \tag{13}$$

$$\mathbf{C}\_{13} = \sqrt{\left(\mathbf{C}\_{11} + \mathbf{C}\_{55} - 2\rho \mathbf{V}\_{13/13}^2\right) \left(\mathbf{C}\_{33} + \mathbf{C}\_{55} - 2\rho \mathbf{V}\_{13/13}^2\right)} - \mathbf{C}\_{55} \tag{14}$$

$$\mathbf{C}\_{33} = \rho \mathbf{V}\_{3/3}^2 \tag{15}$$

For measuring the velocities, the standard ultrasonic pulse-echo ASTM E797/E797-M-15 was accomplished according to Reference [26]. In this study, the immersion procedure via water between probe and sample was utilized and the effect of pressure derived from a hand placed into the probe was decreased. In addition, the longitudinal frequency probe was 5.4 MHz, and to decrease the extension and depreciation of waves, especially the more energetic waves, a lens with a specific curvature was utilized according to the standard of Reference [28]. To measure the value of C66, a high amplitude of curve was carried out; therefore, changing the rotation angles was useful. To create the transverse waves, a probe of 2.3 MHz and a high viscos interface material (honey) were used and finally the pressure on the sample was adjusted by obtaining the better and smoother curve [29]. Taking into account each point of the samples, a three-dimensional axis, such as X1, X2 and X3, can be performed. In this case, according to Figure S5, X1 is the radial coordinate, X2 is the circumferential coordinate and X3 is the axial coordinate. Vi/j is the denoted velocity of an ultrasound wave that can be propagated in the Xi direction with particle displacements in the Xj direction simultaneously. Vi/j, with the same i and j, is longitudinal and with i = j being transverse waves. For measuring quasi-longitudinal or quasi-transverse velocity (Vij/ij), specimens should be cut (bezel) on the edges toward surfaces of perpendicular X directions (Figure S5). The obtained values of velocities are shown in Table S2. On the other hand, C11 is in good agreement with the longitudinal distortion and longitudinal compression/tension; thus, C11 can be introduced as the hardness. Moreover, the transverse distortion depends on the C12, and C12 comes from the transverse expansion related to the Poisson's ratio. Additionally, C44 corresponds to the shear modulus, and C44 is in the settlement with C11 and C12 [22,30]. Accordingly, the shear modulus is proportional to the Burgers vector and the Young's modulus; in addition, dislocation density is in agreement with the Young's modulus [31,32].

In the ultrasonic method, longitudinal and transverse waves were utilized for measuring the Young's modulus value [33,34]. According to this method (Equation (16)), based on the velocity of ultrasound waves and density of sample, the Young's modulus value was determined.

$$\mathcal{E} = \frac{\rho c\_1^2 \left[3\left(\frac{c\_1}{c\_t}\right)^2 - 4\right]}{\left(\frac{c\_1}{c\_t}\right)^2 - 1} \tag{16}$$

In Equation (16), ρ, cl and ct are density, velocity of longitudinal and transverse ultrasound waves tandemly. Furthermore, according to Equation (17), the velocity of longitudinal and transverse waves can be registered by determining the length of specimen and the differences between two echoes (t = t2 − t1) in the signals [35].

$$\mathbf{c} = \frac{2\mathbf{L}}{\mathbf{t}} \tag{17}$$

Here, L is the length of the sample and t is the difference between two echoes, and the density of the sample can be detected by measuring the mass and volume of the sample [36]. Additionally, with the substitution of Equations (16) and (17), the main equation for the calculation of Young's modulus is Equation (18).

$$\mathbf{E} = \frac{4\rho \left(\frac{1}{\mathbf{f}\_s}\right)^2 \left(3\mathbf{t}\_s^2 - 4\mathbf{t}\_1^2\right)}{\mathbf{t}\_s^2 - \mathbf{t}\_1^2} \tag{18}$$

ts and tl are differences between two echoes in longitudinal and transverse waves separately [37]. The results of the theoretical calculation, and the experimental measurements of elastic stiffness constant values of hydroxyapatite (from the literature and this study), are shown in Table 2. In addition, taking into account Equations (10)–(15), the elastic compliance values of this study were calculated and are recorded in Table 3. As a result, the theoretical values were in good agreement with the theoretical values extracted by References [38,39]. The replicated values (five times) for measuring the time of transverse and longitudinal waves of the samples are presented in Table 4.

**Table 2.** Theoretical and experimental values of the Elastic constant of hydroxyapatite.


**Table 3.** Theoretical and experimental values of the Elastic compliance of hydroxyapatite.


**Table 4.** Resulted values of time for transverse and longitudinal waves of samples.


### *2.4. The Young's Modulus versus Planar Density of Unit Cell and Super Cells (2* × *2* × *2) of Hydroxyapatite*

X-ray diffraction has provided data on diffracted planes and the location of atoms in each plane. In the previous section, the planar density of each diffracted plane was calculated and it could play an important role in the mechanical properties of each plane. The Young's modulus of each plane (Ehkl) of a hydroxyapatite lattice can be calculated as Equation (19) [5,43]. In this equation, h, k and l are the plane indices, a and c are the lattice parameters, C and S are the elastic stiffness constant and elastic compliance, respectively. The values of Young's modulus of 32 diffracted planes of hydroxyapatite in unit cells and super cells (2 × 2 × 2),E(h1k1l1), E(h2k2l2), E(h3k3l3) ............ ., E(h32k32l32), related to the literature and the present study, are reported in Table S3.

$$\mathbf{E}\_{\mathbf{h}\mathbf{k}l} = \frac{\left[\mathbf{h}^2 + \frac{(\mathbf{h} + 2\mathbf{k})^2}{3} + \left(\frac{\mathbf{a}\mathbf{l}}{\mathbf{c}}\right)^2\right]^2}{\mathbf{S}\_{11}\left(\mathbf{h}^2 + \frac{(\mathbf{h} + 2\mathbf{k})^2}{3}\right)^2 + \mathbf{S}\_{33}\left(\frac{\mathbf{a}\mathbf{l}}{\mathbf{c}}\right)^4 + (2\mathbf{S}\_{13} + \mathbf{S}\_{44})\left(\mathbf{h}^2 + \frac{(\mathbf{h} + 2\mathbf{k})^2}{3}\right)\left(\frac{\mathbf{a}\mathbf{l}}{\mathbf{c}}\right)^2} \tag{19}$$

By using the least squares method between the Young's modulus and the planar density of diffracted planes (based on our proposed method), the Young's modulus value of hydroxyapatite was determined with high precision. Consequently, the calculation of C and S parameters for crystalline materials is essential for the application of this method. To show the feasibility and accuracy of our proposed method for determining the Young's modulus, the values of Young's modulus of each plane versus the planar density of the unit cell and supercells (2 × 2 × 2) are depicted in Figure 4. The Young's modulus values (intercept) of the unit cell and super cells (2 × 2 × 2) of hydroxyapatite, obtained from the least squares method, are tabulated in Table 5.

**Figure 4.** Young's modulus of each plane (**a**) unit-cell (**b**) super cells (2 × 2 × 2) of hydroxyapatite extracted by XRD patterns and planar density.


**Table 5.** Young's modulus values of unit cell and super cell lattices of hydroxyapatite.

According to the uncertainty measurement (the measurements were replicated five times (Table 4)) and Equation (18), the Young's modulus value gained 113.08 ± 0.14 GPa by ultrasonic measurement. This value is in good agreement with the reported values of our study in Table 5. In this study, the difference between theory and experiment values for both unit cell and super cells (2 × 2 × 2) are identical. This difference for unit cell and super cells (2 × 2 × 2) are 10.47 GPa and 10.57 GPa, respectively. This means that the theoretical calculation is valid and, by reducing it by about ~10 GPa, the experimental values can be obtained.

### **3. Result and Discussion**
