DOI:
10.1039/C9NR08209C
(Paper)
Nanoscale, 2020,
12, 14641477
Temperatureinduced nanostructural evolution of hydrogenrich voids in amorphous silicon: a firstprinciples study†
Received
24th September 2019
, Accepted 14th November 2019
First published on 14th November 2019
The paper presents an ab initio study of temperatureinduced nanostructural evolution of hydrogenrich voids in amorphous silicon. By using large aSi models, obtained from classical moleculardynamics simulations, with a realistic voidvolume density of 0.2%, the dynamics of Si and H atoms on the surface of the nanometersize cavities were studied and their effects on the shape and size of the voids were examined using firstprinciples densityfunctional simulations. The results from ab initio calculations were compared with those obtained from using the modified Stillinger–Weber potential. The temperatureinduced nanostructural evolution of the voids was examined by analyzing the threedimensional distribution of Si and H atoms on/near void surfaces using the convexhull approximation, and computing the radius of gyration of the corresponding convex hulls. A comparison of the results with those from the simulated values of the intensity in smallangle Xray scattering of aSi/aSi:H in the Guinier approximation is also provided, along with a discussion on the dynamics of bonded and nonbonded hydrogen in the vicinity of voids.
I. Introduction
Amorphous silicon and its hydrogenated counterpart have a wide range of applications, from photovoltaics to thinfilm technology.^{1,2} Thin films of hydrogenated amorphous silicon (aSi:H) are extensively used for passivation of crystalline silicon (cSi) surfaces, which is essential to produce a high opencircuit voltage in siliconbased heterojunction solar cells.^{3,4} The passivation of Si dangling bonds by H atoms in the vicinity of aSi:H/cSi interfaces plays a key role in improving the properties of silicon heterostructures.^{5} Postdeposition annealing is routinely used for structural relaxation of aSi:H samples and incorporation of hydrogen at the aSi:H/cSi interfaces. The presence of hydrogen to an extent thus is often preferred^{6} for the preparation of aSi:H samples with good interface properties. Experimental studies, such as spectroscopic ellipsometry^{7,8} and Raman spectroscopy,^{9} indicate that, although the presence of nanometersize voids can degrade the electronic quality of the samples, a small concentration of voids can improve the degree of surface passivation. This is generally believed to be due to the presence of mobile H atoms in a voidrich environment, where hydrogen can reach aSi:H/cSi interfaces, via diffusion or other mechanisms, in comparison to dense aSi:H samples.
Heat treatment or annealing is an effective tool for nanostructural relaxation of laboratorygrown samples. Annealing at low to medium temperature (400–700 K) can considerably improve the quality of asdeposited samples, by removing unwanted impurities and reducing imperfections (e.g., defects) in the samples.^{10} For amorphous materials, annealing also increases the degree of local ordering in the amorphous environment. Earlier studies by smallangle Xray scattering (SAXS) reported the presence of columnarlike geometric structure^{11} and blister formation in annealed samples of aSi:H.^{12} Likewise, experiments on aSi/aSi:H, using positronannihilation lifetime spectroscopy, have indicated that annealing above 673 K can lead to the formation of voids via vacancy clustering, which affects the nanostructural properties of the network.^{13,14} It has been suggested that at higher temperature, near 800 K, a considerable restructuring can take place in the vicinity of the void boundary that can modify the shape and size of the voids.^{15} Although tilting SAXS^{15,16} can provide some information on the geometrical shape of the voids, a direct experimental determination of annealing effects on the morphology of voids in the amorphous environment is highly nontrivial, as scattering experiments generally provide integrated or scalar information on the size and shape of the voids. By contrast, computational modeling can yield reliable structural information on the morphology of voids,^{17,18} provided that highquality large models of aSi/aSi:H – with a linear dimension of several nanometers – are readily available and accurate totalenergy functionals to describe the atomic dynamics of Si and H atoms are in place. However, despite numerous studies on computational modeling of aSi and aSi:H over the past several decades, there exist only a few computational studies^{17–21} that truly attempted to address the structure and morphology of voids on the nanometer length scale in recent years.
The purpose of the present study is to conduct ab initio simulations of hydrogenrich voids in large models of aSi at low and high annealing temperatures. We are particularly interested in the structural reconstruction of the void surfaces, induced by annealing, and the resulting effects on the intensity in smallangle Xray scattering (SAXS). Further, we intend to study the evolution of the shape of voids with temperature in the presence of hydrogen inside the voids. To this end, we develop significantly large models of aSi, which are characterized by realistic distributions of voids as far as the results from SAXS,^{16,22} nuclear magnetic resonance (NMR),^{23} infrared (IR),^{24,25} and positronannihilation^{26} measurements are concerned. The term realistic here refers to the fact that the models must be able to accommodate nanometersize voids with a voidvolume fraction of 0.2–0.3%, as observed in SAXS, IR, and NMR measurements. The threedimensional shapes of the voids are studied by using the convexhull approximation^{17} of the boundary region of the voids, upon annealing at 300 K and 800 K, and the effects of the shape and size of the voids on the simulated SAXS intensity from the models are examined. The presence of molecular hydrogen and the statistics of various siliconhydrogen bonding configurations near the void surfaces are also examined by introducing H atoms within the voids, with a number density of hydrogen (inside cavities) as observed in infrared spectroscopy and NMR measurements.^{24,27}
The rest of the paper is as follows. In section II, we briefly describe the simulation of large models of aSi with a linear dimension of a few tens of nanometers and the generation of nanometersize voids in the resulting models in order to incorporate a realistic distribution of voids, with a voidvolume fraction as observed in SAXS and IR measurements. Section III discusses the results of our work, with particular emphasis on the threedimensional shape and size of the voids upon annealing at 300 K and 800 K. The correlation between the threedimensional shape and size of the voids and the intensity of SAXS from simulations is discussed here. A comparison of the results with those from positronannihilation lifetime (PAL) spectroscopy^{13} and Dopplerbroadening positronannihilation spectroscopy (DBPAS)^{28} are presented. This is followed by a discussion on the atomic dynamics of hydrogen and the time evolution of monohydride siliconhydrogen (Si–H) bond formation in the vicinity of the voids, using ab initio moleculardynamics (AIMD) simulations over several picoseconds. Section IV presents the conclusions of our work.
II. Computational methods
A. Generation of large models of aSi
We began by generating a large model of amorphous silicon in order to be able to study the intensity due to SAXS in the smallangle region of scattering wavevector k ≤ 1.0 Å^{−1}. Since the key purpose of the study is to examine the structural evolution of the voids in aSi on the nanometer length scale, by joint analyses of simulated SAXS intensities and threedimensional distributions of atoms near void surfaces, it is necessary to simulate a sufficiently large model that can accommodate a realistic distribution of nanometersize voids with the correct voidvolume fraction, in the range from 0.01% to 0.3%, as observed from SAXS and IR measurements.^{15,16} To this end, a model aSi network, consisting of 262400 atoms, was generated using classical moleculardynamics (MD) simulations, followed by geometry relaxation, by employing the modified Stillinger–Weber potential.^{29,30} The initial configuration was obtained by randomly generating atomic positions in a cubic box of length 176.123 Å, which corresponds to the experimental mass density, 2.24 g cm^{−3}, of aSi.^{31} Thereafter, constanttemperature MD simulations were performed using a chain of Nosé–Hoover thermostats^{32,33} and a time step of Δt = 1 fs. During the MD runs, the system was heated to a temperature of 1800 K for 10 ps, and then it was cooled from 1800 K to 300 K at the rate of 10 K ps^{−1}. The heating–cooling cycle was carried out for at least 20 iterations. For each cycle, structural properties of the system, such as the coordinationnumber statistics, the bondangle distribution, and the structure factor of the evolving configurations were examined. The procedure continues until it satisfies the convergence criteria that entail imposing a restriction on the upper limit of the rootmeansquare (RMS) width of the bondangle distribution (Δθ_{max} ≤ 10°) and the coordinationdefect density (of 3% or less) in the network. It may be noted that these variables are in conflict with each other, meaning that one can be reduced at the expense of increasing the other, and vice versa. Thus, care was taken to ensure that these criteria were simultaneously satisfied. Once the structural properties satisfied the aforementioned convergence criteria, the MD run was terminated and the system was relaxed by minimizing the total energy of the system with respect to the atomic positions. The resulting configuration was then used as a structural basis for the generation of additional models of aSi by introducing voids in the network. The details of the simulation procedure and an analysis of the structural and electronic properties of these large aSi models can be found in ref. 34.
B. Void geometry and distribution in aSi
Having obtained a large model of aSi, we proceed to generate a realistic distribution of voids. Since computational modeling of void formation in aSi, by mimicking the actual deposition and growth processes, is highly nontrivial, we generated a void distribution using experimental information on the shape, size, and the number density of voids. Experimental data from an array of measurements, such as smallangle scattering (SAS),^{16,22,35,36} infrared (IR) spectroscopy,^{24} nuclear magnetic resonance (NMR),^{23} positronannihilation lifetime (PAL) spectroscopy,^{26} and calorimetry measurements,^{27,37,38} suggest that the linear size of the voids in aSi typically lies between 10 Å and 40 Å and that the voidvolume fraction can range from 0.01% to 0.3% of the total volume, depending upon the growth process, deposition rate, and the consequent electronic quality of the samples. Following these observations, we constructed a void distribution in the large model with the following properties: (a) the voids are spherical in shape and have a diameter of 12 Å; (b) the number density of the voids corresponds to a voidvolume fraction of 0.2%; (c) the voids are sparsely and randomly distributed, so that the distance between any two void centers is significantly larger than the dimension of the voids.
Fig. 1 shows an example of a void embedded in a twodimensional amorphous network. The void region is indicated by the gray circle of radius R_{v}, which is surrounded by a layer of Si atoms (yellow) with a thickness d. Silicon atoms in the layer constitute the surface/wall of the void, which we shall refer to as the voidsurface atoms. The remaining Si atoms, within the cubic box (of length l), will be referred to as bulk Si atoms and are indicated as black and blue circles. To study the effect of hydrogen on the evolution of void surfaces, H atoms were introduced inside the voids so that no two H atoms were at a distance less than 1 Å from each other. The presence of H atoms in the network requires that the problem must be treated at the quantummechanical level using, for example, firstprinciples densityfunctional theory. However, given the size of the large model, such a task is hopelessly difficult and computationally infeasible, and an approximation of some sort is necessary. Here, we employed the local approach to electronicstructure calculations by invoking the principle of nearsightedness of an equilibrium system, as proposed by Kohn,^{39} and address the problem by studying an appropriately large subsystem, defined by a cubical box of length l surrounding each void (see Fig. 1). Following an (approximate) invariance theorem on electronic structure of disordered alloys, due to Heine and coworkers,^{40,41} the quantum local density of states near the voids is essentially independent of the boundary conditions at a sufficiently large distance,^{42} for example, the presence of a fixed layer of Si atoms (blue) on the edge of the box. Referring to Fig. 1, we chose R_{v} = 6 Å, d = 3 Å, and l = 36 Å. To ensure that the voidvolume fraction is consistent with the experimentally observed value of 0.01%–0.3%, we created 12 voids, reflecting a voidvolume fraction of 0.2%, which were randomly distributed in the large aSi model of linear dimension 176.123 Å. An outer layer of Si atoms (blue) of thickness 3 Å was held fixed at the edge of the box of length l for each void during firstprinciples simulations and totalenergy relaxation of the subsystem so that the latter can be put back into the original large model for the calculation of the SAXS intensity, upon completion of annealing and geometry relaxation without boundary perturbations. To ensure that the SAXS intensity from the large aSi sample with voids is not affected by the atoms on the fixed boundary layer (shown in blue color in Fig. 1) of the subsystems or supercells, the (undercoordinated) atoms on the boundary layer were left undisturbed during simulations. All other atoms were allowed to move during classical and quantummechanical simulations. Each subsystem, consisting of a cubical box of length 36 Å with a central void region of radius 6 Å, contained about 2200 atoms. The size of the cubical box was chosen carefully so that the local structural and electronic properties of the atoms in the vicinity of void surfaces were not affected by the boundary layer on the edge of the box, at a distance of 12 Å from the void surfaces. Care was also taken to ensure that the centers of the voids were sparsely distributed in the original large model so that the individual supercells did not overlap with each other. Specifically, the positions of the center of the voids were generated randomly in such a way that the centertocenter distance between any two neighboring voids in the large model is at least 36 Å or more. To study the effect of hydrogen on void surfaces, each void was loaded with 30 H atoms, as observed in IR and NMR measurements.^{43} The initial distributions of the H atoms were generated to satisfy the following two criteria: (a) the H atoms were randomly distributed so that no two H atoms were separated by a distance less than 1 Å; (b) the distance between any H atoms inside a void and a Si atom on the void surface was greater than or equal to 2 Å. The latter was achieved by confining the H atoms within a sphere of radius 4 Å from the center of the void. These criteria prevent the system from developing a sudden large force on light H atoms and an instantaneous formation of siliconhydrogen bonds on void surfaces or shooting off H atoms from the void region at the beginning of AIMD simulations.

 Fig. 1 A schematic representation of a void (gray region) in two dimensions. The annular region with Si atoms (yellow circles) indicates the void boundary, with a few H atoms (red circles) inside the void. The bulk Si atoms are shown as black and blue circles – the latter indicate a layer of fixed Si atoms for the purpose of local AIMD simulations, as discussed in the text.  
