Morphology on Reaction Mechanism Dependency for Twin Polymerization

An in-depth knowledge of the structure formation process and the resulting dependency of the morphology on the reaction mechanism is a key requirement in order to design application-oriented materials. For twin polymerization, the basic idea of the reaction process is established, and important structural properties of the final nanoporous hybrid materials are known. However, the effects of changing the reaction mechanism parameters on the final morphology is still an open issue. In this work, the dependence of the morphology on the reaction mechanism is investigated based on a previously introduced lattice-based Monte Carlo method, the reactive bond fluctuation model. We analyze the effects of the model parameters, such as movability, attraction, or reaction probabilities on structural properties, like the specific surface area, the radial distribution function, the local porosity distribution, or the total fraction of percolating elements. From these examinations, we can identify key factors to adapt structural properties to fulfill desired requirements for possible applications. Hereby, we point out which implications theses parameter changes have on the underlying chemical structure.


Introduction
In the last years, the design and development of simple process pathways for application-oriented materials has came up as an important field of research [1][2][3][4][5][6] in material science, chemistry, polymer engineering, and polymer technology. Especially in the topic of adapting materials to specific requirements of applications, the analysis of the morphology and an in-depth knowledge of the occurring structure formation processes induced by the reaction mechanism or the manufacturing process are key features to understand and afterwards to design materials for desired requirements.
One recently developed simple synthesis route to create nanostructured hybrid materials is twin polymerization [6,7]. Twin polymerization starts with twin monomers consisting of at least two different components, an organic and an inorganic one, that are connected via two differently behaving bonds. An exemplary twin monomer is shown in Figure 1a, where the different components and bonds are highlighted in different colors. Due to the complex reaction mechanism these twin monomers react acid-, base-, or thermally induced into two separated but highly interweaved organic and inorganic networks with typical domain sizes of 0.5 to 3 nm [8,9]. By simple post-processing [10], a microporous organic phenolic resin or mesoporous inorganic silica network can be obtained for various applications. It appeared that these materials are of special interest in the fields of anodes for accumulators and gas adsorption systems [11] or for lithium ion batteries [12].   Figure 2. Reaction centers with a solid square are bonded, with a dashed square are non-bonded and with darker filling are blocked. Note that to make it easier to differentiate between the reaction centers that are all placed in the center of the beads, here they are plotted next to each other. a) b) c) d) Figure 2. Exemplary snapshots of a reduced 3D cubical lattice of reduced size 24 × 24 × 24 at (a) t ini and (b) t fin for the parameter combination (I) (see Equation (9)) and two corresponding cutouts of size 6 × 6 × 6 (c): t ini ; (d): t fin ) that are magnified to show the actual existing bond vectors. The colors are chosen in analogy to Figure 1. In (c) only twin monomer structures are observed (see Figure 1b). As over reaction time bonds may cleave and form, a final structure emerge as depict in (d), where a small organic (gray gray) and inorganic (red bonds) network emerge. Note that some initial bonds (green) may survive till the end. The origin of both sub cubes is at position {6, 6, 6} within the large 3D cubes. Note that to make it easier to differentiate between the bond vector types that are connected to the center of the beads, here they are plotted next to each other.
To modify the obtained twin polymer structure, a wide class of different twin monomers has been introduced and various approaches and several experimental studies have been conducted [6,7,13]. However, a detailed understanding of the morphology to reaction mechanism dependency is still missing. Important questions like: "Which elements of the twin monomer structure need to be changed to obtain an nanoporous hybrid materials with a desired pore size distribution?" or "Where to tune the reaction mechanism to generate a specific application based material with desired specific surface area?" are still difficult to answer. To bridge this gap, in the last years several theoretical approaches [9,[14][15][16][17][18][19] at different levels of detail have been established. From these investigations, we gain important insights into the reaction mechanism of the thermal, acid-and base-induced twin polymerization. Furthermore, a reactive bond fluctuation model (rBFM) [9,18,19] has been derived which connects the structure formation process of the reaction mechanism with the macroscopic structural properties of the synthesized twin polymers. It showed that important structural properties of the twin polymers, i.e., the pore size distribution [9], are reproduce in very good agreement with experiments.
The rBFM [9,18,19] is an extended lattice-based bond fluctuation model [20][21][22]. It opens up the possibility to systematically investigate the effects of all reaction relevant parameters and to evaluate their impact on the final morphology. With this knowledge we can deduce a) the parameters that need to be adapted to achieve a specific change in the final structure and b) their implications on the underlying chemical structure. The main focus in this work will be the extraction of major influence factors on structural properties such as bond formation, phase separation, bulk and local porosity, specific surface area, radial distribution function and percolation of the final material. This will be done for the organic-inorganic hybrid materials, as well as for the organic and the inorganic networks theirselves. With this knowledge one will be able to derive chemical characteristics that needs to be fulfilled to synthesize specific morphologies, or specific applications.
The paper is structured as follows: First, we propose the rBFM and its components, second, the twin polymerization and its mapping to the rBFM, and third, all reaction relevant and structural properties that will be investigated. We will introduce and explain all examined model parameters. Then, we will present the obtained results and discuss them. During the discussion, we will identify the major influence factors to modify the final morphology. At the end, a summary of the overall results and possible perspectives for further research are given.

