Mesoscale Modeling of Phase Separation Controlled by Hydrosilylation in Polyhydromethylsiloxane (PHMS)-Containing Blends

Controlling morphology of polysiloxane blends crosslinked by the hydrosilylation reaction followed by pyrolysis constitutes a robust strategy to fabricate polymer-derived ceramics (PDCs) for a number of applications, from water purification to hydrogen storage. Herein, we introduce a dissipative particle dynamics (DPD) approach that captures the phase separation in binary and ternary polymer blends undergoing hydrosilylation. Linear polyhydromethylsiloxane (PHMS) chains are chosen as preceramic precursors and linear vinyl-terminated polydimethylsiloxane (v-PDMS) chains constitute the reactive sacrificial component. Hydrosilylation of carbon–carbon unsaturated double bonds results in the formation of carbon–silicon bonds and is widely utilized in the synthesis of organosilicons. We characterize the dynamics of binary PHMS/v-PDMS blends undergoing hydrosilylation and ternary blends in which a fraction of the reactive sacrificial component (v-PDMS) is replaced with the non-reactive sacrificial component (methyl-terminated PDMS (m-PDMS), polyacrylonitrile (PAN), or poly(methyl methacrylate) (PMMA)). Our results clearly demonstrate that the morphology of the sacrificial domains in the nanostructured polymer network formed can be tailored by tunning the composition, chemical nature, and the degree of polymerization of the sacrificial component. We also show that the addition of a non-reactive sacrificial component introduces facile means to control the self-assembly and morphology of these nanostructured materials by varying the fraction, degree of polymerization, or the chemical nature of this component.


Introduction
Controlling reaction kinetics along with the dynamics of phase separation plays an important role in polymer processing, allowing one to tailor a range of material properties. In particular, controlling the phase separation process in Si-based polymers followed by pyrolysis constitutes a robust strategy to fabricate porous polymer-derived ceramics (PDCs) for a number of applications from water purification [1] and microfiltration [2] to hydrogen storage [3]. The pore size of PDCs fabricated using phase separation can range from a few nanometers [4,5] to several hundred microns [6]. The processing of porous PDCs typically involves blending a sacrificial component and preceramic precursors; crosslinking is typically required prior to the pyrolysis [7][8][9]. Polysiloxanes (PSO) are versatile materials with superior chemical, physical, and electrical properties [10]. Ceramic materials derived from PSO exhibit excellent mechanical strength and high thermal stability while also having relatively low processing temperatures [7,11]. The crosslinking of preceramic polymers prior to pyrolysis prevents the porous structures from collapsing and results in high yields of ceramics [11]. Crosslinking of linear PSO, such as linear polyhydromethylsiloxane (PHMS), Nanomaterials 2022, 12, 3117 2 of 15 results in the formation of a ladder-type network structure [12,13]. The hydrosilylation reaction between the silicon-hydrogen group of PHMS and the vinyl groups of an olefin, such as divinylbenzene (DVB) or vinyl-terminated polydimethylsiloxane (v-PDMS), provides a facile method of crosslinking and is widely used in synthesis of organosilicons [14][15][16]. Due to the complete decomposition during pyrolysis, PDMS is often selected as the pore former (or the sacrificial component) to fabricate porous PDCs [9,17,18].
The morphology of phase-separated precursors significantly affects the structure of the PDCs. It had been shown that the obtained macroporous ceramics closely resemble the microstructure of the crosslinked precursors [19]. While the total porosity depends on the fraction of the sacrificial components used [1,2,6], the pore size and pore connectivity of PDCs can be controlled via regulating the dynamics of the phase separation process. For example, choosing a sacrificial component of higher viscosity [20] or a higher fraction of a sacrificial component [2,19,20] results in larger pores. In the system of PHMS crosslinked with v-PDMS and a crosslinking enhancer, the size of the generated mesopores was shown to increase with an increase in the molecular weight of v-PDMS [21]. The miscibility between the polymeric components directly affects the phase separation and subsequent microstructure of PDCs [18,22].
During the hydrosilylation reaction, silicon-hydrogen (Si-H) bond is added across the unsaturated carbon-carbon double bond (C=C) of an olefin, resulting in the formation of the silicon-carbon (Si-C) bond [23]. This reaction takes place in the presence of a catalyst and is often accompanied by various side reactions (such as dehydrogenative silylation, hydrogenation of olefins, isomerization and oligomerization of olefins) [23]. Extensive efforts were directed towards designing efficient and selective catalysts for this reaction; notably, hydrosilylation is regarded as one of the most important industrial applications of homogeneous catalysis [24]. Platinum-based (Pt-based) catalysts are primarily used in various industrial applications [24] and are shown to significantly reduce the side products. For example, a Pt-N-heterocyclic carbene (Pt-NHC) catalyst is reported to reduce the side products to a 1% yield, selectively forming primary alkylsilanes [23]. In our DPD framework, we neglect the formation of any side products and only consider the formation of silanes since this reaction is shown to be a dominant reaction of the Ptcatalyzed hydrosilylation [23,24].
Herein, we focus on modeling the phase separation coupled with crosslinking reactions in binary and ternary polymer blends containing linear preceramic precursors (PHMS) and sacrificial linear polymer chains. We characterize the phase separation process hindered by the crosslinking and focus on the effects of chemical nature and degree of polymerization of the sacrificial polymer on the evolution of the characteristic length scale in the system. We use the dissipative particle dynamics (DPD) approach [25][26][27], which is an efficient mesoscale approach utilizing soft repulsion potentials. DPD has been used to model a variety of polymer systems [28][29][30][31][32][33][34][35] and has been adapted to model various reactions within polymer solutions, melts, and polymer networks [36][37][38][39][40][41][42][43][44]. In what follows we first introduce our modeling approach and model parameters. We then utilize this approach to model the phase separation and crosslinking due to the hydrosilylation reaction. Specifically, we first focus on the phase separation in the binary system containing PHMS blended with v-PDMS. Prior experiments have shown [21] that v-PDMS can serve both as a solvent and as a crosslinker and thereby is effectively acting as templating agent controlling the pore size. We characterize the effects of the degree of polymerization and composition on the phase separation and resulting characteristic length scale in the systems. We then consider ternary systems, in which we introduce a second sacrificial non-reactive component. By considering methyl-terminated PDMS (m-PDMS), polyacrylonitrile (PAN), and poly(methyl methacrylate) (PMMA), we probe the effect of the chemical nature of this third non-reactive component on the characteristic features of the resulting morphologies formed.