C. Ab initio simulations of voids in aSi
In the preceding section, we have generated a number of models (supercells/subsystems), containing as many as 2200 atoms and a central void of radius 6 Å, using the modified Stillinger–Weber (SW) potential. However, since the SW potential may not accurately describe the chemistry of silicon surfaces, it is necessary to treat these models using a firstprinciples totalenergy functional from density functional theory (DFT) in order to reduce artifacts, such as the number of excessive dangling bonds or defects, associated with the SWgenerated void surfaces. The standard DFT protocol here is to compute the selfconsistentfield (SCF) solution of the Kohn–Sham (KS) equation, and the presence of voids or inhomogeneities in the system suggests that the generalized gradient approximation (GGA) should be employed in order to deal with the atomic density fluctuations near void surfaces. However, it is simply not feasible to solve the KS equations selfconsistently for a problem of this size (i.e., about 2200 atoms) and complexity using the existing planewavebased/localbasis DFT methodology. In particular, the computational cost for conducting SCF AIMD runs for several picoseconds is prohibitively large, even for a single configuration. To reduce the computational complexity of the problem, it is thus necessary to employ a suite of approximations for solving the KS equation. Here, we invoke the nonselfconsistent Harrisfunctional approach^{44,45} in the local density approximation (LDA), using a minimal set of basis functions for Si atoms. While these approximations may affect the accurate determination of some aspects of the electronic structure of void surfaces, the results from test calculations on 1000atom models with a void of radius 5 Å, obtained via the full SCF procedure for solving the KS equation in the LDA and GGA, indicate that the structural properties of the void surfaces in aSi are minimally affected by these approximations. In particular, the shape and size of the voids, as well as the values of the average bond angle and bond length of the silicon atoms on the reconstructed void surface, exhibit a very little to no changes with the varying level of DFT approximations (see ESI† for results from test calculations). A discussion on the accuracy and computational efficiency of the Harrisfunctional approach to treat bulk aSi and aSi:H can be found in ref. 46.
Ab initio calculations and totalenergy optimization were conducted using the firstprinciples densityfunctional code SIESTA.^{47} The latter employs a localized basis set and normconserving Troullier–Martins pseudopotentials,^{48} which are factorized in the Kleinman–Bylander^{49} form to remove the effect of core electrons. Electronic and exchange correlations between electrons are handled using the LDA, by employing the Perdew–Zunger parameterization of the exchange–correlation functional.^{50} As stated earlier, we employed the Harrisfunctional approach^{44,45} that involves the linearization of the Kohn–Sham equations for structural relaxation and AIMD simulations. The latter were conducted over a period of 10 ps, with a time step of 1 fs. Singlezeta (SZ) and doublezetapolarized (DZP) basis functions were employed for Si and H, respectively. Here, SZ functions correspond to one s and three p orbitals for Si atoms, whereas DZP functions refer to two s and three p orbitals for H atoms. The size of the simulation box or supercell is found to be large enough in this work so that the Brillouin zone collapses onto the Γ point, which was used for the Brillouinzone sampling.
D. Simulation of SAXS intensity and voidsize determination
The SAXS intensity for a binary system can be expressed with the aid of the partial structure factors and the atomic form factors of the constituent atoms. For a binary system, the scattering intensity, I(k), can be expressed as,^{51} 
 (1) 
where 
 (2) 
