Porous Si Partially Filled with Water Molecules—Crystal Structure, Energy Bands and Optical Properties from First Principles

The paper reports the results on first-principles investigation of energy band spectrum and optical properties of bulk and nanoporous silicon. We present the evolution of energy band-gap, refractive indices and extinction coefficients going from the bulk Si of cubic symmetry to porous Si with periodically ordered square-shaped pores of 7.34, 11.26 and 15.40 Å width. We consider two natural processes observed in practice, the hydroxylation of Si pores (introduction of OH groups into pores) and the penetration of water molecules into Si pores, as well as their impact on the electronic spectrum and optical properties of Si superstructures. The penetration of OH groups into the pores of the smallest 7.34 Å width causes a disintegration of hydroxyl groups and forms non-bonded protons which might be a reason for proton conductivity of porous Si. The porosity of silicon increases the extinction coefficient, k, in the visible range of the spectrum. The water structuring in pores of various diameters is analysed in detail. By using the bond valence sum approach we demonstrate that the types and geometry of most of hydrogen bonds created within the pores manifest a structural evolution from distorted hydrogen bonds inherent to small pores (∼7 Å) to typical hydrogen bonds observed by us in larger pores (∼15 Å) which are consistent with those observed in a wide database of inorganic crystals.


Introduction
Silicon, Si, has been one of the most important technological materials for the past decades. The main physical properties of bulk silicon are well known being widely used in technology. The bulk silicon is a diamond type crystal of cubic symmetry. Another morphological form of silicon, porous silicon, characterized by a set of elongated pores of various sizes and shapes is a macroscopically anisotropic material. The range of pore widths is very broad comprising the micropores (<2 nm), mesopores (2-50 nm) and macropores (>50 nm) depending on bulk Si substrate Table 1. Crystalline and computational parameters of the bulk and nanoporous Si structures. n uc is the number of atoms in the unit cell; n(H 2 O) corresponds to the number of water molecules deposited in the unit cell; E g is the band-gap energy (in eV), notations dir and ind correspond to direct and indirect energy band-gaps, respectively; f max is the maximal force acting on each atom (in eV/Å). Experimental value of E g [45] was obtained at 302 K. Structural optimization was done within the Broyden-Fletcher-Goldfarb-Shanno algorithm [46]. In the worst case of the largest 6aSi supercell, the maximal forces acting on each atom were reduced to the value lower than 1.3 × 10 −4 eV/Å, whereas for all other structures considered in our work these forces were reduced to much lower values (see Table 1). For the cut-off energy, E cut , we used the same value of 1632 eV with the cut-off smearing of 13.6 eV for all structures considered here.

