Constructing bilayers with tuneable ring statistics and topologies

ABSTRACT A computationally tractable method is developed and described to generate two-dimensional networks with the aim of producing configurations for thin films (bilayers) of SiO and related materials. The method developed allows ideal (defect-free) networks of any given shape to be grown from seeds with both tuneable ring statistics (ring distributions) and topologies, the latter characterised by the Aboav–Weaire parameter, α. The method developed is demonstrated by growing networks which differ in their ring distributions and topologies as controlled by a combination of the choice of the ‘allowed’ rings and the effective growth ‘temperature’. Configurations are generated with Aboav–Weaire parameters commensurate with those obtained from an analysis of experimental configurations, improving significantly on previous methods for generating these networks (which systematically underestimate this parameter). The ability to efficiently grow configurations allows us to explore the structural basis of Lemaitre's law (which couples the underlying network ring distribution second moment with the fraction of six-membered rings, ) in terms of maximum entropy. A rationale for the commonly observed value of is presented as a balance between entropic and enthalpic contributions to the free energy. The deviations of the respective ring areas from the ideal are discussed and the relative insensitivity of the ring area (in systems of this type) to relatively strong distortions is highlighted. GRAPHICAL ABSTRACT


I. Introduction
Thin films of silicon dioxide, SiO 2 , have been prepared experimentally by chemical vapour deposition on metal and graphene supports [1,2]. Consisting of two symmetric layers of corner sharing [SiO 4 ] tetrahedra, these silica bilayers can be viewed topologically as two-dimensional (2D) networks of percolating rings ( Figure 1). Indeed, scanning tunnelling microscopy (STM) has been used to directly visualise their atomistic structure, revealing both crystalline and glassy arrangements and even the interface between the two [3,4].
The novelty of these effective 2D networks has led to multiple complementary computational studies on silica bilayers. These have included both ab initio methods [5,6] and molecular dynamics studies using classical forcefields of varying levels of theory [7,8]. In order to perform these simulations, it is necessary to have a starting atomistic configuration. This can be achieved in multiple ways. The most straightforward is to take one of the existing experimental distributions. These are however limited in size and number, and can contain defects or areas which could not be fully imaged. Computational techniques are therefore often preferable. These include construction from (potentially topologically related) amorphous graphene distributions [7] and bond-switching algorithms where crystalline graphene sheets are amorphised to match the properties of experimental samples [9].
Generation of configurations with all the required properties has proved surprisingly difficult. Networks consist of rings of varying size, n, defined by the number of constituent vertices. For silica bilayers, experimental configurations indicate that n = 4−10. The ring structure of a network can be characterised by two measures. The first is the ring distribution, p n (also known as the ring statistics), which gives the proportion of rings in the network of size n (see Figure 2(a)). A key constraint that arises from Euler's formula for connected planar graphs is that the mean ring size must be equal to six, n = 6. The second is the ring correlations, the tendency of rings of different sizes to be adjacent to one another. The ring correlations can be described by a single parameter, α, using the empirical Aboav-Weaire law [10,11], which states: where m n is the mean ring size about a ring of size n and μ = n 2 − n 2 , the second moment about the mean. The parameter α is simply determined from a linear fit, as shown in Figure 2(b). A larger value of α indicates an increased tendency of small rings to be adjacent to large rings. For experimental silica bilayers, the value of α > 0.3 is higher than in many natural systems [12], suggesting a strong tendency for large rings to be positioned next to small rings. Whilst bond-switching and molecular dynamics can generate configurations with similar ring distributions to experiment, replicating α, and therefore exact network topology, has proved elusive [8,13]. Both the experimental ring distribution and correlations have been achieved by targeted optimisation of the dual lattice, but this is a purely statistical process and does not consider the energetics of the system [14]. The motivation of this work was therefore to develop a construction algorithm to generate samples of silica bilayers which can capture the full 2D network topology; both the ring distribution and correlations. The model should be able to explore all phases from crystalline to amorphous yet computationally efficient enough to generate configurations suitable for further high throughput calculations. To achieve this, we have developed a grow-from-seed Monte Carlo algorithm, where rings are individually added to build up a sample. This allows for 'organic' formation, where growth cannot be influenced by enforced periodicity. This model is in a way similar to the first hand-built models [15], which were superseded by computational techniques designed to generate periodic configurations. However, the recent development of techniques to simulate aperiodic samples, such as sliding boundary conditions for molecular dynamics [16], makes this constraint no longer essential, and we may benefit from the added freedom of an aperiodic model.
The structure of this paper is as follows. We will first explain our model for silica bilayers and the details of the construction algorithm. We will then provide results showing the capabilities of the method and compare to experimental data and existing computational models. The results will then be considered in the context of maximum entropy. Finally, we discuss the applicability of an area law to atomic 2D networks.