is the partial structure factor associated with the atoms of type i and j such that S_{ij} = S_{ji}, and g_{ij}(r) is the partial pairdistribution function. Here, and f_{i} are the atomic fraction and the atomic form factor of the atoms of type i, respectively. In eqn (2), the partial pairdistribution function, g_{ij}(r), is defined as g_{ij}(r) = n_{ij}/n_{j}, where n_{ij} is the average number density of atoms of type j with an atom of type i at the center, n_{j} = x_{j}n_{0}, and n_{0} = N/V. Once the scattering intensity is available, either numerically or experimentally, the size of the voids, or any extended inhomogeneities, can be estimated by invoking the Guinier approximation^{52} in the smallangle limit, kR_{G} ≤ 1, 
 (3) 
where R_{G} is the characteristic size of the inhomogeneities in the dilute concentration limit. While eqn (3) provides a reasonable estimate of void sizes from experimental data, difficulties arise in extracting the R_{G} value from a lnI–k^{2} plot in simulations, owing to the finite size of the models. Since the effective smallangle region of k is approximately given by 4π/L < k < 1/R_{G}, where L is the linear dimension of the sample, a sufficiently large model is needed in order to obtain the mean value of R_{G}, via averaging over several independent configurations. However, it has been observed^{53} that, even for very large models, the lower limit of k cannot be reduced arbitrarily by increasing L due to the presence of noise in rg(r) beyond r ≈ 25–30 Å. The term r[g(r) − 1] in eqn (2) for r > 25 Å introduces artifacts in S(k) for small values of k, which make it difficult to accurately compute R_{G} from the simulated intensity data. This necessitates the construction of alternative measures of void sizes and shapes directly from the distribution of atoms in the vicinity of voids. A simple but effective measure is to compute the gyrational radius from the distribution of the atoms on void surfaces between radii R_{v} and R_{v} + d (=R′_{v}), where R_{v} is the void radius and d is the width of the void surface, as discussed in section IIB. This is schematically illustrated in two dimensions in Fig. 2, where the positions of the voidsurface atoms before and after annealing are shown in yellow and red colors, respectively. Assuming that the voidsurface atoms (yellow) can displace, in the meansquare sense, from their initial position by x during annealing, the radius of gyration, R_{g}, for the assembly of void atoms (red) can be expressed as, 
 (4) 

 Fig. 2 An illustration of the (re)construction of the shape and size of a void formed by an assembly of voidsurface atoms in two dimensions. The distributions of the atoms before and after annealing are shown in yellow () and red () colors, respectively. The convex polygon (red line) approximates the void region in two dimensions after annealing.  
Here, r_{i} and r_{cm} are the position of the i^{th} atom and the center of mass of n_{s} voidsurface atoms, respectively. The atoms on the reconstructed surface are distributed within a spherical region of radius R′_{v} + x and H[y] is the Heaviside step function that asserts that only the atoms with y ≥ 0 contribute to R_{g}. Here, we chose a value of x = 1.5 Å for Si atoms in order to define the boundary of a reconstructed void surface.^{53} Since annealing at high temperature can introduce considerable restructuring of void surfaces, it is often convenient to invoke a suitable convex approximation to estimate the size and shape of the voids. To this end, we employ the convexhull approximation that entails constructing the minimal convex polyhedron that includes all the voidsurface atoms on the boundary. The size of the convex region can be expressed as the radius of gyration of the convex polyhedron, or convex hull,

 (5) 
In eqn (5), n_{H} is the number of atoms (or vertices) that defines the convex polyhedral surface, and _{h} is the center of mass of the polyhedron. To determine the degree of deviation of the shape of a (convex) void from an ideal sphere, we employ the convexhull volume (V_{H}) and surface area (A_{H}), and obtain the sphericity parameter,^{54}

 (6) 
The definition above is frequently used to measure the shape and roundness of sedimentary quartz particles, and it expresses the ratio of the surface area of a sphere to that of a nonspherical particle having an identical volume.^{54} It may be noted that a highly distorted void may not be accurately represented by a convex polyhedron and that the convexhull volume provides an upper bound of the actual void volume, leading to R_{H} ≥ R_{g}. Likewise, as observed by Porod,^{55} the R_{G} value obtained from the Guinier approximation may not be accurate for very anisometric voids.
III. Results and discussion
Before addressing the temperatureinduced nanostructural changes in voids, we briefly examine the large model of aSi, developed in section IIA, in order to validate its structural properties. Fig. 3 shows the simulated values of the Xray scattering intensity, along with the experimental values of the intensity for asimplanted (green) and annealed (red) samples of pure aSi, reported by Laaziri et al.^{56} The simulated values closely match with the experimental data in Fig. 3, especially for the annealed samples, in the wideangle region from 1 Å^{−1} to 10 Å^{−1}. Some characteristic structural properties of the large aSi model are also listed in Table 1. Together with the Xray intensity, the average bond angle and its deviation, 109.23° ± 9.26°, and the number of fourfoldcoordinated atoms, 97.6%, suggest that the structural properties of the model obtained from classical moleculardynamics simulations in section IIA are consistent with experimental results. It may be reasonably assumed that the presence of a small amount of coordination defects (≈2.4%), which involve a length scale of 2–3 Å and are sparsely distributed in the amorphous environment, would not affect the scattering intensity appreciably in the smallangle region of the scattering wavevector.^{57} A minor deviation in the height of the first peak in the simulated data in Fig. 3 from experimental values can be attributed in part to the use of the classical SW potential and partly to the large size of the model. A discussion on the validation of the structural and electronic properties of these large models can be found in ref. 34.

 Fig. 3 The Xray scattering intensity for aSi from experiments and classical moleculardynamics simulations. The experimental data^{56} correspond to the results from asimplanted () and annealed () samples, while the simulated values () of the intensity correspond to the modified Stillinger–Weber model of aSi.  
Table 1 Structural properties of the large model of aSi. N, L, ρ, C_{4}, r, θ, and Δθ are the number of atoms, simulation box length, mass density, fourfold coordination (in %), average bond length, average bond angle (in degree), and the rootmeansquare width of the bond angles, respectively
N

L (Å) 
ρ (g cm^{−3}) 
C
_{4} (%) 
r (Å) 
θ 
Δθ 
262400 
176.123 
2.24 
97.6 
2.39 
109.23° 
9.26° 
A. Structural evolution of voids in aSi at 300 K and 800 K
In this section, we discuss the threedimensional shape and size of the voids in pure aSi obtained from classical and quantummechanical annealing of the subsystems at 300 K and 800 K, followed by totalenergy optimization. Starting with an examination of the classical and ab initiogenerated void surfaces, we analyze the computed SAXS spectrum obtained from the large model, which is embedded with the supercells containing voids annealed at 300 K and 800 K.
Fig. 4 shows the structure of two representative voids (V3 and V16) after thermalization at 300 K for 10 ps, using classical and ab initio MD simulations. The void surfaces shown in Fig. 4 corresponds to Gaussianconvoluted isosurfaces from the Xcrysden package.^{58} The surface has been generated by placing a threedimensional Gaussian function at the center of each atom and ensuring that the function has a value of 2.0 and 1.0 at the center and on the surface of the atom, respectively. By choosing an appropriate isovalue, between 0 to 2, and a suitable radius for Si atom, from standard crystallographic databases (e.g., Cambridge Structural Database), it is possible to examine the morphological changes of voids, associated with the structural relaxation of the voidsurface atoms at different temperature. A close examination of Fig. 4(a) and (b) suggests that the resultant structures are very similar to each other as far as the void surfaces are concerned, except for small changes as indicated by white rings. A similar observation can be made for V16 in Fig. 4(c) and (d). This observation was found to be true for the rest of the voids as well, which showed very little changes on the void surfaces. These results are consistent with the variation of the simulated intensity values with wavevectors shown in Fig. 5. A somewhat more pronounced scattering from the classicallytreated voids (CMDR) in comparison to the ab initiotreated voids (AIMDR) can be attributed to a minor expansion of the voids during classical simulations. This is clearly evident from Fig. 6(a) and (b), where the radius of gyration and the volume of the voids, respectively, are obtained from the convexhull approximation of the void region. The polyhedral radii (R_{H}) for the voids from classical simulations were found to be consistently larger than the corresponding radii from ab initio simulations in Fig. 6(a). A similar conclusion applies to polyhedral volumes (V_{H}), which are plotted in Fig. 6(b). Once the linear size, area, and volume of the voids are available in the convex approximation, the sphericity of voids, Φ_{S}, can be determined from eqn (6). A further measure of the size of the voids follows from the variation of the SAXS intensity, I(k), with the wavevector, k, in the smallangle region. By invoking the Guinier approximation (see eqn (3)), the Guinier radius, R_{G}, can be obtained from ln I(k) vs. k^{2} plots. Fig. 6(c) shows the Guinier fits in the wavevector range, 0.15 Å ≤ k ≤ 0.3 Å, which resulted in an R_{G} value of 6.38 Å for the CMDR model and 8.07 Å for the AIMDR model. It may be noted that a somewhat smaller value of R_{G} was obtained for classicallytreated voids compared to their AIMD counterparts (cf.Table 2), despite increased scattering by the former in the smallangle region (see Fig. 5). This minor discrepancy between the two R_{G} values originated from the difficulty in computing the Guinier radius from the SAXS data in the wavevector region of 0.15–0.3 Å^{−1}. For computational consistency, we employed an identical wavevector range, namely 0.15–0.3 Å^{−1}, for the evaluation of R_{G}. However, this range may not necessarily be the same for CMDR and AIMDR cases and can affect the R_{G} value obtained from the Guinier's fit.

 Fig. 4 An isosurface representation of two voids (V3 and V16), obtained from the Gaussianconvolution of atomic positions, after thermalization at 300 K from classical and AIMD simulations: (a) classical V3; (b) AIMD V3; (c) classical V16; and (d) AIMD V16. The white rings in (a) and (b) indicate minor structural differences between the two surfaces in V3. See section IIIA for details.  

 Fig. 5 The SAXS intensity of the large model of aSi from classical (CMDR) and ab initio (AIMDR) MD simulations of voids at 300 K, followed by geometry relaxation.  

 Fig. 6 (a) The radius of gyration of the voids from the convexhull approximation at 300 K. The results for the CMDR () and AIMDR () models are presented here. (b) Estimated void volumes in the convex approximation. (c) Guinier plots for the CMDR and AIMDR models at 300 K. The Guinier radii from the plots correspond to a value of 6.38 Å (CMDR) and 8.07 Å (AIMDR).  