Reactive Bond Fluctuation Model
The reactive bond fluctuation model (rBFM) is a Monte Carlo method that is based on the bond fluctuation model (BFM) known from literature [20][21][22]. Within in BFM chemical structures, like polymers and macromolecules, are represented not by their atoms and interacting forces, but by coarse grained elements that can move with certain probabilities on a (typically cubic) lattice. In our extension, more complex connections between coarse grained elements as well as reactions, i.e., bond formation and cleavage, are included additionally.
Chemical structures that do not alter there steric structure during the reaction process are merged to cubic beads in the framework of the rBFM. Each bead can contain several reaction centers in the center of the beads, which are the connection points between two beads and which are accessible from all sides, i.e., orientation is nonrelevant. A reaction center represents a typical chemical element of the original chemical structure where a bond may form or cleave. A reaction center can take different states like bonded, non-bonded or blocked. The bonds are represented by bond vector types, which are the shortest distance vector between two reaction centers. A bond vector type is defined via the two reaction centers it connects and each bond vector type has an allowed set of different finite lengths, depending on the underlying chemical structure [21,23,24]. In Figure 1 an exemplary structure of a typical twin monomer 2,2 -spiro[4H-1,3,2-benzodioxasiline] (1) and its mapping to possible bead types, reaction centers, and bond vector types is shown. The corresponding 3D representation of such a structure before and after an exemplary reaction is shown in Figure 2. Note that in Figures 1 and 2 for a better differentiation the reaction centers, and bond vectors are plotted next to each other instead of in the center of the beads.
The beads can move on the sites of cubic lattice L of size L × L × L, with an edge length a (lattice constant) between two neighboring sites. Here, we use the lattice notation of Shaffer [21], i.e., the beads move on the even grid lattice nodes, whereas on the odd grid lattice nodes the bond crossings are counted. For structure determination only the even grid lattice nodes are taken into account, which will be called reduced lattice. Typically, periodic boundary conditions are defined, but others are also possible. Each bead occupies one lattice site. Double occupation of lattice sites as well as crossing of bonds is prohibited. In order to include possible effects like phase separation, as it occurs in immiscible polymer blends [25,26] or systems with different hydrophobicities [27], we introduce an attraction parameter λ to allow different pairwise interaction energies between identical (αα) and different (αβ) bead types. Therefore we define a rate constant for the movability as k move = k 0 e β ∆E with k 0 as prefactor, β is the inverse temperature 1/(k B T) and ∆E as energy difference between the system before and after the potential move, to move a bead one allowed position on the lattice. Transferring this approach to the rBFM we obtain with the attraction parameter λ ∈ R representing β and the energy difference ∆E in units of k B T given by the number of neighboring beads of identical or different bead type n αα and n αβ before (t − 1) and after (t) the nMC step. From this we can derive the metropolis probability to accept the move as Note that in the case of λ = 0, all neighboring lattice sites are chosen randomly from a uniform distribution with equal probability for the nMC step. Furthermore, we do not differentiate for λ between αα, αβ, βα, and ββ in order to keep the model as simple as possible but as complex as necessary. Our choice means that beads of the same type attract each other equally strong for both types and that beads of different bead types shows a repulsion strength of the same size. The interaction energy can be attractive or repulsive depending on the local environment as indicated by Equation (2). The simplest version of attraction and repulsion is also studied in [28] .
Parallel to the nMC steps the beads can also perform reactive Monte Carlo (rMC) steps. Therefore, first, an initial bead, second, a possible reaction for this bead, and third, a reaction partner (second bead) within its surrounding is chosen randomly. If • the reaction is allowed and • no crossing of bonds occur after the rMC step, the rMC steps is accepted with a metropolis rate p reac , which will be called reaction probability in the following.The ratio of non-reactive to reactive MC steps can be varied via the ratio m.