II. Methods
As with the bulk glass, the basic building blocks of silica bilayers are vertex-sharing [SiO 4 ] tetrahedra. These are arranged such that three of the vertices are connected to tetrahedra in the same layer, with the final vertex being shared between layers acting as a 'bridge'. A consequence of these bridging oxygen atoms is to enforce a symmetry plane in the middle of the upper and lower layers, which has been observed experimentally and in simulation [17]. The two layers are therefore exact mirror images and their corresponding ring structure necessarily identical.
Owing to this relationship, it is possible to capture the full topology of silica bilayers with a simplified representation. Discarding one of the layers, with the bridging oxygens, and projecting the remaining atoms onto the horizontal plane leaves a 2D network of vertex-sharing [SiO 3 ] triangles (see Figure 1). This allows the system to be conveniently mapped onto two dimensions whilst retaining all the key information regarding ring statistics and ring connectivities. Additionally, the bilayer structure can be easily regenerated by the reverse process.
The focus of this work is on generating a large number of samples with varying ring statistics, to be used as a base for further calculations. Therefore, we propose that working with this reduced representation is sufficient, as it provides a computationally efficient way to produce networks with the required topology. The precise geometry of the bilayer can be refined with advanced optimisation techniques if required.
In order to simulate bilayer systems in two dimensions, a suitable potential model is needed which captures the essential physics of the system: the local triangular environment of the [SiO 3 ] units and the varying relative energies of rings of different sizes. The model used here is modified from a relatively simple potential used in allatom bilayer calculations [7,18]. Each [SiO 3 ] unit has a harmonic potential acting between all three Si-O pairs, and the three nearest-neighbour O-O pairs, given by: where k is a constant, d is the interatomic separation and d 0 AB the equilibrium interatomic separation between species of elements A,B. The equilibrium separations are set such that d 0 OO = √ 3d 0 SiO , driving the system towards a set of ideal [SiO 3 ] triangles.
The Si-O-Si angle, which determines the strain associated with different ring sizes, is controlled by a shifted and cut 24-12 potential of the form: where is a constant and r is the Si-Si separation between atoms in adjacent triangles. It is the value of r 0 which sets the Si-O-Si angle at which strain begins to be felt, denoted φ 0 . The angle φ 0 can in turn be related to the minimum angle in the ring, θ 0 , by φ 0 = 60 • + θ 0 (see Figure 3). The separation r 0 can then be calculated using the equation: We take the ideal hexagonal lattice as being the zero in energy, such that θ 0 = 120 • and r 0 = √ 2d SiO . Rings in which the angle deviates increasingly from θ 0 will therefore incur an increasingly energetic penalty. . The strain in rings of different sizes can be related to the angle, θ, between adjacent triangles. For rings with n ≤ 6 (panels (a), (b)), this is related to the interior angle of the ring, but for n > 6 it is related to the exterior angle (panel (c)).
Our primary aim here is to generate topologies suitable for later investigation using more detailed (and hence more accurate but more computationallydemanding) potential models. As a result, the harmonic springs simply control the local (triangular) geometries whilst the 24-12 potential controls the repulsion between these local polyhedra. These functions are chosen as deliberately simple to improve computational efficiency and achieve high throughput of idealised networks. Furthermore, the parameters k and need have no direct physical meaning, simply controlling the meaning of the system 'temperature' as discussed below. Our only requirement is that they generate energies of the same magnitude to allow for efficient structural evolution.
Using the model detailed above, a Monte Carlo construction algorithm has been developed which allows two-dimensional networks to be built ring by ring in the shape of a specified function. The main steps of the algorithm are outlined below: (1) Take a starting seed, such as a single ring or experimental configuration (2) Select triangles on which to build the next ring (i) Overlay a function on the network (e.g. circle, square) (ii) Check for atoms with dangling bonds lying inside the function region (iii) If no such atoms exist, systematically increase the function size until an atom is found (iv) Find the next nearest atoms which also have a dangling bonds (v) Choose the two triangles that correspond to the largest starting ring size (see Figure 4) (3) Determine the probability of constructing rings of different sizes (i) Build trial rings in the range n min to n max (see Figure 5) (ii) Geometry optimise the local structure and calculate minimised potential energy  (iii) Calculate the probabilities of each ring occurring, P n , using Equation (5) (4) Accept single trial ring according to the probability distribution (5) Repeat steps (2) → (4) until the target number of rings is reached The probability of a ring of size n being accepted, P n , is given by the equation: where E n and E 0 correspond to the energy of the trial structure and lowest energy of all trial structures respectively, and T is a 'temperature'. The parameter T controls how easily the potential energy landscape can be explored, and therefore how accessible strained rings become. In the low T limit, the acceptance probabilities Table 1. Variation of acceptance probabilities with temperature for the configurations in Figure 5. are dominated by the energy term, and the rings which are selected will be those with the lowest energy. Note that this is not necessarily the 6-ring, but rather is dependent on the local environment. On the other hand, in the high T limit, the acceptance probabilities are approximately equal, and rings are selected on a more random basis. This is demonstrated in Table 1, using the example configurations from Figure 5. The 'temperature' parameter is therefore the primary method for controlling the distribution of ring sizes in constructed networks.