Table 2 Estimated volumes (Å^{3}), linear sizes (Å), and the average sphericities (Φ_{S}) of the voids in aSi with and without hydrogen inside voids. R_{X} indicates the radius of the void obtained from the convexhull (X = H) and Guinier approximations (X = G), and the gyrational radius of the voidsurface atoms (X = g)
Models 
R
_{g}

R
_{G}

R
_{H}

V
_{H}

Φ
_{S}

aSi (without hydrogen in voids)

CMDR300 
7.60 
6.38 
8.34 
2100 
0.66 
AIMDR300 
7.55 
8.07 
8.28 
2040 
0.67 
CMDR800 
7.56 
6.57 
8.31 
2080 
0.63 
AIMDR800 
6.84 
9.48 
8.67 
2250 
0.68 

aSi (with hydrogen in voids)

AIMDR300 
6.66 
8.66 
8.36 
2120 
0.56 
AIMDR800 
6.61 
10.01 
8.76 
2290 
0.71 
Having studied the microstructure of voids at 300 K, we now examine the structure of voids at the high temperature of 800 K. The results for the SAXS intensity obtained from the large model at 800 K are presented in Fig. 7. It is evident from Fig. 7 that the high temperature annealing at 800 K via AIMD simulations produces notable changes in the SAXS spectrum by reducing intensity values in the smallangle region. This reduction in the SAXS intensity can be understood by examining the shape of the individual voids, obtained from AIMDR models at 800 K. Fig. 8 shows the shape of a representative void (V3) in aSi after annealing at 800 K from classical and ab initio MD simulations. A comparison between the two structures in Fig. 8(a) and (b) shows the regions on the two surfaces with a varying degree of reconstruction. The shape of the respective void surfaces obtained from using the convexhull approximation is also shown in Fig. 8(c) and (d). Table 2 lists some of the parameters that characterize the void geometry in terms of the convexhull radius (R_{H}), the radius of gyration of the voidsurface atoms (R_{g}), the Guinier radius (R_{G}), the hull volume (V_{H}), and the average sphericity of the voids (Φ_{S}) after annealing at 300 K and 800 K. Although these scalar parameters do not vary much in the convex approximation, it is evident from the Gaussian isosurface representation of the void in Fig. 8(a) and (b), and the respective convex hulls in Fig. 8(c) and (d), that classical results noticeably differ from ab initio results, indicating the limitation of the classical force field in describing the restructuring of void surfaces at high temperature, namely at 800 K. This observation is also reflected in Fig. 9a–c, where the convexhull radii, the hull volumes, and the intensity of SAXS at 800 K, respectively, are presented.

 Fig. 7 The SAXS intensity of the large model of aSi from classical (CMDR) and ab initio (AIMDR) MD simulations of voids at 800 K, followed by geometry relaxation.  

 Fig. 8 An isosurface representation of a void (V3) in aSi obtained from: (a) classical; (b) ab initio simulations upon annealing at 800 K, followed by geometry relaxation. The corresponding shapes of the voids in the convexhull approximation are shown in (c) and (d), respectively.  

 Fig. 9 (a) The convexhull radius for the voids from classical CMDR () and AIMDR () simulations at 800 K, followed by geometry relaxation. (b) The corresponding volume of voids from the convexhull approximation. (c) Guinier plots for the models at 800 K. The Guinier radii from the plots correspond to a value of about 6.57 Å (CMDR) and 9.48 Å (AIMDR).  
B. Structural evolution of hydrogenrich voids in aSi
Smallangle Xray scattering measurements on hotwire chemically vapor deposited (HWCVD) films indicate that the presence of preexisting H clusters can affect the nanostructure of voids due to H_{2} bubble pressure in the cavities, which, upon annealing, can lead to geometric changes, as observed via tilting SAXS and surface transmission electron microscopy (STEM).^{15} Although the degree of such nanostructural changes may depend on the growth mechanism and the rate of film growth, it is instructive to address the problem from a computational viewpoint using CRN models of aSi with nanometersize voids, which are characterized by an experimentally consistent voidvolume fraction and clusters of H atoms inside the voids.
Fig. 10 shows the structure of a hydrogenrich nanovoid surface (of V3) obtained after annealing at 300 K and 800 K, followed by totalenergy minimization. The evolution of Si and H atoms during annealing was described using ab initio densityfunctional forces and energies, from the localbasis DFT package SIESTA.^{47} The procedure to conduct such AIMD simulations of voids, which are embedded in a large model of aSi, has been described in section IIC. A comparison of Fig. 10(a) (with H atoms inside the void at 300 K) with Fig. 4(b) (with no H atoms at 300 K) reveals that the thermalization of the void (V3) and its neighboring region at 300 K does not introduce notable changes on the void surface, with the exception of the distribution of H atoms. The initial random distribution of H atoms, which was confined within a radius of 4 Å from the void center, evolved to produce several monohydrides (SiH) and a few dihydrides (SiH_{2}), along with a few H_{2} molecules inside the cavity. The majority of H atoms continue to stay near the void at 300 K. By contrast, a considerable number of Si and H atoms have been found to be thermally driven out of the void region V3 at 800 K. An analysis of the distribution of Si atoms on the void surface, shown as yellow blob in Fig. 10, and H atoms (red) inside the void reveals that approximately 12.4% of the original Si atoms and 23.3% of the H atoms left the void region V3 of radius 10.5 Å after annealing at 800 K. This is apparent from the isosurface plot in Fig. 10(b), where a considerable number of Si atoms can be seen to move away from the voidsurface region, creating a somewhat diffused or scattered isosurface. This observation is found to be true for most of the other voids. On average, over 12 independent voids, approximately 8.3% of the Si atoms and 16.1% of the H atoms were found to leave the void regions at 800 K. By contrast, the corresponding average values for the same at 300 K for Si and H atoms are 0% and 2.6%, respectively. It thus follows that at 800 K the restructuring of a void surface can be considerably affected via thermal and Hinduced motion of Si atoms, as well as through the formation of various siliconhydrogen bonding configurations on void surfaces. These changes can be also observed in the convexhull approximation of the void surface (V3) in the presence and absence of H atoms at 800 K. The convex polyhedral region in Fig. 11b, with H atoms present inside the cavity, consists of 78 triangular faces and 41 vertices. These values are noticeably higher than the corresponding values of 62 faces and 32 vertices for the same cavity with no H atoms in Fig. 11a. The structural changes on the void surface have been reflected in the number of triangular faces, which are needed to construct the surface in the convexhull approximation. Table 3 presents the statistics of various siliconhydrogen bonding configurations in the vicinity of each void, which is defined by a radius of 10.5 Å.

 Fig. 10 The structure of a representative void (V3), in isosurface representation, from AIMD simulations in the presence of H atoms inside the cavity at: (a) 300 K; (b) 800 K. The yellow and red blobs represent Si and H atoms, respectively. For visualization and comparison on the same footing, an identical set of isosurface parameters was used to generate the surfaces.  

 Fig. 11 The shape of a void (V3) obtained from the convexhull approximation in AIMDR simulations at 800 K in the presence and absence of H atoms in the void. (a) The convex surface with no H atoms. (b) The convex surface with 30 H atoms inside the void.  