Space Morphology and Stability Conditions
Bulk silicon, Si, has a cubic crystal structure of diamond type (Fd3m space group, #227, a cubic = 5.4331 Å). In the present work, we investigate porous Si with periodically arranged square shaped through pores. In this case, the cubic symmetry of bulk silicon is lowered to the tetrahedral one (P4m2, #115, a = b = 3.8403, c = 5.4310 Å) which is a subgroup of Fd3m space group. Three tetrahedral supercells of porous silicon with periodic boundary conditions are considered in our paper, that is, 4a*4a*c (hereafter abbreviated as 4aSi), 5a*5a*c (5aSi) and 6a*6a*c (6aSi) with the following corresponding size of pores 7.34, 11.26 and 15.40 Å, respectively (see Figure 1). Note that a and c correspond to lattice parameters of bulk Si in the tetrahedral setting. The real microporous Si (pores < 2 nm) consists of the net of irregular pores of various shapes, sizes and inter-pore walls. In order to simulate the electronic properties of this material by first-principles methods that use periodic boundary conditions, we need to consider the periodic sequence of the same pores. To simplify the transition from the cubic structure inherent to the bulk Si to tetrahedrally coordinated pores, we decided to use the square-shaped voids. Placing the additional silicon ions into the corners of square-shaped pores of the smallest 7.34 Å diameter led to a stairs-like pore with the shape far from both circle and square. To preserve the uniformity of our consideration, we operated with square-shaped voids also for all higher dimensionality supercells. Note that the square-shaped void of 5 * 5 µm is a quite typical geometry for macroporous Si (>50 nm) used for the fabrication of photonic crystals [47]. Recently, the square-shaped 400 * 400 nm Si pores were successfully realized for nucleation of sub-micrometer protein crystals in a silicon mesh [48]. In our opinion, the shape of a pore plays a relatively minor role compared with pore diameter, supercell dimension and porosity of porous Si.
Due to the high chemical activity of hydroxyl OH groups, the pores of porous Si are normally filled with OH radicals that are present in the air at ambient conditions. We took into consideration this phenomenon in our calculations. Figure 2 presents the optimized crystal structure of three various supercells of porous Si decorated with hydroxyl groups. As seen in Figure 2b,c, in two supercells 5aSi and 6aSi, OH groups create a certain covalent bonding with the Si atoms deposited on the surfaces of pores. However, in the case of the smallest 4aSi supercell, the OH radicals turn out to be destroyed forming some amount of almost free hydrogen atoms shifted to the center of nanopore (Figure 2a). Both oxygen atoms coordinate with two Si atoms: the first has a linear coordination environment (R(O-Si) = 1.66 Å), the second lies in position with angular neighbor location (R(O-Si) = 1.74 Å, ∠O-Si-O -115.0 • ). Optimized interatomic O-H distances are in the range of 2.24-3.01 Å whereas the H-H distances turn out to be even smaller, i.e., 1.7 and 1.66 Å. For comparison, the inter-hydrogen distance in molecular hydrogen is ∼0.74 Å and in the molecular ion H + 2 is 1.06 Å [49]. From the crystal-chemical point of view, such a coordination geometry of hydrogen atoms seems to be rather problematic for practical implementation, though it explains a metal-like conductive behavior of this 4aSi structure that will be discussed below.
1.7 and 1.66Å. For comparison, the inter-hydrogen distance in molecular hydrogen is ∼0.74Å and in the molecular ion H + 2 is 1.06Å [48]. From the crystal-chemical point of view, such a coordination geometry of hydrogen atoms seems to be rather problematic for practical implementation, though it explains a metal-like conductive behavior of this 4aSi structure that will be discussed below. The next level of our modeling was embedding the water molecules inside the pores already decorated with hydroxyl OH groups. We dealt with a comparatively small number of H 2 O molecules in pores. Choosing a particular number of water molecules, we aimed at keeping a reasonable balance between, on the one hand, the reliability of the physical phenomenon of water penetration into the Si pores due to the air humidity and the available computing resources, on the other hand. The water has a definite impact on the energy band spectrum of porous Si. One may come to this conclusion by comparing the band-gap energies for the same 5aSi supercell filled with the different number of water molecules, namely E g (5aSi+8H 2 O)= 0.809 eV and E g (5aSi+16H 2 O)= 0.766 eV (see Table 1). However, the precise quantitative expressions showing how a certain number of water molecules influences the electronic band spectrum for certain Si supercells could not be realized within our research because of the considerable computing time needed to The next level of our modeling was embedding the water molecules inside the pores already decorated with hydroxyl OH groups. We dealt with a comparatively small number of H 2 O molecules in pores. Choosing a particular number of water molecules, we aimed at keeping a reasonable balance between, on the one hand, the reliability of the physical phenomenon of water penetration into the Si pores due to the air humidity and the available computing resources, on the other hand. The water has a definite impact on the energy band spectrum of porous Si. One may come to this conclusion by comparing the band-gap energies for the same 5aSi supercell filled with the different number of water molecules, namely E g (5aSi+8H 2 O) = 0.809 eV and E g (5aSi+16H 2 O) = 0.766 eV (see Table 1). However, the precise quantitative expressions showing how a certain number of water molecules influences the electronic band spectrum for certain Si supercells could not be realized within our research because of the considerable computing time needed to accomplish this task. Figure 3 depicts four arrangements of water molecules in porous Si selected by us for further consideration. The first and fourth cases, Figure 3a,d, correspond to the smallest and largest pores, 4aSi and 6aSi, filled with 8 and 24 H 2 O molecules, respectively. The second and third arrangements, Figure 3b,c, represent the middle pore structure, 5aSi, filled with 8 and 16 molecules of water, respectively. As seen in this Figure, an increase of water molecules in the same supercell, 5aSi, causes a complete reconstruction of water arrangement. Consequently, the whole net of hydrogen bonds, undergoes a significant rearrangement which will be discussed below. Note that tetrahedral symmetry is inherent to the whole set of hydroxyl groups, water molecules and hydrogen bonds of all superstructures considered in our work.
It turns out that the geometry of H 2 O molecules greatly depends on the size of the pores. In Table 2 Figure 3a) placed in the smallest 4aSi pores (pore size is 7.34 Å). Such a comparably large deviation may be explained by a rather peculiar crystal field created by a set of oxygen and hydrogen atoms confined within the tight pore of 7.34 Å. With increasing pore size to 11.26 and 15.40 Å of 5aSi and 6aSi structures, respectively, the geometry of water molecules tends to its typical average form (see Table 2). In 6aSi*H 2   A complex network of hydrogen bonds created in Si pores is another peculiar feature of porous Si that deserves to be noted here. In Figure 3 we depicted only hydrogen bonds of the lengths smaller than 2.6 Å due to the baffling complexity of the total net of O-H· · ·O bonds. According to the paper by Zhang and Xue [51] based on the analysis of 404 hydrogen bonding inorganic systems from 128 structure types, most of single hydrogen bonds and the majority of multi-furcated hydrogen bonds satisfy the criterion R(O-H· · ·O) ≤ 2.35 Å and ∠O-H· · ·O ≥ 140 • . Most of the hydrogen bonds listed in Table 2 match these conditions. More precisely, all R(O-H) distances calculated in the present work lie within the range 0.97-1.03 Å. Nearly 75% of hydrogen bonding systems from Zhang and Xue's database [51] manifest similar R(O-H) values. Almost all R(H· · ·O) distances and ∠O-H· · ·O angles calculated in our work are placed in the 1.57-2.26 Å and 137-179 • regions, respectively, which are consistent with ∼90% of the corresponding experimentally observed values of single bonds and major components of multi-furcated hydrogen bonds [51].
We would like to emphasize some structural peculiarities of 4aSi+H 2 O supercell. Almost all oxygen atoms are involved in strong covalent O-Si or donor O-H chemical bonds, with the only exception of O 2 atom which is involved in acceptor H 3 · · ·O 2 , H 1 · · ·O 2 and H 2 · · ·O 2 bonds (Figure 3a). H 2 atom is placed in the middle between two neighbor O 2 atoms creating thus some unexpected kind of O 2 · · ·H 2 · · ·O 2 hydrogen bond. A similar spatial structure is observed in HF − 2 ion [52]. However, the presence of symmetric hydrogen bonds has not been yet established among oxygen-containing compounds by diffraction methods. Meanwhile, this structural feature resembles the structural peculiarity observed in some hydrogen bonding crystals, for example, of KH 2 PO 4 type in which the center location of H atoms in hydrogen bonds describes the static time-averaged picture of high-symmetry center-symmetric phase. On the dynamical scale of very short times, the H atoms make a flip-flop motion between two off-center positions in the hydrogen bond [31][32][33][34][35]. The involvement of H 2 atom in O 2 · · ·H 2 · · ·O 2 bond eliminates the metal-like conductivity inherent to 4aSi+OH supercell creating a definite energy gap E g = 0.752 eV in 4aSi+H 2 O supercell (see Figure 4c). Therefore, one may come to a rather paradoxical conclusion that embedding the water molecules into 4aSi+OH pores changes its conductive properties from metal-like to insulator-like ones.   In order to check the crystal chemistry correctness of our structural modeling, we use the bond valence sum (BVS) concept [53]. According to this approach, the V H valence of a hydrogen atom is the sum of the individual bond valences v i of O-H or O· · ·H bonds surrounding a certain H hydrogen atom, that is, where R i is an observed length of i-th bond, R 0 = 0.831 Å and b = 0.565 Å are the fitting constants evaluated for most of the inorganic crystals from the database analyzed in Reference [53] that satisfied the V H ≈ 1 condition. Calculating the BVS for certain hydrogen atoms of all structures considered in our paper we took into account only those oxygen atoms that satisfy the geometrical criterion R(O· · ·H) ≤ 3.04 Å, that is, a sum of van der Waal's radii of oxygen. As seen in Table 2