Dissipative Particle Dynamics Approach
We first briefly introduce the DPD approach we use to model the phase separation and crosslinking reaction; more details of the DPD methodology can be found in the original publications [25][26][27]. The beads in DPD represent groups of atoms; the motion of these beads is governed by Newton's equations of motion [27]: where r i , v i , and p i = mv i are the position, velocity, and momentum of a particle i of mass m, respectively, and F ij is the pairwise additive force exerted on the particle i by the particle j. The summation is taken over all the particles within the chosen cut-off distance r c , which effectively introduces a characteristic length scale of the model. The total force acting between the non-bonded beads consists of three contributions [27]: ij is a purely repulsive conservative force, F D ij is dissipative, and F R ij is a random force. We chose the conservative force based on the soft repulsion potential [27]: where a ij is the maximum repulsion between the beads i and j, r ij = r ij is the distance between the centers of these beads, r ij = r i − r j , and e ij = r ij /r ij . The repulsion parameter between the beads of the same type, a ii , is typically calculated based on the degree of coarsegraining [27,45]. Herein, we chose a ii = 78 k B T/r c ; this parametrization is derived based on the compressibility of water and corresponds to coarse-graining of three water molecules into a single bead [45], so that the volume of a single bead is approximately 90 Å 3 . The repulsion parameter for the dissimilar beads, a ij , is chosen based on the affinity between the moieties these beads represent as detailed below. Remaining contributions to the total force read [27] F D ij = −γω D r ij e ij ·v ij e ij and F R ij = σω R r ij ζ ij ∆t −1/2 e ij , respectively, where γ and σ are the strengths of the dissipative and random forces, v ij = v i − v j is the relative bead velocity, and ∆t is the simulation time step. The ζ ij is a symmetric Gaussian distributed random variable with zero mean and unit variance, < ζ ij (t) >= 0, ζ ij (t)ζ ij (t ) = (δ ik δ jl + δ il δ jk )δ(t − t ). The F D ij and F R ij are coupled through the fluctuation-dissipation theorem [26] for which the following two conditions are imposed [27]: ω D r ij = ω R r ij 2 , and σ 2 = 2γk B T/m, where k B is the Boltzmann constant and T is temperature. We chose the weight function [27] ω R r ij = 1 − r ij /r c for r ij ≤ r c . Within the polymer strands, the beads are connected by harmonic bonds with the interaction potential U bond = K bond r ij − r b 2 /2, where K bond = 10 3 k B T/r 2 c and r b = 0.7 r c are chosen [46] as the spring constant and the equilibrium bond length, respectively.
To minimize unphysical topological crossings of bonded polymer chains (a known limitation of DPD), we adapted the modified Segmental Repulsive Potential (mSRP) [47][48][49] formulation. Notably, the mSRP DPD formulation captures the effects of entanglements in polymer melts [47]. Within the mSRP DPD framework, pseudo beads are introduced at the center of each bond for all the bonds initially present in the system; the pseudo beads interact with other pseudo beads with a soft repulsion potential defined as [47]: where b defines the strength of the interaction between the pseudo beads, d ij = d ij is the distance between these beads located at the centres of the bonds, e S ij = d ij /d ij , and d c is the cut-off distance for the mSRP interactions. Herein, the mSRP framework is modified [48] to create additional pseudo beads for the bonds created during crosslinking. For the mSRP potential, we set [47] b = 80 k B T/r 2 c and d c = 0.8 r c ; these parameters are shown to effectively minimize the bond crossings [47,49]. The strength of the random force is chosen as σ = 3, the number density in simulations is ρ = 3/r c 3 , and the time step is ∆t = 0.02 τ, where τ is the reduced unit of time [27]. The r c , temperature, and the mass of a DPD bead are set at 1.0 in reduced DPD units [27,45]. With the above choice of the degree of coarse-graining, the dimensionless units of length and time correspond to [45] r c ≈ 0.65 nm and τ ≈ 88 ps, respectively. All the results below are provided in the reduced DPD units. PHMS is represented by a linear chain that consists of 30 DPD beads (N PHMS = 30) and we vary the number of beads comprising the v-PDMS chains, N v−PDMS . The size of the simulation box is L × L × L, where L = 40 in reduced DPD units as provided above, periodic boundary conditions are applied in all directions. With the chosen density of 3, the total number of beads in the system is 1.92 × 10 5 . Prior to the production runs, the blends are equilibrated for 6 × 10 6 simulation timesteps with the same value of the repulsion coefficient (a ij = 78) chosen between all the beads within each system; no reactions occur during the equilibration. At the initial timestep of the production run, the repulsion parameters are modified between the beads using the values of a ij provided below. The equations of motion are integrated using the DPD module as implemented in the LAMMPS simulation package [50,51] with mSRP code [47] and the most recent modification of the mSRP code allowing one to use mSRP approach simultaneously with breaking or creating polymer bonds [48,49]. All simulation snapshots were visualized using the Visual Molecular Dynamics (VMD) software [52], and Wolfram Mathematica software [53] was used to generate all the plots.