Table 3 Statistics of various siliconhydrogen bonding configurations near hydrogenrich voids at 300 K and 800 K from AIMD simulations. R_{H} and Φ_{S} indicate the convexhull radius and sphericity of the voids, whereas asterisked entries in column 1 indicate at least one H atom left the void region of radial size 10.5 Å
Void 
H_{iso} 
H_{2} 
SiH 
SiH_{2} 
SiH_{3} 
R
_{H}

Φ
_{S}

300 K

V1* 
2 
1 
20 
2 
0 
8.34 
0.51 
V2 
0 
3 
22 
1 
0 
8.33 
0.56 
V3 
0 
3 
17 
3 
0 
8.5 
0.59 
V4 
0 
0 
26 
1 
0 
8.4 
0.58 
V5 
0 
1 
25 
1 
0 
8.2 
0.62 
V6* 
1 
2 
22 
1 
0 
8.34 
0.54 
V7 
0 
1 
26 
1 
0 
8.31 
0.54 
V8 
1 
1 
23 
2 
0 
8.5 
0.57 
V9 
0 
3 
20 
0 
1 
8.35 
0.51 
V10* 
1 
3 
18 
2 
0 
8.28 
0.60 
V11 
2 
4 
19 
0 
0 
8.34 
0.59 
V12 
2 
2 
23 
0 
0 
8.52 
0.56 

800 K

V1 
0 
0 
20 
2 
0 
8.36 
0.90 
V2 
0 
2 
20 
2 
0 
8.84 
0.66 
V3* 
1 
2 
16 
1 
0 
9.07 
0.79 
V4 
0 
1 
23 
1 
0 
8.64 
0.62 
V5 
0 
0 
20 
2 
0 
8.95 
0.70 
V6 
0 
2 
14 
1 
1 
8.67 
0.83 
V7 
0 
3 
14 
3 
0 
8.94 
0.64 
V8* 
1 
1 
14 
3 
0 
8.95 
0.73 
V9 
0 
1 
22 
1 
0 
9.07 
0.62 
V10 
0 
1 
21 
2 
0 
8.71 
0.62 
V11 
2 
3 
15 
1 
0 
8.38 
0.81 
V12 
1 
0 
21 
2 
0 
8.56 
0.56 
To obtain the linear size of the voids, we computed the radius of gyration, R_{g}, using eqn (4). Likewise, the convexhull approximation of a void region, defined by a set of voidsurface atoms, provides the radius of the convex polyhedron, R_{H}, from eqn (5). Fig. 12 shows the values of R_{H} and the corresponding hull volumes, V_{H}, of 12 voids, each of which has been computed upon annealing and geometry relaxation. Similarly, the R_{g} values obtained from the realspace distribution of the void atoms are listed in Table 2. The apparent deviation of R_{H} from R_{g} can be attributed to the fact that the computation of R_{g} involves all the atoms on the surface and interior regions of the voids, whereas the convexhull approximation includes only those atoms on the voidsurface region that define the minimal convex polyhedral volume. Thus, R_{H} provides an upper bound of the radial size of the voids. In order to compare R_{H} and R_{g} values with the linear size of the voids from the simulated intensity plots in SAXS, we have invoked the Guinier approximation to estimate the Guinier radius, R_{G}, from lnI(k)–k^{2} plots, as shown in Fig. 12(c). Table 2 suggests that the R_{G} values obtained from the simulated intensity plots reasonably match with the corresponding average values of R_{H} at 300 K and 800 K. It may be noted that, while the intensity plot in Fig. 12(c) is not particularly sensitive to the shape of the voids, the slightly pronounced scattering in the region k^{2} ≤ 0.05 Å^{−2} is possibly indicative of the expansion of the void volume at 800 K, as observed in the convexhull volume in Fig. 12(b). The exact cause of this increased scattering intensity is difficult to determine but it appears that a combination of the pressure due to H_{2} bubbles, surface roughening of the voids, and the displacement of the atoms on the void surfaces at high temperature could lead to this additional scattering. This observation is corroborated by the experimental results from a smallangle Xray scattering study of nanovoids in HWCVD amorphous Si films that indicated an increase of individual void volumes upon annealing at 813 K.^{15}

 Fig. 12 (a) The radius of gyration of hydrogenrich voids from AIMD simulations at 300 K () and 800 K () in the convexhull approximation. (b) Estimated convexhull volumes at 300 K and 800 K. (c) Guinier fits for the SAXS intensities in AIMDR models of aSi:H at 300 K and 800 K.  
C. Hydrogen dynamics near voids in aSi
In this section, we shall make a few observations on the dynamics of hydrogen atoms in the vicinity of void surfaces from ab initio moleculardynamics simulations. Since the mass of H atoms is significantly smaller than that of Si atoms, the motion of H atoms is more susceptible to the annealing temperature and, more importantly, to the resulting temperatureinduced structural changes on the void surfaces. This is particularly noticeable at high temperature. Having established earlier that the structural changes on the void surfaces are more pronounced at 800 K, we shall now examine the character of hydrogen dynamics inside (and near the surface of) the voids at 800 K. Despite the limited timescale of the present AIMD simulations, spanning 10 ps only, it is reasonable to assume that, at high temperature, the dynamical character of hydrogen motion inside the voids would be reflected in the first several picoseconds. It may be noted that, unlike structural properties, the dynamics of hydrogen in aSi may depend on the accuracy of total energies and forces from DFT calculations. The use of the Harrisfunctional approximation in the present work, in an effort to address the structure of voids on the nanometer length scale, suggests that the results on hydrogen dynamics presented here are somewhat approximate in nature. A full selfconsistentfield DFT calculation, over a period of a few tens of picoseconds, using planewave or even local basis sets for 2200atom models is computationally overkill and outside scope of the present work on nanostructural evolution.
Fig. 13(a) depicts the time evolution of hydrogen atoms at 800 K within two voids, V12 and V7. Initially, all hydrogen atoms/molecules were uniformly distributed inside the voids within a radius of 4 Å from the center of the voids (see Fig. 1). The meansquare displacements (MSD), 〈R^{2}(t)〉, of hydrogen atoms – averaged over all H atoms in the respective void – are shown in Fig. 13(a) from 0 to 10 ps. The corresponding realspace distributions of Si and H atoms on the void surfaces and their immediate vicinity are also shown in Fig. 13(b) and (c) for V12 and V7, respectively. The earlytime behavior of the dynamics, from 0 to 0.5 ps, can be roughly understood from elementary considerations. Given that the majority of H atoms are away from the void surface at the beginning of simulation, by at least 2 Å or more, the earlytime behavior of the dynamics would be primarily determined by the temperature of the system for a reasonably welldefined initial H distribution. Since the rootmeansquare (RMS) speed of H atoms at 800 K for the MaxwellBoltzmann distribution is about 0.045 Å fs^{−1}, it would be reasonable to assume that the initial behavior can last for about a fraction of a picosecond, after which H atoms become affected by the presence of the void surface. It is the marked difference in the value of 〈R^{2}(t)〉 in Fig. 13(a), from 1 ps to 8 ps, that we find most interesting. This notable difference between the two sets of MSD values beyond 0.5 ps can be understood by taking into account the realspace structure of the void surfaces. In addition to the trivial thermal and configurational fluctuations – caused by temperature and the disorder in H distributions, respectively – the motion of H atoms is also affected by the varying degree of restructuring of Si atoms and the disorder of the void surfaces after 0.5 ps. The surface structure of V7 can be seen to be more open or scattered compared to its counterpart in V12, which is more compact or welldefined, as shown in Fig. 13(c) and (b), respectively. This has been also confirmed by analyzing the distributions of voidsurface atoms in the interior of V12 and V7. The presence of Si atoms inside V7 may reduce the effective diffusion rate of H atoms, due to additional scattering from interior Si atoms, during intermediatetime evolution from 1 ps onward. By contrast, the more welldefined void surface of V12, with few Si atoms inside, provides less resistance to H atoms during their intermediatetime evolution. Thus, hydrogen diffusion in V12 proceeds relatively uninterruptedly, in the presence of few Si atoms from void surfaces, leading to a notable difference in their MSD values during the first several picoseconds of simulations. This observation has been found to be true for other pairs of voids as well, for which a notable difference in the MSD values exists. Additional results on the time evolution of H atoms in voids V5, V9 and V11, which exhibit a similar MSD behavior, are provided as ESI.† In view of our observation that the intermediatetime behavior (i.e., from 1 ps to 7–8 ps) of the H dynamics can be influenced by the roughness of void surfaces, it might be necessary to take into account the appropriate size and the accurate surface structure of the voids in the calculations. However, a full SCF study of hydrogen dynamics inside nanometersize voids in aSi for several tens of picoseconds is computationally prohibitive and outside the scope of the present study.

 Fig. 13 (a) The time evolution of the meansquare displacements (MSD) of H inside two voids, V12 () and V7 (), at 800 K. The difference between the two sets of MSD values, from 1 ps to 8 ps, can be attributed to the degree of voidsurface restructuring, as discussed in section IIIC. (b) The compact or smooth structure of the V12 surface, obtained from an isosurface representation of the void, compared to a relatively diffused or scattered distribution of Si atoms in (c) V7 at 800 K.  