Twin Polymerization
We start our systematical investigation of the structure formation process by briefly reviewing the reaction mechanism of twin polymerization for 2,2 -spiro[4H-1,3,2-benzodioxasiline] and its mapping to the rBFM.
The rBFM is intentionally established to investigate the acid catalyzed twin polymerization of the typical twin monomer 2,2 -spiro[4H-1,3,2-benzodioxasiline]. As shown in Figure 3 1 twin polymerizes to a phenolic resin (2) and silica network (3). From experimental [6][7][8]10,29,30] and theoretical [8,[14][15][16][17][18] investigations it is known, that the acid catalyzed reaction mechanism is a reaction in a melt [8]. It can be characterized by three main reaction steps. First, the methylene bonds open (ring opening process), second, the organic network formation via organic bond formations starts, and third, somehow later the inorganic network formation starts by opening the aryl bonds and forming the siloxane bonds at the same time. For further details to the reaction mechanism we refer the reader to [17][18][19]. Mapping the twin polymerization of 1 to the rBFM (see Figure 3), leads to two different types of beads A and B that can be identified with the organic and inorganic component of 1, respectivly. The organic bead type A consists of three reaction centers, C, O, R, whereas the inorganic bead type B has two times two different reaction centers, Si x , O x , as indicated by the index x = {1, 2}. Note that all reaction centers despite O have two states, bonded and non-bonded, whereas O can also be blocked, as depict in Figure 1 .
Due to the reaction mechanism of 1, we introduce four different bond vector types, i.e., methlyene bond (OSi x ), aryl bond (CO x ), organic bond (CR), and siloxane bond (Si x O x ), which can take the bond vector lengths {1, √ 2, √ 3}. We allow the following reactions during the rMC step: • cleavage of OSi x , • formation of CR, • cleavage and formation of CO x , and For our further analysis we keep the basic structure of the original twin monomer 1, as it includes all important features of a twin monomer, i.e., an inorganic and organic components that are connected via two differently behaving bonds, as well as the reaction mechanism itself. However, we varied all reaction probabilities, the ratio m of non-reactive to reactive MC steps, and the attraction parameter λ. Within Table 1 the varied model parameters, their symbols and the range of variation are specified. Hereby it is important to note that the variation of the model parameters can be identified with the following changes in the structure of the twin monomer or the real experimental setup. Modifying the reaction probabilities means a change of the bonding energies, which could be achieved by changing the underlying chemical structure, i.e., replacing atoms. Adjusting m is used as a basis of discussion how diffusive effects influence the structure formation process. The parameter λ does also represent a change of the chemical structure of beads, so that the single components are more or less attractive among each other.

Simulation Details
For a typical simulation run, we place N beads, i.e., N A = 2 N/3 and N B = N/3 beads, randomly on the cubic lattice, where always three beads are arranged as shown in Figure 1b. We choose N as N = L 3Φ with the volume fractionΦ = 0.5. This volume fraction represents a dense melt in the bond fluctuation model [20]. An exemplary filled cubic lattice of reduced size 24 × 24 × 24 is given in Figure 2a and a corresponding cutout of size 6 × 6 × 6 in Figure 2c.We checked the influence of the system size by performing the simulations for L = {48, 72} with a = 1 (lattice notation as Shaffer [21]). We observed no qualitative and quantitative changes of the quantities under investigation.
The reaction time is measured in Monte Carlo Cycles (MCCs). In one MCC each bead can perform randomly a reactive or non-reactive MC step. The ratio m can be varied but it is fixed per simulation run. For the first 1000 time-steps the rMC steps are switched off, thus the system can relax. Afterwards, the rMC steps are switched on and the reaction process starts. Note that we always use the time step before the rMC steps are switched on as initial time t ini . The simulation ends when the relative changes of number of reaction centers in each state are not larger than 5 × 10 −5 over 10 successive time steps on the logarithmic time scale. The final time of each simulation is t fin . An exemplary final cubic lattice of reduced size 24 × 24 × 24 is given in Figure 2b and a corresponding cutout of size 6 × 6 × 6 in Figure 2d.
During the simulation the number of reaction centers in each state and the positions of all beads are tracked over time. This allows us to count the number of the non-bonded reactions centers C-, O-, R-, O x -, Si x -and of the bond vectors CO x , OSi x , CR, Si x O x . In following these reaction centers and bonds are indicated by ξ. Additionally we also determine the number of all nearest-neighbors (n αα (t), n αβ (t)) and next-nearest-neighbors (n αα (t), n αβ (t)) per time step, where αα represents counts between identical bead types and αβ between different bead types.

Process to Structure Analysis
In previous work [9,18,19] we showed that with a proper choice of parameters the rBFM is capable to reproduce the typical three-phase-reaction behavior of the twin polymerization of 1 [18,19], as one observes the initial ring opening phase, the organic, and the inorganic network formation phase, as supposed by previous findings [8,14]. And it can reproduce qualitatively the correct morphology of the final nanoporous hybrid material, as for instance shown for the pore size distribution [9].
In order to qualitatively understand the effects of the model parameters of the rBFM on the final morphology we will perform a systematic analysis of all effects and a discussion of the corresponding consequences in the following. Interesting questions that we want to answer are: Which parameter lead to which changes in the morphology? Is there a different effect on the full nanoporous hybrid material, than on the organic or inorganic components itself? Does the same parameter combination lead to different effects for the different structural properties? Therefore we define chemical properties that quantify the outcome of the reaction mechanism and structural properties that describe the morphology of the final material. A comparison of these quantities can give us the morphology to reaction mechanism dependency.
In order to characterize the chemical properties of the reaction process we focus on two quantities. First, we define the bond fraction BF ξ of reaction center in a specific state ξ or of bond vector types ξ. BF ξ is the ratio of the actual number and the maximal possible number of reaction centers in the specific state or of bond vector types. Second, we determine the ratio of the mean number of actual to initial nearest-neighbor contacts between beads of identical type ( n αα (t fin ) / n αα (t ini ) ) or different types ( n αβ (t fin ) / n αα (t ini ) ). Both quantities reflects the connectivity and network formation between the beads during the reaction process. Note in the following α and β are general replacements for different bead types. It can be either A, B, or A∪B.
To characterize the structural properties of the final rBFM morphologies, we determine the following structural properties  Table 2.
The bulk porosityΦ α and specific surface area S V α , i.e., surface area to volume ratio, are typical structural properties of porous materials. They give an insight whether the material has a large internal surface or volume, which is an important information for applications of the material for filter system, catalyst or gas storage. They are determined as where V = L 3 is the total volume of L, V α is the bulk volume of the α beads, and A α is corresponding the surface area. Counting the number of α beads on the lattice gives V α and A α is obtained by counting the number of free surfaces of α beads, i.e., surfaces where the neighboring lattice site is not occupied by an α bead. As V of L and N are fixed per simulation run, the bulk porosities are known in advance. For the bead type AΦ A = 1/6, for the bead type BΦ B = 1/3, and for the overall material A∪B we haveΦ A∪B = 1/2. To characterize the pore structure of the nanoporous hybrid material experimentally the pore size distribution is determined. In terms of numerical simulations typically the radial distribution function (RDF) g βα and the local porosity distribution µ α [31] are evaluated, which can be related to the pore size distribution [9]. g βα (r) characterizes typical radial distances of clusters of a bead type α to a centered β bead. This can be obtained by counting number of α beads dn βα in a shell of width dr at a radial distance r from a centered β bead. The RDF g βα is then given as where the index βα = {AA, BB}. Note that if all α beads are distributed homogenous g βα (r) = 1. The local porosity distribution µ α (Φ α , K) is an empirical probability density distribution that characterizes fluctuations of the porosityΦ α in the material. It can be derived by defining measurement cells K( x, K) at position x of size K × K × K with K < L and a corresponding total volumeV. These measurement cells are placed at k possible sites on the lattice L and then the local porosity distribution can be determined as with the local porosity Φ α ,V α being the volume of α beads inV, and the δ-distribution.Note that we place K on lattice sites that have a distance of at least K/2 lattice sites from the boundary. Additionally, we investigate the connectivity of the beads by analyzing fluctuations of the percolation probability ϑ α,d and the total percolation fraction θ α,d in direction d (see Table 2). A material is percolating, if there is a closed path of α beads from one side of K to the another side via direction d. The percolation probability is defined as ϑ α,d represents the fraction of measurement cells K of side length K with local porosity Φ α that are percolating in direction d. If d = x, y, z the corresponding space directions if d = 3 in all directions at the same time, if d = c in at least one direction, and if d = 0 no percolation at all is indicated.From this one can determine the total fraction of percolating cells as which is an important quantity to design a network representation of the morphology, as θ α,d gives the fraction of network elements (bonds, sites, . . . ) that have to be permeable.

Results and Discussion
The model parameters under investigation and the analyzed valuesare specified in Table 1. We chose the parameter values in order to sample a broad parameter space close to a log-scale around the control parameter set, which reproduced the twin polymerization as shown in [9]. Thus, there are 288 different parameter combinations per system size L. In the following we discuss the influence of the different parameters and the occurring effects on the reaction mechanism, the structure formation process, and the morphology. During the analysis process, we first investigated and compared the single parameter combinations with each other. Thereby, typical values for the quantities under investigation are obtained for several parameter combinations and a clustering of the results suggested itself in order to reduce the complexity the results. We defined groups of parameter combinations, where certain parameter values are fixed, representing the typical clusters found.In Table 3 the notation and the corresponding fixed parameter value of the resulting 17 different parameter groups G(i) are given. For the analysis of the local porosity and the percolation we choose K = {4, 6,8,12,24, 48} for the system size L = 108 and smaller K for the other system sizes, respectively. Table 3. 17 parameter groups G(i) are defined, where typically one parameter value is fixed. The corresponding parameter values for each group are given in first column and the short notation in the second.

Model Parameter Notation Model Parameter Notation
average over all

Bond Fraction
First, we evaluate the influence of the model parameters on the time development of the overall reaction mechanism. With respect to the bond fraction BF ξ we find a huge variety of time behaviors for the different reaction phases. In Figure 4, the time development of BF ξ for all bond vector types and the non-bonded reaction centers are given over logarithmic time t for two exemplarily parameter combinations An exemplary initial and final 3D structure of a simulation run for the parameter combination (I) is shown in Figure 2.
For the combination (I) (Figure 4a,b) all CO x bond vectors cleave at first and then subsequent most of the OSi x bond vectors. The available non-bonded reaction centers allow then to first form a nearly complete network of A-beads via CR bond vector formation and then to some extend the B-bead network via Si x O x bond vectors grow. For the combination (II) (Figure 4c,d) first the bond vectors of type CO x and OSi x cleave completely. In this case first the B-bead network forms and second the A-bead network, contrary to (I) (compare Figure 4a,c). This has a drastical influence on BF O-, BF O x -, and BF Si x -(compare Figure 4b,d).
Overall, we find a complex behavior of the reaction mechanism due to the different reaction steps that influence each other for the 288 parameter combinations. In detail it can be observed that:

•
With decreasing m all reactions are shifted to later times. With increasing p c CO x the reactions, where C, R, CR, and CO x are related, are shifted to significantly shorter times, whereas the other processes are only influenced in a minor way. • With decreasing p f CR the reaction, where C-, R-, and CR participate, are shifted to significantly later times. The other process steps are not affected. • The parameter group G OSi x influences reactions in various ways, which are connected with O-, O x -, Si x -, OSi x , and Si x O x . In principle these parameter combinations change the duration of the reaction process, so that one can order it by increasing reaction duration. Ordering the results leads Despite the influence of p c CO x , all effects are reflected by the parameter combination (I) and (II) shown in Figure 4.
From the above list, we can sum up that all parameters influence the time behavior of the reaction mechanism, but it seems that only some effects survive as a change in the final values of the bond fractions till the end of the reaction process and thus influence the final morphology. So we examine these changes in the bond fractions BF ξ for t fin depending on parameter groups G(i) (see Table 3). For each parameter combination we get one value per bond fraction per reaction center. leading to a histogram per bond fraction per reaction center containing 288 values. In order to visualize the effects of the different parameters combinations, we highlighted all histogram counts that are counted among a sub group i of a special parameter group G x (i) in a different color. Additionally, we determined the averaged bond fraction BF ξ per reaction center ξ for each parameter group G(i) as the arithmetic mean of all corresponding values. Note that also for the following quantities under investigation the arithmetic mean is used to determine the corresponding average. In Figure 5 the resulting histograms for the BF ξ for ξ = {O-, Si x -, OSi x } for the parameter groups G λ (left column) and G OSi x (right column) are given. Note that the vertical dashed lines indicate the averaged bond fraction BF ξ of the corresponding parameter groups G(i). For all three cases ξ = {O-, Si x -, OSi x } each distribution exhibits one bigger (SD1) and two smaller subdistributions (SD2, SD3). For BF O-SD1 is at (0.95, 1.0], SD2 af (0.75, 0.8), and SD3 at (0.7, 0.75), which can be identified analog for BF Si xand BF OSi x , whereas for BF OSi x the order is reversed. The two smaller subdistributions (SD2, SD3) are fully characterized by G OSi x (2), where the values of the bond fraction are significantly smaller for O-and Si x -(larger for OSi x ) than the values of SD1. SD1 is characterized by G OSi x (1), G OSi x (3), and G OSi x (4), where one observes an internal ordering with the smallest (largest) bond fraction values for G OSi x (1) and the largest (smallest) bond fraction values for G OSi x (3) for O and Si x (OSi x ). Coloring the same bond fraction distribution with respect to G λ , we find for BF O-and BF Si xthat SD1 has contributions of all G λ (i), with increasing (decreasing) order of i, i.e., G λ (1) < G λ (2) < G λ (3). The same order is found for SD2 and SD3, however there SD2 is completely described by G λ (3) and SD2 by G λ (1) < G λ (2). Note that for BF OSi x the same effect but reversed order is found. The same behavior is also reflected by the averaged bond fractions BF ξ over the parameter groups G(i) as shown in Figure 6. G(1) = G all represents the averaged bond fraction over all parameter combinations. We observe variations of BF ξ over G(i) with i > 1 from the overall BF ξ of G all for ξ = {O-, Si x -, OSi x } for two parameter groups G OSi x and G λ , where the effect of OSi x is inverse to O-and Si x -. There is a significant variation for G OSi x (i), where we can order the parameter groups with decreasing (increasing) order of BF ξ of O-and Si x -(OSi x ) as G OSi (2), which is also derived from the distributions above. Furthermore, we find the same minor increase (decrease) of BF ξ for the parameter combinations G λ (1) < G λ (2) < G λ (3) for O-and Si x -(OSi x ) as from the distributions itself.  Table 3.

Phase Separation
In the next step we investigate the phase separation behavior of the reaction mechanism depending on G(i). Again, for each parameter combination we obtain one value respectively for the mean nearest-neighbor contacts n αβ (t fin ) / n αβ (t ini ) with αβ being the bead type combinations AA, AB, and BB. We determine the averaged mean nearest-neighbor contacts n αβ (t fin ) / n αβ (t ini ) as the arithmetic mean of all values that belong to the parameter combination G(i). In Figure 7 n αβ (t fin ) / n αβ (t ini ) over the parameter groups G(i) (see Table 3) are depict. We observe a more complex behavior for n αβ (t fin ) / n αβ (t ini ) than for BF ξ (compare Figure 6). Again a significant change of n αβ (t fin ) / n αβ (t ini ) for G OSi x and G λ is found. Additionally, deviations for G m , G CR , and G Si x O x for BB in major and for ABin minor occurrence can be identified. We find that with increasing m decreasing p f CR and p f Si x O x the relative number mean nearest-neighbor contacts increases, whereas the effect is inverted for AB. This means that the reduced movability and low network formation reactions hinder the formation of clusters, indicated by a high number of nearest neighbor contact.
From the above analysis of the time development of the reaction process we can conclude that if highly interweaved material with high internal connectivity and no phase separation is need, the twin monomers need to fulfill two important features: A and B should be as miscible as possible (λ low) and the movability should be low, while one of the connecting bonds has to have a high activation energy (G OSi x (2)) to compensate phase separation and the corresponding formation of larger cluster.

Analysis of the Model Parameters on the Structure
In the following, we present the results for the structure properties of the overall material A∪B, as well as for the subsystems A and B. We analyze the specific surface area S V α , the local porosity distribution µ α , the percolation probability ϑ α,d , and the percolation fraction θ α,d with α = {A,B,A∪B}. and the RDF g αβ with αβ = {AA, BB}. As in the previous section, next to the quantities itself we also look at the arithmetic means (· ), i.e., called averaged values in the following, over G(i), if possible.

Bulk Porosity and Specific Surface Area
In Section 2.4 the bulk porosities for our system are given asΦ A = 1/6,Φ B = 1/3, andΦ A∪B = 1/2. Looking at the specific surface areas S V α for α = {A, B, A∪B} we find significant changes of the averaged relative specific surface area S V α (t fin )/S V α (t ini ). In Figure 8 the averaged relative specific surface areas S V α (t fin )/S V α (t ini ) are given over the parameter groups G(i). In the case of A∪B the strongest influence factor is the attraction parameter λ, where one observes a decrease of S V A∪B (t fin )/S V A∪B (t ini ) with increasing larger λ. For G OSi x and G Si x O x a minimal change in S V A∪B (t fin )/S V A∪B (t ini ) is observable. These effects are magnified for the subsystems A and B. We find an increasing S V α (t fin )/S V α (t ini ) for the parameter groups in the order G OSi x (3) < G OSi x (4) < G OSi x (1) < G OSi x (2), where the difference between G OSi x (2) and the other three groups is in same order as between G λ (1) and G λ (3). Additionally, we find a secondary effect for the parameter groups G m , G CR , and G Si x O x . There an increase of S V α (t fin )/S V α (t ini ) for α = {A, B} for decreasing m and increasing values of p f CR and p f Si x O x is recognized. Similar to the bond fraction, and nearest-neighbor analysis we find that nanoporous material with a high specific surface area are obtained, when the twin monomer is designed of components A and B that are are miscible (low λ) and the energy barrier of the final connecting bond is high (G OSi x (2)). Furthermore, a low movability (low m) and a prefered and fast bond formation between identical bead types (high p f CR and p f Si x O x ) support a high specific surface area of the material. Figure 8. The averaged relative specific surface area S V α (t fin )/S V α (t ini ) is given for the parameter groups G(i) (see Table 3).

Radial Distribution Function
The analysis of the RDF g βα (r), i.e., per each parameter combination there is one resulting RDF, exhibited that the influence of p c CO x is neglectable and an increase of p f CR leads to a minor increase of g AA and p f Si x O x of g BB . The main effects are observed for the parameter groups G λ (a)-b)), G OSi x (c)-d)), and G m (e)-f)), as depicted in Figure 9. Here, we present the arithmetic mean RDF g βα of all RDFs per parameter combinations that belong to the G(i) (thick colored line), as well as the minimum (lowest thin colored line) and maximum (highest thin colored line) g βα (r) out of the corresponding parameter group G x (i). The area between the minimum and the maximum curve is covered by the various parameter combination within each parameter group and colored in the same color. In the left column the results for βα =AA and in the right column for BB are presented. Note, that the corresponding subgroups i of each parameter group G x (i) are highlighted in different colors. From the plots we find two parameter groups with major variations and one with minor one. The major effects are again observed for λ and for G OSi x . With increasing λ g AA (r) and g BB (r) gets higher and wider by a factor of around 1.5-2 in height and width. For BB additionally we observe that for small values of λ the maximum is truncated. For G OSi x a major difference is observed between G OSi x (2) and G OSi x (i) with i = 1, 3, 4. The distributions for G OSi x (2) are narrow and spiky with lower maximal values, whereas the distributions of the other three groups are broader with higher values, depending on λ as well. Note that for G OSi x the width increases up to a factor of 5, whereas the height increases by factor of 1.5-2. Minor difference is displayed by G m , where find that with increasing m the RDF g βα (r) primary gets wider (factor around 1.5) and only a minimal increase in height. Figure 9. The arithmetic mean, i.e., averaged (thick colored line), the minimum (lowest thin colored line) and the maximum (highest thin colored line) radial distribution function g βα (r) over radial distance r are given for ba = fAA, BBg in (a,c,e) and (b,d,f) and for the parameter groups G λ , G OSi x , and G m in (a,b), (c,d), and (e,f). Each subgroup i is highlighted with a different color.