Defining Repulsion Parameters a ij between Dissimilar Beads
The repulsion parameter for the dissimilar beads is chosen based on the affinity between the moieties these beads represent as [27] a ij = a ii + 3.27χ ij for the chosen density, where χ ij is the Flory-Huggins interaction parameter. The values of χ ij can be estimated using the respective solubility parameters δ i and δ j as is the molar volume of solvent and R is the gas constant. Based on direct comparison with a number of experimental studies, a modified empirical expression for χ ij is often found to be more reliable for polymer solutions [54,55]:  Table S1 of Supplementary Materials. The solubility parameter of PHMS was estimated using the group molar attraction constants, F z , as [55][56][57] density, M is the molar weight of the repeat unit, and the summation is taken over all groups of the repeat unit. The relevant values of F z are listed in Table S2 of Supplementary Materials. The same repulsion parameters are used for all the beads within each polymer chain, it is also assumed that changing a bead state from reactive to non-reactive does not affect the beads' affinity. Finally, the values of χ ij calculated based on the solubility coefficients and the respective repulsion coefficients a ij are provided in Tables S3 and  S4 of Supplementary Materials, respectively. We take the temperature T = 298.15 K at which the reference polymer system of interest in this study was crosslinked in the prior experiments [21]. Notably, with the estimated value of a PDMS−PHMS = 79.12, single beads of PHMS are miscible with PDMS, while polymer chains of PHMS and PDMS of the chosen length are immiscible. Hence this choice of the interaction parameter is consistent with experimental observations, showing both initial compatibility between the PDMS and PHMS molecules, and the phase separation between these moieties for sufficiently high molecular weight of PHMS [19,58].