Electronic Properties
The electronic spectrum of porous silicon along with the density of electronic states is depicted in Figures 4 and 5, respectively. The energy levels of the valence band are drawn in Figure 4 in blue color whereas the conduction band is depicted in red color. Energy band-gap of each silicon porous structure is presented in Table 1.
Transformation of the space morphology from the bulk to porous form causes a splitting of the valence band into two valence sub-bands separated by a distinct energy gap. Several features merit special attention. The most peculiar feature of electronic spectra of porous Si is the absence of an energy gap between the valence and conductive bands in the spectrum of the smallest 4aSi pores containing hydroxyl OH groups (see Figure 4b). This 4aSi+OH superstructure reveals metal-like conductive properties. The probable reason is the presence of four unbounded protons in the supercell. As seen from the representative partial DOS of hydrogen, oxygen and silicon atoms of 4aSi+OH hydroxylated structure (Figure 6), the main contribution to DOS near Fermi level comes from hydrogen atoms (zero set to Fermi level). Oxygen and silicon atoms give a certain but much smaller contribution compared to hydrogen.
Four free hydrogen ions are also inherent to 4aSi+H 2 O superstructure (4aSi+OH supercell + 8 molecules of H 2 O, see Figure 3a). However, this superstructure manifests a quite distinct band gap of 0.752 eV (Figures 4c and 5c). The most probable reason of this difference is the hydrogen O 2 · · ·H 3 · · ·O 2 bonds which effectively tie together four H 3 protons in 4aSi+H 2 O supercell.  Another peculiarity of the energy band spectra of all nanoporous Si materials calculated by us is the lower energy band gap compared to the experimental E g value observed in the bulk Si. As indicated in Table 1, the calculated E bulk g = 0.575 eV of the bulk silicon is almost two times smaller than the experimental value of 1.14 eV [45]. The reason of such underestimation is the well-known DFT band-gap problem [54]. However, as follows from Table 1, the energy band gap tends to increase with an increasing pore diameter of d. In Figure 7, we plot the band-gap values as a function of the pore spacing ratio, d a superc , where a superc is a unit cell constant of supercell. In this figure, we used the average E g value of all corresponding values calculated for a certain supercell, that is, of 4aSi, 5aSi and 6aSi. As seen in this figure, the energy band gap obeys non-linear law, The increase of E g with an increasing porosity does not seem strange at all because the total volume of the empty space increases with an increase of the pore diameter. Note that the porosity changes the type of the forbidden energy gap, thus changing it from an indirect type inherent to bulk homogeneous Si to a direct type detected by us in porous Si with 6aSi supercell (see Table 1).

