Back to the Basics: Probing the Role of Surfaces in the Experimentally Observed Morphological Evolution of ZnO

Although the physics and chemistry of materials are driven by exposed surfaces in the morphology, they are fleeting, making them inherently challenging to study experimentally. The rational design of their morphology and delivery in a synthesis process remains complex because of the numerous kinetic parameters that involve the effective shocks of atoms or clusters, which end up leading to the formation of different morphologies. Herein, we combined functional density theory calculations of the surface energies of ZnO and the Wulff construction to develop a simple computational model capable of predicting its available morphologies in an attempt to guide the search for images obtained by field-emission scanning electron microscopy (FE-SEM). The figures in this morphology map agree with the experimental FE-SEM images. The mechanism of this computational model is as follows: when the model is used, a reaction pathway is designed to find a given morphology and the ideal step height in the whole morphology map in the practical experiment. This concept article provides a practical tool to understand, at the atomic level, the routes for the morphological evolution observed in experiments as well as their correlation with changes in the properties of materials based solely on theoretical calculations. The findings presented herein not only explain the occurrence of changes during the synthesis (with targeted reaction characteristics that underpin an essential structure–function relationship) but also offer deep insights into how to enhance the efficiency of other metal-oxide-based materials via matching.


Introduction
The determination of the surface-dependent properties of materials is essential for the structure-property relationship and the rational design for their high performance. Surface properties (i.e., surface energy, atomic structures, electronic structures, etc.) make a large difference to the stability and performance of materials. Thus, it is highly desirable that these properties be tuned in order to maximize performance, since the efficiency of these systems depends on the ability to control electronic levels on surfaces and at interfaces. Nevertheless, developing an adequate design of reaction conditions for the synthesis of a desirable morphology is a complex and difficult process. In principle, it is possible to predict equilibrium morphologies once the specific surface energies of exposed crystal surfaces become available. However, to determine the precise surface structures and energies and the morphological evolution of a given material, it is necessary to carry out multiple measurements, which is a time-consuming and resource-intensive task. Quantum mechanical simulations based on first-principles calculations have been commonly used to illuminate these phenomena at a fundamental level. Combining density functional theory (DFT) with ab initio atomistic thermodynamics can be a good strategy to overcome such experimental drawbacks and allow the investigation of exposed surfaces and morphological

Theoretical Methods
The theoretical morphologies of ZnO were studied using the Wulff construction obtained through the surface energy (E sur f ) values of the 1120 , 1010 , 1011 , 0001 , 1012 , 1122 , and 1121 planes reported by Na and Park [51]. The authors used the Vienna ab initio simulation package, adopting the LDA + U and PAW schemes. For the surface energy calculations, they employed periodically repeated slab geometry, which was separated by a vacuum layer of proper thickness [51]. To achieve the theoretical morphologies, the methodology proposed by our research group [33,52,53] was applied to obtain the available set of morphologies of ZnO. According to this methodology, the crystal morphology depends on the ratios between the surface energy and the crystal symmetry and structure [33,54].
The polyhedron energy (E pol ) and the percentage of contribution of each surface in the ZnO morphology are calculated. E pol is obtained by the following equation: where C i is the contribution of the surface area to the total surface area of the polyhedron (C i = A i /A pol ) and E i sur f is the surface energy value of the corresponding surface i [37]. On the other hand, the well-known Wulff construction is a convenient method to evaluate the formation of a macroscopic surface B of orientation (h 2 k 2 l 2 ) on a surface A of orientation (h 1 k 1 l 1 ). The relative energy, ∆E, can be calculated by the following expression: where E A sur f is the surface energy (per unit area) of surface A (of orientation (h 1 k 1 l 1 )), E B sur f is the surface energy (per unit area) of surface B (of orientation (h 2 k 2 l 2 )), θ is the angle between surfaces A and B, and the cosθ factor corresponds to the change in surface area when facets are formed [55]. According to this expression, if ∆E is negative, surface B can grow stably on surface A, i.e., the growth process takes place along the surface with lowest surface energy. This is the classical growth mechanism of Ostwald ripening, which describes the growth of smaller crystals into larger ones through diffusion in order to reduce the total surface energy [56][57][58].