Introducing Hydrosilylation Reaction within DPD Framework
We model the crosslinking between PHMS and v-PDMS chains due to hydrosilylation [21,58] as a stochastic process, similar to the approach previously used within the DPD framework to model polymerization of linear chains [36,38] and various reactions within polymer networks [37,39,40]. The mapping of the group of atoms within the polymer chains into the coarse-grained DPD beads is shown in Figure 1a. In this study, all initial components are linear chains (Figure 1b) composed of N X beads, where X indicates the type of chemical species. Correspondingly, N X represents the degree of polymerization of the individual chains prior to crosslinking. For example, a PHMS chain with N PHMS = 30 consists of 28 beads of type 2 and two beads of type 6 ( Figure 1a,b). The end vinyl group of v-PDMS and the repeat unit of PHMS are represented by a single DPD bead each, as shown in Figure 1a. The effective reactivity of these beads is set by selecting their states (reactive or non-reactive). The reactive moieties are classified as type 1 and type 2 beads for v-PDMS and PHMS, respectively (Figure 1a). Once a bond forms between the type 1 and type 2 beads, they become non-reactive and their type is changed to type 3 and type 4, respectively ( Figure 1d). For simplicity of representation, we do not highlight the state of reactive beads in the snapshots below. The crosslinking reaction is attempted at each reaction time step τ r , which is ten times larger than the DPD time step (τ r = 10∆t), same as in previous DPD simulations of reactive systems [36][37][38][39][40].  Table S4 of Supplementary Materials.
During each reaction step, for each reactive bead all complementary reactive beads that are within the reaction interaction distance [36][37][38][39][40], r int = 0.7, are identified as possible bond partners. Further, within this list of possible bond partners, the closest bead is identified as the sole bond partner for the particular bead. The reaction probability is only evaluated for the pairs in which both beads are identified as respective sole partners. A random number between 0 and 1 is generated; the bond is formed if this number is smaller than the chosen probability of reaction, P c . Each successful reaction results in an irreversible bond formation, with the bond energy given by U bond defined above; the state of the two beads involved in the crosslinking is correspondingly changed to non-reactive. We assume that changing the bead state from reactive to non-reactive does not affect the affinity of this bead. We set P c = 10 −3 as the reference value of the reaction probability. As we show below, this choice of the reference value of P c results in the first order reaction kinetics of the consumption of the vinyl groups provided that there is an abundant number of PHMS with respect to the vinyl groups in the system. Notably, prior experimental studies showed the consumption of silicon-hydrogen groups to obey first-order kinetics [59][60][61]. In experiments, an increase in the catalyst concentration or reaction temperature speeds up the hydrosilylation reactions [61]. Finally, prior to all the production runs in this study, the systems are equilibrated as described above; an example of an equilibrated system with f v−PDMS = 0.5, N PHMS = 30, and N v−PDMS = 100 is shown in Figure 1c.