Local Porosity Distribution
Alternatively to the RDF, the local porosity distribution µ α is used to characterize the pore structure of the material. In Figure 10 the results for the local porosity distributions µ α (Φ α , K) for α = {A,B,A∪B} are shown exemplary for K = {8, 12, 24} at t ini in (a,b), and t fin in (c,d) for the two parameter combinations (I) and (II) given in Equations (9) and (10). Analyzing the distributions depending on K, we observe that in all cases the distributions get localized aroundΦ α . This is accepted, as µ α represents the fluctuations of Φ α for K and when K → L follows that Φ α →Φ α [31]. Furthermore, we find that the initial distributions for α = {A,B,A∪B} for (I) and (II) are similar in shape independed of the parameters. This reflects the random initialization of the beads on the lattice. However for t fin the distributions are quite different. For (I) the final distributions is more narrow and peaked (see Figure 10c), whereas for (II) they get wider and flatter (see Figure 10d). As the changing width is an important feature of these distributions we will also analyze the variance (µ α (Φ α , K) − µ α (Φ α , K) ) 2 of the local porosity distributions, next to their mean values µ α (Φ α , K) ( Figure 11).    The averaged mean values of the local porosity distributions µ α (Φ α,K ) are nearly constant at values of the corresponding bulk porosityΦ α (see Figure 11a). The variances however, display the previous findings, i.e., changes depending on the parameter groups G m , G λ , G OSi x , G CR , and G Si x O x . In Figure 10b the averaged variances (µ α (Φ α , K) − µ α (Φ α , K) ) 2 are given over the parameter groups G(i) (see Table 3). We observe for smaller measurement cells larger the fluctuations. For K > 6 in the case of α =A∪B and for K > 12 for the subsystems the variances goes down to 0, thus only results for K ≤ 12 are plotted. The observed changes in the variance fit very well with the previous findings. The major impact is obtained by varying G λ and G OSi x . For all α the variance increases for G λ (i) with increasing i and we observe the smallest variance for G OSi x (2) and the largest for G OSi x (3). This means the more attractive the beads the wider gets the distribution and the more clusters of different sizes are found. This effect is enhanced if the final connecting bond has a low energy barrier and beads separate easily. Furthermore, the subsystems A and B show a minor dependency on G m , G CR , and G Si x O x . We find a decrease of (µ α (Φ α , K) − µ α (Φ α , K) ) 2 for decreasing m or for increasing p f CR and p f Si x O x , i.e., the less the beads move and the faster the networks are formed the more homogenous are the cluster sizes in the final material. Figure 11. In a) the averaged mean values of the local porosity distribution µ α (Φ α , K) and in b) the averaged variances (µ α (Φ α , K) − µ α (Φ α , K) ) 2 for different lengths K = {4, 6, 8, 12} of the measurement cells for α = {A, B, A∪B} are shown over the parameter groups G(i). Note that for K > 12 all variances are 0 and thus are not shown here.