A computational Road to Morphology
This paper provides an alternative approach for the efficient generation of the available morphologies (morphology map) of a given material, in addition to a quantitative structurereactivity relationship based on quantum chemistry and the Wulff construction. The first step of this investigation involves the study of the bulk (unit cell).
The crystallographic unit cell of ZnO is shown in Figure 1. The ZnO structure is fully determined by the lattice parameters a = b and c, belongs to the space group P6 3 mc, and is formed by a two-unit formula per cell (Z = 2). In the ZnO structure, the Zn cations have a coordination number of four, which means that they are surrounded by a tetrahedron of O 2− anions. Therefore, the ZnO structure has [ZnO 4 ] clusters as building blocks, but with some local disorder, as illustrated in Figure 1. A note on terminology: to avoid confusion, the term cluster will be used exclusively to denote the local coordination of the Zn cation corresponding to the number of neighboring oxygen anions both in the bulk and on the exposed surfaces in the morphology.
The next step consists of the investigation of surfaces that can be cut through the bulk. Previous studies have applied several different functionals to calculate the surface energy values of ZnO, as shown in Table 1. The next step consists of the investigation of surfaces that can be cut through the bulk. Previous studies have applied several different functionals to calculate the surface energy values of ZnO, as shown in Table 1. By applying our methodology [33,52,53] and combining the surface energy values reported by Na and Park [51] and the Wulff construction, we were able to obtain a map of the available morphologies of ZnO, as illustrated in Figure 2. In the center of this figure, it is possible to see the starting morphology using the surface energy values calculated by Na and Park [51]. From this morphology, we obtained the available morphologies by decreasing the surface energy values of one (or two) surfaces using different synthesis methods.
This map becomes a powerful tool for experimentalists during the morphological characterization of materials, since it allows matching the experimental morphologies to the theoretical ones. At this point, it is possible to note that the variation in relative surface energy values is more important than their calculation, which avoids the technical drawback resulting from the fact that different functionals provide different calculated surface energy values, as observed in Table 1.

Where Will This Road Take Us?
This work provides insights to better understand the underlying mechanisms of crystal growth in the synthesis process. To this end, we delineated the reaction pathways that connect the different morphologies.  By applying our methodology [33,52,53] and combining the surface energy values reported by Na and Park [51] and the Wulff construction, we were able to obtain a map of the available morphologies of ZnO, as illustrated in Figure 2. In the center of this figure, it is possible to see the starting morphology using the surface energy values calculated by Na and Park [51]. From this morphology, we obtained the available morphologies by decreasing the surface energy values of one (or two) surfaces using different synthesis methods.
This map becomes a powerful tool for experimentalists during the morphological characterization of materials, since it allows matching the experimental morphologies to the theoretical ones. At this point, it is possible to note that the variation in relative surface energy values is more important than their calculation, which avoids the technical drawback resulting from the fact that different functionals provide different calculated surface energy values, as observed in Table 1.