Hydrosilylation Reaction Arrests Domain Growth in Blends Containing PHMS
In the first series of simulations, we characterized the time evolution of the 50/50 reactive blend (the fraction of v-PDMS, f v−PDMS = 0.5) containing 3200 PHMS chains with the degree of polymerization N PHMS = 30 and 960 v-PDMS chains with N v−PDMS = 100. Representative simulation snapshots are shown in Figure 2a; the time steps corresponding to these snapshots are marked by the red dots on the arrow below the images. Initially, the chains are well intermixed and equilibrated as detailed in Materials and Methods section. The phase separation between the PHMS and v-PDMS chains results in the aggregation of the chains of the same type with time as evident in these snapshots and in the corresponding increase in the characteristic domain size quantified below. The sacrificial phase is shown in green (bead type 5) and blue (bead types 1 and 3), and the matrix phase (PHMS chains) is shown in dark pink (bead types 2 and 4) and black (bead type 6), respectively (please refer to Figure 1a for definition of bead type). Note that in the absence of the hydrosilylation reaction, macroscopic phase separation is observed in the same system with two large domains forming upon equilibration ( Figure S1a in Supplementary Materials). The hydrosilylation reaction reduces the mobility of polymer chains and arrests the macroscopic phase separation, correspondingly limiting the characteristic domain size.
To quantitatively characterize evolution in this system, we track the time evolution of the characteristic size within the system, l(t), and a fraction of vinyl groups that remain reactive (beads of type 1 in the scheme given in Figure 1a) with respect to the initial number of the vinyl groups available. In the multi-component system, it is convenient to use partial radial distribution functions (RDFs) g αβ (r, t) to assess the structural properties during the evolution and upon reaching an equilibrium, where the indexes α and β refer to the types of beads in the multicomponent system. The g αβ (r, t) for the beads of types α and β is computed via the beads' coordinates as [62]: where ρ is the total number density of the system of N beads; c α,β = N α,β /N denotes the fractions of beads of types α and β in the multicomponent system, r i,j are the coordinates of beads α and β, respectively; δ(. . .) is the Dirac δ-function, and . . . denotes an ensemble average at the given time instant t. It is instructive to group all the v-PDMS beads in our systems (bead types 1, 3, and 5), and all the beads comprising PHMS chains (bead types 2, 4 and 6). We will refer to these groups of beads as sacrificial (denoted by the subscript "s") and as matrix beads (denoted by the subscript "m"), respectively. To characterize the domains formed by the sacrificial beads, we track the time evolution of the partial RDFs of all the sacrificial beads g ss (r, t) in the simulations below. Notably, the total RDF in this system can be calculated as g(r) = c 2 ss g ss (r, t) + 2c ss c mm g sm (r, t) + c 2 mm g mm (r, t), where g mm (r, t) is the partial RDF for all the beads comprising PHMS chains. The inset shows evolution at early times.
To quantitatively characterize evolution in this system, we track the time evolution of the characteristic size within the system, ( ), and a fraction of vinyl groups that remain reactive (beads of type 1 in the scheme given in Figure 1a) with respect to the initial number of the vinyl groups available. In the multi-component system, it is convenient to use partial radial distribution functions (RDFs) ( , ) to assess the structural properties during the evolution and upon reaching an equilibrium, where the indexes and refer to the types of beads in the multicomponent system. The ( , ) for the beads of types and is computed via the beads' coordinates as [62]: where is the total number density of the system of beads; , = , ⁄ denotes the fractions of beads of types and in the multicomponent system, , are the coordinates of beads and , respectively; (… ) is the Dirac -function, and 〈… 〉 denotes an ensemble average at the given time instant . It is instructive to group all the v-PDMS beads in our systems (bead types 1, 3, and 5), and all the beads comprising PHMS chains (bead types 2, 4 and 6). We will refer to these groups of beads as sacrificial (denoted by the subscript " ") and as matrix beads (denoted by the subscript " "), respectively. To characterize the domains formed by the sacrificial beads, we track the time evolution of the partial RDFs of all the sacrificial beads ( , ) in the simulations below. Notably, the total RDF in this The shifted partial RDFs of the sacrificial component, g ss (r, t) − 1, at the time instants corresponding to the morphologies in the snapshots in Figure 2a, are plotted in Figure 2b. All g ss (r, t) − 1 curves exhibit primary peaks at the separation r ≈ 0.7 corresponding to the bond length between the beads of the polymer chains set by the choice of the harmonic bond potential (see Materials and Methods section above). There are also additional smaller amplitude oscillations in g ss (r, t) at relatively short distances, which can be attributed to the packing of the nearest-neighbor and next-nearest-neighbor [63] beads. At all distances beyond this small distance, g ss (r, t) − 1 attains approximately zero values at initial times (t = 0, solid black line in Figure 2b), indicating that the sacrificial beads are essentially uniformly distributed at these length scales.
The magnitude of g ss (r, t) − 1 increases with time until the system reaches an equilibrium, confirming the aggregation of the sacrificial beads into distinct domains during the phase separation arrested by the hydrosilylation reaction. The characteristic length scale within this system can be defined in a straightforward manner via zero crossing [64][65][66] of g ss (r, t) − 1. Specifically, herein we quantify the characteristic length scale, l(t), via first value of r at which the function g ss (r, t) − 1 crosses zero from a positive to a negative value (Figure 2b). Alternatively, the characteristic length scale can also be calculated from the first minimum [66] of g ss (r, t) at distances exceeding the short distance correlations discussed above. It had been previously shown that the characteristic length scales calculated from zero crossings of g αα (r, t) − 1 obey expected scaling during spinodal decomposition [64,66]. Further, it had also been shown that the ratio between the characteristic length scales corresponding to the minima of RDF and those calculated from the zero crossing of the shifted RDF remains approximately constant in various systems undergoing phase separation [64,66]. This is also the case for the simulation in Figure 2a, as shown in Figure S2 of Supplementary Materials.
The characteristic length scale l(t) (in blue, left axis) and the number of the reactive vinyl groups (beads of type 1), N I (t), with respect to the initial number of these beads, N I (0) (in red, right axis), are given in Figure 2c. The value of N I (0) is defined by the number of v-PDMS chains in the system and can be expressed as: where ρ and L are the dimensionless beads number density and the linear size of the simulation box as defined in Models section. The open circles on this plot mark the time instants that are chosen in the snapshots in Figure 2a, the inset shows the time evolution of both l(t) and N I (t)/N I (0) at early times. These results clearly show that the reactive vinyl beads are consumed significantly faster than the time required for the phase separation process to reach equilibrium and correspondingly for l(t) to reach its value at saturation, l. The observed ratio N I (t)/N I (0) decays exponentially; however, the apparent reaction rate calculated via exponential fitting is lower due to the contribution from the diffusion processes than that in the case of an abundant number of PHMS ( Figure S3a). In the latter case ( f v−PDMS = 0.1), the consumption of the vinyl groups obeys the first order reaction kinetics with the reaction rate ≈ P c /τ r due to the abundant number of reactive PHMS beads available for each reactive vinal group bead ( Figure S3a). Notably, prior experimental studies showed that the consumption of vinyl groups follows the first-order kinetics [59][60][61]67]. The effect of the crosslinking reaction kinetics on the equilibrium morphology is probed by reducing the reaction probability from the reference value P c = 10 −3 to 5 × 10 −4 ( Figure S3b,c). The first order reaction kinetics with the reaction rate P c /τ r in the limit of an abundant number of PHMS ( f v−PDMS = 0.1) also holds for the lower value of P c probed. The evolution of the characteristic length scale ( Figure S3c) at these two probabilities shows that a slower crosslinking reaction results in somewhat faster increase in l at early times, but approximately the same equilibrium domain size is reached ( Figure S3c). More significant decrease in P c would result in an increase of the characteristic size, with macro phase separation observed in the limit of zero reaction rates, as shown in Figure S1.