Optical Properties
Crystalline silicon is an optically isotropic material with the same refractive index in any direction (Fd3m cubic symmetry). Porous silicon loses the cubic symmetry inherent to the bulk Si and becomes an optically anisotropic material. We have chosen the tetragonal subgroup P4m2 of cubic space group for theoretical treatment of porous Si. This symmetry admits two different refractive indices (optical anisotropy) which we hereafter denote as n 1 = n 2 and n 3 . Calculating the high-frequency dispersion of the complex dielectric ε * function and using the following relations one may easily obtain the dispersion of refractive indices n 1 and n 3 and extinction coefficients k 1 and k 3 .
In Figure 8 we present the dispersion of refractive indices and extinction coefficients of various Si space morphologies. As seen in Figure 8a,b, there is a very good agreement in the visible range between the calculated refractive index and the extinction coefficient of the bulk cubic Si and the experimental data [55]. Compared to the bulk silicon, the porosity of Si generally decreases the refractive indices within the whole calculated spectral range and increases the extinction coefficient k in the visible optical region (Figure 8a,b). The same tendency of increasing the extinction coefficient in the visible range at transition from bulk three-dimensional material to the ordered sequence of one-dimensional nano-rods was formerly observed in PbHPO 4 crystal [56]. Moreover, the Si porosity increases the optical anisotropy (n 1 > n 3 ) especially in ultraviolet region and in the range higher than the near ultraviolet one (above 3.0 eV). For representative purpose in Figure 8c,d we depict only the data of 5aSi porous superstructure. Filling the pores with water molecules causes a decrease of the n 1 refractive index to the value comparable with n 3 (see Figure 8c), thus reducing the optical birefringence.
It is quite instructive to compare the results of our calculations with experimental data. Formerly, the model dielectric function (MDF) approach [19] was successfully used for characterization of ellipsometric reflectivity spectra measured on porous sponge-like Si thin films with different nanocrystal (wall) sizes, i.e., 12, 6 and 3 nm. Based on the fitted parameters, the authors of Reference [19] established a high-frequency dielectric dispersion of porous Si thin films with various nano-grains, which are reproduced in Figure 9 by symbols. In this figure, we also present the dispersion of real and imaginary parts, Re(ε 1 ) and Im(ε 1 ), of dielectric function ε * 1 of the same Si space morphologies as those depicted in Figure 8a,b. All the data simulated by us are drawn in Figure 9 by lines. As seen in this Figure, for two marginal cases, that is, for the bulk Si and for the smallest 4aSi supercell, there is a good qualitative agreement between our data and those published by Petrik et al. [19]. Real part, Re(ε 1 ), obtained within MDF approach for 12 nm nano-sized Si crystals almost coincides in the visible optical range with our results calculated for bulk Si crystal. Imaginary part, Im(ε 1 ), shows a worse correlation since our data present the sequence of oscillation peaks whereas the Petrik et al's data manifest two broad peaks. There is a good quantitative correlation for real part Re(ε 1 ) between our 4aSi data (supercell dimension is 1.536 nm) and 3 nm result of Reference [19] (see Figure 9a). The quantitative correlation for Im(ε 1 ) is observed below ∼3.5 eV, showing only a similar qualitative energy dependence above 3.5 eV. As seen in Figure 9, Re(ε 1 ) and Im(ε 1 ) of 6 nm sized Si thin films manifest an intermediate energy dispersion between two marginal cases, 3 nm and 12 nm. The 6 nm Si supercell was not considered in our study. Summarizing this comparison with Petrik et al.'s data, one may draw the following conclusion. First-principles calculation of dielectric properties of nano-porous Si feels the dimensionality effects up to ∼10 nm supercell. Larger Si supercells give the dielectric response similar to the bulk Si modification. Smallest Si supercells of an order of ∼1.5 to 3 nm, demonstrate practically the same dielectric response below 4 eV irrespective of the calculation method, that is, ab initio or MDF model.    Figure 9. Dispersion of dielectric function, Re(ε 1 ) (a) and Im(ε 1 ) (b), of bulk and porous Si of 4aSi and 5aSi supercells. The data represented by symbols are reproduced from Reference [19] for porous Si films with nano-grains of different size. All lines correspond to data calculated in our paper.
Optical anisotropy of nanoporous Si composites in the optical transparency region (E < 1 eV) deserves a special discussion. Isotropic media with uniaxially aligned cylindrical nanopores are characterized by optical anisotropy [57] also referred to as geometric birefringence ∆n. Figure 10 shows the calculated refractive indices and respective ∆n in the optical transparency region for two porous Si supercells, 4aSi and 5aSi, decorated with hydroxyl or hydroxyl-water surface layers, that is, 4aSi+OH+8H 2 O, 5aSi+OH and 5aSi+OH+16H 2 O. For pure porous Si of porosity P = 0.25, that is, for the same porosity as in the case of our 4aSi structure, the generalized Bruggeman equation (for evaluation methodology see Reference [57]) gives the refractive indices n 1 = 2.753 and n 3 = 3.035 and thus a positive birefringence (∆n = n 3 − n 1 >0). Note that the effective medium approach provides here a positive birefringence in the entire range of porosities (0 < P < 1). Our ab initio calculations, by contrast, predict a weak but negative birefringence, ∆n < 0, for porous Si with hydroxyl 5aSi+OH or hydroxyl-water surface layers 5aSi+OH+16H 2 O, see Figure 10. Probably, such a difference of optical properties may be attributed to the strong optical anisotropy of the ordered OH and OH+H 2 O layers which may result in a considerable negative contribution to the optical birefringence. In other words, the effective (averaged) birefringence, ∆n, may be represented as superposition of the positive geometrical birefringence of the porous Si matrix alone, ∆n Si > 0, and the negative birefringence of the anisotropic hydroxyl or hydroxyl-water layer(s) covering the pore walls, ∆n OH+H 2 O < 0. Ab initio simulation shows that the molecular arrangement in the "adsorbed" OH and H 2 O layers should play a crucial role in the effective optic anisotropy of nanoporous Si+OH+H 2 composite materials. A prominent example here may be the 4aSi+OH+8H 2 O composite, which in contrast to the 5aSi+OH and 5aSi+OH+16H 2 O ones, is characterized by a weak but positive birefringence. Apparently, this means that ∆n Si + ∆n OH+H 2 O = ∆n < 0 for 5aSi+OH and 5aSi+OH+16H 2  ∆n > 0 Figure 10. Refractive indices, n 1 and n 3 , and respective birefringence ∆n of two supercells, 4aSi and 5aSi, filled with hydroxyl groups and water molecules. All full symbols correspond to n 1 index whereas all empty symbols correspond to n 3 index.

