PorosityPlus: characterisation of defective, nanoporous and amorphous materials

As both porous and amorphous semiconductors have different advantages the challenge becomes knowing how to select one over the other, and knowing how to anticipate the degree of crystallinity or the fraction of voids as a function of a controllable feature such as the density. These sorts of relationships can be modelled computationally but unambiguous characterisation of the porosity of complex and disordered structures requires specialist tools. In this paper we demonstrate the use of PorosityPlus to investigate porosity in the vacancy induced amorphisation of defective crystals of carbon, silicon and germanium. The PorosityPlus software allows for the identification of vacancies, twin vacancies and larger pores, along with their relative locations and their respective populations. Void migration and coalescence along with the associated density changes can also be calculated. We show that, with increasing initial vacancies (reduced density) carbon (diamond), silicon and germanium exhibit characteristic density-dependent porosity profiles, coupled with simultaneous amorphisation. This is an ideal tool for integration into advanced computational workflows, such as creating fingerprints for topological data analysis or machine learning, since the porosity profile for each configuration is unique.

Porosity is defined as the fraction of void space and can be determined easily by weight measurements [14], but results are strongly correlated to the density which can be influenced by other factors. Amorphisation is a competing process that can also explain changes in density, or occur simultaneously to produce amorphous materials with pores. This is not necessarily a bad thing. Amorphous carbon and hydrogenated amorphous carbon exhibit some advantageous physical and chemical properties [28], such as high mechanical hardness, high optical transparency in the visible and near infrared, high thermal conductivity, low friction coefficient, and chemical inertness to some corrosive agents [29,30], making them suitable as protective tribological coating, radiation protection, electron field emitters, and in coatings in biomedicine [31][32][33][34][35]. Amorphous silicon and hydrogenated amorphous silicon are important in a range of applications from thin-film transistors, active matrix displays, image-sensor arrays, multi-junction solar cells to multilayer colour and thin-film position detectors [36]. Amorphous germanium exhibits interesting optical properties [37][38][39] useful for high-quality flexible displays and solar cells [40]. Photovoltaic applications often employ hydrogenated amorphous silicon [41,42] and while hydrogenation is not included in this work, it may be investigated in future applications.
As both porous and amorphous semiconductors have different advantages the challenge becomes knowing how to select one over the other, and knowing how to anticipate the degree of crystallinity or the fraction of voids as a function of a controllable feature such the density. These sorts of relationships can be modelled computationally but unambiguous characterisation of the porosity of complex and disordered structures requires specialist tools. In this paper we model the vacancy induced formation of porous carbon, silicon and germanium, while simultaneously considering the issue of amorphisation. This is done with new software that characterizes the pore geometry within an atomic configuration by calculating useful properties such as the porosity, internal surface area and pore size distribution (PSD) using a variety of Monte Carlo (MC) based integration methods. The results for carbon (diamond), silicon and germanium reveal characteristic densitydependent porosity, coupled with simultaneous amorphisation. At low density each material eventually undergoes complete amorphisation, with carbon retaining the greatest porosity.