Effect of Degree of Polymerization of the Sacrificial Chains, N v−PDMS
In the next series of simulations, we focus on the effect of the degree of polymerization of v-PDMS, N v−PDMS , on the average characteristic size upon reaching an equilibrium, l. We vary N v−PDMS within the following range: N v−PDMS = [10, 30, 48, 60, 80, 100, 120]. Representative snapshots upon reaching the equilibrium with selected values of N v−PDMS are shown in Figure 3a. The shifted partial RDFs, g ss (r) − 1, corresponding to the snapshots in Figure 3a are provided in Figure 3b for the same selected values of N v−PDMS . This plot shows that at the lowest degree of polymerization of the sacrificial component considered, N v−PDMS = 10, g ss (r, t) − 1 remains approximately zero beyond the initial small distances, indicating that the blend remains well intermixed on these length scales. An increase in N v−PDMS leads to the distinct increase in l. To understand this effect, recall that at a constant fraction of the sacrificial components ( f v−PDMS = 0.5 in this series of simulations), higher values of N v−PDMS correspond to a smaller number of chains and, respectively, to a lower number of reactive vinyl groups available, N I (0) Equation (4). Hence, an increase in N v−PDMS results in a decrease in the number of crosslinks formed in the system upon equilibration, in turn providing less constraints on the phase separation and resulting in formation of larger sacrificial domains. This tendency is clearly evident from the shift of zero crossing of g ss (r) − 1 to the higher values of r for larger N v−PDMS (Figure 3b). ered, = 10, ( , ) − 1 remains approximately zero beyond the initial small distances, indicating that the blend remains well intermixed on these length scales. An increase in leads to the distinct increase in ̅ . To understand this effect, recall that at a constant fraction of the sacrificial components ( = 0.5 in this series of simulations), higher values of correspond to a smaller number of chains and, respectively, to a lower number of reactive vinyl groups available, (0) (eqn. 4). Hence, an increase in results in a decrease in the number of crosslinks formed in the system upon equilibration, in turn providing less constraints on the phase separation and resulting in formation of larger sacrificial domains. This tendency is clearly evident from the shift of zero crossing of ( ) − 1 to the higher values of for larger (Figure 3b). Here and below, ̅ is averaged over the data taken within the last five frames (upon reaching an equilibrium) for four independent simulation runs; the error bars correspond to the standard deviation.
To characterize the average equilibrium length scale of the sacrificial domains, ̅ , herein and in the subsequent series of simulations, we output the RDF data every 5 × 10 time steps and average the value of calculated within the last five frames (upon reaching an equilibrium) for four independent simulation runs. An increase in with an increase in is shown in Figure 3c, with error bars corresponding to standard To characterize the average equilibrium length scale of the sacrificial domains, l, herein and in the subsequent series of simulations, we output the RDF data every 5 × 10 4 time steps and average the value of l calculated within the last five frames (upon reaching an equilibrium) for four independent simulation runs. An increase in l with an increase in N v−PDMS is shown in Figure 3c, with error bars corresponding to standard deviation. Hence our mesoscale simulations demonstrate formation of larger sacrificial domains with an increase in N v−PDMS ; similar trend showing an increase in the pore sizes with an increase in the molecular weight of v-PDMS was observed in a previous experimental study [21].