Where Will This Road Take Us?
This work provides insights to better understand the underlying mechanisms of crystal growth in the synthesis process. To this end, we delineated the reaction pathways that connect the different morphologies.
Some of the morphologies displayed in Figure 2 were reported in the literature. From these morphologies, we were able to calculate the E pol values and build a reaction pathway for use in the synthesis process, as shown in Figure 3. It is important to note that the selected images of the morphology in Figure 3 correspond to static or steady-state values. Some of the morphologies displayed in Figure 2 were reported in the literature. From these morphologies, we were able to calculate the values and build a reaction pathway for use in the synthesis process, as shown in Figure 3. It is important to note that the selected images of the morphology in Figure 3 correspond to static or steady-state values.
The starting morphology was obtained by Debroye et al. [9] and synthesized as described by Kiomarsipour and Shoja Razavi [66] by a simple hydrothermal process at low temperature without any additional surfactant, organic solvent, or catalytic agent. The elongated hexagonal morphology (a) was obtained by Amin et al. by the hydrothermal method using different experimental parameters such as pH, precursor concentration, growth time, and temperature [67]. This morphology was also reported by some of us through the doping of ZnO with Ni and Fe to enhance its photocatalytic activity [68]. These results are clear-cut examples of how the reaction conditions and the synthesis methods can modulate the final morphology.
The elongated octahedral morphology (b) was also obtained by Wu et al. by the hydrothermal method, but using water and methanol during the synthesis [69]. Zhang et al. observed a lance-shaped morphology (c) in a flower-like architecture by using different conditions in a controlled hydrothermal process (water/ethanol as a solvent and different ratio of precursors) [70].
Through the solvothermal method, Liu et al. modulated the reaction conditions, i.e., reaction time and additive (tetramethylammonium hydroxide, TMAH) concentration, to obtain a wide range of morphologies of the as-synthetized ZnO samples [71]. By using the values of , we were able to calculate the reaction pathway that connected the morphologies obtained by Liu et al., as illustrated in Figure 4. The starting morphology was obtained by Debroye et al. [9] and synthesized as described by Kiomarsipour and Shoja Razavi [66] by a simple hydrothermal process at low temperature without any additional surfactant, organic solvent, or catalytic agent. The elongated hexagonal morphology (a) was obtained by Amin et al. by the hydrothermal method using different experimental parameters such as pH, precursor concentration, growth time, and temperature [67]. This morphology was also reported by some of us through the doping of ZnO with Ni and Fe to enhance its photocatalytic activity [68]. These results are clear-cut examples of how the reaction conditions and the synthesis methods can modulate the final morphology.
The elongated octahedral morphology (b) was also obtained by Wu et al. by the hydrothermal method, but using water and methanol during the synthesis [69]. Zhang et al. observed a lance-shaped morphology (c) in a flower-like architecture by using different conditions in a controlled hydrothermal process (water/ethanol as a solvent and different ratio of precursors) [70].
Through the solvothermal method, Liu et al. modulated the reaction conditions, i.e., reaction time and additive (tetramethylammonium hydroxide, TMAH) concentration, to obtain a wide range of morphologies of the as-synthetized ZnO samples [71]. By using the values of E polyhedron , we were able to calculate the reaction pathway that connected the morphologies obtained by Liu et al., as illustrated in Figure 4. ls 2023, 13, x FOR PEER REVIEW 6 of 15 Figure 3. Calculated values and the reaction pathway used in the synthesis process to obtain the most common experimental morphologies (inset): starting [9], (a) elongated hexagonal [67], (b) elongated octahedral [69], and (c) lance-shaped [70]. The percentages of each surface area are also provided for comparison purposes in a pie chart. Reprinted with permission from [9], under the terms of the Creative Commons CC-BY license. Reprinted with permission from [64], under the terms of the Creative Commons CC license. Reprinted with permission from [66]; Copyright 2020, Elsevier. Reprinted with permission from [67]; Copyright 2011, John Wiley and Sons.
As it can be seen, the morphology after the "initial stage" (b) has a higher value of than the starting morphology (a), i.e., 0.87 J/m 2 vs. 0.84 J/m 2 . From this point, two alternative routes can be opened as functions of synthesis process and time, resulting in two different situations to be analyzed: changes in time and TMAH concentration during synthesis. In the first case, the surface composition changes as a function of time. The contribution of the (101 0) surface decreases from 58% to 33%, whereas the contribution of (101 1) increases from 14% to 48%. However, the value does not vary, remaining 0.87 J/m 2 (see morphology (c)). On the other hand, when the amount of TMAH is increased, the value of morphology (d) is found to be lower (0.81 J/m 2 ). The contributions of (101 0) and (0001 ) decrease to 15% and 14%, respectively, while the contribution of (101 1) increases from 14% to 72%. This means that the stabilization of the (101 1) surface renders a more stable crystal morphology. These results can rationalize those reported by Liu et al. and explain why the (101 1) surface is more favorable for photocatalysis than the (101 0) plane [71].