D. Siliconhydrogen bonding configurations on void surfaces
Experimental data from positronannihilation lifetime (PAL) spectroscopy,^{13} Rutherford backscattering spectrometry (RBS),^{59} and Fouriertransform infrared spectroscopyattenuated total reflection (FTIRATR)^{59} indicate that bonded and nonbonded hydrogens play an important role in characterizing the structural and optical properties of aSi:H. Sekimoto et al.^{13,59} recently demonstrated that the presence of nonbonded hydrogens (NBH) at concentrations beyond 2.8 at% in amorphous silicon networks led to changes in the vacancysize distribution and induced formation of nanometersize voids in the network via relaxation of internal stress. The formation of H_{2} molecules and their evolution can be probed by an analysis of the lowtemperature (LT) peak of the hydrogeneffusion profile^{13,60,61} near 750 K, whereas the number density of H_{2} molecules can be estimated from the collisioninduced weak IR absorption in the frequency region of 4100 cm^{−1} to 5500 cm^{−1}. Chabal and Patel^{24} reported a value of 10^{21} H_{2} per cm^{3} by analyzing data obtained from IR measurements at room temperature. Although the present study, based on pure aSi networks, does not permit us to address the crucial role of siliconhydrogen bonding configurations in aSi:H in general, it does provide new insights into the role of bonded and nonbonded hydrogens in the vicinity of the walls of Hrich voids. Tables 3 and 4 list the statistics and the average concentration of various bonded and nonbonded hydride configurations inside the voids, respectively. Nonbonded hydrogens chiefly comprise H_{2} molecules, along with one or two isolated H atoms near the voids. The latter possibly originate from the short annealing time of 10 ps of AIMD simulations. The number of H_{2} molecules per void ranges from 0 to 4 and 0 to 3 at 300 K and 800 K, respectively, which translates into an average density of 1.5–2.2 × 10^{21} H_{2} per cm^{3} for nanometersize voids. This value matches very well with the experimental value of 10^{21} H_{2} per cm^{3} from IR measurements by Chabal and Patel,^{24} as mentioned earlier. These results are also consistent with the reported data from RBS measurements.^{13} On the other hand, the great majority of bonded hydrogens appeared in monohydride (SiH) configurations. An examination of the local bonding environment of silicon atoms near the void surfaces showed that approximately 72.5 at% and 61.1 at% of total H atoms formed Si–H bonds at 300 K and 800 K, respectively (see Table 4). Of the remaining H atoms, about 8.3% (11.7%) were found to be bonded with Si as SiH_{2} configurations at 300 K (800 K). Apart from a few isolated H atoms, the rest of the hydrogen has left the spherical void region of radius 10.5 Å. A graphic distribution of the presence of H_{2} molecules and the bonded SiH/SiH_{2} configurations on the wall of V3 at 300 K is depicted in Fig. 14. Several SiH and a few SiH_{2} configurations can be seen to have formed on the surface of the void (V3). This observation is consistent with experimental results from IR measurements^{24} and the recent computational studies based on informationdriven atomistic modeling of aSi:H.^{17,21}Fig. 15 shows the concentration of hydrogen associated with various siliconhydrogen bonding configurations. The results clearly establish that a considerable number of H_{2} molecules, of about 9–13% of the total H, can form in the vicinity of voids, depending upon the annealing temperature, concentration of H atoms, the source of hydrogen, and the preparation methods of aSi samples.

 Fig. 14 Bonded and nonbonded hydrogens inside a void (V3) from AIMD simulations at 300 K. (a) The interior wall of the void can be seen to be decorated with several SiH (yellowred) and three SiH_{2} (bluered) configurations, as well as three H_{2} molecules (white). (b) A magnified view, showing protruding monohydrides (yellowred) and dihydrides (bluered) on the walls, and a hydrogen molecule (white) inside the void.  

 Fig. 15 Average statistics of siliconhydrogen bonding configurations near void surfaces in aSi:H after annealing at 300 K and 800 K. The results were obtained by averaging over 12 independent voids, distributed over a distance of several nanometers in the large aSi model.  
Table 4 The average concentration of hydrogen (in % of total H atoms) in siliconhydrogen bonding configurations, H_{2} molecules and isolated H near voids (within a radius of 10.5 Å). Bonded and nonbonded hydrogen that left the voids during annealing are listed under the category Exvoid
Temp. 
H_{iso} 
H_{2} 
SiH 
SiH_{2} 
SiH_{3} 
Exvoid 
300 K 
2.5 
13.3 
72.5 
8.3 
0.8 
2.6 
800 K 
1.4 
8.9 
61.1 
11.7 
0.8 
16.1 
We conclude this section with the following observation on siliconhydrogen bond formation during AIMD simulations. The results from the preceding paragraph suggest that the internal surfaces of hydrogenrich voids in aSi are principally decorated with monohydride Si–H bonds. It is thus instructive to examine the time evolution of the number of Si–H bonds during annealing via MD simulations at low and high temperatures. Since a statistically significant number of H atoms participate in forming Si–H bonds on void surfaces, we restrict ourselves to address the formation of monohydride Si–H configurations only. Fig. 16 shows the evolution of the average hydrogen content in SiH bonds in the vicinity of voids at 300 K (blue) and 800 K (red). The results were obtained by sampling atomic configurations (of voids) at an interval of every 10 fs during the course of AIMD simulations, for a period of 7–8 ps, which were then averaged over 12 individual voids. It is apparent from Fig. 16 that the formation of monohydride Si–H bonds takes place very rapidly at 800 K, which can saturate the void surfaces within the first two picoseconds. By contrast, it takes several picoseconds for hydrogen atoms inside the cavities to passivate the void surfaces at 300 K. A comparison of the saturated values of the SiH bonds in Fig. 16, obtained from the MDannealed models at 300 K and 800 K, with those from the corresponding final AIMDR models in Table 4 suggests that the results are consistent with each other. A small difference between the saturated values of the number of SiH bonds at 300 K and 800 K (in Fig. 16) with those in Table 4 can be attributed to the totalenergy relaxation of the final MDannealed models, which were used in the calculation of siliconhydrogen bonding statistics in Table 4.

 Fig. 16 The time evolution of monohydride SiH bonds during AIMD simulations at 300 K (blue) and 800 K (red). The vertical axis indicates the amount of hydrogen (in percent of total H atoms) that contribute to SiH bond formation. The results were obtained via configurational averaging over 12 voids.  