Effect of the Fraction of the Sacrificial Component, f v−PDMS
In the next series of simulations, we focus on the effect of the volume fraction of the sacrificial component,  Figure 4a. The corresponding snapshots of the sacrificial domains only are provided in Figure S4, where beads belonging to the distinct clusters are marked using the same color. These snapshots show that the sacrificial beads form large, interconnected clusters ( Figure S4). the sacrificial domains only are provided in Figure S4, where beads belonging to the distinct clusters are marked using the same color. These snapshots show that the sacrificial beads form large, interconnected clusters ( Figure S4). A decrease in results in a clear shift of the ( ) − 1 zero crossing towards higher values of (Figure 4b), corresponding to the formation of larger sacrificial domains. Secondary shallow peaks in ( ) are distinguishable at lower values of (inset in Figure 4b), indicating an average distance between the smaller domains. The characteristic average length scale in equilibrium, ̅ , decreases with an increase in for both values of probed in these series of simulations (data corresponding to = 30 and = 60 are shown in black and blue, respectively). This decrease is most pronounced for the highest , while the value of ̅ approximately saturates at lower . The vinyl groups in all the simulations are completely consumed via the crosslinking reactions, resulting in a smaller number of crosslinks in systems with lower . Hence, larger sacrificial domains are formed in the systems with the smaller number of crosslinks (Figure 4c).

Ternary Systems Incorporating Reactive and Non-Reactive Sacrificial Components
To further understand the dependence of the equilibrium characteristic length scale on the properties of the sacrificial components, in the remaining series of simulations we consider ternary systems encompassing a second sacrificial component in addition to the reactive v-PDMS chains. In the first example of the ternary system considered, we replace a fraction of v-PDMS chains in the system with methyl-terminated PDMS (m-PDMS) linear chains, which also serve as the sacrificial component [18,20,58]. We assume that m-PDMS chains only consist of the beads of type 5 and do not have any reactive beads. In this series of simulations, we keep the fraction of the matrix at f PHMS = 0.  (cases II, IV, and V) is found to be significantly more pronounced than a moderate increase in l with an increase in N m−PDMS (cases I-III). These trends are consistent with the trends observed in the binary systems (Figures 3c and 4c). Further, a direct comparison between the values of l in cases III and V indicates that using longer reactive chains and shorter non-reactive sacrificial chains (N v−PDMS = 100, N m−PDMS = 60, case V) results in significantly more pronounced increase in the characteristic length scale l than using longer non-reactive and shorter reactive chains (N v−PDMS = 60, N m−PDMS = 100, case III).
for v-PDMS, and (0.5 − ) and for m-PDMS. The details of each system are listed in Figure 5a and the corresponding snapshots for each system in equilibrium are shown in Figure 5b. The corresponding equilibrium characteristic length scales of the sacrificial domains ̅ (incorporating all the beads comprising m-PDMS and v-PDMS chains) for each system are compiled in Figure 5c. The effects of the degree of polymerization of the non-reactive component, , (cases I-III) and the degree of polymerization of the reactive component, , (cases II, IV, and V) are probed. The increase in the average value of ̅ with an increase in (cases II, IV, and V) is found to be significantly more pronounced than a moderate increase in ̅ with an increase in (cases I-III). These trends are consistent with the trends observed in the binary systems (Figures 3c and 4c). Further, a direct comparison between the values of ̅ in cases III and V indicates that using longer reactive chains and shorter non-reactive sacrificial chains ( = 100, = 60, case V) results in significantly more pronounced increase in the characteristic length scale ̅ than using longer non-reactive and shorter reactive chains ( = 60, = 100, case III). Specifically, we chose f v−PDMS = 0.3 , 0.4, and 0.5 in cases VI, II, and VII, respectively. Recall that the total fraction of the sacrificial component is fixed in these simulations ( f v−PDMS + f m−PDMS = 0.5), so that the case VII corresponds to the binary system considered above and is provided in Figure 5 for comparison. These results quantify the decrease in the characteristic length scale l with a relative increase in the fraction of the reactive sacrificial component, f v−PDMS .
Finally, we probe the effect of the chemical nature of the non-reactive sacrificial component on the phase-separated morphology and the characteristic length scale within the system. In the final series of simulations, PAN and PMMA linear chains are used as non-reactive sacrificial components in addition to the v-PDMS reactive chains. In effect, these series of simulations correspond to the case II described in Figure 5, but with m-PDMS chains replaced with PAN or PMMA chains with the same number of beads. The repulsion parameters a ij used between all the moieties are provided in Table S4 of Supplementary Materials. Note that PAN has stronger repulsion with both PHMS and v-PDMS compared to that of PMMA. Due to the miscibility and non-reactivity, m-PDMS chains mix and spread within the simulation box along with v-PDMS chains (bottom row in Figure 6a, third column). In contrast, due to relatively low affinity between PAN (and PMMA) with both PHMS or v-PDMS, we observe that PAN and PMMA chains aggregate and form isolated approximately spherical clusters dispersed in the blends (Figure 6a, bottom row, first two columns). These results show that increasing the repulsion between the nonreactive sacrificial component and matrix (or, similarly, between the non-reactive sacrificial component and reactive sacrificial component) can lead to the pronounced decrease in the characteristic length scale l (Figure 6b). Note that l in Figure 6b is calculated based on the RDFs of all the sacrificial beads ( Figure S5). To summarize, our results point out that introducing non-reactive sacrificial polymer chains could provide a facile method to tailor the morphology of the pores.
system. In the final series of simulations, PAN and PMMA linear chains are used as nonreactive sacrificial components in addition to the v-PDMS reactive chains. In effect, these series of simulations correspond to the case II described in Figure 5, but with m-PDMS chains replaced with PAN or PMMA chains with the same number of beads. The repulsion parameters used between all the moieties are provided in Table S4 of Supplementary Materials. Note that PAN has stronger repulsion with both PHMS and v-PDMS compared to that of PMMA. Due to the miscibility and non-reactivity, m-PDMS chains mix and spread within the simulation box along with v-PDMS chains (bottom row in Figure 6a, third column). In contrast, due to relatively low affinity between PAN (and PMMA) with both PHMS or v-PDMS, we observe that PAN and PMMA chains aggregate and form isolated approximately spherical clusters dispersed in the blends (Figure 6a, bottom row, first two columns). These results show that increasing the repulsion between the non-reactive sacrificial component and matrix (or, similarly, between the non-reactive sacrificial component and reactive sacrificial component) can lead to the pronounced decrease in the characteristic length scale ̅ (Figure 6b). Note that ̅ in Figure 6b is calculated based on the RDFs of all the sacrificial beads ( Figures S5). To summarize, our results point out that introducing non-reactive sacrificial polymer chains could provide a facile method to tailor the morphology of the pores.