Which Way Does the Morphology Go?
As observed, the hydrothermal method is one of the most frequently used to synthetize ZnO, with the reaction conditions being responsible for the changes in surface stabil- . Calculated E pol values and the reaction pathway used in the synthesis process to obtain the most common experimental morphologies (inset): starting [9], (a) elongated hexagonal [67], (b) elongated octahedral [69], and (c) lance-shaped [70]. The percentages of each surface area are also provided for comparison purposes in a pie chart. Reprinted with permission from [9], under the terms of the Creative Commons CC-BY license. Reprinted with permission from [64], under the terms of the Creative Commons CC license. Reprinted with permission from [66]; Copyright 2020, Elsevier. Reprinted with permission from [67]; Copyright 2011, John Wiley and Sons.
As it can be seen, the morphology after the "initial stage" (b) has a higher value of E pol than the starting morphology (a), i.e., 0.87 J/m 2 vs. 0.84 J/m 2 . From this point, two alternative routes can be opened as functions of synthesis process and time, resulting in two different situations to be analyzed: changes in time and TMAH concentration during synthesis. In the first case, the surface composition changes as a function of time. The contribution of the 1010 surface decreases from 58% to 33%, whereas the contribution of 1011 increases from 14% to 48%. However, the E polyhedron value does not vary, remaining 0.87 J/m 2 (see morphology (c)). On the other hand, when the amount of TMAH is increased, the E polyhedron value of morphology (d) is found to be lower (0.81 J/m 2 ). The contributions of 1010 and 0001 decrease to 15% and 14%, respectively, while the contribution of 1011 increases from 14% to 72%. This means that the stabilization of the 1011 surface renders a more stable crystal morphology. These results can rationalize those reported by Liu et al. and explain why the 1011 surface is more favorable for photocatalysis than the 1010 plane [71].
ity. The favorable growth directions of ZnO can be disclosed by using the Wulff construction to evaluate the formation of a macroscopic surface B of orientation (h2k2l2) on a surface A of orientation (h1k1l1).  Table 2.

Which Way Does the Morphology Go?
As observed, the hydrothermal method is one of the most frequently used to synthetize ZnO, with the reaction conditions being responsible for the changes in surface stability. The favorable growth directions of ZnO can be disclosed by using the Wulff construction to evaluate the formation of a macroscopic surface B of orientation (h 2 k 2 l 2 ) on a surface A of orientation (h 1 k 1 l 1 ).
From the surface energy values used in the construction of the morphology map, it was possible to calculate the ∆E values for all possible surfaces B on all surfaces A. These values are presented in Table 2.
By using the values of ∆E, it is possible to predict the preferential crystal growth direction of ZnO. According to the literature [10,[72][73][74], the growth of ZnO crystals along the [0001] direction is the most reported. Nonetheless, Cho et al. observed that when the triethyl citrate is used as a surfactant, a lateral growth of each spine along the six symmetric directions can be noted [75]. Therefore, by combining the ∆E values listed in Table 2, we could thermodynamically estimate the most favorable side surface growth along these directions, as shown in Table 3. As it can be seen in Table 3, the combination of 1120 / 1010 , 1010 / 1011 , and 1120 / 1011 surfaces results in a favorable crystal growth process along the [0001] direction. For the [1010] direction, the combination of 1120 / 1011 and 0001 / 1012 also indicates a favorable crystal growth process. To further explore this behavior, it is necessary to calculate the ∆E values corresponding to the combination of surface B of orientation (h 2 k 2 l 2 ) on surface A of orientation (h 1 k 1 l 1 ) using the surface energy values of the morphologies depicted in Figures 3 and 4. These values are presented in Tables 4 and 5, respectively. A detailed analysis of the results in Table 4 shows that in the starting morphology,  Table 2). This explains the lance-shaped morphology depicted in (c), where the contribution of the 1010 surface corresponds to 82% against 18% for the 1011 plane. However, for morphology (b), these values are 35% and 65%, respectively. Table 4. Calculated values of relative energy (∆E) of surface B on surface A (J/m 2 ) and angle (θ, degree) between surfaces A and B for the crystal morphologies reported in Figure 3.  Table 5. Calculated values of relative energy (∆E) of surface B on surface A (J/m 2 ) and angle (θ, degree) between planes A and B for all crystal shapes reported in Figure 4. According to the results presented in Table 5 for the morphologies reported by Liu et al. [71] in Figure 4, the reaction conditions change the ZnO morphologies, resulting in a stabilization of the 1011 surface in morphology (b). Initially, the formation energy in morphology (a) stabilizes the growth of 1010 on 0001 (∆E = −1.22 J/m 2 ). After the initial stage, as the time and additive concentration increase, the growth of 1011 is also stabilized, as confirmed by the negative ∆E values for the combination of 0001 \ 1011 surfaces. Another important fact is that the percentage contribution of the 1010 plane decreases, whereas the contribution of 1011 increases, as seen in Figure 4.