IV. Conclusions
In this paper, we studied the temperatureinduced nanostructural changes of the voids in aSi at low and high annealing temperature in the range 300–800 K using classical and quantummechanical simulations, in the absence and presence of hydrogen near voids. This was achieved by generating a large model of aSi, consisting of more than 262000 atoms, in order to be able to produce a realistic distribution of nanometersize voids in the amorphous environment of Si atoms, as observed in SAXS, PAL, IR, and hydrogeneffusion measurements. An examination of the distribution of the atoms near the voids reveals that the reconstruction of void surfaces in pure aSi at 300 K led to minimal changes in the shape and size of the voids, which can be readily understood from a classical treatment of the problem. By contrast, the hightemperature annealing at 800 K caused significant changes in the shape and size of the voids, with noticeable structural differences between classical and quantummechanical results. This observation appears to indicate that classical potentials, such as the modified SW potential, might not be particularly suitable for an accurate description of the dynamics of Si atoms near voids and the resulting reconstruction of the void surfaces at high temperatures.
An important outcome of the present study is that the dynamics of bonded and nonbonded hydrogens near the voids and the degree of voidsurface reconstruction are found to be intrinsically related to each other. Ab initio annealing of aSi networks with a voidvolume fraction of 0.2% at 300 K and 800 K suggests that the presence of hydrogen within voids can facilitate surface reconstruction through the formation of silicon monohydride and dihydride configurations, as well as the displacement of Si atoms on the void surfaces. The rearrangement of atoms on void surfaces affects the diffusion of nonbonded hydrogens, which in turn produce surface bumps and changes the shape and size of the original void. This observation is reflected in the time evolution of hydrogen during annealing. The presence of Si atoms inside a heavily reconstructed rough void surface can reduce the effective diffusion rate of hydrogen, due to additional scattering from the interior Si atoms, in their intermediatetime evolution from 1 ps onward. In contrast, a compact void surface provides less resistance to hydrogen during their evolution. Thus, hydrogen diffusion inside a smooth or welldefined void proceeds with little or no hindrance, leading to a notable difference in their meansquaredisplacement (MSD) values during their intermediatetime evolution. Finally, the results from the study show the presence of bonded hydrogens (BH), mostly SiH, SiH_{2}, and nonbonded hydrogens (NBH) in the form of H_{2} molecules. The concentration of BH and NBH found in our study is consistent with the experimental values obtained from infrared (IR), Rutherford back scattering (RBS), and hydrogenforward scattering (HFS) measurements. An examination of the siliconhydrogen bonding configurations near the voids suggests that the interior walls of the voids are decorated with SiH_{2}, which is supported by experimental results from infrared and ellipsometric studies. Likewise, the concentration of H_{2} molecules obtained from the firstprinciples densityfunctional simulations in the present study is found to be consistent with the experimental value estimated from the collisioninduced weak IR absorption by H_{2} molecules in the frequency region of 4100–5500 cm^{−1}. A preliminary study using more accurate DFT calculations on 1000atom models, via the SCF solutions of the Kohn–Sham equation in the generalizedgradient approximation, with 20–40 H atoms inside the voids indicates that more SiH_{2} configurations tend to form at high concentration of hydrogen in the cavities. However, further studies are needed for an accurate determination of SiH, SiH_{2} and H_{2} molecules with varying concentrations of H atoms inside the voids.
Conflicts of interest
There are no conflicts of interest to declare.
Acknowledgements
This work was partially supported by the U.S. National Science Foundation under Grants No. DMR 1507166 and No. DMR 1507118. We acknowledge the Texas Advanced Computing Center at the University of Texas at Austin for providing HPC resources that have contributed to the results reported in this work. One of us (P. B.) thanks Prof. David Drabold (Ohio University) and Mr Dil Limbu (USM) for discussions.
References
 M. J. Thompson and H. C. Tuan, Amorphous Silicon Electronic Devices and their Applications, IEEE Aero. Electron. Sys. Mag., 1987, 2, 25–29 Search PubMed.

R. A. Street, Technology and Applications of Amorphous Silicon, SpringerVerlag, 2000, vol. 37 Search PubMed.
 S. De Wolf, A. Descoeudres, Z. Holman and C. Ballif, Highefficiency silicon heterojunction solar cells: A review, Green, 2012, 2, 7–24 CAS.
 A. Kanevce and W. K. Metzger, The role of amorphous silicon and tunneling in heterojunction with intrinsic thin layer (HIT) solar cells, J. Appl. Phys., 2009, 105, 094507 CrossRef.

M. Zeman and D. Zhang, Heterojunction Silicon Based Solar Cells, Springer Berlin Heidelberg, 2012, pp. 13–43 Search PubMed.
 J. Ge, Z. P. Ling, J. Wong, T. Mueller and A. Aberle, Optimisation of Intrinsic aSi:H Passivation Layers in Crystallineamorphous Silicon Heterojunction Solar Cells, Energy Procedia, 2012, 15, 107–117 CrossRef CAS.
 H. Fujiwara, M. Kondo and A. Matsuda, Realtime spectroscopic ellipsometry studies of the nucleation and grain growth processes in microcrystalline silicon thin films, Phys. Rev. B: Condens. Matter Mater. Phys., 2001, 63, 115306 CrossRef.
 R. Collins, A. Ferlauto, G. Ferreira, C. Chen, J. Koh, R. Koval, Y. Lee, J. Pearce and C. Wronski, Evolution of microstructure and phase in amorphous, protocrystalline, and microcrystalline silicon studied by real time spectroscopic ellipsometry, Sol. Energy Mater. Sol. Cells, 2003, 78, 143–180 CrossRef CAS.
 H. Daxing, J. D. Lorentzen, J. WeinbergWolf, L. E. McNeil and Q. Wang, Raman study of thin films of amorphoustomicrocrystalline silicon prepared by hotwire chemical vapor deposition, J. Appl. Phys., 2003, 94, 2930–2936 CrossRef.
 B. Macco, J. Melskens, N. Podraza, K. Arts, C. Pugh, O. Thomas and W. Kessels, Correlating the silicon surface passivation to the nanostructure of lowtemperature aSi:H after rapid thermal annealing, J. Appl. Phys., 2017, 122, 035302 CrossRef.