Percolation
Finally, we investigate the total fraction of percolating measurement cells along different directions d given in Table 2 for different cell sizes K. We find that the obtained structures are fully percolating or not at all in nearly all cases, i.e., θ α,d (K) = 1 for d = {x, y, z, 3, c} and 0 for d = 0. There are two exceptions: α =A∪B and d = {x, y, z, 3, c} at small K = {6, 8} and α = {A, B} and d = 0. These two cases are shown in Figure 12a,b. As the effects for K = {6, 8} are similar, but the magnitude of the effect for K = 6 is larger, we present here only the results for K = 6.  Table 3). The results are shown in (a) for α=A∪B, K = 6, and d = {x, y, z, 3, c}, whereas in (b) for α = {A, B} and d = 0.
In Figure 12a the averaged total fraction of percolating measurement cells θ α,d (K) is shown for d = 0 and α =A∪B at K = 6. We find an overall high value of θ A∪B,d (6). The highest value is observed for d = c, i.e., there is nearly always a percolating direction within the material. The second highest θ A∪B,d (6) are found for any single percolating direction (d = x, y, z). Here all three directions are nearly the same, thus there is no biased direction within the material. The lowest value is obtained for a percolation all three directions at the same time, i.e., d = c. Independent of the percolating direction, we find an decreasing θ α,d (6) for increasing λ, as the beads form packed clusters that are not necessarily connected to both sides of the measurement cell and thus can not percolate. We find that an enhanced network formation (high value of p f CR , p f Si x O x ) and a low movability of the beads (low value of m), increases θ A∪B,d (K). This can be accounted to the high connectivity of identical beads due to the fast network formation or due to the homogeneous distribution of beads, as they can not move too far away from the initial homogeneous distribution. Interestingly, G OSi x shows here only a minor influence on the percolation for A∪B. In Figure 12b θ α,0 (K) we show for A and B for all K.
Here, the case of no percolation in any direction is analyzed. As mentioned above, for θ A∪B,0 (K) = 0, but for the subsystems it is not. Analyzing the general behavior of θ α,0 (K) over K we find that for the bead type A that percolation is restored for system size K > 12, whereas for B it already restored for K > 8. For both subsystems we find changes in θ α,0 (K), depending on G λ , G OSi x , G CR , G Si x O x , and G λ , similar to the previous results. With increasing λ θ α,0 (K) decreases, as the probability of percolation increases. There is a drastic drop of θ α,0 (K) for G OSi x (2), as the high connectivity of the beads for this parameter combination typically leads to percolation. This effect is also indicated by the decrease of θ α,0 (K) for increasing p f CR , p f Si x O x . Finally the decrease of θ α,0 (K) with decreasing m is due to the low movability that keeps the beads close to their initial random starting position on the lattice, which enhances the chance of percolation.