Conclusions
We performed a detailed computational study of porous silicon filled with water molecules. We established the impact of nanoconfinement on the water structuration and hydrogen bonding within the pores of 7.34, 11.26 and 15.40 Å. Hydroxyl OH groups confined within the Si pores of 7.34 Å diameter appear to be disintegrated, thus forming free non-bonded protons which might be the reason for proton conductivity. However, the embedding of water molecules in these pores of the smallest diameter destroys the proton conductivity, thus retrieving the small but distinct energy gap in the electronic spectrum.
Energy band-gap of porous Si filled with water molecules increases with an increasing pore diameter according to nearly linear law. According to our calculations, the energy band-gap of bulk and porous Si of small diameter d 11.26 Å (4aSi and 5aSi supercells) is of indirect type and the water filling of the Si pores does not transform the type of E ind . However, with an increase of the pore diameter to 15.40 Å (6aSi supercell), the indirect band-gap E ind changes to the direct E dir one. Embedding the water to 6aSi supercell enlarges the E dir band-gap.
The calculated dispersion of refractive indices and the extinction coefficient of bulk Si nicely agree with the experimental data proving the reliability of our calculations. Si porosity augments the optical anisotropy (n 1 > n 3 ), especially in ultraviolet region and enlarges the extinction coefficients k of porous Si in the visible spectral range.
We described the structural evolution of water molecules and hydrogen bonds at an increased Si pore diameter and found that the geometry of hydrogen bonds deviates from the typical values due to the space confinement in the smallest 7.34 Å pores but tends to the normal values with enlarging the pore diameter to 15.40 Å.
Optical properties of porous Si decorated with hydroxyl OH and water H 2 O layers strongly depend on Si porosity and on OH+H 2 O morphology inside the pores. The first-principles modeling, in this respect, may be considered as an efficient tool in the technology of nanomaterials, for preliminary evaluation, prediction and/or tailoring of their optical properties.
Author Contributions: Conceptualization, methodology, calculation, data analysis and discussion, Y.S.; crystal structure visualization, hydrogen bonding analysis, O.P.; project administration, funding acquisition and discussion, A.S.A.; discussion, computing resources, S.V.; data analysis, discussion, computing resources and funding acquisition, A.V.K. All authors have read and agreed to the published version of the manuscript.

Conflicts of Interest:
The authors declare no conflict of interest.