D. L. Williamson, S. J. Jones and Y. Chen, SmallAngle xray Scattering Studies of Microvoids in AmorphousSiliconBased Semiconductors, 1994, p. 6395, NREL/TP4516395.UC Category 271. DE94006931 Search PubMed.
 M. Serényi, C. Frigeri, Z. Szekrényes, K. Kamarás, L. Nasi, A. Csik and N. Q. Khánh, On the formation of blisters in annealed hydrogenated aSi layers, Nanoscale Res. Lett., 2013, 8, 84 CrossRef.
 T. Sekimoto, M. Matsumoto, A. Sagara, M. Hishida and A. Terakawa, Changes in the vacancy size distribution induced by nonbonded hydrogens in hydrogenated amorphous silicon, J. NonCryst. Solids, 2016, 447, 207–211 CrossRef CAS.
 D. T. Britton, A. Hempel, M. Härting, G. Kögel, P. Sperr, W. Triftshäuser, C. Arendse and D. Knoesen, Annealing and recrystallization of hydrogenated amorphous silicon, Phys. Rev. B: Condens. Matter Mater. Phys., 2001, 64, 075403 CrossRef.
 D. L. Young, P. Stradins, Y. Xu, L. M. Gedvilas, E. Iwaniczko, Y. Yan, H. M. Branz, Q. Wang and D. L. Williamson, Nanostructure evolution in hydrogenated amorphous silicon during hydrogen effusion and crystallization, Appl. Phys. Lett., 2007, 90, 081923 CrossRef.
 D. L. Williamson, A. H. Mahan, B. P. Nelson and R. S. Crandall, The observation of microvoids in device quality hydrogenated amorphous silicon, J. NonCryst. Solids, 1989, 114, 226–228 CrossRef CAS.
 P. Biswas and S. R. Elliott, Nanoscale structure of microvoids in aSi:H: a firstprinciples study, J. Phys.: Condens. Matter, 2015, 27, 435201 CrossRef.
 P. Biswas, D. Paudel, R. AttaFynn, D. A. Drabold and S. R. Elliott, Morphology and Number Density of Voids in Hydrogenated Amorphous Silicon: An Ab Initio Study, Phys. Rev. Appl., 2017, 7, 024013 CrossRef.
 R. Biswas, I. Kwon, A. M. Bouchard, C. M. Soukoulis and G. S. Grest, Intense small wavevector scattering from voids in amorphous silicon: A theoretical simulation, Phys. Rev. B: Condens. Matter Mater. Phys., 1989, 39, 5101–5106 CrossRef CAS.
 R. Brahim and A. Chehaidar, SmallAngle XRay Scattering of Amorphous Germanium: Numerical Modeling, Adv. Mater. Phys. Chem., 2013, 3, 19–30 CrossRef.
 P. Biswas, D. A. Drabold and R. AttaFynn, Microstructure from joint analysis of experimental data and ab initio interactions: Hydrogenated amorphous silicon, J. Appl. Phys., 2014, 116, 244305 CrossRef.
 S. Muramatsu, S. Matsubara, T. Watanabe, T. Shimada, T. Kamiyama and K. Suzuki, Smallangle Xray scattering studies of aSi_{1−x}Ge_{x}: H alloys, J. NonCryst. Solids, 1992, 150, 163–166 CrossRef CAS.
 J. B. Boyce and M. Stutzmann, Orientational Ordering and Melting of Molecular H_{2} in an aSi Matrix: NMR Studies, Phys. Rev. Lett., 1985, 54, 562–565 CrossRef CAS PubMed.
 Y. J. Chabal and C. K. N. Patel, Infrared Absorption in $a$Si: H: First Observation of Gaseous Molecular ${\mathrm{H}}_{2}$ and SiH Overtone, Phys. Rev. Lett., 1984, 53, 210–213 CrossRef CAS.
 A. von Keudell and J. R. Abelson, Evidence for atomic H insertion into strained SiSi bonds in the amorphous hydrogenated silicon subsurface from in situ infrared spectroscopy, Appl. Phys. Lett., 1997, 71, 3832–3834 CrossRef CAS.
 Y. J. He, M. Hasegawa, R. Lee, S. Berko, D. Adler and A.L. Jung, Positronannihilation study of voids in aSi and aSi:H, Phys. Rev. B: Condens. Matter Mater. Phys., 1986, 33, 5924–5927 CrossRef CAS PubMed.
 J. E. Graebner, B. Golding, L. C. Allen, D. K. Biegelsen and M. Stutzmann, Solid Hydrogen in Hydrogenated Amorphous Silicon, Phys. Rev. Lett., 1984, 52, 553–556 CrossRef CAS.
 J. Melskens, S. W. H. Eijt, M. Schouten, A. S. Vullers, A. Mannheim, H. Schut, B. Macco, M. Zeman and A. H. M. Smets, Migration of Open Volume Deficiencies in Hydrogenated Amorphous Silicon During Annealing, IEEE J. Photovoltaics, 2017, 7, 421–429 Search PubMed.
 R. L. C. Vink, G. T. Barkema, W. F. van der Weg and N. Mousseau, Fitting the StillingerWeber potential to amorphous silicon, J. NonCryst. Solids, 2001, 282, 248–255 CrossRef CAS.
 F. H. Stillinger and T. A. Weber, Computer simulation of local order in condensed phases of silicon, Phys. Rev. B: Condens. Matter Mater. Phys., 1985, 31, 5262–5271 CrossRef CAS.
 J. S. Custer, M. O. Thompson, D. C. Jacobson, J. M. Poate, S. Roorda, W. C. Sinke and F. Spaepen, Density of amorphous Si, Appl. Phys. Lett., 1994, 64, 437–439 CrossRef CAS.
 S. Nosé, A unified formulation of the constant temperature molecular dynamics methods, J. Chem. Phys., 1984, 81, 511–519 CrossRef.
 W. G. Hoover, Canonical dynamics: Equilibrium phasespace distributions, Phys. Rev. A, 1985, 31, 1695–1697 CrossRef.
 R. AttaFynn and P. Biswas, Nearly defectfree dynamical models of disordered solids: The case of amorphous silicon, J. Chem. Phys., 2018, 148, 204503 CrossRef.
 P. D'Antonio and J. H. Konnert, SmallAngleScattering Evidence of Voids in Hydrogenated Amorphous Silicon, Phys. Rev. Lett., 1979, 43, 1161–1163 CrossRef.
 A. C. Wright, A. C. Hannon, R. N. Sinclair, T. M. Brunier, C. A. Guy, R. J. Stewart, M. B. Strobel and F. Jansen, Neutron scattering studies of hydrogenated, deuterated and fluorinated amorphous silicon, J. Phys.: Condens. Matter, 2007, 19, 415109 CrossRef.
 M. S. Conradi and R. E. Norberg, Molecular h_{2}: Nuclearspinrelaxation centers for protons in aSi: H, Phys. Rev. B: Condens. Matter Mater. Phys., 1981, 24, 2285–2288 CrossRef CAS.
 H. V. Löhneysen, H. J. Schink and W. Beyer, Direct Experimental Evidence for Molecular Hydrogen in Amorphous Si: H, Phys. Rev. Lett., 1984, 52, 549–552 CrossRef.
 W. Kohn, Density Functional and Density Matrix Method Scaling Linearly with the Number of Atoms, Phys. Rev. Lett., 1996, 76, 3168–3171 CrossRef CAS.

H. Ehrenreich, F. Seitz and D. Turnbull, Solid State Physics, Elsevier Science, 1980 Search PubMed.
 R. Haydock, V. Heine and M. Kelly, Electronic structure based on the local atomic environment for tightbinding bands. II, J. Phys. C: Solid State Phys., 1975, 8, 2591–2605 CrossRef.
 Strictly speaking, the boundary must be located at a distance greater than the characteristic decay length of the density matrix of the system. The presence of an electronic gap in aSi asserts that its density matrix 62 can be characterized by a welldefined decay length.
 The rationale behind adding 30 H atoms inside the voids is based on calorimetry and infrared measurements, 24 which suggest the presence of about 12 to 15 “bulk” H_{2} molecules inside nanometersize voids.
 J. Harris, Simplified method for calculating the energy of weakly interacting fragments, Phys. Rev. B: Condens. Matter Mater. Phys., 1985, 31, 1770–1779 CrossRef CAS.
 E. Zaremba, Extremal properties of the Harris energy functional, J. Phys.: Condens. Matter, 1990, 2, 2479–2486 CrossRef.
 R. AttaFynn, P. Biswas, P. Ordejón and D. A. Drabold, Systematic study of electron localization in an amorphous semiconductor, Phys. Rev. B: Condens. Matter Mater. Phys., 2004, 69, 085207 CrossRef.
 J. M. Soler, E. Artacho, J. D. Gale, A. Garciá, J. Junquera, P. Ordejón and D. SánchezPortal, The SIESTA method for ab initio order N materials simulation, J. Phys.: Condens. Matter, 2002, 14, 2745 CrossRef CAS.
 N. Troullier and J. L. Martins, Efficient pseudopotentials for planewave calculations, Phys. Rev. B: Condens. Matter Mater. Phys., 1991, 43, 1993–2006 CrossRef CAS.
 L. Kleinman and D. M. Bylander, Efficacious Form for Model Pseudopotentials, Phys. Rev. Lett., 1982, 48, 1425–1428 CrossRef CAS.
 J. P. Perdew and A. Zunger, Selfinteraction correction to densityfunctional approximations for manyelectron systems, Phys. Rev. B: Condens. Matter Mater. Phys., 1981, 23, 5048–5079 CrossRef CAS.

N. E. Cusack, The Physics of Structurally Disordered Matter, Adam Hilger, 1987 Search PubMed.

A. Guinier, G. Fournet and K. L. Yudowitch, Smallangle scattering of Xrays, Wiley, New York, 1995 Search PubMed.
 D. Paudel, R. AttaFynn, D. A. Drabold, S. R. Elliott and P. Biswas, Smallangle Xray scattering in amorphous silicon: A computational study, Phys. Rev. B, 2018, 97, 184202 CrossRef CAS.
 H. Wadell, Volume, Shape, and Roundness of Quartz Particles, J. Geol., 1935, 43, 250–280 CrossRef.

O. G. Porod, Small Angle Xray Scattering, ed. O. Glatter and O. Kratky, Academic Press, New York, 1982, pp. 17–51 Search PubMed.
 K. Laaziri, S. Kycia, S. Roorda, M. Chicoine, J. L. Robertson, J. Wang and S. C. Moss, Highenergy xray diffraction study of pure amorphous silicon, Phys. Rev. B: Condens. Matter Mater. Phys., 1999, 60, 13520–13533 CrossRef CAS.
 The observation rests on the assumptions that the variation of the local mass density near the defect sites and the corresponding length scale of 2–3 Å, associated with a sparse distribution of a few isolated defects, would not be too large to induce additional scattering in the smallangle region of the scattering wavevector. Suffice it to mention that, for clustered defects, this may not be necessarily true. These assumptions have been verified computationally by creating artificial defects in 100% 4foldcoordinated aSi networks here.
 A. Kokalj, XCrySDen¿a new program for displaying crystalline structures and electron densities, J. Mol. Graphics Modell., 1999, 17, 176–179 CrossRef CAS.
 T. Sekimoto, M. Matsumoto and A. Terakawa, Dense restructuring of amorphous silicon network induced by nonbonded hydrogen, Jpn. J. Appl. Phys., 2018, 57, 08RB07 CrossRef.
 W. Beyer, Hydrogen effusion: a probe for surface desorption
and diffusion, Phys. B, 1991, 170, 105–114 CrossRef CAS.
 W. Beyer, Diffusion and evolution of hydrogen in hydrogenated amorphous and microcrystalline silicon, Sol. Energy Mater. Sol. Cells, 2003, 78, 235–267 CrossRef CAS.
 S. Goedecker, Linear scaling electronic structure methods, Rev. Mod. Phys., 1999, 71, 1085–1123 CrossRef CAS.
Footnote 
† Electronic supplementary information (ESI) available. See DOI: 10.1039/C9NR08209C 

This journal is © The Royal Society of Chemistry 2020 