Summary
The focus of this work was the derivation of a structure on reaction process dependency for twin polymerization and the identification of influencing factors to modify the final morphology to allow an application-oriented adaptation of the final twin polymers. Therefore, we utilized the previously introduced reactive bond fluctuation model for twin polymerization of the typical twin monomer 2,2 -spiro[4H-1,3,2-benzodioxasiline]. We kept the reaction mechanism and the basic composition of the twin monomer for our analysis as both can be recovered for many other twin monomers as well. However, we varied all reaction relevant parameters and analyzed their effects on the final morphology. In this context, we investigate various chemical and structural properties, that are directly or indirectly connected with the reaction mechanism, as the bond fractions, the specific surface area, the bulk and the local porosity, percolation effects and the radial distribution function for a wide range of possible parameter combinations.
Our investigations showed that there are four main components that affect the final morphologies; most of all: That are the attraction between the two chemical components (beads), the strength, i.e., activation energy, of the final connecting bond between two different bead types, the probability of identical beads to form bonds and the movability of the beads. The more attractive identical bead types are the more often the beads are bonded and thus the larger and more heterogenous are the clusters sizes of identical beads. Consequently, the specific surface area decreases, the cluster size and number of beads per cluster increases, and the percolation probability decreases. The higher activation energy of the final connecting bond, the longer are the beads of different type connected, and thus, the smaller are the cluster sizes and beads per cluster leading to more homogenous materials. This is also reflected in the percolation property of the obtained materials, as a high connectivity increases the probability to find percolation. The described effects of the strength of the final connecting bond can be enhanced by increasing the probability to form networks between identical bead types. Additionally, we found that the movability influences the structural properties of the material as well. The higher the movability, the larger the cluster sizes and number of beads per clusters, i.e., the material is more heterogenous due to the different cluster sizes and this decreases the percolation probability, however, note that the bond fractions have not been affected by the movability.
Hence, we found a direct dependency of the reaction relevant parameters to structural properties of the final morphology. We were able to deduce direct proportionalities on how to change the underlying components to design certain effects, i.e., a high specific surface area or homogeneous cluster sizes, within the material. At this point it is now of special interest to analyze the already known class of twin monomers with respect to the above-mentioned main components that affect the morphology and to compare the resulting findings with the properties of the used 2,2 -spiro[4H-1,3,2-benzodioxasiline]. So we will be able to classify the twin monomers within our parameter groups and to check our prediction on the real materials.
Funding: This research was funded by the German Research Foundation/DFG in the framework of the DFG Forschergruppe (FOR) 1497, grant number PR 1580/1-1 and the publication costs of this article were funded by the German Research Foundation/DFG -392676956 and the Technische Universität Chemnitz in the funding program Open Access Publishing.

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

Abbreviations
The following abbreviations are used in this manuscript: BFM bond fluctuation model rBFM reactive bond fluctuation model RDF radial distribution function nMC, rMC non-reactive Monte Carlo, reactive Monte Carlo