A. Network growth
We will demonstrate that our method is robust and controllable, and is able to generate configurations which have the same topological properties as those visualised from experiment. As such, results will largely focus on the system where n = 4−10, denoted {4, 10}. For reference, six example configurations are given in Figure 6, which are generated with a range of temperatures and growth geometries. Panels (a) to (d) provide a good qualitative analysis of the effect of temperature on the ring structure. At low temperature a phase boundary can be seen separating crystalline and amorphous regions, as seen in experimental silica bilayers [4]. In these samples although the proportion of small and large rings is low, their positions are highly correlated and chain structures of alternating rings sizes are clearly present. These motifs are reminiscent of defects found in a wide range of materials, including amorphous graphene and thin silicon and germanium oxides [17,[19][20][21]. The increase in temperature is coupled with the emergence of rings of more extreme sizes and regions which could be viewed as nano-crystalline are dispersed. The high temperature limit reveals an amorphous structure. Panels (e) and (f) give examples of the diverse geometries in which samples may be constructed. It is interesting to note that even 'difficult' shapes, such as those containing concave regions and cusps, do not prevent growth. Although the shape does not affect the network topology and is in a sense arbitrary, certain calculations may benefit from the different configurational shapes. For instance, molecular dynamics with sliding boundary conditions requires fitting of a smooth function to the sample perimeter, which is facilitated by having a near-circular form.
The quantitative relationship between temperature and ring structure was investigated for three systems of varying ring size ranges; {5, 7}, {4, 8} and {4, 10}. For each system, 100 samples consisting of 1000 rings were grown at temperatures between T = 10 −4.5 → 10 −1.5 . The evolution of the combined ring statistics with temperature is presented in Figure 7. Panels  of just three ring sizes, Euler's formula demands that p 5 = p 7 = (1 − p 6 )/2 and so the system is relatively well defined. Hence as the 5 and 7-rings are more strained than the 6-ring, p 5 and p 7 show a systematic increase with temperature. Furthermore, the uniform equilibrium distribution can only satisfy Euler's formula when the ring size range is symmetric about 6, as is observed for {5, 7} and {4, 8}. The form of the ring statistics at intermediate temperatures and for {4, 10} will be discussed later in terms of the maximum entropy distribution.
The ring distribution for {4, 10} is also shown as a function of temperature in Figure 7 panel (d), along with the value of the Aboav-Weaire parameter, α, allowing for more facile comparison with experiment. The temperature which gives the best agreement between our model and amorphous experimental samples is highlighted by the vertical dashed line. The values of p 7 and α are provided in Table 2, alongside results from two experimental samples. We see that our model can be successfully tuned to match the topology of the experimental system. Not only are the ring distributions in very good accordance, but also the ring correlations, which have until now proved difficult to capture. We are therefore confident that our simplified but physically motivated model is able to reproduce the behaviour of real systems. Table 2 also lists the ring statistics obtained from previous computational studies which used both Monte Carlo [13] and molecular dynamics [8] methods. Roy et al. applied an effective pair potential coupled with molecular dynamics (analogous to a method employed previously to study amorphous graphene [24]) to generate monolayers of stoichiometry M 2 X 3 . Such methods offer the potential for generating realistic configurations but are difficult to control as the cooling rates which must be  applied are necessarily huge compared to experimental rates. A potential artefact of the high cooling rates is the effectively freezing in of defect states, either in terms of local coordination environments [24] or highly-strained (three-membered) rings [8]. A key difference in the ring distributions obtained in the present work, compared to those obtained previously, lies in the final Aboav-Weaire parameter, α, which describes the ring connectivity. Previous work appears to systematically underestimate these values, indicative of too little structural ordering.
As an additional check that the model behaves physically, the angle distribution between adjacent [SiO 3 ] units, f (θ ), was calculated for the {4, 10} system across the range of temperatures studied. The results are summarised in Figure 8. The angle distributions are necessarily symmetric about 120 • , as each triangle pair contributes two opposing angles. At lower temperatures the distribution is dominated by angles close to 120 • , as a consequence of the large proportion of near strainless six-membered rings. Furthermore, at the temperature corresponding to the amorphous experimental region, T = 10 −3 , the distribution has a similar extent to the angle distribution found in experimental samples (see for example Figure 7 ref. [8]). However, as the temperature increases, the form of f (θ ) does not simply broaden as might be expected, but becomes bimodal. This can be rationalised by considering the angles that would be present in regular polygons of different sizes, marked by vertical lines in Figure 8. These ideal angles are clustered away from the mean value of 120 • , and hence increasing the diversity of ring sizes through temperature acts to shift the most commonly observed angles from the central value of 120 • . It is therefore interesting to note that increasing structure in the angle distribution does not necessarily translate to increased order in the atomic configurations.

B. Maximum entropy and Lemaitre's law
The question remains as to the form of the ring distribution in this amorphous phase. Many studies fit a lognormal curve to observed distributions, as good agreement was originally noted in vitreous silica models by Shackelford [25]; but this does not explain why such a functional form arises. This question is however a generic problem of 2D networks, and an answer based on maximum entropy was provided by Lemaitre et al., whilst studying the assembly of hard discs on an air-table [23]. The information entropy of a discrete distribution can be written, As expected, the entropy of an unconstrained system will be maximised by a uniform distribution. We know however that 2D networks are constrained and Lemaitre et al. argued for the following, n np n = 6, where Equations (7) and (8) are just statements of normalisation and Euler's law, respectively. The final Equation (9) is an empirical constraint based on an observed property of the system. Lemaitre originally chose f n = 1/n as a result of measurements of the area of rings in 2D mosaics, although the exact form of the function does not greatly affect the following conclusions [26]. Applying Lagrange's method, the ring distribution can be written, where p n satisfies the above constraints. This can be used to yield the maximum entropy distribution for a given the value of p 6 . These results are summarised in Lemaitre's law, which displays the second moment, μ, against p 6 for the maximum entropy distribution, as in Figure 9(a). Many systems throughout biology, chemistry and physics have been shown to follow Lemaitre's law, including hard discs, epithelial cells and geographical regions [27][28][29].
The maximum entropy distribution is found to closely match the theoretical and experimental ring statistics for amorphous silica bilayers, as highlighted in Table 2 and Figure 2(a). Furthermore, as demonstrated in Figure 9, our computational growth model also follows Lemaitre's law across the entire temperature range. In this figure, the results from the individual 1000 ring samples are presented, coloured by temperature. Panels (a) and (b) compare the observed μ and S of the generated configurations to those expected from Lemaitre's law. This law provides a good fit, with only a small deviation observed for p 6 > 0.5. Panel (c) plots the geometry optimised potential energy of the samples against p 6 , which increases as the ring sizes become more diverse. The curve is split into two regimes, with gradual increase in energy from p 6 = 1.0 → 0.4 followed by exponential increase for p 6 < 0.4. This is consistent with the information in Figure 7(d) which shows that below p 6 ≈ 0.4, not only does the number of extreme ring sizes increase rapidly, but they become less correlated with a lower α, decreasing the number of favourable small-large ring pairings. We can now propose why it is that the experimental amorphous distributions are found with a value of p 6 ≈ 0.4. The system aims to maximise entropy by obtaining a ring distribution along the Lemaitre curve with the minimum p 6 possible. However, for p 6 < 0.4 the energetic cost becomes prohibitively large, as higher entropy distributions can only be achieved by increasing the proportion of extreme ring sizes at the expense of relatively low strain 5-and 7-rings.

C. Ring areas
Inspection of amorphous experimental samples reveals that the rings appear highly regular in shape. This can be quantified by determining the average dimensionless area for each ring size, A n , and comparing it to the area of the corresponding regular polygon, A 0 n , where: A 0 n = n 4 tan(π/n) .
As the regular polygon has the maximum achievable area for a given ring size, the ratio A n /A 0 n is expected to lie in the range 0 → 1, with a lower value corresponding to increased deviation from regularity, and assuming d 0 SiSi to be fixed.
The study by Kumar et al. found that whereas for experimental configurations, A n /A 0 n ≈ 1, configurations generated using bond-switching generally displayed ratios much less than unity [24], indicative of large distortions in the ring structure. For larger rings, a value of A n /A 0 n > 1 was also found, which can only be achieved if there is appreciable bond stretching (see Equations (11), (12)).
The analogous results for the method presented in this work can be found in Figure 10(a), for T = 10 −3 . This figure demonstrates that there is now good agreement between experimental and computational results. We note that in both cases the deviation from regularity increases with ring size, as the flexibility of the rings increases. We propose that the difference between current and previous methods could be due to the lack of enforced periodicity on the system. By allowing the network to grow relatively freely, the system can avoid a build up of strain associated with maintaining periodic boundaries.
Even with this analysis, an argument can be made that by visual inspection the rings in the experimental configurations are still more regular than those generated from computational samples. We therefore consider if deformation of a ring should be expected to lead to significant reduction in area. This can be explored by considering the distortion of a circle to an ellipse. The degree of distortion can be described by the eccentricity of the ellipse, = (1 − b 2 /a 2 ) 1/2 , where a, b are the major and minor axis radii respectively. If the circle and the ellipse have the same circumference, it can be shown that, such that 0 ≤ A n /A 0 n ≤ 1. This function is shown in Figure 10 panel (b). We see that a large degree of Figure 10. Panel (a) compares the regularity of rings in computational and experimental amorphous configurations, with points indicating the mean value of A n /A 0 n and bars corresponding to the standard deviation. Experimental data is taken from ref. [13]. Panel (b) shows the effect on the area when distorting a circle to an ellipse whilst maintaining a constant perimeter length. The solid line gives the exact result whilst the dashed line approximates the perimeter of the ellipse as √ 2π(a 2 + b 2 ) 1/2 . eccentricity is needed for a significant change in the observable area. For example, if a = 1.5b, the area is still 92.3% of the area of the corresponding circle. For silica networks, the Si-Si distances lie in a narrow range because of the covalent nature of the atomic bonding and the near-linear Si-O-Si bridges which join the two layers. Hence we would expect similar behaviour to occur, with ring areas relatively invariant to distortions in the ring shape (this same analysis would not be expected to hold for foams for example, where the length of the boundary is much more flexible). This suggests that the ring area is not the most suitable metric for quantifying the regularity of rings in systems such as this, and could explain any disagreement between the seemingly near ideal ring areas and the visual evidence. As previously stated, although the potential model used is physically motivated, it is lightweight in order to facilitate generation of a large number of configurations with the correct network topology. In future, it would be informative to see if the required regularity can be achieved by geometry optimising the resulting bilayer configurations with a more accurate potential, such as the TS potential which includes potentially significant electrostatic interactions including many-body polarisation effects [30].

IV. Conclusion
In this paper, we have developed a method for the effective growth of two dimensional networks from a given seed, allowing for control over the ring size distributions and the system topologies. The latter is often characterised by the Aboav-Weaire parameter, α, and the values obtained here are more commensurate with those obtained from experimental imaging compared with previously constructed configurations. The high throughput method has allowed a detailed analysis of Lemaitre's law and has highlighted why the fraction of six-membered rings observed in real systems is often ∼ 0.4. Finally, a consideration of the ring areas shows our configurations to contain more regular polyhedra than a number of previous configurations. However, the area itself is shown to be a relatively poor measure of a deviation from ideality for systems of this type.