Methods
Initial diamond lattices for carbon, silicon and germanium were prepared with the simulation cell parameters shown in table 1.
The density of diamond-structured carbon, silicon and was varied by the removal of random atoms. At each density ten MC simulations were performed to obtain more accurate statistics for testing the thermal stability of the defective lattices. Each simulation ran for 100 million attempted steps with a different random seed to ensure that different atoms are removed from the initially perfect lattice. The defective lattices were simulated at a temperature below the experimental melting points, beginning at 3500 K, 1500 K, 1100 K for carbon, silicon and germanium, respectively. For 25 million steps the simulations were annealed at the starting temperatures before being quenched linearly to 300 K over 50 million steps. Finally, the system was annealed at 300 K for 25 million steps. This scheme is illustrated in figure 1. The environment dependent interatomic potential was used for carbon [43] and silicon [44] while a Stillinger-Weber potential was used for germanium [45]. The MC simulations were performed using the hybrid reverse Monte Carlo method described elsewhere [46][47][48].
Porosity in itself is a macroscopic parameter that does not yield any information regarding the microstructure of the material, and so it is conventionally described in terms of the pore size and PSD within the sample. Given these parameters porous materials are divided into three categories based on the size of its pores: macroporous, mesoporous, and microporous. In the present study the pore size and PSD was probed via a MC integration sampling approach outlined by Gelb et al [49], using our software PorosityPlus [50]. PorosityPlus allows for the characterization of pore space using numerous MC technicals to determine quantities such as porosity, PSD and internal surface areas. To determine the porosity a probe particle of some radius is inserted into the cell at a random position and is checked to ensure that it does not overlap any atoms. If it does not, the insertion is accepted and the fraction of accepted moves to total attempted moves is the value of porosity or volume fraction. While Porosityplus is here applied to the study of atomic systems, being a gemetrical characterization technique, it can also be used to characterize larger scale systems such as spherical nanoparticle aggregates and mechanical particle stacking. The calculation of the PSD also begins with random trial point insertions. If there is no overlap, a determination of the maximum sphere size which both encompasses the original insertion point at the centre of the sphere and which does not overlap any atoms is found. The software finds this maximum sphere diameter by a random walk away from the original insertion point and once found, the point is binned into the sphere diameter bin at the maximum discovered value. This is repeated over many insertion points to determine the PSD. The atomic radius for the carbon, silicon and germanium atoms was set to 0.075 nm, 0.11 nm and 0.12 nm respectively. The radius for the probe particle was set at 0.08 nm for the carbon, 0.13 nm for the silicon and 1.3 nm for the germanium models. The atomic radii are approximately the experimental covalent radii while the probe particle radii were minimized to the chosen values with the requirement that in a perfect lattice, no insertions are possible. Thus, the choice of the probe radius is selected to a value just large enough to avoid detection of spacing between atoms in a perfect lattice. Since we track the porosity change, defined as the difference in the porosities of the pre and post thermal ramping configurations, as a function of density for a given material, the conclusions are not dependent upon the exact values of these parameters.
The level of amorphisation of a system is quantified by looking at the height of the third peak in the radial distribution function, g(r). This peak is seen to disappear as the system moves from a fully crystalline diamond lattice into an amorphous state. A similar criterion was used by Fairchild et al [51].
A flowchart of the computational workflow is provided in figure 2.

Results
The porosity was characterized using the pore volume and the PSD for each density, averaged over the 10 simulations to acquire the statistics. Examples of final results at different defect densities are provided in figure 3 for carbon (green), silicon (red) and germanium (blue); these colours are used consistently for these elements throughout this article. In these examples the initial and final configurations for each material at different densities are shown; the density for carbon 2.60 g cm −3 , for silicon 2.14 g cm −3 , and for germanium 4.55 g cm −3 . These figures are an illustration of the maximum spheres detected by PorosityPlus, and are therefore a visualization of the spheres in the histogram making up the PSD, producing the arbitrary pore geometry. In each case the pores are superimposed with the atomic structure, shaded to illustrate the material and pore diameter; pores > 0.4 nm exhibit the brightest shading, 0.35-0.40 nm pores exhibit medium/bright shading, 0.30-0.35 nm pores exhibit medium/dark shading, and 0.25-0.30 nm pores exhibit the darkest shading. Note that carbon does not show single missing atoms in the initial structure due to the carbon lattice parameters being smaller than silicon and germanium. The visualisation methodology used here is described in the manual for PorosityPlus.   Since randomly removing atoms (introducing vacancies) in itself changes the density, this characterisation procedure was applied to both the initial and final states and the change in the pore volume and the PSD was recorded. This effectively corrects for the artificial base line, and the use of different probe sizes for each material. The initial and final PSD for carbon, silicon and germanium are shown in figure 4. In the initial structures we can clearly see the impact of single vacancies and larger clusters of vacancies accumulated through random deletions at low densities. In each case the averaged PSD is linearly dependent on the density. In the final structures we can see that significant relaxation of the lattice has occurred to lower the energy while accommodating the lower densities. In general, regardless of the material this relaxation results in lower PSD and larger pores as atoms relax and small pores combine.
Although the density is conserved over the simulation cell at each average, the porosity (defined as the volume of empty space divided by the total volume of the cell) is not conserved, since the total combined volume of the pores can vary due to local changes in the bonding. Indeed, the decrease in density eventually leads to a complete rearrangement of bonds through transition to an amorphous phase. Different bond angles and bond lengths can facilitate a local lattice contraction or expansion that means the area under the curve in figure 4 will not be the same. For this reason the change in the dimensionless porosity is also important and is included in