Conclusions
We utilized modified the Segmental Repulsive Potential (mSRP) Dissipative Particle Dynamics (DPD) approach to simulate the phase separation coupled with hydrosilylation in binary and ternary blends. Our approach introduces the first approach to capture hydrosilylation reaction along with the phase separation of preceramic precursors on the mesoscale. Our results shows that the phase separation between polyhydromethylsiloxane (PHMS) and the vinyl-terminated polydimethylsiloxane (v-PDMS) is arrested due to the hydrosilylation reaction. We characterized the dynamics of binary PHMS/v-PDMS blends undergoing hydrosilylation and ternary blends in which a fraction of the reactive sacrificial component (v-PDMS) is replaced with the non-reactive sacrificial component. As the sacrificial non-reactive component in these ternary blends, we chose methyl-terminated PDMS (m-PDMS), polyacrylonitrile (PAN), and poly(methyl methacrylate) (PMMA). Our results clearly demonstrate that the morphology of the sacrificial domains in the phase-separated network can be tailored by tuning the composition, chemical nature, and the degree of polymerization of the sacrificial component. Our results show a decrease in the characteristic length scale of the sacrificial domains with either an increase in the fraction of v-PDMS chains or a decrease of the degree of polymerization of these chains. We also show that an addition of a non-reactive sacrificial component introduces facile means to control the morphology by varying the fraction, degree of polymerization, or the chemical nature of this component. In particular, a decrease in affinity between the non-reactive sacrificial component and the preceramic precursors results in a decrease in the characteristic length scale of the sacrificial phase in ternary blends.

Supplementary Materials:
The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/nano12183117/s1, Table S1: Solubility parameters and molar volumes of chosen monomers; Table S2: The group molar attraction parameters; Table S3: Flory-Huggins interaction parameters χ ij (T = 25 • C); Table S4: Repulsion parameters α ij between all the types of beads; Figure S1: Phase separation in the non-reactive blend with f v−PDMS = 0.5 and N v−PDMS = 100 (no hydrosilylation reactions); Figure S2: Comparison between characteristic length scales defined by l and l min in the system shown in Figure 2 of the main text ( f v−PDMS = 0.5 and N v−PDMS = 100); Figure S3: Effects of reaction probability, P c ; Figure S4: Snapshots of sacrificial domains in the blends with N v−PDMS = 60 and various f v−PDMS ; Figure S5: Shifted RDFs in ternary systems corresponding to Figure 6 at equilibrium. References [52,57,[68][69][70][71][72][73][74] are cited in the supplementary materials.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.