A\B
Several works have used experimental characterization techniques such as XRD and TEM as valuable tools to investigate the growth direction of crystals [76][77][78][79]. For instance, Chang and Waclawik controlled morphological transformations by varying the reaction temperature and molar ratio (benzylamine/Zn 2+ concentration from 1 to 10) [20] and obtained ZnO with a nano-bullet-like morphology exhibiting exposed 1011 and 1010 surfaces (Figure 5a). An analysis of the results previously reported in Table 2 shows that the combination of both surfaces provokes a thermodynamically favorable crystal growth process with ∆E < 0. As shown in Table 2, the growth in the [0001] direction is favored when these surfaces are combined, as observed in the elongated nano-bullet-like morphology (see (a) in Figure 5). By decreasing the synthesis temperature from 210 to 170 • C, a hexagonal cone-like morphology with an exposed 1010 surface can be obtained (see (b) in Figure 5), while an increase in the benzylamine/Zn 2+ molar ratio results in a 2D plate-like morphology with an exposed (0001) surface (see (c) in Figure 5). Ahmed et al. prepared ZnO nanocrystals by the hydrothermal method. Single ZnO nanorods were transformed into sharp sword-like tips by increasing the reaction time to 30 min. After 60 min of reaction, a single ZnO semi-hollow nanorod (pyramid-like morphology) was obtained [80]. This can be supported by the HRTEM images, which depict a lattice spacing corresponding to the distance between the (002) planes in the obtained ZnO structures growing along the [0001] direction. According to the results in Table 2, the growth process along the [0001] direction is favored when different surfaces are combined. In the synthesis process, the final morphology is dependent on the growth velocity, which in turn is affected by the nonuniformity and variability of the precursor solution throughout the reaction time.
Interesting reaction pathways were proposed by Liu and Liu by employing the pulsed-laser deposition technique to obtain a given morphology of ZnO [81]. The authors demonstrated that laser-induced crystal growth is a practical tool to tune the morphology of nanomaterials in a precise and effective manner. The growth directions of the ZnO crystallization in a hydrothermal reaction were selected by adjusting the laser irradiation conditions (power and time) to control the appearance of the final morphology (hexagonal versus pyramid-like, corresponding to the kinetic and thermodynamic reaction pathways, respectively). The above results can be rationalized by tuning the surface energy values of (0001 ), (101 0), and ( 101 1 [20] through the combination of (0001), 1010 , and 1011 surfaces. Ahmed et al. prepared ZnO nanocrystals by the hydrothermal method. Single ZnO nanorods were transformed into sharp sword-like tips by increasing the reaction time to 30 min. After 60 min of reaction, a single ZnO semi-hollow nanorod (pyramid-like morphology) was obtained [80]. This can be supported by the HRTEM images, which depict a lattice spacing corresponding to the distance between the (002) planes in the obtained ZnO structures growing along the [0001] direction. According to the results in Table 2, the growth process along the [0001] direction is favored when different surfaces are combined. In the synthesis process, the final morphology is dependent on the growth velocity, which in turn is affected by the nonuniformity and variability of the precursor solution throughout the reaction time.
Interesting reaction pathways were proposed by Liu and Liu by employing the pulsedlaser deposition technique to obtain a given morphology of ZnO [81]. The authors demonstrated that laser-induced crystal growth is a practical tool to tune the morphology of nanomaterials in a precise and effective manner. The growth directions of the ZnO crystallization in a hydrothermal reaction were selected by adjusting the laser irradiation conditions (power and time) to control the appearance of the final morphology (hexagonal versus pyramid-like, corresponding to the kinetic and thermodynamic reaction pathways, respectively). The above results can be rationalized by tuning the surface energy values of 0001 , 1010 , and 1011 . From the kinetically controlled growth product with surface energy values of 1.20, 0.84, and 1.73 J/m 2 , respectively, a need arises: to overcome an energy barrier of 0.11 J/m 2 in order to obtain an intermediate morphology. This pathway is achieved by increasing or decreasing the surface energy value of the 0001 and 1011 surfaces to 1.55 and 1.05 J/m 2 , respectively. From this intermediate morphology, it is possible to obtain a thermodynamically controlled growth product by increasing or decreasing the surface energy values of 0001 and 1011 to 2.20 and 0.87 J/m 2 , respectively. This results in energy barriers from the intermediate to the kinetically and thermodynamically controlled morphologies reported by Liu and Liu of 0.11 and 0.22 J/m 2 , respectively. A schematic illustration of the reaction pathways is presented in Figure 6. tively. This results in energy barriers from the intermediate to the kinetically and thermodynamically controlled morphologies reported by Liu and Liu of 0.11 and 0.22 J/m 2 , respectively. A schematic illustration of the reaction pathways is presented in Figure 6.  [81].

Conclusions
In this study, we selected ZnO to evaluate the importance of our computational method created to reciprocate findings on morphologies and crystal growth processes from experiments. Based on the surface energy values and the Wulff construction, this strategy was found to be very useful for unraveling the morphologies that are challenging to characterize experimentally. This work provides a theoretical framework that requires as input data the surface energy values to obtain the available morphologies of a given semiconductor, its polyhedron energy, and the reaction pathways involved in a synthetic road to achieve a certain morphology. Important information can also be taken from the results presented herein for further research on the effects of morphological control on the synthesis of semiconductors. The figures of merit in this morphology map agree with the experimental images obtained by field-emission scanning electron microscopy. The high degree of coincidence between the theory and the experiments makes us believe that the model might have a more general scope of application, which, to the best of our knowledge, has not been discussed in previous literature. Additional studies are in progress to standardize this computational procedure by proving its efficiency in replicating experimental data, as well as the usefulness of techniques employed to predict structural, physical, chemical, and dynamic properties of materials. We also expect that these new insights will guide researchers aiming to apply the potential of computational methods to illustrate minute details of various types of materials.

Conclusions
In this study, we selected ZnO to evaluate the importance of our computational method created to reciprocate findings on morphologies and crystal growth processes from experiments. Based on the surface energy values and the Wulff construction, this strategy was found to be very useful for unraveling the morphologies that are challenging to characterize experimentally. This work provides a theoretical framework that requires as input data the surface energy values to obtain the available morphologies of a given semiconductor, its polyhedron energy, and the reaction pathways involved in a synthetic road to achieve a certain morphology. Important information can also be taken from the results presented herein for further research on the effects of morphological control on the synthesis of semiconductors. The figures of merit in this morphology map agree with the experimental images obtained by field-emission scanning electron microscopy. The high degree of coincidence between the theory and the experiments makes us believe that the model might have a more general scope of application, which, to the best of our knowledge, has not been discussed in previous literature. Additional studies are in progress to standardize this computational procedure by proving its efficiency in replicating experimental data, as well as the usefulness of techniques employed to predict structural, physical, chemical, and dynamic properties of materials. We also expect that these new insights will guide researchers aiming to apply the potential of computational methods to illustrate minute details of various types of materials.