Discussion
We can see from these results that the porosity change profile is unique to each semiconductor, as carbon, silicon and germanium have both a different tolerance for porosity and a different propensity for alternative bonding configurations and local disorder.
In the case of carbon we can see a bimodal change in the porosity indicating two competing mechanisms that influence the final structures ( figure 4(a)). When the density it close to the ideal diamond density of 3.52 g cm −3 the lattice relaxes to convert some of the single vacancies to larger pores (micropores). At lower densities either small 1.5-2 Å pores (vacancies or divacancies) or 2.4-2.9 Å pores are energetically preferred, with a distribution that reflects the initial sampling with minor reconstructions around each pore as the carbon-carbon bonds rehybridise from sp 3 to sp 2 . As the density decreases further, pores at intermediate and larger sizes are accommodated (mesopores). At a density of 2.90 g cm −3 the material undergoes a significant restructuring and involves the coalescence of vacancies into much larger pores indicating that a nanoporous carbon material may be stable. Diamond-like carbon remains microporous, but amorphous carbon becomes mesoporous. It should be noted that an outlier occurs in this dataset at 2.7 g cm −3 where one of the simulations did not result in an mesoporous structure. In this case the distribution of the vacancies introduced during the sampling was uneven and resulted in a metastable defective diamond lattice with local reconstructions but only micropores. In general the density-dependent pore diameters are consistent with observations of activated porous carbon [52], ordered porous carbon [53], carbon fibres [54] and supercapacitive and functionalized porous carbon derived from plants [55,56].
An alternative explanation for the bimodal change in the porosity involves the influence of the competing process of amorphisation. As mentioned above, PorosityPlus also characterizes crystallinity via the third neighbour peak of the g(r), as shown in figure 6. The results for carbon show an almost linear reduction in the number of third nearest neighbours that is remarkably consistent among the statistical samples at each density. This trend correlates with the restructuring around the two types of pores identified by comparing the pore size and PSD, and the density-dependent change in the porosity shown in figure 5(a). At 2.9 g cm −3 the trend changes and a more rapid decline in the third peak height is observed, indicating the collapse of the lattice between the pores and the formation of amorphous carbon. This is consistent with observations of mesoporous carbon materials [2].
In the case of silicon we can see a dramatic change in the structure within a very narrow density range ( figure 5(b), noting the difference in the density increments compared to carbon). The pore diameter range is similar to carbon, and we observe almost no change in the porosity as the density decreases from 2.33 to 2.21 g cm −3 ; in each of the statistical samples silicon resists the formation of pores and preserves the distributed vacancies (see figure 3(b)). When the density shifts below 2.20 g cm −3 an exponential increase in porosity is observed, changing from microporous to mesoporous, and maximising at 2.14 g cm −3 . This range is characterised by large variances as the results vary considerably with the initial configuration of vacancies, but also the range correlates with the disappearance of the third peak in the g(r) indicating that silicon lattice has amorphised (see figure 6(b)), and that the pores reside in amorphous silicon. At densities below 2.14 g cm −3 the porosity declines proportionally with the density, indicating that the amorphous silicon is accommodating the reduction in density without requiring more, or larger, pores. This is consistent with observations that Si exhibits swelling as well as a density reduction upon amorphisation [57,58].
Germanium, in contrast, presents a different porosity profile than carbon and silicon ( figure 5(c)). Germanium exhibits a linear increase in porosity as the density decreases until 4.55 g cm −3 , at which point the porosity diminishes and is indistinguishable (within uncertainties) from the bulk germanium state when the density reaches 4.25 g cm −3 . The onset of amorphisation is rapid and consistent across this range indicating that these mechanisms are acting concurrently, and that mesopores in crystalline germanium are never stable. The porosity change in amorphous germanium is small compared to carbon and silicon and has a greater variance, particularly when the density is low.

Conclusions
The PorosityPlus programme was applied to the characterization of porosity in the vacancy induced amorphisation of defective crystals of carbon, silicon and germanium. With increasing initial vacancies (reduced density) it was found that carbon undergoes two transformations, representing the formation of nanoporous and amorphous solid, respectively; characterised by a bimodal profile of the porosity change. Silicon resists the formation of pores until a threshold (reduced) density is reached upon which it rapidly transforms to a porous amorphous material. With further decreasing density the swelling of amorphous silicon restricts the size of pores. In germanium the onset of porosity is strongly coupled to the formation of amorphous germanium as the density is reduced, and a reduction in porosity is predicted beyond a threshold (reduced) density as the amorphisation is completed.
The PorosityPlus software allows for the identification of vacancies, twin vacancies and larger pores, along with their locations and respective populations. Void migration and coalescence along with the associated density changes can be calculated. Applications are not limited to atomic systems and can be applied to nanoparticle aggregates and even larger macroscopic porous media. This is an ideal tool to integrate into advanced computational workflows, such as creating fingerprints for topological data analysis or machine learning [59], since the porosity profile for each configuration is unique.